How to make a random draw from a mgcv GAM model Hi there statistical wizards!
I have a (maybe) minor problem, where I would like to use the predicted fit (fit) and standard error (se.fit) from a GAM model to take a random draw. Ideally I would like to know how to convert this standard error to a standard deviation so as to use rnorm(mean, sd) on each fitted value.
To compute the standard deviation from the standard error we need an N value, but I'm not sure how to compute this.
So far my code looks like this:
library(mgcv)
library(ggplot2)


eReefs.chl = c(0.24596960, 0.15762042, 0.10942775, 0.12170488, 0.08990742, 0.07395359, 0.08712383, 0.08182129, 0.08367566)
SON = c(-1.7, -1.1,  0.2, -0.2,  0.4,  2.4, -0.7, -0.7,  0.7)
years = 2010:2018

dat = data.frame(cbind(eReefs.chl, SON, years))

gamfit = mgcv::gam(eReefs.chl~s(SON, bs = "cs", k=4), data = dat)
mgcv::gam.check(gamfit)
plot(gamfit, residuals = T, pch=19, cex=1, xlab = "SON Nino3.4", ylab = "Median Chlorophyll Concentration", 
     all.terms = T, shade=T, shift = coef(gamfit)[1])

newdat = data.frame(SON = c(0.3,  0.8, -0.2,  0.1,  0.7, -1.0, -0.4,  2.3, -1.4, -1.3, -0.6, -0.3,  1.2,  0.3,  0.7, -0.3,
                            0.7, -1.4, -0.4,  1.0, -1.7, -1.1,  0.2, -0.2,  0.4,  2.4, -0.7, -0.7, 0.7), years = 1990:2018)

pred = predict(gamfit, newdata = newdat, se.fit = T)

# Now I want to draw from a particular year to run stochastic simulations, 
# so I want to create  a standard deviation from the SE, but I don't know N

#stddev = pred$se.fit*n
#rnorm(1, pred$fit[1], stdvev[1])

I would really appreciate some help on this, or an alternative approach. For context I am fitting about 1500 of these gam's so thats why I want to use the fit and standard error to draw from later after the model has been fit.
 A: That standard error is the uncertainty in the estimate for the expected/fitted value. For the variance of the posterior distribution, you want the estimated value of the variance (or scale parameter) for the model you fitted. You fitted something like
$$y_i = \beta_0 + f(x_{1i}) + \varepsilon_i$$
where $\varepsilon_i \sim N(0, \sigma^2)$, where the variance of the residuals (data) is estimated $\hat{\sigma}^2$. Another way of writing this is
$$Y \sim N(\hat{\mu}_i, \hat{\sigma}^2)$$
$$\hat{\mu}_i = \beta_0 + f(x_{1i})$$
Notice that $\hat{\sigma}^2$ is separate from the predicted value $\hat{\mu}_i$, and is an estimate of the variance of the data not the uncertainty in the predicted value. It is this $\hat{\sigma}^2$ that you want to get.
One way to do this for the Gaussian models you fitted requires three things:


*

*the expected value, which will be the value you pass to mu

*the estimated variance/scale parameter,

*observation weights, which we presume to be 1 here.


You already have mu as this is the fitted value for the new observation(s). The estimated variance/scale parameter can be extracted from the model as component sig2 of the fitted model. (This may not exist in all GAMs fitted by mgcv, so we can try to extract it as
summary(model)[["dispersion"]]

but for these Gaussian GAMs, model$sig2 will suffice, where model is your fitted GAM model.) Let's call this variance/scale parameter sig2.
In the case of weights, the standard deviation for use in rnorm() is
sqrt(sig2/wt)

where wt is the weight for that observation. As I presume you will want a weight of 1, then you can just use
sqrt(sig2)

This process extends to all of the models the gam() is capable of fitting — although the detail will be different.
You can get a function that given mu, sig2 and wt will generate draws from the correct distribution for a GAM using the fix.family.rd() function in mgcv. You pass it the family function for the fitted model & extract the newly added rd component.
rdfun <- fix.family.rd(family(model))$rd

For Gaussian models this is:
> fix.family.rd(gaussian())$rd

function (mu, wt, scale) 
{
    rnorm(mu, mean = mu, sd = sqrt(scale/wt))
}
<bytecode: 0xd0e81b8>
<environment: 0x950aa40>

Which is where we see how this is implemented/computed.
As for the detail for other models, if you are unsure, consult the rd functions mgcv adds. For example, for the Poisson:
> fix.family.rd(poisson())$rd
function (mu, wt, scale) 
{
    rpois(length(mu), mu)
}
<bytecode: 0xd0e95a0>
<environment: 0xcd55138>

we see that we only need mu, the expected value to make draws from the posterior distribution, but for the negative binomial we need something extra, the $\theta$ parameter of the (particular parameterisation of the) negative binomial distribution
> fix.family.rd(nb())$rd
function (mu, wt, scale) 
{
    Theta <- exp(get(".Theta"))
    rnbinom(mu, size = Theta, mu = mu)
}
<environment: 0xe4f9c80>

