模擬一些資料
n <- 50
beta0 <- 2
beta1 <- 0.32
x <- runif(n=n, min=0, max=5)
mu <- exp(beta0 beta1 * x)
y <- rpois(n=n, lambda=mu)
data <- data.frame(x=x, y=y)
plot(data$x, data$y)
應用泊松回歸
glm.pois <- glm(formula= y ~ x, data=data, family=poisson(link="log"))
summary(glm.pois)
計算組態檔 95% CI
(pCI <- confint(glm.pois))
nice.xs <- sort(data$x)
pred.pCI <- apply(t(pCI), 1, FUN=function(x){exp( x[1] x[2]*nice.xs)})
使用 SE 計算 95% CI
new.dat <- data.frame(phat = predict(glm.pois, type="response"),
x = data$x)
new.dat <- new.dat[with(new.dat, order(x)), ]
ilink <- poisson()$linkinv
## Add fit and se.fit on the link scale
new.dat <- bind_cols(new.dat,
setNames(as_tibble(predict(glm.pois, new.dat,
se.fit = TRUE)[1:2]),
c('fit_link','se_link')))
new.dat <- mutate(new.dat,
fit = ilink(fit_link),
right_upr = ilink(fit_link (2 * se_link)),
right_lwr = ilink(fit_link - (2 * se_link)))
陰謀
ggplot(data=new.dat, aes(x=x, y= phat))
geom_point()
geom_ribbon(aes(x=x, ymin =right_lwr, ymax = right_upr),
fill = "darkgreen", alpha=0.3)
geom_ribbon(aes(x=nice.xs, ymin = pred.pCI[,1], ymax = pred.pCI[,2]),
fill = "darkorange", alpha=0.3)
theme_classic(base_size=15)
我的問題是為什么輪廓 95% CI（橙色絲帶）這么大？我想我計算它們是錯誤的，但我不知道我的錯誤在哪里。
uj5u.com熱心網友回復：
截距和斜率之間有很強的相關性：
cov2cor(vcov(glm.pois))
(Intercept) x
(Intercept) 1.0000000 -0.9147526
x -0.9147526 1.0000000
您的組態檔 CI 正在計算 (logistic(intercept_lwr slope_lwr*x)和logistic(intercept_upr slope_upr*x)，忽略了這種相關性。
計算派生值（即引數本身以外的數量）的剖面置信區間并不容易。
轉載請註明出處，本文鏈接：https://www.uj5u.com/gongcheng/354135.html
標籤：r
上一篇：從服務類呼叫方法時，SpringBoot@AutowiredNullPointerException
下一篇：一種熱編碼僅R配方中的因子變數