0

我试图了解在 survreg 模型上运行残差函数产生的“ld”残差是什么？

例如

library(survival)
mod <- survreg(Surv(time, status -1) ~  age , data = lung)
residuals( mod , "ldcase")
residuals( mod, "ldshape")
residuals( mod , "ldresp")

残差函数的文档说明如下：

Escobar 和 Meeker 在一篇文章中讨论了基于这些数量的诊断。主要是对案例权重 (ldcase)、响应值 (ldresp) 和形状的扰动的似然位移残差。

参考

Escobar, LA 和 Meeker, WQ (1992)。使用审查数据评估回归分析中的影响。生物识别技术 48, 507-528。

以案例权重“ldcase”为例，我从参考论文中的理解是，这些残差代表了原始模型与通过将主体 i 的权重设置为 2 来拟合的相同模型之间的对数似然差异的两倍的估计。

但是，当我尝试自己手动编写代码时，我的派生值似乎与残差函数产生的值完全没有关系（下面是完全可重现的示例）：

library(survival)
library(ggplot2)

mod <- survreg( Surv(time, status -1) ~  age ,    data = lung)

get_ld <- function(i, mod){
    weight <- rep(1 , nrow(lung))
    weight[i] <- 2
    modw <- survreg( 
        Surv(time, status -1) ~  age  , 
        data = lung , 
        weights = weight
    )
    2 * as.numeric(logLik(mod) - logLik(modw))
}

dat <- data.frame(
    ld =  sapply( 1:nrow(lung), get_ld , mod = mod),
    ld_est = residuals(mod , "ldcase")
)

ggplot( data = dat , aes( x = ld_est , y = ld)) + geom_point()

在此处输入图像描述

此外，从论文来看，这些残差应该以 2 * chisq ( p + 2 ) 分布分布，在这种情况下，p = 1 给出 15.62 的单边 95% 截止点，这意味着我手动导出的残差是至少在正确的规模上，这让我对“ldcase”返回的残差实际上是什么感到非常困惑？

4

0 回答 0