Generalized linear mixed effect model problem

February 16, 2010

(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.