72

Rでは、predict.lmは線形回帰の結果に基づいて予測を計算し、これらの予測の信頼区間を計算することもできます。マニュアルによると、これらの間隔はフィッティングの誤差分散に基づいていますが、係数の誤差間隔には基づいていません。

一方、ロジスティック回帰とポアソン回帰に基づいて予測を計算するpredict.glmには、信頼区間のオプションがありません。そして、ポアソンとロジスティック回帰に意味のある洞察を提供するために、そのような信頼区間をどのように計算できるかを想像するのにも苦労しています。

そのような予測に信頼区間を提供することが意味のある場合はありますか？それらはどのように解釈できますか？そして、これらの場合の仮定は何ですか？

4

2 に答える 2

88

通常の方法は、線形予測子のスケールで信頼区間を計算することです。ここで、物事はより正規（ガウス）になり、リンク関数の逆関数を適用して、線形予測子スケールから応答スケールに信頼区間をマッピングします。

これを行うには、2つのことが必要です。

  1. で呼び出すpredict()、type = "link"および
  2. predict()で呼び出しますse.fit = TRUE。

1つ目は線形予測子のスケールで予測を生成し、2つ目は予測の標準誤差を返します。擬似コードで

## foo <- mtcars[,c("mpg","vs")]; names(foo) <- c("x","y") ## Working example data
mod <- glm(y ~ x, data = foo, family = binomial)
preddata <- with(foo, data.frame(x = seq(min(x), max(x), length = 100)))
preds <- predict(mod, newdata = preddata, type = "link", se.fit = TRUE)

predsfit次に、コンポーネントとのリストse.fitです。

線形予測子の信頼区間は次のようになります

critval <- 1.96 ## approx 95% CI
upr <- preds$fit + (critval * preds$se.fit)
lwr <- preds$fit - (critval * preds$se.fit)
fit <- preds$fit

critval必要に応じてtまたはz（正規）分布から選択され（どのタイプのGLMにどちらを使用するか、およびプロパティは何であるかを正確に忘れています）、必要なカバレッジがあります。これ1.96は、95％のカバレッジを与えるガウス分布の値です。

> qnorm(0.975) ## 0.975 as this is upper tail, 2.5% also in lower tail
[1] 1.959964

ここでfit、について、リンク関数の逆関数をそれらに適用する必要がありますupr。lwr

fit2 <- mod$family$linkinv(fit)
upr2 <- mod$family$linkinv(upr)
lwr2 <- mod$family$linkinv(lwr)

これで、3つすべてとデータをプロットできます。

preddata$lwr <- lwr2 
preddata$upr <- upr2 
ggplot(data=foo, mapping=aes(x=x,y=y)) + geom_point() +         
   stat_smooth(method="glm", method.args=list(family=binomial)) + 
   geom_line(data=preddata, mapping=aes(x=x, y=upr), col="red") + 
   geom_line(data=preddata, mapping=aes(x=x, y=lwr), col="red")

ここに画像の説明を入力してください

于 2013-01-20T12:15:07.413 に答える
2

ブートストラップまたはシミュレーションアプローチを使用してポアソン推定の問題を解決するLiuWenSuiの方法に出くわしました。

著者からの例

pkgs <- c('doParallel', 'foreach')
lapply(pkgs, require, character.only = T)
registerDoParallel(cores = 4)
 
data(AutoCollision, package = "insuranceData")
df <- rbind(AutoCollision, AutoCollision)
mdl <- glm(Claim_Count ~ Age + Vehicle_Use, data = df, family = poisson(link = "log"))
new_fake <- df[1:5, 1:2]

boot_pi <- function(model, pdata, n, p) {
  odata <- model$data
  lp <- (1 - p) / 2
  up <- 1 - lp
  set.seed(2016)
  seeds <- round(runif(n, 1, 1000), 0)
  boot_y <- foreach(i = 1:n, .combine = rbind) %dopar% {
    set.seed(seeds[i])
    bdata <- odata[sample(seq(nrow(odata)), size = nrow(odata), replace = TRUE), ]
    bpred <- predict(update(model, data = bdata), type = "response", newdata = pdata)
    rpois(length(bpred), lambda = bpred)
  }
  boot_ci <- t(apply(boot_y, 2, quantile, c(lp, up)))
  return(data.frame(pred = predict(model, newdata = pdata, type = "response"), lower = boot_ci[, 1], upper = boot_ci[, 2]))
}
 
boot_pi(mdl, new_fake, 1000, 0.95)

sim_pi <- function(model, pdata, n, p) {
  odata <- model$data
  yhat <- predict(model, type = "response")
  lp <- (1 - p) / 2
  up <- 1 - lp
  set.seed(2016)
  seeds <- round(runif(n, 1, 1000), 0)
  sim_y <- foreach(i = 1:n, .combine = rbind) %dopar% {
    set.seed(seeds[i])
    sim_y <- rpois(length(yhat), lambda = yhat)
    sdata <- data.frame(y = sim_y, odata[names(model$x)])
    refit <- glm(y ~ ., data = sdata, family = poisson)
    bpred <- predict(refit, type = "response", newdata = pdata)
    rpois(length(bpred),lambda = bpred)
  }
  sim_ci <- t(apply(sim_y, 2, quantile, c(lp, up)))
  return(data.frame(pred = predict(model, newdata = pdata, type = "response"), lower = sim_ci[, 1], upper = sim_ci[, 2]))
}
 
sim_pi(mdl, new_fake, 1000, 0.95)
于 2021-04-19T16:19:09.667 に答える