为什么使用 R 中的生存包进行区间回归有时会给出 0 的可疑标准误差？

如何解决为什么使用 R 中的生存包进行区间回归有时会给出 0 的可疑标准误差？

我正在 R 中使用区间数据的生存包运行区间回归模型。在标准误差为 0 且 p 值无穷大的情况下，我不断得到估计值，而这只是没有通过嗅探测试。我还将我的结果与 Stata 中相同数据的区间回归进行了比较，有些结果相同，有些则不同。

我尝试阅读代码以确定发生了什么，但我一直无法弄清楚。数据是否打破了某些假设？数值优化在 R 中做了什么奇怪的事情吗？这是一个错误吗？

library(survival)

data <- structure(list(lb = c(80000,30000,50000,20000,1e+05,60000,130000,160000),ub = c(1e+05,40000,80000,160000,2e+05),group = c(0,1,0)),row.names = c(NA,-24L),class = "data.frame")

所以这个，除了最后一个观察之外的所有内容，都可以正常工作。

Y <- with(data[1:23,],Surv(lb,ub,event=rep(3,nrow(data[1:23,])),type="interval"))
summary(survreg(Y~group,data=data[1:23,dist="gaussian",robust=T))

Call:
survreg(formula = Y ~ group,data = data[1:23,dist = "gaussian",robust = T)
               Value Std. Err (Naive SE)     z       p
(Intercept) 4.56e+04 9.28e+03   9.42e+03  4.92 8.7e-07
group       5.29e+03 1.36e+04   1.36e+04  0.39     0.7
Log(scale)  1.04e+01 1.87e-01   1.54e-01 55.40 < 2e-16

Scale= 32232 

Gaussian distribution
Loglik(model)= -52.4   Loglik(intercept only)= -52.4
    Chisq= 0.15 on 1 degrees of freedom,p= 0.7 
(Loglikelihood assumes independent observations)
Number of Newton-Raphson Iterations: 2 
n= 23

现在，添加最后的观察结果。

Y <- with(data,nrow(data)),data=data,data = data,robust = T)
                Value  Std. Err (Naive SE)     z       p
(Intercept) 55475.198  8221.817   8328.581  6.75 1.5e-11
group       -4364.673     0.000      0.000  -Inf < 2e-16
Log(scale)     10.609     0.202      0.150 52.65 < 2e-16

Scale= 40490 

Gaussian distribution
Loglik(model)= -59   Loglik(intercept only)= -59.1
    Chisq= 0.07 on 1 degrees of freedom,p= 0.79 
(Loglikelihood assumes independent observations)
Number of Newton-Raphson Iterations: 2 
n= 24 

在 24 个观测数据集中添加一个观测值会产生 0 的标准误差，这似乎很疯狂。所以我将结果与 Stata 中的区间回归进行了比较。第一个区间回归在 R 和 Stata 中产生相同的结果。然而，第二个在 Stata 中完全不同。

. intreg lb ub group

Fitting constant-only model:

Iteration 0:   log likelihood = -59.128989  
Iteration 1:   log likelihood = -59.054728  
Iteration 2:   log likelihood = -59.054631  
Iteration 3:   log likelihood = -59.054631  

Fitting full model:

Iteration 0:   log likelihood =  -59.15747  
Iteration 1:   log likelihood = -59.019687  
Iteration 2:   log likelihood = -59.019345  
Iteration 3:   log likelihood = -59.019345  

Interval regression                             Number of obs     =         24
                                                LR chi2(1)        =       0.07
Log likelihood = -59.019345                     Prob > chi2       =     0.7905

------------------------------------------------------------------------------
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       group |  -4441.439   16708.64    -0.27   0.790    -37189.77    28306.89
       _cons |   55510.53   11336.39     4.90   0.000     33291.61    77729.45
-------------+----------------------------------------------------------------
    /lnsigma |   10.60881   .1501972    70.63   0.000     10.31443    10.90319
-------------+----------------------------------------------------------------
       sigma |   40490.08   6081.497                       30164.8    54349.64
------------------------------------------------------------------------------
             0  left-censored observations
             0     uncensored observations
             0 right-censored observations
            24       interval observations

如有任何想法，我们将不胜感激。

解决方法

我通过电子邮件向软件包作者发送了电子邮件，以寻求这方面的帮助。他不确定这种行为的确切原因，但确认存在某种错误。

他提出了将所有值缩小到 0 到 200 之间的简单临时解决方案。当将两个边界除以 1000 并将系数解释为千位时，问题确实解决了。

,

当我跑步时：

 summary(survreg(Surv(lb,ub,event=rep(3,nrow(data[1:24,])),type="interval")~group,data=data,dist="gaussian",robust=T))

...我没有发现错误。这不是我第一次遇到通过不在全局环境中构造 Surv 对象来解决模糊错误情况的情况，而是将它们放在公式中以确保它们从传递给生存函数的数据对象中提取.然而，当我回去尝试复制错误时我没有成功，所以我们现在需要更多的细节来了解为什么你的体验与我的不同。这是我的 sessionInfo()：

Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.5 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/libopenblasp-r0.2.20.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8    LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] survival_3.2-11

loaded via a namespace (and not attached):
[1] compiler_4.0.4  Matrix_1.3-4    tools_4.0.4     splines_4.0.4   grid_4.0.4      lattice_0.20-44

我希望您不会使用更新的 R 和匹配的生存包运行新会话并重现您的错误。然后，您需要将该示例作为错误报告发送给 Terry Therneau。