Impossible d'ajuster la régression binomiale négative dans R (tentative de réplication des résultats publiés)

Tenter de reproduire les résultats de l'article récemment publié,

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

(Les données et le code stata sont disponibles sur http://www.aeaweb.org/aer/data/feb2013/20100973_data.zip ).

Je n'ai aucun problème à recréer les 5 premières régressions dans R (en utilisant les méthodes OLS et poisson), mais je suis tout simplement incapable de recréer leurs résultats de régression binomiale négative dans R, tandis qu'en stata, la régression fonctionne bien.

Plus précisément, voici le code R que j'ai écrit, qui ne parvient pas à exécuter une régression binomiale négative sur les données:

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))

En exécutant ce qui précède dans R, j'obtiens la sortie suivante:

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

J'ai essayé d'utiliser un certain nombre de valeurs initiales différentes pour thêta, ainsi que de faire varier le nombre maximal d'itérations sans succès. Le code stata fourni par les auteurs fonctionne bien, mais je n'arrive toujours pas à contraindre R à faire fonctionner le modèle. Existe-t-il des méthodes d'ajustement alternatives pour glm.nb () qui pourraient être plus robustes au problème que je rencontre?

Il existe de nombreuses publications sur le paramétrage stable des modèles non linéaires. Pour une raison quelconque, cela semble être largement ignoré dans R. Dans ce cas, la «matrice de conception» pour le prédicteur linéaire bénéficie de certains travaux. Laisser M$M$$M$ être la matrice de conception et $p$être les paramètres du modèle. Le prédicteur linéaire des moyennes$\mu$ est donné par

$$\mu=\exp(Mp)$$ Le reparameterization est accompli par le gramme-Schmidt modifié qui produit une matrice carrée $\Sigma$ tel que $$M=O\Sigma$$ où les colonnes de $O$sont orthonormés. (En fait, certaines des colonnes dans ce cas sont 0, de sorte que la méthode doit être légèrement modifiée pour y faire face.) Ensuite $$Mp=O\Sigma p$$ Laissez les nouveaux paramètres $q$ satisfaire $$p=\Sigma q$$ Pour que l'équation des moyens devienne $$\mu=\exp(Oq)$$ Cette paramétrisation est beaucoup plus stable et on peut adapter le modèle pour $q$ puis résoudre le $p$. J'ai utilisé cette technique pour ajuster le modèle avec AD Model Builder, mais cela pourrait fonctionner avec R. En tout état de cause, après avoir ajusté le modèle, on devrait regarder les "résidus" qui sont la différence au carré entre chaque observation et sa moyenne divisée par le estimation de la variance. Comme cela semble courant pour ce type de modèle, il existe d'énormes résidus. Je pense que ces points devraient être examinés avant que les résultats du document ne soient pris au sérieux.