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.