Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It only takes a minute to sign up.
Sign up to join this community
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!
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[[1]]$coeff
for (i in 2:reps){
bhat=bhat+results[[i]]$coeff
}
bhat = bhat/reps
W=results[[1]]$Vp
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)
}
B=B/(reps-1)
Vb = W+(1+1/reps)*B
dfr=results[[1]]$df.residual
for (i in 2:reps){
dfr=dfr+results[[i]]$df.residual
}
dfr = dfr/reps
MI = results[[1]]
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).
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.
