No se puede ajustar la regresión binomial negativa en R (intentando replicar los resultados publicados)

Intentando replicar los resultados del artículo publicado recientemente,

Aghion, Philippe, John Van Reenen y Luigi Zingales. 2013. "Innovación y propiedad institucional". American Economic Review, 103 (1): 277-304.

(Los datos y el código de estado están disponibles en http://www.aeaweb.org/aer/data/feb2013/20100973_data.zip ).

No tengo problemas para recrear las primeras 5 regresiones en R (usando métodos OLS y Poisson), pero simplemente no puedo recrear sus resultados de regresión binomial negativa en R, mientras que en el estado la regresión funciona bien.

Específicamente, aquí está el código R que he escrito, que no puede ejecutar una regresión binomial negativa en los datos:

library(foreign)
library(MASS)
data.AVRZ <- read.dta("results_data2011.dta",
                  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,
   link = log,
   control=glm.control(trace=10,maxit=100))

Al ejecutar lo anterior en R, obtengo el siguiente resultado:

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'
In addition: Warning message:
step size truncated due to divergence

Intenté usar una cantidad de valores iniciales diferentes para theta, así como variar la cantidad máxima de iteraciones sin suerte. El código stata suministrado por los autores funciona bien, pero parece que todavía no puedo obligar a R a hacer que el modelo funcione. ¿Existen métodos de ajuste alternativos para glm.nb () que puedan ser más sólidos para el problema que estoy encontrando?

Existe mucha literatura sobre la parametrización estable de modelos no lineales. Por alguna razón, esto parece ignorarse en gran medida en R. En este caso, la "matriz de diseño" para el predictor lineal se beneficia de algún trabajo. Dejar METRO$M$$M$ ser la matriz de diseño y $p$Ser los parámetros del modelo. El predictor lineal de las medias.$\mu$ es dado por

$$\mu=\exp(Mp)$$ La reparametrización se realiza mediante gram-Schmidt modificado que produce una matriz cuadrada $\Sigma$ tal que $$M=O\Sigma$$ donde las columnas de $O$son ortonormales (De hecho, algunas de las columnas en este caso son 0, por lo que el método debe modificarse ligeramente para tratar esto). $$Mp=O\Sigma p$$ Deja que los nuevos parámetros $q$ satisfacer $$p=\Sigma q$$ Para que la ecuación para las medias se convierta $$\mu=\exp(Oq)$$ Esta parametrización es mucho más estable y se puede ajustar el modelo para $q$ y luego resolver el $p$. Utilicé esta técnica para ajustar el modelo con AD Model Builder, pero podría funcionar con R. En cualquier caso, después de ajustar el modelo, uno debe mirar los "residuos", que son la diferencia al cuadrado entre cada observación y su media dividida por el Estimación de la varianza. Como parece ser común para este tipo de modelo, hay algunos residuos enormes. Creo que deberían examinarse antes de que se tomen en serio los resultados del trabajo.