I'm trying to replicate the results of the first model of this article:
Hultman, Lisa, Jacob Kathman, and Megan Shannon. 2013. “United Nations Peacekeeping and Civilian Protection in Civil War.” American Journal of Political Science 57(4): 875–91.
Replication material can be found here: http://thedata.harvard.edu/dvn/dv/ajps/faces/study/StudyPage.xhtml?studyId=87987&tab=files
The Stata code provided in the do file runs without problems. However, when I try to replicate the model in R. The error strongly resembles the one in this CV question. Here's the R code I've written to replicate the results:
library(MASS)
library(foreign)
pko <- read.dta("HKS_AJPS_2013.dta")
pko_model1 <- glm.nb(osvAll ~ troopLag +
policeLag +
militaryobserversLag +
brv_AllLag +
osvAllLagDum +
incomp +
epduration +
lntpop, data = pko, link = log)
This produces the following error:
Error in glm.fitter(x = X, y = Y, w = w, etastart = eta, offset = offset, :
NA/NaN/Inf in 'x'
If I include the control=glm.control(trace=10,maxit=100) option in glm.nb it produces the following output:
Deviance = 75787029 Iterations - 1
Deviance = 28247900 Iterations - 2
Deviance = 11043902 Iterations - 3
Deviance = 4952253 Iterations - 4
Deviance = 2896062 Iterations - 5
Deviance = 2286069 Iterations - 6
Deviance = 2152722 Iterations - 7
Deviance = 2135621 Iterations - 8
Deviance = 2134804 Iterations - 9
Deviance = 2134801 Iterations - 10
Deviance = 2134801 Iterations - 11
theta.ml: iter 0 'theta = 0.000609'
theta.ml: iter1 theta =0.00120256
theta.ml: iter2 theta =0.00234778
theta.ml: iter3 theta =0.00449211
theta.ml: iter4 theta =0.0082914
theta.ml: iter5 theta =0.0143798
theta.ml: iter6 theta =0.0225089
theta.ml: iter7 theta =0.0302781
theta.ml: iter8 theta =0.0342821
theta.ml: iter9 theta =0.0349533
theta.ml: iter10 theta =0.0349683
Initial value for 'theta': 0.034968
Deviance = 3634.951 Iterations - 1
Deviance = 1160161 Iterations - 2
Error in glm.fitter(x = X, y = Y, w = w, etastart = eta, offset = offset, :
NA/NaN/Inf in 'x'
If I exclude the epduration variable or the incomp variable, the error disappears and I can roughly replicate the results from the article, but the parameter estimates of course vary because I don't include all variables in the model (and I don't use robust, clustered standard errors in R).
Two questions:
- Why does this run in Stata without any complaint, but not in R?
- How can I make this work in R? The answers to this question suggest a possible solution by first estimating a Poisson regression and then feeding the results as starting parameter values into an MLE estimation. Yet, I haven't been able to make this work in R.
I realize that this might be a duplicate to the existing CrossValidated question but since I'm having a similar problem with different data, this might be a more general problem.
Best Answer
This answer did not seem to be too popular since it was not be done in R. For better or worse a lot of people seem to be restricted to R. I recall that I helped produce a package, glmmADMB based on AD Model builder, to do among other things negative binomial mixed models in R. I was curious to see if the software could handle the special case of zero random effects i.e not a mixed model. The results look pretty good so I think that one could use glmmADMB to fit negative binomial (or ZINB) models in R even when there are no random effects. Of source this does not change the fact that the model does not appear to fit the data very well. The R script is
It appears to get the solution without any problem. The output is
I also fit the alternative variance model nbinom1
GLMM's in R powered by AD Model Builder:
The AD Model Builder solution for the standard parameterization of the overdispersion is
For theta in glm.nb the estimate is 1/exp(log_b)=0.059604. So that explains how the r code got the answer. It was started at the correct value of theta. For ADMB I did not use any special value for the equivalent parameter.
Now glm.nb uses the paramaterization 1+mu/theta for the overdispersion. This is done mostly because it fits in with the glm approach, not because it is the only way to parameterize the overdispersion. Another way is to use tau for the overdispersion ie the variance is tau * mu rather than mu * (1+mu/theta). I also fit this model. The parameter estimates are
The objective function is the negative log-likelihood so that the log-likelihood difference is over 300 better.
The parameter estimates appear quite different. At this point I think one should investigate the scaled squared residuals. I think it is remarkable that this paper was accepted for publication without any discussion over the lack of fit when the model appears to be so inadequate. These are the worst fitting observations. Note that they depend somehwat on how the overdispersion is parameterized
For the R parameterization of the overdispersion the biggest residuals are
Feedback welcome.