如何获得 smooth.spline 的置信区间？

Posted 2023-03-12

【中文标题】如何获得 smooth.spline 的置信区间？

【英文标题】：How to get confidence interval for smooth.spline?

【发布时间】：2014-07-14 04:29:43

【问题描述】：

我使用smooth.spline 为我的数据估计三次样条。但是当我使用方程计算 90% 的逐点置信区间时，结果似乎有点偏离。如果我做错了，有人可以告诉我吗？我只是想知道是否有一个函数可以自动计算与smooth.spline 函数关联的逐点区间。

boneMaleSmooth = smooth.spline( bone[males,"age"], bone[males,"spnbmd"], cv=FALSE)
error90_male = qnorm(.95)*sd(boneMaleSmooth$x)/sqrt(length(boneMaleSmooth$x))

plot(boneMaleSmooth, ylim=c(-0.5,0.5), col="blue", lwd=3, type="l", xlab="Age", 
     ylab="Relative Change in Spinal BMD")
points(bone[males,c(2,4)], col="blue", pch=20)
lines(boneMaleSmooth$x,boneMaleSmooth$y+error90_male, col="purple",lty=3,lwd=3)
lines(boneMaleSmooth$x,boneMaleSmooth$y-error90_male, col="purple",lty=3,lwd=3)

因为我不确定我是否做得正确，所以我使用了mgcv 包中的gam() 函数。

它立即给出了一个置信区间，但我不确定它是 90% 还是 95% CI 还是其他什么。如果有人能解释一下就太好了。

males=gam(bone[males,c(2,4)]$spnbmd ~s(bone[males,c(2,4)]$age), method = "GCV.Cp")
plot(males,xlab="Age",ylab="Relative Change in Spinal BMD")

【参考方案1】：

我不确定smooth.spline 的置信区间是否像lowess 那样具有“不错的”置信区间。但我从CMU Data Analysis course 中找到了一个代码示例，用于制作贝叶斯引导置信区间。

以下是使用的函数和示例。主要功能是spline.cis，其中第一个参数是一个数据框，其中第一列是x 值，第二列是y 值。另一个重要参数是B，它表示要执行的引导复制次数。 （有关完整详细信息，请参阅上面链接的 PDF。）

# Helper functions
resampler <- function(data) 
    n <- nrow(data)
    resample.rows <- sample(1:n,size=n,replace=TRUE)
    return(data[resample.rows,])


spline.estimator <- function(data,m=300) 
    fit <- smooth.spline(x=data[,1],y=data[,2],cv=TRUE)
    eval.grid <- seq(from=min(data[,1]),to=max(data[,1]),length.out=m)
    return(predict(fit,x=eval.grid)$y) # We only want the predicted values


spline.cis <- function(data,B,alpha=0.05,m=300) 
    spline.main <- spline.estimator(data,m=m)
    spline.boots <- replicate(B,spline.estimator(resampler(data),m=m))
    cis.lower <- 2*spline.main - apply(spline.boots,1,quantile,probs=1-alpha/2)
    cis.upper <- 2*spline.main - apply(spline.boots,1,quantile,probs=alpha/2)
    return(list(main.curve=spline.main,lower.ci=cis.lower,upper.ci=cis.upper,
    x=seq(from=min(data[,1]),to=max(data[,1]),length.out=m)))


#sample data
data<-data.frame(x=rnorm(100), y=rnorm(100))

#run and plot
sp.cis <- spline.cis(data, B=1000,alpha=0.05)
plot(data[,1],data[,2])
lines(x=sp.cis$x,y=sp.cis$main.curve)
lines(x=sp.cis$x,y=sp.cis$lower.ci, lty=2)
lines(x=sp.cis$x,y=sp.cis$upper.ci, lty=2)

这给出了类似的东西

实际上，使用折刀残差计算置信区间似乎可能有一种更参数化的方法。此代码来自S+ help page for smooth.spline

  fit <- smooth.spline(data$x, data$y)      # smooth.spline fit
  res <- (fit$yin - fit$y)/(1-fit$lev)      # jackknife residuals
sigma <- sqrt(var(res))                     # estimate sd

upper <- fit$y + 2.0*sigma*sqrt(fit$lev)   # upper 95% conf. band
lower <- fit$y - 2.0*sigma*sqrt(fit$lev)   # lower 95% conf. band
matplot(fit$x, cbind(upper, fit$y, lower), type="plp", pch=".")

结果是

就gam 置信区间而言，如果您阅读print.gam 帮助文件，则有一个se= 参数默认为TRUE，文档说

当 TRUE（默认值）时，在 2 个标准误差上方和下方的 1-d 绘图中添加上线和下线，而对于 2-d 绘图，表面处于 +1 和 -1 标准误差被轮廓化并覆盖在轮廓图上以进行估计。如果提供了一个正数，则在计算标准误差曲线或曲面时，该数字乘以标准误差。另请参阅下面的阴影。

所以你可以通过调整这个参数来调整置信区间。 （这将在print() 调用中。）

【讨论】：

第一种使用引导方法的方法是有意义的，但与我从gam() 得到的结果相比，结果给出了完全不同的模式。使用折刀残差的那个也给出了一个非常明显的模式，并且这个图中的 CI 非常颠簸。

@YuDeng 好吧，不同的方法给出不同的结果。我想你必须决定你是常客还是贝叶斯主义者。您可能希望咨询统计学家，了解哪种方法最适合您的数据。

【参考方案2】：

R 包mgcv 计算平滑样条曲线和贝叶斯“置信区间”。这些不是通常（频率论者）意义上的置信区间，但数值模拟表明几乎没有差异；请参阅 mgcv 帮助文件中 Marra 和 Wood 的链接文件。

library(SemiPar)
data(lidar)
require(mgcv)

fit=gam(range~s(logratio), data = lidar)
plot(fit)
with(lidar, points(logratio, range-mean(range)))