如何为smooth.spline获取置信区间？

Question

我使用了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%的置信区间或其他什么。如果有人能解释一下就太好了。

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

- Yu Deng

Answer 1

我不确定smooth.spline的置信区间像lowess那样有"好"的置信区间。但我从CMU数据分析课程中找到了一个代码示例，可以使用贝叶斯自助法来计算置信区间。

以下是使用的函数和示例。主要函数是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)

这就得到了类似于以下的内容：

实际上，使用jackknife残差可能有更参数化的方法来计算置信区间。此代码来自于smooth.spline的S+帮助页面

  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帮助文件，您会发现有一个默认值为TRUE的se=参数，并且文档中写道：

当为TRUE（默认）时，在正在绘制的光滑估计值的上方和下方2个标准误差处添加上下线到1-d图中，对于2-d图，对于+1和-1标准误差的曲面轮廓在估计值的轮廓图上被叠加。如果提供了一个正数，则在计算标准误差曲线或曲面时，将该数字乘以标准误差。另请参见下面的shade。

因此，您可以通过调整此参数来调整置信区间。（这将在print()调用中进行。）

- MrFlick

---

使用Bootstrap方法的第一种方法是有意义的，但与我从gam()中得到的结果完全不同。使用jackknife剩余项的方法也给出了非常独特的图案，并且这个绘图中的CIs非常崎岖。 - Yu Deng

@YuDeng 好的，不同的方法会得出不同的结果。我想你需要决定你是频率学派还是贝叶斯学派。你可能需要咨询一位统计学家，以确定哪种方法最适合你的数据。 - MrFlick

@MrFlick这是一篇旧帖子，关于Jackknife残差方法，您还记得重新乘以杠杆fit $ lev来计算区间宽度背后的逻辑吗？ - foreignvol

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

- Peter Straka