ggplot2를 사용하여 견고한 표준 오류로 예측 된 값 플롯

Question

ggplot2를 사용하여 음수 이항 회귀의 예측 값을 그리는 중입니다. 이진 변수가 켜져 있고 다른 하나는 꺼져 있습니다. 비교할 수있는 두 개의 플롯이 있습니다.

링크 here은 페이지 하단에서이를 수행하는 방법을 보여 주지만 견고한 표준 오류를 사용하여 예측 된 값의 플롯 주변에 음영을 만들 수 있기를 원합니다. 나는 이것을 predict() 함수에서 얻는 방법을 모르겠습니다. 이 코드 예제에서 주위에 그늘진 선에 강한 표준 오류를 얻으려는 작업이 있습니까?

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

그리고 내가 사용하고있는 데이터의 샘플이 있습니다 :

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 

내가 사용하고있는 모델은 이것이다 :

나는 강력한 표준 오차를 생성하기 위해 여기 this site의 코드를 사용 :

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)) 

Answer 1

제공하는 링크는 모델을 만들고, ne 예측자는 전체 범위에 따라 달라지며 모델 및 합성 데이터 집합을 전달() 한 다음 결과 예측을 표시합니다. 다르게 수행해야하는 유일한 작업은 강력한 표준 데이터를 데이터 프레임에 넣어 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") 

(합성 dataframe 모델에 필요한 모든 변수가 될 때까지 위의 계속)

#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)