2

我正在尝试使用 ggplot2 绘制负二项式回归的预测值，一个打开二进制变量，另一个关闭它。所以会有两个可以比较的两个地块。

此处的链接演示了如何在页面底部执行此操作，但我希望能够使用稳健的标准误差在预测值的图周围创建阴影。我看不出如何从 predict() 函数中得到这个。此代码示例是否有任何解决方法来获得强大的标准错误以在绘制的线周围进行阴影处理？

我使用此站点的代码来生成可靠的标准错误：

require(sandwich)
cov.nb1 <- vcovHC(nb1, type = "HC0")
std.err <- sqrt(diag(cov.nb1))
r.est <- cbind(Estimate = coef(nb1), `Robust SE` = std.err, `Pr(>|z|)` = 2 *
    pnorm(abs(coef(nb1)/std.err), lower.tail = FALSE), LL = coef(nb1) - 1.96 *
    std.err, UL = coef(nb1) + 1.96 * std.err)

r.est

我使用的模型是这样的：

nb1 <- glm.nb(citecount ~ expbin*novcr + expbin*I(novcr^2) + disease + length +
as.factor(year), data = nov4d.dt)

我正在使用的数据样本是这样的：

nov4d.dt  <-
    structure(list(PMID = c(1279136L, 1279186L, 1279186L, 1279187L, 
    1279187L, 1279190L, 1279257L, 1279317L, 1279332L, 1279523L), 
        min = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), max = c(32L, 
        32L, 32L, 32L, 32L, 32L, 32L, 32L, 32L, 32L), mean = c(11L, 
        13L, 13L, 19L, 19L, 16L, 24L, 15L, 8L, 19L), length = c(45L, 
        120L, 120L, 78L, 78L, 136L, 45L, 36L, 171L, 78L), threslength = c(13L, 
        20L, 20L, 7L, 7L, 26L, 4L, 6L, 77L, 14L), novlength = c(5L, 
        6L, 6L, 3L, 3L, 6L, 3L, 3L, 36L, 5L), novind = c("TRUE", 
        "TRUE", "TRUE", "TRUE", "TRUE", "TRUE", "TRUE", "TRUE", "TRUE", 
        "TRUE"), novcr = c(0.111111, 0.05, 0.05, 0.0384615, 0.0384615, 
        0.0441176, 0.0666667, 0.0833333, 0.210526, 0.0641026), novcrt = c(0.288889, 
        0.166667, 0.166667, 0.0897436, 0.0897436, 0.191176, 0.0888889, 
        0.166667, 0.450292, 0.179487), year = c(1991L, 1991L, 1992L, 
        1992L, 1992L, 1992L, 1992L, 1992L, 1991L, 1992L), disease = structure(c(1L, 
        4L, 2L, 4L, 2L, 1L, 4L, 4L, 2L, 4L), .Label = c("alz", "bc", 
        "cl", "lc"), class = "factor"), citecount = c(5L, 8L, 8L, 
        12L, 12L, 0L, 1L, 0L, 92L, 0L), novind2 = c(TRUE, TRUE, TRUE, 
        TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE), rad = c(FALSE, 
        FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE
        ), exp = c(260, 351, 351, 65, 65, 480, 104, 273, 223, 0), 
        novind4 = c(FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, 
        FALSE, TRUE, FALSE), novind5 = c(FALSE, FALSE, FALSE, FALSE, 
        FALSE, FALSE, FALSE, FALSE, TRUE, FALSE), novind6 = c(FALSE, 
        FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE
        ), expbin = c(TRUE, TRUE, TRUE, FALSE, FALSE, TRUE, FALSE, 
        TRUE, TRUE, FALSE), expbin2 = c(TRUE, TRUE, TRUE, FALSE, 
        FALSE, TRUE, FALSE, TRUE, TRUE, FALSE)), .Names = c("PMID", 
    "min", "max", "mean", "length", "threslength", "novlength", "novind", 
    "novcr", "novcrt", "year", "disease", "citecount", "novind2", 
    "rad", "exp", "novind4", "novind5", "novind6", "expbin", "expbin2"
    ), sorted = "PMID", class = c("data.frame"), row.names = c(NA, 
    -10L))
4

1 回答 1

3

您提供的链接创建一个模型，创建一个合成数据集，其中一个预测变量在其整个范围内变化，将模型和合成数据集传递给 predict()，然后绘制结果预测。您需要做的唯一重要的事情是将强大的 std.err 放入数据框中以计算 CI。

#look at how model thinks citecount ~ novcr for two values of expbin 
#make synthetic data with a range of range(df$novcr)
#include logical predictor variable expbin
#such that each level of expbin has all the novcr values

newdata2 <- data.frame(novcr = rep(seq(from = min(nov4d.dt$novcr), 
    to = max(nov4d.dt$novcr), length.out = 100), 2), 
    expbin  = rep(0:1, each = 100))

#convert expbin type to logical
newdata2$expbin <- as.logical(newdata2$expbin)

# add in the mean or default values of other predictors
# because I assume predict() needs vals for all parameters in the model
newdata2$length <- mean(nov4d.dt$length,na.rm=T)
newdata2$disease <- factor("alz")
newdata2$year <- factor("1992")

（继续上述操作，直到合成数据框具有模型所需的所有变量）

#make predict and add it to synthetic data
newdata2$fit <- predict(nb1, newdata2, type = "response")

# include CIs based on your robust se
newdata2$LL <- newdata2$fit - 1.96 * std.err["novcr"]
newdata2$UL <- newdata2$fit + 1.96 * std.err["novcr"]

#plot
ggplot(newdata2, aes(novcr, fit)) + 
    geom_ribbon(aes(ymin = LL, ymax = UL, fill = expbin), 
    alpha = 0.25) + geom_line(aes(colour = expbin), size = 2)
于 2012-11-06T23:57:22.470 回答