I am very new to R and not very experienced in statistics. I have this general question regarding applying Generalized Additive Models (GAM) in multiple imputation dataset.

I used R package mice for multiple imputation and mgcv for GAM. But I have no idea how to combine the results from each imputed data set.

Does anyone have any good idea about this?

Thanks a lot!

  • $\begingroup$ I have never used mice but I used a lot the GAM package. If you decide to move to the gam package maybe I can help you $\endgroup$
    – Donbeo
    May 27, 2014 at 20:05
  • $\begingroup$ Yea, I'm actually trying both. I'm all ears! $\endgroup$ May 27, 2014 at 20:08
  • $\begingroup$ Can you write in the question what did you try with the package GAM. Maybe I can help you $\endgroup$
    – Donbeo
    May 27, 2014 at 20:08
  • $\begingroup$ e.g. This is the model I wanna run fit=gam(y ~ s(x,5),data=df). I have imputed 5 data sets. I wanna run the same model in each of the data set. say the imputed dataset is stored in dat.imp. Then the code would be, with(dat.imp, gam(y~s(x,5)). When it returns, it shows you the results from each dataset. But I need to think of a way to combine them. $\endgroup$ May 27, 2014 at 20:11
  • $\begingroup$ sorry no clear to me. you have different datasets df1,df2 ,.. ? in each dataset you have an input value x and an ouput value y? $\endgroup$
    – Donbeo
    May 27, 2014 at 20:13

2 Answers 2


First, GAM's are basically a combination of GLMs, splines, and ridge regression (loosely speaking). You might be well advised to work on your understanding of these three things before attempting to work with their combination.

Second, I've done multiple imputation with GAMs, but not with mice.

In general, you're interested in fitting a model: $$ y = f(\mathbf{X}) + \epsilon $$ (in the gaussian, identity-linked case). $\mathbf{X}$ has holes in it -- missing data. Without going into the detail about how multiple imputation works (Gelman and Hill have a very good and less-technical chapter explaining missing data imputation), what it does is create many different $\mathbf{X}$'s, call them $\mathbf{X}_m$ where $m$ indexes imputed datasets from 1 to the number of imputations, with the missing values filled in with plausible guesses about what the missing values might be -- based on correlations in your dataset. You then fit your model $$ y = f(\mathbf{X}_m) + \epsilon $$ to each of those $m$ partially-imputed datasets. Since you're fitting $m$ models to $m$ different datasets, you're going to get $m$ different vectors of regression coefficients. These are combined according to "Rubin's Rules" (after the statistician who invented MI). Basically you average the coefficients; $$\hat\beta = \frac{1}{M}\displaystyle\sum_{m=1}^M \hat\beta_m$$ The variance-covariance matrix $\hat V_\beta$ of the estimated parameters is calculated by first averaging variance-covariance matrices, and then adding a correction to account for variation between imputation models: \begin{equation} \hat V_{\hat\beta} = W + \left(1 + \frac{1}{M} \right)B \end{equation} where $W = \frac{1}{M}\displaystyle\sum_{m=1}^M \widehat{VCV}_m$, $\widehat{VCV}$ is the estimated variance-covariance matrix of the estimated parameters, and $B = \frac{1}{M-1}\displaystyle\sum_{m=1}^M \left(\hat\beta_m - \hat{\bar\beta}\right)\left(\hat\beta_m - \hat{\bar\beta}\right)^T$. This procedure inflates the standard errors on coefficients about which the imputation model is relatively less certain, either due to a lot of missing data, or due to a poorly-informative imputation model.

You can then use the vector of coefficients and the VCV as if they were gotten from a single complete-case model.

OP's question in a comment makes me realize that the below isn't quite right! Read EDIT2 below to see why

So what is different about doing this with a GAM? Only a couple of things. Your variables represented by smooth functions are associated with a number of coefficients (mgcv's default is 10). During estimation they are subject to ridge penalties, which smooths the estimated function. But otherwise they can be slotted into Rubin's rules just like the parametric coefficients. Each of the different $m$ models will have different estimated smoothing parameters, as they are estimated from the data. I don't think there is a problem with this, though I'd be interested to hear if someone has another perspective. You do want to take care that you pre-specify the location of the knots of each spline so that they are uniform across the $m$ models -- otherwise your coefficients won't be comparable, and combining them won't be appropriate.

Once you have calculated the combined coefficients and VCV, you can simply stuff them into a GAM object (one of your $m$ models). Functions that summarize data and make plots draw from those two objects, and will thereby use your imputation estimates rather than those from the $m$th model. There is code online illustrating the whole process from this paper.

EDIT I just saw the comment about how you've got your models fitted to your imputed datasets already. If you can coerce those into a list object with the 5 GAM models, then you can run something like the following to combine them:

    for (i in 2:reps){
bhat = bhat/reps

    for (i in 2:reps){
    	W = W+results[[i]]$Vp
    W = W/reps
    B= (results[[1]]$coeff-bhat) %*% t(results[[1]]$coeff-bhat)
    for (i in 2:reps){
    	B = B+(results[[i]]$coeff-bhat) %*% t(results[[i]]$coeff-bhat)

Vb = W+(1+1/reps)*B

    for (i in 2:reps){
dfr = dfr/reps

MI = results[[1]]

    MI$Vp = Vb
MI$df.residual = dfr

EDIT2 OP's question in a comment about p-values makes me realize that the above isn't quite right. It should get you appropriate coefficient and VCV estimates. But p-values are based on a reduced-rank Wald statistic that relies on the model matrix for its calculation. Obviously the model matrix will differ between different imputations. So if you take a single one of many gam objects (which will all have their own model frames) and stuff the VCV and coefficient vectors into it, you won't get the same result as you would if you chose a different one of your models fit to a different imputed dataset. I'm not sure how Rubin's Rules would be generalized to combine imputations here! Maybe a short-term hack would be to also average the model R matrices from each model? This should be about as valid as averaging the effective degrees of freedom (which may also be a hack).

  • $\begingroup$ This is very helpful information! Thank you very much! $\endgroup$ May 27, 2014 at 20:46
  • $\begingroup$ I have another question. In this case, how do you calculate the p-value for the smoothed term? $\endgroup$ May 27, 2014 at 22:38
  • 1
    $\begingroup$ Good question! I was about to post that you can simply rely on what summary.gam gives you after you replace the adjusted coefficient vector and VCV into the gam object. That is what I have done previously. But as I think about it, it's not right! The details of p-values for smooth functions are in Simon Wood's 2013 Biometrika paper. The reduced-rank wald statistic that summary.gam calculates needs the model R matrix (from a QR decomposition of the X matrix), as well as the VCV, the beta vector, and the residual degrees of freedom. Obviously the X matrix will vary between imputations. $\endgroup$ May 28, 2014 at 2:04

See https://stat.ethz.ch/R-manual/R-devel/library/mgcv/html/missing.data.html
Or enter help(mgcv::missing.data) in an R session. An approach that can be effective, with sample code on the help page, is then to

   substitute a simple random effects model in which the by variable 
   mechanism is used to set s(x) to zero for any missing x, while a 
   Gaussian random effect is then substituted for the ‘missing’ s(x).

Factors are required, one for each variable that has missing values, for use as missing value indicators, in each case with as many levels as there are missing values. The NAs are replaced in each case, by the mean for the relevant variable.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.