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.
3 comments:
Hmm that's quiet interessting but frankly i have a hard time visualizing it... wonder what others have to say..
Approvingly your article helped me terribly much in my college assignment. Hats high to you enter, will look ahead for more interdependent articles promptly as its sole of my favourite issue to read.
I was getting your "Hessian non-positive definite. No variance!" warning as well, and puzzled by what was going on.
I poked around the source code (here: https://github.com/cran/glmmML/blob/master/src/glmmml.c) and it seems like that warning is _always_ produced when `fix.sigma==FALSE`.
It went away when I set `fix.sigma==TRUE`.
Not sure why on earth the option's there if it's guaranteed to give a warning like that...
Post a Comment