负二项回归的预测区间(R)

new9mtju 于 2023-09-27 发布在其他

关注(0)|答案(1)|浏览(107)

我试图从负二项回归模型中估计 * 预测 * 区间（不是置信区间）。
我知道我不能简单地使用predict()，因为predict.glm()不像predict.lm()，不允许指定interval = "prediction"-所以它会返回均值的置信区间，而不是预测区间。

head(cars)
#>   speed dist
#> 1     4    2
#> 2     4   10
#> 3     7    4
#> 4     7   22
#> 5     8   16
#> 6     9   10

我已经看到你可以使用ciTools::add_pi()来实现这个目的。这在基本场景中工作正常：

fit <- cars |>
    MASS::glm.nb(dist ~ speed, data = _)

cars |>
    ciTools::add_pi(fit) |>
    head()
#>   speed dist     pred LPB0.025 UPB0.975
#> 1     4    2 11.42897        3       23
#> 2     4   10 11.42897        3       24
#> 3     7    4 15.63866        5       30
#> 4     7   22 15.63866        5       31
#> 5     8   16 17.36189        6       34
#> 6     9   10 19.27501        7       37

但是，我也希望预测一组新的数据（例如：speed= c(58, 63, 121)）。当我尝试这样做时，我得到一个长度不匹配错误。如果我在线性模型上使用predict()，就不会发生这种情况。

cars |>
    tibble::add_row(speed = c(58, 63, 121)) |>
    ciTools::add_pi(fit)
#> Warning in mu * rgamma(k, theta): longer object length is not a multiple of
#> shorter object length
#> Error in rpois(k, (mu * rgamma(k, theta))/theta): dims [product 50] do not match the length of object [53]

有什么想法吗？
(EDIT：在创建一组新数据时，我最初错误地写了“dist”，而我的意思是“速度”-感谢Dave在答案中指出这一点）。

r

来源：https://stackoverflow.com/questions/77039248/prediction-intervals-for-negative-binomial-regression-r

1条答案

按热度按时间

gfttwv5a1#

生成置信区间的函数使用所提供的数据来创建模型矩阵，该模型矩阵以列表方式删除输入数据-从预测中删除假设行。
您可以找到实际执行工作的函数的“核心”，并将该函数弯曲以满足您自己的需要。当你调用add_pi()时，它会调用add_pi.negbin()，add_pi.negbin()又会调用sim_pi_negbin()，而sim_pi_negbin()又会调用函数get_sim_response_nb()来完成繁重的工作。此参数接受数据、模型对象和模拟次数（作为nSims）。以下是ciTools:::get_sim_response_nb()的外观：

function (df, fit, nSims) 
{
    nPreds <- NROW(df)
    modmat <- model.matrix(fit, data = df)
    response_distr <- fit$family$family
    inverselink <- fit$family$linkinv
    sims <- get_negbin_sims(fit, nSims)
    sim_response <- matrix(0, ncol = nSims, nrow = nPreds)
    for (i in 1:nSims) {
        yhat <- inverselink(modmat %*% sims@coef[i, ])
        sim_response[, i] <- rnegbin(n = nPreds, mu = yhat, theta = fit$theta)
    }
    sim_response
}

尝试使用此函数时，重要的是您需要能够提供模型矩阵，而不是使用自动计算的矩阵。然后，您可以更改df参数，使其期望模型矩阵。您可能还希望该函数不生成预测区间。您可以添加一个alpha参数来标识要设置间隔的级别。

fit <- cars |>
  MASS::glm.nb(dist ~ speed, data = _)

my_sim_response_nb <- function (modmat, fit, nSims, alpha=.05) 
{
  nPreds <- nrow(modmat)
  response_distr <- fit$family$family
  inverselink <- fit$family$linkinv
  sims <- ciTools:::get_negbin_sims(fit, nSims)
  sim_response <- matrix(0, ncol = nSims, nrow = nPreds)
  for (i in 1:nSims) {
    yhat <- inverselink(modmat %*% sims@coef[i, ])
    sim_response[, i] <- MASS::rnegbin(n = nPreds, mu = yhat, theta = fit$theta)
  }
  t(round(apply(sim_response, 1, function(x)quantile(x, probs=c(alpha/2, 1-alpha/2)))))
}

然后，您只需使用适当的模型矩阵调用该函数，另一端的结果就是您想要的预测区间。

set.seed(519)
newdat <- data.frame(speed = c(12,15,19))
my_sim_response_nb(modmat = cbind(1, newdat$speed), 
                   fit = fit, 
                   nSims = 2500, 
                   alpha = .05)
#>      2.5% 97.5%
#> [1,]   10    50
#> [2,]   14    68
#> [3,]   23   101

创建于2023-09-04使用reprex v2.0.2

赞(0）回复(0）举报 2023-09-27

我来回答

负二项回归的预测区间(R)

1条答案

相关问题

热门标签

最新问答