r - 为什么使用 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
任何想法将不胜感激。
解决方案
我给包作者发了电子邮件寻求帮助。他不确定这种行为的确切原因,但确认存在某种错误。
他提出了一个简单的临时解决方案,将所有值缩小到 0 到 200 之间。当将两个边界除以 1000 并将系数解释为数千时,这个问题确实得到了解决。
推荐阅读
- api - Mule Anypoint Studio - 测试连接失败 - 进行连接测试时出现未知错误 - 自定义 API
- python - 如何从 Python 中的偏态正态分布中抽取随机样本?
- android-studio - 无法在 android studio 3.6.1 中构建 APK 获取 app:compileReleaseAidl 错误
- javascript - 我想用 JavaScript 交换 2 个相同的 HTML 元素
- laravel - 在 laravel 5.7 中调用数组上的成员函数 links()
- c++ - 使用 set/multiset 查找以 's' 开头的单词
- python - 从某个类中获取所有文本,仅当它是另一个类的孩子时
- r - 如何将点与R中散点图中的一条线连接起来
- python - RuntimeError:无法在需要 grad 的变量上调用 numpy()。改用 var.detach().numpy()
- javascript - 一个处理完成后打开多个jQuery确认框