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'

在尝试估计托比特模型之后。

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

*编辑：

通过删除一个协变量(我尝试了X3，i和X4，i，)并按照BondedDust的建议将下限审查设置为-0.001，它工作得很好，我甚至将重复次数推到了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")) )

}

请注意，上面我用等效的pmax替换了笨重且可能低效的sapply( ... max)。

> 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

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