# Unable to fit negative binomial regression in R (attempting to replicate published results)

Attempting to replicate the results from the recently published article,

Aghion, Philippe, John Van Reenen, and Luigi Zingales. 2013. "Innovation and Institutional Ownership." American Economic Review, 103(1): 277-304.

(Data and stata code is available at http://www.aeaweb.org/aer/data/feb2013/20100973_data.zip).

Am having no problem recreating the first 5 regressions in R (using OLS and poisson methods), but am simply unable to recreate their negative binomial regression results in R, while in stata the regression works fine.

Specifically, here's the R code I've written, which fails to run a negative binomial regression on the data:

library(foreign)
library(MASS)
convert.underscore=TRUE)
sicDummies <- grep("Isic4", names(data.AVRZ), value=TRUE)
yearDummies <- grep("Iyear", names(data.AVRZ), value=TRUE)
data.column.6 <- subset(data.AVRZ, select = c("cites",
"instit.percown",
"lk.l",
"lsal",
sicDummies,
yearDummies))
data.column.6 <- na.omit(data.column.6)

glm.nb(cites ~ .,
data = data.column.6,
control=glm.control(trace=10,maxit=100))


Running the above in R, I get the following output:

Initial fit:
Deviance = 1137144 Iterations - 1
Deviance = 775272.3 Iterations - 2
Deviance = 725150.7 Iterations - 3
Deviance = 722911.3 Iterations - 4
Deviance = 722883.9 Iterations - 5
Deviance = 722883.3 Iterations - 6
Deviance = 722883.3 Iterations - 7
theta.ml: iter 0 'theta = 0.000040'
theta.ml: iter1 theta =7.99248e-05
Initial value for 'theta': 0.000080
Deviance = 24931694 Iterations - 1
Deviance = NaN Iterations - 2
Step halved: new deviance = 491946.5
Error in glm.fitter(x = X, y = Y, w = w, etastart = eta, offset = offset,  :
NA/NaN/Inf in 'x'
step size truncated due to divergence


Have tried using a number of different initial values for theta, as well as varying the maximum number of iterations with no luck. The authors' supplied stata code works fine, but I still can't seem to coerce R into making the model work. Are there alternative fitting methods for glm.nb() that may be more robust to the problem I'm encountering?

• It will converge for fewer variables but not all of them - there are a lot of variables and a very ugly outcome variable. Perhaps try emailing r-help@r-project.org for help with this if no replies in a day or two. Mar 17, 2013 at 23:35
• Right, I've had no problems getting convergence with the first two covariates (instit.percown, lk.l), but when I add the sales variable, lnsal, convergence fails (not to mention the >100 SIC and state level dummies). Thanks for the tip on the email, curious though why stata has little problem with this model but I just can't get it to work in R.
– jayb
Mar 18, 2013 at 0:32
• Eventually was able to estimate this by running a poisson regression to get starting parameter values, then feeding these into a MLE on the log likelihood function. Will post the solution to this soon.
– jayb
Apr 1, 2013 at 1:07

There is a lot of literature about the stable parameterization of nonlinear models. For some reason this seems to be largely ignored in R. In this case the "design matrix" for the linear predictor benefits from some work. Let $M$ be the design matrix and $p$ be the parameters of the model. The linear predictor for the means $\mu$ is given by $$\mu=\exp(Mp)$$ The reparameterization is accomplished by modified gram-Schmidt which produces a square matrix $\Sigma$ such that $$M=O\Sigma$$ where the columns of $O$ are orthonormal. (In fact some of the columns in this case are 0 so that the method must be modified slightly to deal with this.) Then $$Mp=O\Sigma p$$ Let the new parameters $q$ satisfy $$p=\Sigma q$$ So that the equation for the means becomes $$\mu=\exp(Oq)$$ This parametrization is much more stable and one can fit the model for $q$ and then solve for the $p$. I used this technique to fit the model with AD Model Builder, but it might work with R. In any event, having fit the model one should look at the "residuals" which are the squared difference between each observation and its mean divided by the estimate for the variance. As seems to be common for this type of model there are some huge residuals. I think these should be examined before the results of the paper are taken seriously.