无法在 R 中拟合负二项式回归(试图复制已发布的结果)

机器算法验证 r 回归 负二项分布
2022-03-22 05:32:18

试图复制最近发表的文章的结果,

Aghion、Philippe、John Van Reenen 和 Luigi Zingales。2013.“创新和制度所有权”。美国经济评论,103(1):277-304。

(数据和 stata 代码可在http://www.aeaweb.org/aer/data/feb2013/20100973_data.zip获得)。

在 R 中重新创建前 5 个回归(使用 OLS 和泊松方法)没有问题,但我根本无法在 R 中重新创建它们的负二项式回归结果,而在 stata 中回归工作正常。

具体来说,这是我编写的 R 代码,它无法对数据进行负二项式回归:

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

在 R 中运行上述内容,我得到以下输出:

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

已经尝试使用许多不同的 theta 初始值,以及在没有运气的情况下改变最大迭代次数。作者提供的 stata 代码工作正常,但我似乎仍然无法强迫 R 使模型工作。glm.nb() 是否有替代拟合方法可能对我遇到的问题更稳健?

1个回答

有很多关于非线性模型的稳定参数化的文献。由于某种原因,这似乎在 R 中很大程度上被忽略了。在这种情况下,线性预测器的“设计矩阵”受益于一些工作。M是设计矩阵和p是模型的参数。均值的线性预测器μ是(谁)给的

μ=exp(Mp)
重新参数化是通过修改后的 gram-Schmidt 完成的,它产生一个方阵Σ这样
M=OΣ
其中的列O是正交的。(实际上这种情况下的某些列是 0,因此必须稍微修改方法来处理这个问题。)然后
Mp=OΣp
让新参数q满足
p=Σq
因此,均值方程变为
μ=exp(Oq)
这种参数化更加稳定,可以拟合模型q 然后解决p. 我使用这种技术将模型与 AD 模型生成器一起拟合,但它可能与 R 一起使用。无论如何,拟合模型后,应该查看“残差”,即每个观察值与其均值之间的平方差除以估计方差。这种类型的模型似乎很常见,存在一些巨大的残差。我认为这些应该在论文的结果被认真对待之前进行检查。