R-helpers:
I am using R 2.5 on Windows XP, packages all up to date. I have run
into an issue with the MIcombine function of the mitools package that I
hoped some of you might be able to help with. I will work through a
reproducible example to demonstrate the issue.
First, make a dataset from the pbc dataset in the survival package
---------------
# Make a dataset
library(survival)
d <-
pbc[,c('time','status','age','sex','hepmeg','platelet',
'trt', 'trig')]
d[d==-9] <- NA
d[,c(4,5,7)] <- lapply(d[,c(4,5,7)], FUN=as.factor)
str(d)
summary(d)
---------------
Second, since there is missing data for several (but not all) of the
variables, investigate the patterns.
---------------
library(Hmisc)
na.pattern(d)
clus <- naclus(d, method='complete')
par(mfrow=c(2,2))
naplot(clus, which='all')
print(clus)
detach(package: Hmisc)
---------------
After examining the missing data patterns, impute 10 datasets using the
amelia function from the Amelia package. Check the densities of the
continuous variables to make sure they make sense.
---------------
library(Amelia)
am.imp <- amelia(d, m=10, p2s=1, startvals=1, write.out=F,
noms=c('status','sex','hepmeg','trt'))
compare.density(data=d, output=am.imp, var='trig')
compare.density(data=d, output=am.imp, var='platelet')
---------------
Since everything looks ok, fit Cox models to each of the 10 imputed
datasets using the functions of the mitools package.
---------------
library(mitools)
miset <- imputationList(list(am.imp[[1]],am.imp[[2]],am.imp[[3]],
am.imp[[4]],am.imp[[5]],am.imp[[6]],am.imp[[7]],am.imp[[8]],
am.imp[[9]],am.imp[[10]]))
mifit <- with(miset, coxph(Surv(time, status) ~ age + sex +
hepmeg + platelet + trt + trig))
mifit
---------------
The "mifit" object shows the individual coxph fits for each imputed
dataset. Now we have to pool these fitted models to obtain an overall
model. The code and result are shown below.
---------------
fit.am <- MIcombine(mifit)
summary(fit.am)
Multiple imputation results:
with.imputationList(miset, coxph(Surv(time, status) ~ age + sex +
hepmeg + platelet + trt + trig))
MIcombine.default(mifit)
results se (lower upper) missInfo
age 0.035548792 0.0082506946 0.019373545 0.0517240397 4 %
sex1 -0.070760613 0.2563372831 -0.580309741 0.4387885156 34 %
hepmeg1 0.932378808 0.2026274576 0.532555416 1.3322021993 23 %
platelet -0.001757899 0.0009480636 -0.003620446 0.0001046469 14 %
trt2 0.137413162 0.1924230007 -0.243815288 0.5186416117 29 %
trig 0.003979287 0.0012010053 0.001607266 0.0063513078 25 %
---------------
The question I have concerns the meaning of this result. The missInfo
column of the summary function suggests that age was missing data, when
in fact it was not. Is there a different interpretation for the missInfo
column?
Thanks in advance for any help.
Brant Inman
