Negative binomiale Regression in R kann nicht angepasst werden (Versuch, veröffentlichte Ergebnisse zu replizieren)

Der Versuch, die Ergebnisse des kürzlich veröffentlichten Artikels zu replizieren,

Aghion, Philippe, John Van Reenen und Luigi Zingales. 2013. "Innovation und institutionelles Eigentum." American Economic Review, 103 (1): 277-304.

(Daten- und Datencode finden Sie unter http://www.aeaweb.org/aer/data/feb2013/20100973_data.zip ).

Ich habe kein Problem damit, die ersten 5 Regressionen in R neu zu erstellen (unter Verwendung von OLS- und Poisson-Methoden), kann aber einfach ihre negativen binomialen Regressionsergebnisse in R nicht wiederherstellen, während in Regata die Regression gut funktioniert.

Hier ist der von mir geschriebene R-Code, der keine negative binomiale Regression für die Daten ausführt:

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

Wenn ich das Obige in R ausführe, erhalte ich die folgende Ausgabe:

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

Ich habe versucht, eine Reihe verschiedener Anfangswerte für Theta zu verwenden und die maximale Anzahl von Iterationen ohne Glück zu variieren. Der von den Autoren bereitgestellte Stata-Code funktioniert einwandfrei, aber ich kann R immer noch nicht dazu zwingen, das Modell zum Laufen zu bringen. Gibt es alternative Anpassungsmethoden für glm.nb (), die möglicherweise robuster für das Problem sind, auf das ich stoße?

Es gibt viel Literatur über die stabile Parametrisierung nichtlinearer Modelle. Aus irgendeinem Grund scheint dies in R weitgehend ignoriert zu werden. In diesem Fall profitiert die "Entwurfsmatrix" für den linearen Prädiktor von einigen Arbeiten. Sei die Entwurfsmatrix und die Parameter des Modells. Der lineare Prädiktor für das Mittel ist gegeben durch Die Neuparametrisierung erfolgt durch modifiziertes Gramm-Schmidt, das eine quadratische Matrix so dass wobei die Spalten von orthonormal sind. (Tatsächlich sind einige der Spalten in diesem Fall 0, so dass die Methode geringfügig geändert werden muss, um dies zu bewältigen.) Dann ist M$M$$M$$p$$\mu$

$$\mu=\exp(Mp)$$ $\Sigma$ $$M=O\Sigma$$ $O$ $$Mp=O\Sigma p$$ Lassen Sie die neuen Parameter erfüllen, so dass die Gleichung für die Mittelwerte Diese Parametrisierung ist viel stabiler und man kann das Modell für und dann nach . Ich habe diese Technik verwendet, um das Modell mit AD Model Builder anzupassen, aber es könnte mit R funktionieren. In jedem Fall sollte man nach dem Anpassen des Modells die "Residuen" betrachten, die die quadratische Differenz zwischen jeder Beobachtung und ihrem Mittelwert geteilt durch die sind Schätzung für die Varianz. Wie es für diesen Modelltyp üblich zu sein scheint, gibt es einige große Residuen. Ich denke, diese sollten geprüft werden, bevor die Ergebnisse des Papiers ernst genommen werden.$q$ $$p=\Sigma q$$ $$\mu=\exp(Oq)$$ $q$$p$