Tobit：lm.fit 中的错误（X.vlm，y = z.vlm，...）：'x' 中的 NA/NaN/Inf（模拟数据）

Question

我正在进行一项蒙特卡洛研究。我有一个具有异方差性的线性模型，因变量左删失为 0。删失率的平均值为 25.9。

我得到了错误

Error in lm.fit(X.vlm, y = z.vlm, ...) : NA/NaN/Inf in 'x'

在尝试估计一个 tobit 模型之后。

vglm(y[i,]~X[1,i,]+X[2,i,]+X[3,i,]+X[4,i,],family=tobit(Lower=0))   

我的数据是根据标准分布模拟的，所以问题不应该来自奇数变量。

我发现了另外两个与真实数据有相同问题的问题： lm() NA/NaN/Inf error , lm() NA/NaN/Inf error 但似乎没有任何令人满意的答案。此外，我的数据很容易重现，因此它应该有助于识别问题

代码如下：

library(VGAM)
set.seed(12345)
nobs=100
nsim=100
b=c(2,-2,-3,3)
g=c(1,0.2)
y=matrix(rep(0,nobs*nsim),ncol=nobs,nrow=nsim)
X=array(0,dim=c(4,nsim,nobs))
res=matrix(rep(0,nobs*nsim),ncol=nobs,nrow=nsim)
tobit=vector(mode="list",length=nsim) 
for(i in 1:nsim){
# generate covariates :
X[1,i,]=rlnorm(n=nobs)           
X[2,i,]=runif(n=nobs)<=.75   
X[3,i,]=rnorm(mean = 3,n=nobs)      
X[4,i,]=runif(n=nobs,min=0,max=10)  
res[i,]=(g[1]+g[2]*X[4,i,])*rnorm(n=nobs)   
# generate censored dependent variable
y[i,]=b[1]*X[1,i,]+b[2]*X[2,i,]+b[3]*X[3,i,]+b[4]*X[4,i,]+res[i,]
y[i,]=sapply(y[i,],FUN=function(x){max(0,x)}) #apply censoring
tobit[[i]]<-vglm(y[i,]~X[1,i,]+X[2,i,]+X[3,i,]+X[4,i,],
                              family = tobit(Lower=0)) 
}

这是回溯

traceback()
5: lm.fit(X.vlm, y = z.vlm, ...)
4: vlm.wfit(xmat = X.vlm.save, z, Hlist = NULL, U = U, matrix.out =FALSE, 
   is.vlmX = TRUE, qr = qr.arg, xij = NULL)
3: vglm.fitter(x = x, y = y, w = w, offset = offset, Xm2 = Xm2, 
   Ym2 = Ym2, etastart = etastart, mustart = mustart, coefstart =coefstart, 
   family = family, control = control, constraints = constraints, 
   criterion = control$criterion, extra = extra, qr.arg = qr.arg, 
   Terms = mt, function.name = function.name, ...)
2: vglm(y[1, ] ~ X[1, 1, ] + X[2, i, ] + X[3, i, ] + X[4, i, ], 
   family = tobit(Lower = 0))
1: traceback(vglm(y[1, ] ~ X[1, 1, ] + X[2, i, ] + X[3, i, ] + X[4, 
   i, ], family = tobit(Lower = 0)))

*** 编辑：

通过删除一个协变量（我尝试使用 X[3,i,] 和 X[4,i,]）并将较低的审查设置为 -0.001，正如 BondedDust 建议的那样，它工作正常，我什至推动了复制次数到 1000 没有大问题。

通过将较低的删失设置为 -0.001，并保留所有协变量，我在 100 次迭代中得到了两个错误。值得注意的是，现在的错误是

Error in lm.fit(X.vlm, y = z.vlm, ...) : NA/NaN/Inf in 'y'

除了我收到这些警告

In vglm.fitter(x = x, y = y, w = w, offset = offset, Xm2 = Xm2,  ... :
iterations terminated because half-step sizes are very small

Answer 1

我注意到这在 i=1 处重复失败，因此我认为 vglm 调用本身可能存在问题。查看?tobit 中的示例，我添加了一些与删失分布相关的参数，并开始进行一些额外的迭代。然后我尝试缩小审查范围，并且只有 10% 的失败率获得了更大的成功。所以我最后添加了一个 try() 包装器，让循环在不停止计算的情况下进行迭代，并获得了大部分成功运行：

for(i in 1:nsim){

X[1,i,]=rlnorm(n=nobs)           
X[2,i,]=runif(n=nobs)<=.75   
X[3,i,]=rnorm(mean = 3,n=nobs)      
X[4,i,]=runif(n=nobs,min=0,max=10)  
res[i,]=(g[1]+g[2]*X[4,i,])*rnorm(n=nobs)   

y[i,]=b[1]*X[1,i,]+b[2]*X[2,i,]+b[3]*X[3,i,]+b[4]*X[4,i,]+res[i,]
y[i,]=pmax(0,y[i,])
tobit[[i]]<-try( vglm(y[i,]~X[1,i,]+X[2,i,]+X[3,i,]+X[4,i,], crit = "coeff", 
                              family = tobit(Lower=-.001, Upper=30, type.f = "cens")) )

}

请注意，我将您笨拙且可能效率低下的sapply( ... max) 替换为等效的pmax。

> table( sapply(tobit, class))

try-error      vglm 
       12        88 

您可以通过以下方式循环成功返回：

sapply( tobit[ sapply(tobit, class) == "vglm"],  coefficients)

结果顶部：

                    [,1]      [,2]      [,3]       [,4]       [,5]       [,6]
(Intercept):1  2.8460081  1.910137  1.672237  1.2888827  2.4970536  1.0006290
(Intercept):2  0.9183935  1.042424  1.094658  0.9767228  0.9263946  0.9250609
X[1, i, ]      1.7777788  1.880506  1.662835  1.6204394  1.4412304  1.6275208
X[2, i, ]     -3.0847792 -0.453110 -1.152709 -0.9900163 -2.4705355 -0.9651577
X[3, i, ]     -2.4272169 -2.094114 -2.314748 -2.4628501 -1.9001385 -2.1076416
X[4, i, ]      2.6225234  2.245107  2.460182  2.7027493  2.3653673  2.3841989
                    [,7]       [,8]      [,9]      [,10]     [,11]      [,12]
(Intercept):1  0.9520376  1.6319010  1.572563  1.4709517  1.616158  2.4992492
(Intercept):2  0.8698777  0.9005506  1.147485  0.9285724  1.012186  0.9229233
X[1, i, ]      1.6483879  1.6789573  1.718641  1.6544123  1.599116  1.7204001
X[2, i, ]     -0.3718720 -1.8690782 -2.408657 -1.7278915 -1.208939 -2.0037999
X[3, i, ]     -2.2601637 -1.9118288 -2.359274 -1.7828438 -2.257556 -2.3778443
X[4, i, ]      2.5381367  2.3091630  2.583869  2.3582418  2.333988  2.4389336

在获得这种适度的成功后，我尝试将 Lower 设置回 0 并得到所有错误。增加上限值似乎不会影响有限测试的成功率。我无法解释这些发现，但也许可以咨询软件包作者。