在半对数图上绘制二次方程的置信区间

Question

给定一个数据框，Data 形式

     x          y
1  250 1.00000000
2  345 0.03567766
3  290 0.16654457
4  260 0.58363858
5  270 0.38754579
6  280 0.24713065
7  290 0.17142857
8  300 0.11709402
9  310 0.09047619
10 320 0.06439560
11 330 0.05098901

我能够使用

推导和绘制数据拟合

library(ggplot2)
Data$x2<-Data$x^2
quadratic.model <- lm(log(Data$y) ~ Data$x + Data$x2)

fun_quad <- function(x){return(exp(
    quadratic.model$coef[[3]] * x ^ 2 + 
    quadratic.model$coef[[2]] * x    + 
    quadratic.model$coef[[1]]  
    ))}

chartObj <- ggplot() +

    stat_function(
        fun         = fun_quad,
        aes(color   = factor(0)),
        size        = 1.3, 
        linetype    = "dotdash"
        )+  
            
    geom_point(data = Data,
        aes(x = x, y = y, fill = factor(0)),
        color = "black", shape = 22, stroke = 0.7, size = 2.2) +
            
    coord_trans(y = 'log10', 
            limx = c(250,350), limy = c(.025,1))+ 
            
    theme_bw() + 
    
    guides(fill=F,color=F,linetype=F)
    
chartObj

呈现

.

我还尝试使用 confint 和 geom_ribbon 绘制 CI。

ribbon.ymin <- function(x){return(exp(
    confint(quadratic.model)[[3]]*x^2 + 
    confint(quadratic.model)[[2]]*x   + 
    confint(quadratic.model)[[1]] 
    ))}

ribbon.ymax <- function(x){return(exp(
    confint(quadratic.model)[[6]]*x^2 + 
    confint(quadratic.model)[[5]]*x   + 
    confint(quadratic.model)[[4]]
    ))}
            
            
ribbonData <- as.data.frame(cbind(x = seq(250,350,.01)))
attach(ribbonData)
ribbonData$ymin     <- ribbon.ymin(x)
ribbonData$ymax     <- ribbon.ymax(x)
ribbonData$y        <- fun_quad(x)
detach(ribbonData)
head(ribbonData)
            
chartObj <- chartObj + 
            
    geom_ribbon( data = ribbonData,
            aes(x = x, y = 0:0,
                ymin = ymin, ymax = ymax,
                color = factor(0),fill = factor(0)),
            alpha = 0.3)

然而，这呈现如下，再次感觉明显不正确。

那么，如何绘制与 quadratic.model 描述的函数关联的置信区间？

更新

我认为我已经通过使用 predict 命令几乎找到了我正在寻找的东西，具体如下所示，但是这仍然有一些不足之处，尤其是生产色带的边缘。

Data$x2<-Data$x^2
quadratic.model <- lm(log(Data$y) ~ Data$x + Data$x2)

fun_quad <- function(x){return(exp(
    quadratic.model$coef[[3]] * x ^ 2 + 
    quadratic.model$coef[[2]] * x    + 
    quadratic.model$coef[[1]]  
    ))}

ribbonData<-predict(quadratic.model,data.frame(x=Data$x),interval="predict",level=.95)
# "predict" used over "confidence" in this example to show the rough edges better.
ribbonData<-as.data.frame(cbind(x=Data$x,fit=ribbonData[,1],lower=ribbonData[,2],upper=ribbonData[,3]))
ribbonData[,2:4]<-exp(ribbonData[,2:4])

chartObj <- ggplot() +

    geom_ribbon( data = ribbonData,
        aes(x = x, y = fit,
            ymin = lower, ymax = upper,
            color = factor(0),fill = factor(0)),
        alpha = 0.3) +
    
    stat_function(
        fun         = fun_quad,
        aes(color   = factor(0)),
        size        = 1.3, 
        linetype    = "dotdash"
        )+  
            
    geom_point(data = Data,
        aes(x = x, y = y, fill = factor(0)),
        color = "black", shape = 22, stroke = 0.7, size = 2.2) +
            
    coord_trans(y = 'log10', 
            limx = c(250,350), limy = c(.025,1))+ 
            
    theme_bw() + 
    
    guides(fill=F,color=F,linetype=F)

是否有更好的方式来表示上图所呈现的信息？要平滑色带的粗糙边缘？

Answer 1

它可能 "feel obviously incorrect"，但它绘制了它被问到的内容。看不到整个区间，因为设置了limx和limy：

ribbon <- function(x, level = 0.95) {
  data.frame(
    x,
    ymin = exp(
      confint(quadratic.model, level = level)[[3]] * x ^ 2 + 
        confint(quadratic.model, level = level)[[2]] * x + 
        confint(quadratic.model, level = level)[[1]] 
    ),
    ymax = exp(
      confint(quadratic.model, level = level)[[6]]*x^2 + 
        confint(quadratic.model, level = level)[[5]]*x   + 
        confint(quadratic.model, level = level)[[4]]
    )
  )
}

chartObj +
  coord_trans(y = 'log10') +
  geom_ribbon(data = ribbon(seq(250, 350, .01), level = 0.95),
              aes(x = x, ymin = ymin, ymax = ymax,
                  color = factor(0), fill = factor(0)),
              alpha = 0.3)

（注意：我的回答完全是关于使用 ggplot2 进行编程，并没有提及对置信区间求幂的统计有效性）。

根据 OP 的更新问题进行编辑（平滑色带边缘）。

predict()过分多：

quadratic.model <- lm(log(y) ~ x + x2, data = Data)

ribbonData <- data.frame(x = seq(250, 350, 0.01), x2 = seq(250, 350, 0.01) ^ 2)
ribbonData <- cbind(
  ribbonData,
  predict(quadratic.model, ribbonData,
          interval = "prediction", level = 0.95)
)
# "predict" used over "confidence" in this example to show the rough edges better.
ribbonData[, 3:5] <- exp(ribbonData[, 3:5])

ggplot() +

  geom_ribbon( data = ribbonData,
               aes(x = x, y = fit,
                   ymin = lwr, ymax = upr,
                   color = factor(0),fill = factor(0)),
               alpha = 0.3) +

  stat_function(
    fun         = fun_quad,
    aes(color   = factor(0)),
    size        = 1.3, 
    linetype    = "dotdash"
  ) +  

  geom_point(data = Data,
             aes(x = x, y = y, fill = factor(0)),
             color = "black", shape = 22, stroke = 0.7, size = 2.2) +

  coord_trans(y = 'log10', 
              limx = c(250, 350), limy = c(.025, 1)) + 

  theme_bw() + 

  guides(fill = F, color = F, linetype = F)

在半对数图上绘制二次方程的置信区间

Plot the Confidence Interval of a Quadratic on a Semi-Log plot

plot

r

ggplot2

confidence-interval

tidyverse

更新