首页 文章

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

提问于
浏览 次
3

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

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)

enter image description here

因为我不确定我是否正确使用了，所以我使用了 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")

enter image description here

r smoothing confidence-interval splines mgcv

2 回答

  • 8

    我不确定 smooth.spline 的置信区间是否像 lowess 那样具有"nice"置信区间 . 但是我发现了一个来自CMU Data Analysis course的代码样本来制作贝叶斯bootstap置信区间 .

    以下是使用的功能和示例 . 主要功能是 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)

    这给了类似的东西

    bootstrap confidence intervals

    实际上，看起来可能存在使用折刀残差计算置信区间的更具参数的方法 . 此代码来自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=".")

    这导致了

    residual CI estimate

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

    当TRUE（默认）上下线被添加到1-d图时，2个标准误差高于和低于平滑绘图的估计值，而对于2-d图，1和-1标准误差的曲面轮廓和覆盖在轮廓图上以进行估算 . 如果提供正数，则在计算标准误差曲线或曲面时，此数字乘以标准误差 . 另见下面的阴影 .

    因此，您可以通过调整此参数来调整置信区间 . （这将在 print() 电话中 . ）

    回复于
  • 4

    R包 mgcv 计算平滑样条和贝叶斯"confidence intervals."这些不是通常（频率论）意义上的置信区间，但数值模拟表明几乎没有区别;请参阅 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)))

    enter image description here

    回复于

相关问题