[R] Uncertainty propagation

[Ben Bolker]

Maayt <m.lupker <at> hotmail.com> writes:

[snip]
> My first intention was to use a kind of monte carlo routine and run the
> model many times by changing the power law parameters. These power laws were
> obtained by fitting data points under R. I thus have std error associated to
> them: alpha (±da) * WaterHight ^ beta (±db). Is it statistically correct to
> sample alpha and beta for each run by picking them from a normal
> distribution centered on alpha (resp. beta) with a standard deviation of da
> (resp. db) and to perform my statistics (mean and standrad edviation of the
> model result) on the model output?
> It seems to me that da and db are correlated in some way and by doing what I
> entended to, I would overestimate the final error of my model...

  How have you fitted the models?
  Many of the fitting procedures in R give you access not just to
the standard errors of the parameters, but also to their correlations/
covariances.  If you have this information, you can sample the pairs
of parameters from an appropriate multivariate normal distribution.
Typically you could do something like ...

params <- MASS::mvrnorm(1000,mu=coef(modelfit),Sigma=vcov(modelfit))
predictions <- apply(params,1,predictfun)