# How to summarize GAM model result from multiple imputation data in R

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.

Thanks a lot!

• 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 – Donbeo May 27 '14 at 20:05
• Yea, I'm actually trying both. I'm all ears! – user2816711 May 27 '14 at 20:08
• Can you write in the question what did you try with the package GAM. Maybe I can help you – Donbeo May 27 '14 at 20:08
• 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. – user2816711 May 27 '14 at 20:11
• 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? – Donbeo May 27 '14 at 20:13

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:

bhat=results[]$coeff for (i in 2:reps){ bhat=bhat+results[[i]]$coeff
}
bhat = bhat/reps

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

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

dfr=results[]$df.residual for (i in 2:reps){ dfr=dfr+results[[i]]$df.residual
}
dfr = dfr/reps

MI = results[]

MI$coefficients=bhat 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).

• This is very helpful information! Thank you very much! – user2816711 May 27 '14 at 20:46
• I have another question. In this case, how do you calculate the p-value for the smoothed term? – user2816711 May 27 '14 at 22:38
• 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. – generic_user May 28 '14 at 2:04