Generalized linear mixed effect model problem

February 16, 2010
By

(This article was first published on Shige's Research Blog, and kindly contributed to R-bloggers)

I am trying to compare cohort difference in infant mortality using generalized linear mixed model. I first estimated the model in Stata:

xi:xtlogit inftmort i.cohort, i(code)

which converged nicely:


Fitting comparison model:


Iteration 0:   log likelihood = -1754.4476
Iteration 1:   log likelihood = -1749.3366
Iteration 2:   log likelihood = -1749.2491
Iteration 3:   log likelihood = -1749.2491


Fitting full model:


tau =  0.0     log likelihood = -1749.2491
tau =  0.1     log likelihood = -1743.8418
tau =  0.2     log likelihood = -1739.0769
tau =  0.3     log likelihood = -1736.4914
tau =  0.4     log likelihood = -1739.5415


Iteration 0:   log likelihood = -1736.4914  
Iteration 1:   log likelihood = -1722.6629  
Iteration 2:   log likelihood = -1694.9114  
Iteration 3:   log likelihood = -1694.6509  
Iteration 4:   log likelihood =  -1694.649  
Iteration 5:   log likelihood =  -1694.649  


Random-effects logistic regression              Number of obs      =     21694
Group variable: code                            Number of groups   =     10789


Random effects u_i ~ Gaussian                   Obs per group: min =         1
                                                               avg =       2.0
                                                               max =         9


                                                Wald chi2(2)       =      8.05
Log likelihood  =  -1694.649                    Prob > chi2        =    0.0178


------------------------------------------------------------------------------
    inftmort |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
  _Icohort_2 |  -.5246846   .1850328    -2.84   0.005    -.8873422   -.1620269
  _Icohort_3 |  -.1424331    .140369    -1.01   0.310    -.4175513     .132685
       _cons |  -5.214642   .1839703   -28.35   0.000    -5.575217   -4.854067
-------------+----------------------------------------------------------------
    /lnsig2u |   .9232684   .1416214                      .6456956    1.200841
-------------+----------------------------------------------------------------
     sigma_u |   1.586665   .1123528                      1.381055    1.822885
         rho |   .4335015   .0347791                      .3669899    .5024984
------------------------------------------------------------------------------
Likelihood-ratio test of rho=0: chibar2(01) =   109.20 Prob >= chibar2 = 0.000

Then I tried the same model using lme4:

m2 <- glmer(inftmort ~ (1|code) + as.factor(cohort), family=binomial, data=d)

I got:


Warning message:
In mer_finalize(ans) : false convergence (8)

And the results are quite different from what I got from Stata. I tried to estimate the model using "glmmML" and also got into trouble. This time the error message is:


Warning message:
In glmmML.fit(X, Y, weights, cluster.weights, start.coef, start.sigma,  :
  Hessian non-positive definite. No variance!

I then tried MCMCglmm:


prior=list(R=list(V=1, nu=0, fix=1), G=list(G1=list(V=1, nu=0)))
m2 <- MCMCglmm(inftmort ~ 1 + as.factor(cohort), random=~code, family = "categorical", data=d, prior=prior)



It seems to be working and produced estimates that are comparable to what Stata produced (not identical, of course):


Iterations = 3001:12991
Thinning interval = 10 
Number of chains = 1 
Sample size per chain = 1000 


1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:


                      Mean     SD Naive SE Time-series SE
(Intercept)        -5.7145 0.1805 0.005708        0.02295
as.factor(cohort)2 -0.5633 0.1788 0.005653        0.02194
as.factor(cohort)3 -0.1888 0.1471 0.004653        0.01912


2. Quantiles for each variable:


                      2.5%     25%     50%      75%   97.5%
(Intercept)        -6.0464 -5.8324 -5.7251 -5.61120 -5.2817
as.factor(cohort)2 -0.8977 -0.6985 -0.5551 -0.43084 -0.2371
as.factor(cohort)3 -0.4614 -0.2947 -0.1867 -0.09644  0.1232 

This is puzzling.

To leave a comment for the author, please follow the link and comment on his blog: Shige's Research Blog.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Tags: , , , ,

Comments are closed.