I am having problems with the MICE package in R, particularity with pooling the imputed data sets.
I am running a multilevel binomial logistic regression, with Level1 - topic (participant response to 10 questions on different topics, e.g. Darkness, Day) nested within Level2 - individuals.
The model is created using R2MLwiN, the formula is
> fit1 <-runMLwiN( c(probit(T_Darkness, cons), probit(T_Day, cons), probit(T_Light, cons), probit(T_Night, cons), probit(T_Rain, cons), probit(T_Rainbows, cons), probit(T_Snow, cons), probit(T_Storms, cons), probit(T_Waterfalls, cons), probit(T_Waves, cons)) ~ 1, D=c("Mixed", "Binomial", "Binomial","Binomial","Binomial", "Binomial", "Binomial", "Binomial", "Binomial", "Binomial" ,"Binomial"), estoptions = list(EstM = 0), data=data)
Unfortunately, there is missing data in all of the Level1 (topic) responses.
I have been using the mice package ([CRAN][1]) to multiply impute the missing values.
I can fit the model to the imputed datasets, using the formula > fitMI <- (with(MI.Data, runMLwiN( c(probit(T_Darkness, cons), probit(T_Day, cons), probit(T_Light, cons), probit(T_Night, cons), probit(T_Rain, cons), probit(T_Rainbows, cons), probit(T_Snow, cons), probit(T_Storms, cons), probit(T_Waterfalls, cons), probit(T_Waves, cons)) ~ 1, D=c("Mixed", "Binomial", "Binomial","Binomial","Binomial", "Binomial", "Binomial", "Binomial", "Binomial", "Binomial" ,"Binomial"), estoptions = list(EstM = 0), data=data)))
However, when I come to pool the analyses with the call code > pool(fitMI) it fails, with the Error:
Error in pool(with(tempData, runMLwiN(c(probit(T_Darkness, cons), probit(T_Day, :
Object has no coef() method.
I am not sure why it is saying there is no coefficient, as the analyses of the individual MI datasets provide both fixed parts (coefficients) and random parts (covariances)
Any help with what is going wrong would be much appreciated.
I should warn you that this is my first foray into using R and multilevel modelling. Also I know there is a MlwiN package ([REALCOM][2]) that can do this but I don't have the background to use the MLwiN software.
thanks johnny
Update - R reproducible example
Libraries used
library(R2MLwiN)
library(mice)
Subset of data `
T_Darkness <- c(0, 1, 0, 0, 0, 0, 0, 1, 0, 0, NA, 0, 0, 0, NA, 1, 0, NA,NA, 1, 0, 0, 0, 1, 0, 0, 0, NA, 0, 0, 0, NA, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, NA, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, NA, 1, 0)
T_Day <- c(0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, NA, 0, 0, 0, NA, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, NA, 0, 0, 0, 0, NA, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, NA, NA, 0)
T_Light <- c(0, 0, NA, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, NA, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, NA, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, NA, 0, 0)
T_Night <- c(0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, NA, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,NA, 0, NA, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, NA, 0, 0)
T_Rain <- c(1, 0, 0, 1, 1, 0, 0, NA, 0, 1, 0, 0, 1, 0, 0, 0, 0, NA, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, NA, 0, 0, 0, 0, 1, 0, 0, 0, NA, 1, NA, 0, 0, 0, 0, 1, NA, 1, 0, 0, 0, 0, 1, NA, 0, 0)
T_Rainbows <- c(1, 1, 1, 1, 0, 1, 0, 1, 0, 1, NA, 1, 1, 0, 0, 1, 0, NA, 0, 1, 0, NA, 0, 1, 0, 0, 0, 0, 0, NA, 0, 0, 0, NA, 1, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, NA, 1, 0, 1, NA, 0, 0, 1, 0, 1, 1, 1, 0, 1)
T_Snow <- c(0, 0, 1, 0, 0, 0, 1, 1, 0, 1, 0, NA, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, NA, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, NA, 0, 0, 1, NA, 1, 0, 1, 1, 0, 0, 0, 0, 0, NA, 0, 0, 0)
T_Storms <- c(0, 0, 0, 1, 1, 1, 0, 1, 0, 1, NA, 0, 0, 0, 0, 1, 0, NA, 0, 0, 1, 0, 0, NA, 1, 1, NA, 0, 0, NA, 0, 1, 0, NA, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, NA, 1, 0, NA, 0, 0, 0, 1, 1, 0, 1, NA, NA, 1)
T_Waterfalls <- c(0, 0, 0, 0, 0, 0, 0, NA, 0, 0, 0, 0, 0, 0, 0, NA, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, NA, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, NA, 0, 0, 0, 0, 0, NA, 0, 1, 0, NA, 1, 0, 1, 0, 0, 0, NA, 0, 0, 0, NA, NA, 0)
T_Waves <- c(0, 1, 0, 1, 1, 0, 1, NA, 0, 0, NA, 0, 0, 0, NA, 1, 0, 0, 0, 0, 1, 0, NA, 0, NA, 0, 0, NA, 0, 0, 0, 0, 0, 0, NA, 1, 0, 0, 0, 1, 0, 0, NA, 0, 1, 0, 0, 0, 0, 0, 1, 1, NA, 1, 1, NA, 0, 0, 0, NA, 0, 0, 0, NA, 0, 0)
data <- data.frame (T_Darkness, T_Day, T_Light, T_Night, T_Rain, T_Rainbows, T_Snow, T_Storms, T_Waterfalls, T_Waves)
data$cons <- 1
`
Data imputed using mice with
MI.Data <- mice(data,m=5,maxit=50,meth='pmm',seed=500)
This appears to be due to some of the model extraction methods in R2MLwiN not being correctly found, which should have been fixed in the recently released 0.8-2 version of the package. Running your example with this gives me the following results:
> pool(fitMI)
Call: pool(object = fitMI)
Pooled coefficients:
FP_Intercept_T_Darkness FP_Intercept_T_Day FP_Intercept_T_Light FP_Intercept_T_Night FP_Intercept_T_Rain
-0.9687210917 -1.0720602274 -0.9584792256 -1.1816471815 -0.7082406878
FP_Intercept_T_Rainbows FP_Intercept_T_Snow FP_Intercept_T_Storms FP_Intercept_T_Waterfalls FP_Intercept_T_Waves
-0.0455361903 -0.7537600398 -0.3883027434 -1.2365225554 -0.6423609257
RP1_var_bcons_1 RP1_cov_bcons_1_bcons_2 RP1_var_bcons_2 RP1_cov_bcons_1_bcons_3 RP1_cov_bcons_2_bcons_3
1.0000000000 0.0508168936 1.0000000000 0.2744663656 0.1625871509
RP1_var_bcons_3 RP1_cov_bcons_1_bcons_4 RP1_cov_bcons_2_bcons_4 RP1_cov_bcons_3_bcons_4 RP1_var_bcons_4
1.0000000000 0.0013987361 0.0576194786 0.0201622359 1.0000000000
RP1_cov_bcons_1_bcons_5 RP1_cov_bcons_2_bcons_5 RP1_cov_bcons_3_bcons_5 RP1_cov_bcons_4_bcons_5 RP1_var_bcons_5
-0.0220604800 0.1620389074 0.0956511647 -0.0242812764 1.0000000000
RP1_cov_bcons_1_bcons_6 RP1_cov_bcons_2_bcons_6 RP1_cov_bcons_3_bcons_6 RP1_cov_bcons_4_bcons_6 RP1_cov_bcons_5_bcons_6
0.2644620836 0.0555731133 0.1911445856 0.2584619522 0.1523280591
RP1_var_bcons_6 RP1_cov_bcons_1_bcons_7 RP1_cov_bcons_2_bcons_7 RP1_cov_bcons_3_bcons_7 RP1_cov_bcons_4_bcons_7
1.0000000000 0.1877118051 0.0872156173 0.2800109982 0.1433261335
RP1_cov_bcons_5_bcons_7 RP1_cov_bcons_6_bcons_7 RP1_var_bcons_7 RP1_cov_bcons_1_bcons_8 RP1_cov_bcons_2_bcons_8
-0.0006230903 0.1582182944 1.0000000000 -0.0749104023 0.1435756236
RP1_cov_bcons_3_bcons_8 RP1_cov_bcons_4_bcons_8 RP1_cov_bcons_5_bcons_8 RP1_cov_bcons_6_bcons_8 RP1_cov_bcons_7_bcons_8
0.0537744537 0.2291038185 0.2553031743 0.2716509402 0.1914017051
RP1_var_bcons_8 RP1_cov_bcons_1_bcons_9 RP1_cov_bcons_2_bcons_9 RP1_cov_bcons_3_bcons_9 RP1_cov_bcons_4_bcons_9
1.0000000000 0.1936145425 0.2835071683 0.0144172618 0.3326070011
RP1_cov_bcons_5_bcons_9 RP1_cov_bcons_6_bcons_9 RP1_cov_bcons_7_bcons_9 RP1_cov_bcons_8_bcons_9 RP1_var_bcons_9
0.1372590512 0.2854030728 0.0750594735 0.2545967996 1.0000000000
RP1_cov_bcons_1_bcons_10 RP1_cov_bcons_2_bcons_10 RP1_cov_bcons_3_bcons_10 RP1_cov_bcons_4_bcons_10 RP1_cov_bcons_5_bcons_10
0.3137466609 0.3498021364 0.2846792042 0.1126367375 0.2416045219
RP1_cov_bcons_6_bcons_10 RP1_cov_bcons_7_bcons_10 RP1_cov_bcons_8_bcons_10 RP1_cov_bcons_9_bcons_10 RP1_var_bcons_10
0.2137401104 0.1849118918 0.2134640366 0.6101759672 1.0000000000
Fraction of information about the coefficients missing due to nonresponse:
FP_Intercept_T_Darkness FP_Intercept_T_Day FP_Intercept_T_Light FP_Intercept_T_Night FP_Intercept_T_Rain
0.5714367 0.5714367 0.5714367 0.5714367 0.5714367
FP_Intercept_T_Rainbows FP_Intercept_T_Snow FP_Intercept_T_Storms FP_Intercept_T_Waterfalls FP_Intercept_T_Waves
0.5714367 0.5714367 0.5714367 0.5714367 0.5714367
RP1_var_bcons_1 RP1_cov_bcons_1_bcons_2 RP1_var_bcons_2 RP1_cov_bcons_1_bcons_3 RP1_cov_bcons_2_bcons_3
0.5714367 0.5714367 0.5714367 0.5714367 0.5714367
RP1_var_bcons_3 RP1_cov_bcons_1_bcons_4 RP1_cov_bcons_2_bcons_4 RP1_cov_bcons_3_bcons_4 RP1_var_bcons_4
0.5714367 0.5714367 0.5714367 0.5714367 0.5714367
RP1_cov_bcons_1_bcons_5 RP1_cov_bcons_2_bcons_5 RP1_cov_bcons_3_bcons_5 RP1_cov_bcons_4_bcons_5 RP1_var_bcons_5
0.5714367 0.5714367 0.5714367 0.5714367 0.5714367
RP1_cov_bcons_1_bcons_6 RP1_cov_bcons_2_bcons_6 RP1_cov_bcons_3_bcons_6 RP1_cov_bcons_4_bcons_6 RP1_cov_bcons_5_bcons_6
0.5714367 0.5714367 0.5714367 0.5714367 0.5714367
RP1_var_bcons_6 RP1_cov_bcons_1_bcons_7 RP1_cov_bcons_2_bcons_7 RP1_cov_bcons_3_bcons_7 RP1_cov_bcons_4_bcons_7
0.5714367 0.5714367 0.5714367 0.5714367 0.5714367
RP1_cov_bcons_5_bcons_7 RP1_cov_bcons_6_bcons_7 RP1_var_bcons_7 RP1_cov_bcons_1_bcons_8 RP1_cov_bcons_2_bcons_8
0.5714367 0.5714367 0.5714367 0.5714367 0.5714367
RP1_cov_bcons_3_bcons_8 RP1_cov_bcons_4_bcons_8 RP1_cov_bcons_5_bcons_8 RP1_cov_bcons_6_bcons_8 RP1_cov_bcons_7_bcons_8
0.5714367 0.5714367 0.5714367 0.5714367 0.5714367
RP1_var_bcons_8 RP1_cov_bcons_1_bcons_9 RP1_cov_bcons_2_bcons_9 RP1_cov_bcons_3_bcons_9 RP1_cov_bcons_4_bcons_9
0.5714367 0.5714367 0.5714367 0.5714367 0.5714367
RP1_cov_bcons_5_bcons_9 RP1_cov_bcons_6_bcons_9 RP1_cov_bcons_7_bcons_9 RP1_cov_bcons_8_bcons_9 RP1_var_bcons_9
0.5714367 0.5714367 0.5714367 0.5714367 0.5714367
RP1_cov_bcons_1_bcons_10 RP1_cov_bcons_2_bcons_10 RP1_cov_bcons_3_bcons_10 RP1_cov_bcons_4_bcons_10 RP1_cov_bcons_5_bcons_10
0.5714367 0.5714367 0.5714367 0.5714367 0.5714367
RP1_cov_bcons_6_bcons_10 RP1_cov_bcons_7_bcons_10 RP1_cov_bcons_8_bcons_10 RP1_cov_bcons_9_bcons_10 RP1_var_bcons_10
0.5714367 0.5714367 0.5714367 0.5714367 0.5714367