Как построить доверительные интервалы для моделей glm (гамма-семейство)?

Question

Я хотел бы построить линию и 95% доверительный интервал из модели glm (семейная гамма). Для линейных моделей мне ранее удавалось построить доверительные интервалы из прогнозов, так как они включали подгонку, нижний и верхний уровни и с использованием многоугольников, но я не знаю, как это сделать здесь, поскольку прогнозы не включают верхний и нижний уровни , Я также пробовал ggplot, но там, похоже, сглаживание кривой. Заранее спасибо за помощь. См. Код:

library(ggplot2)

# Data
dat <- data.frame(c(45,75,85,2,14,45,45,45,45,45,55,55,65,85,15,15,315,3,40,85,125,115,13,105,
                    145,125,145,125,205,125,155,125,19,17,145,14,85,65,135,45,40,15,14,10,15,10,10,45,37,30),
                  c(1.928607e-01, 3.038813e-01, 8.041174e-02, 0.000000e+00, 1.017541e-02, 1.658876e-01, 2.084661e-01,
                  1.891305e-01, 2.657766e-01, 1.270864e-01, 1.720141e-01, 1.644947e-01, 7.038978e-02, 3.046604e-01,
                  3.111646e-02, 9.443539e-04, 3.590906e-02, 0.000000e+00, 2.384494e-01, 5.955332e-02, 7.703567e-02,
                  5.524471e-02, 9.915716e-04, 1.169936e-01, 1.409448e-01, 1.411809e-01, 1.025096e-01, 2.649503e-01,
                  6.309465e-02, 3.727837e-02, 8.855679e-02, 1.707864e-01, 1.714002e-02, 1.038789e-03, 1.208065e-01,
                  3.541327e-04, 7.492268e-02, 9.633591e-02, 7.414359e-02, 2.235050e-01, 1.489010e-01, 2.478929e-03,
                  2.573364e-03, 5.430035e-04, 1.719905e-02, 1.243006e-02, 6.822957e-03, 1.927544e-01, 1.146918e-01, 9.030385e-03))

colnames(dat) <- c("age", "wood")  

# Model
model<- glm(wood+0.001 ~ log(age) + I(log(age)^2), data=dat, family = Gamma)
summary(model)

p<-predict(model, data.frame(age=1:200), interval="confidence", level=.95)
p.tr <- 1/p # inverse link according to ?glm

# Plot
plot(1:200, p.tr, type="n", ylim = c(0,.4), xlab="Forest age",
     ylab="Proportion",
     main="Wood production", yaxt="n")
axis(2, las=2)
lines(1:200, p.tr, ylim=range(p.tr), lwd=2, col=rgb(0, .4, 1))
# How can I add to this plot the 95% confidence intervals of the model?

# Ggplot
# I use this function because there was a warning of "Ignoring unknown parameters: family" and this solves that
binomial_smooth <- function(...) { 
  geom_smooth(method = "glm", method.args = list(family = "binomial"),formula=y~log(x)+I(log(x)^2),se=FALSE) 
} 

ggplot(dat, aes(x=age, 
                y=wood+0.001)) +
  binomial_smooth() +
  xlab("Forest age") +
  ylab("Proportion") +
  ggtitle("Wood production") +
  xlim(0, 200) +
  ylim(0,0.4) +
  theme_bw() +
  theme (plot.title = element_text(hjust = 0.5), legend.position = "none")

# Why I get this warning (Warning: In eval(family$initialize) : non-integer #successes in a binomial glm!)?
# Why is the curve more smooth here?

Answer 1

Я не математик / статистика, но я думаю, что "family = binomial" дает вам неуместные оценки, так как это неправильное распределение, поскольку ни дерево, ни возраст не являются исчисляемым числом значений.

О доверительных интервалах: я использовал stat_smooth (), см. Ниже. Должен быть таким же, как и geom_smooth ().

dat <- data.frame(c(45,75,85,2,14,45,45,45,45,45,55,55,65,85,15,15,315,3,40,85,125,115,13,105,
                    145,125,145,125,205,125,155,125,19,17,145,14,85,65,135,45,40,15,14,10,15,10,10,45,37,30),
                  c(1.928607e-01, 3.038813e-01, 8.041174e-02, 0.000000e+00, 1.017541e-02, 1.658876e-01, 2.084661e-01,
                    1.891305e-01, 2.657766e-01, 1.270864e-01, 1.720141e-01, 1.644947e-01, 7.038978e-02, 3.046604e-01,
                    3.111646e-02, 9.443539e-04, 3.590906e-02, 0.000000e+00, 2.384494e-01, 5.955332e-02, 7.703567e-02,
                    5.524471e-02, 9.915716e-04, 1.169936e-01, 1.409448e-01, 1.411809e-01, 1.025096e-01, 2.649503e-01,
                    6.309465e-02, 3.727837e-02, 8.855679e-02, 1.707864e-01, 1.714002e-02, 1.038789e-03, 1.208065e-01,
                    3.541327e-04, 7.492268e-02, 9.633591e-02, 7.414359e-02, 2.235050e-01, 1.489010e-01, 2.478929e-03,
                    2.573364e-03, 5.430035e-04, 1.719905e-02, 1.243006e-02, 6.822957e-03, 1.927544e-01, 1.146918e-01,
                    9.030385e-03))

colnames(dat) <- c("age", "wood")

model<- glm(wood+0.001 ~ log(age) + I(log(age)^2), data=dat, family = Gamma)
#summary(model)

p<-predict(model, data.frame(age=1:200), interval="confidence", level=.95)
p.tr <- 1/p # inverse link according to ?glm

prediction <- data.frame(age = as.numeric(names(p)), wood = 1/p)


ggplot(data = dat, aes(x = age, y = wood)) +
  geom_point() +
  geom_line(data= prediction) +
  stat_smooth(data = dat, method = "glm",
              formula = y+0.001 ~ log(x) + I(log(x)^2),
              method.args = c(family = Gamma))