首页 > 解决方案 > 为什么使用 R 中的生存包进行区间回归有时会给出可疑的标准误 0?

问题描述

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

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

library(survival)

data <- structure(list(lb = c(80000, 0, 0, 30000, 30000, 50000, 50000, 
20000, 0, 0, 30000, 30000, 50000, 50000, 20000, 50000, 50000, 
50000, 1e+05, 0, 60000, 130000, 60000, 160000), ub = c(1e+05, 
20000, 20000, 40000, 40000, 60000, 60000, 30000, 20000, 20000, 
40000, 40000, 60000, 60000, 30000, 60000, 60000, 60000, 130000, 
20000, 80000, 160000, 80000, 2e+05), group = c(0, 0, 0, 0, 1, 
1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 0, 1, 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,Surv(lb, ub, event=rep(3, nrow(data)), type="interval"))
summary(survreg(Y~group, data=data, dist="gaussian", robust=T))

Call:
survreg(formula = Y ~ group, data = data, dist = "gaussian", 
    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

任何想法将不胜感激。

标签: rsurvivalstandard-error

解决方案


我给包作者发了电子邮件寻求帮助。他不确定这种行为的确切原因,但确认存在某种错误。

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


推荐阅读