Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

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.