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 に設定することによって適合された同じモデルとの間の対数尤度の差の 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