Next Upcoming Google+ Hangout: Tuesday, August 27 @ 7PM (CST) - To Participate CLICK HERE

Search For Topics/Content

Missing Data/Imputation Discussion > poisson multilevel model after multiple imputation

Hi
I have carried out a multiple imputation in R using Amelia II and have created 5 multiply imputed datasets. The purpose of my research is to fit a Poisson Multilevel Model to the data to estimate numbers of hospital admissions.
Now that I have 5 completed datasets, how do I go about fitting the multilevel model? How do I pool all the 5 datasets to get one final multilevel/mixed model?

Can anyone PLEASE guide me?!

Many thanks and regards...

September 27, 2012 | Unregistered CommenterSwarna

Hi Swarna! First, congrats on coming as far as you have with R and Amelia (and for taking-on an ambitious analysis)! To complete the analysis, I think I'd use the Zelig package, which I believe does poisson and I know integrates well with MI datasets. Here is the link to the Zelig page:

http://cran.r-project.org/web/packages/Zelig/index.html

October 18, 2012 | Registered CommenterJeremy Taylor

Thank you Jeremy...so i have been tirelessly working with zelig and unfortunately i am getting real stuck! here is my R code:

###these are my 5 imputed datasets####

d.1 <- read.table("imp1.csv", header=TRUE,sep=",")
d.2 <- read.table("imp2.csv", header=TRUE,sep=",")
d.3 <- read.table("imp3.csv", header=TRUE,sep=",")
d.4 <- read.table("imp4.csv", header=TRUE,sep=",")
d.5 <- read.table("imp5.csv", header=TRUE,sep=",")

#####my zelig code#####
model3=zelig(DKA.counts~age+as.factor(sex)+bmi_value+hba1c_value,model="poisson",data=mi(d.1,d.2,d.3,d.4,d.5))
summary(model3)


Till so far it is fine... but after this I get separate output for each dataset! I have tried and tried to find the problem but i still cannot get the combined output.
Can anyone please tell me what i'm doing wrong!?!

#####output####
$d.1

Call:
glm(formula = formula, family = poisson(), data = data, weights = weights,
model = F)

Deviance Residuals:
Min 1Q Median 3Q Max
-1.9410 -0.6990 -0.4136 0.1423 10.3703

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.638545 0.097690 -6.536 6.3e-11 ***
age 0.033006 0.004625 7.136 9.6e-13 ***
as.factor(sex)1.5 0.672653 0.289634 2.322 0.0202 *
as.factor(sex)2 0.131818 0.027845 4.734 2.2e-06 ***
bmi_value -0.007784 0.003794 -2.051 0.0402 *
hba1c_value 0.087617 0.006812 12.862 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 3656.9 on 2943 degrees of freedom
Residual deviance: 3299.6 on 2938 degrees of freedom
AIC: 10307

Number of Fisher Scoring iterations: 5


$d.2
:
:
$d.3
:
:
$d.4
:
:
$d.5
Call:
glm(formula = formula, family = poisson(), data = data, weights = weights,
model = F)

Deviance Residuals:
Min 1Q Median 3Q Max
-1.8925 -0.6973 -0.4194 0.1321 10.3741

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.696957 0.097539 -7.145 8.97e-13 ***
age 0.029981 0.004636 6.467 1.00e-10 ***
as.factor(sex)1.5 0.671931 0.289643 2.320 0.0203 *
as.factor(sex)2 0.122598 0.027794 4.411 1.03e-05 ***
bmi_value -0.002243 0.003719 -0.603 0.5465
hba1c_value 0.086084 0.006774 12.709 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 3656.9 on 2943 degrees of freedom
Residual deviance: 3308.9 on 2938 degrees of freedom
AIC: 10316

Number of Fisher Scoring iterations: 5

October 26, 2012 | Unregistered CommenterSwarna

Dear Swarna,

I think I see the problem:

You need to create a single imputation object from our 5 imputed datasets. For example:

d.1 <- read.table("imp1.csv", header=TRUE,sep=",")
d.2 <- read.table("imp2.csv", header=TRUE,sep=",")
d.3 <- read.table("imp3.csv", header=TRUE,sep=",")
d.4 <- read.table("imp4.csv", header=TRUE,sep=",")
d.5 <- read.table("imp5.csv", header=TRUE,sep=",")

IMP.data<- mi(d.1,d.2,d.3,d.4,d.5)

#####my zelig code#####
model3<-zelig(DKA.counts~age+as.factor(sex)+bmi_value+hba1c_value,data=IMP.data, model="poisson")

print(summary(model3))


Let me know if that works for you...

November 3, 2012 | Registered CommenterJeremy Taylor