如何獲得smooth.spline的置信區間？

我已經使用smooth.spline來估計我的數據的三次樣條。但是當我使用方程式計算90％的逐點置信區間時，結果似乎有點偏離。有人可以告訴我，如果我做錯了嗎？我只是想知道是否有一個功能，可以自動計算一個點間隔帶與smooth.spline函數關聯。如何獲得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)

enter image description here

因爲我不知道我是否做了正確的，然後我用gam()功能從mgcv包。

它立即給了一個信心帶，但我不確定它是否是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

來源

2014-05-25 Yu Deng

我不知道爲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)

這給像

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（de故障）上下兩條線被添加到上面和下面的2個標準誤差上面和下面的平滑被繪製的1-d曲線，而對於2-d曲線，在+1和-1標準誤差的曲面被曲線化並覆蓋在估計的等高線圖。如果提供了正數，則在計算標準誤差曲線或曲面時，此數字乘以標準誤差。下面請參閱陰影。

因此，您可以通過調整此參數來調整置信區間。（這將在print()呼叫。）

來源

2014-05-25 06:39:18 MrFlick

使用引導方法第一種方法是有意義但結果給比較到一個我從'GAM得到（）'完全不同的模式。使用摺疊殘差的一個也給出了非常不同的模式，並且在這個圖中的CI非常非常顛簸。 –

@YuDeng好吧，不同的方法會產生不同的結果。我想你必須決定你是一名頻率主義者還是貝葉斯人。您可能希望諮詢統計人員以查看哪種方法最適合您的數據。 – MrFlick

R程序包mgcv計算平滑樣條和貝葉斯「置信區間」。這些不是通常（頻率論者）意義上的置信區間，但數值模擬表明幾乎沒有區別;請參閱Marra和Wood在mgcv幫助文件中的鏈接文件。

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

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

來源

2015-04-24 01:25:39

如何獲得smooth.spline的置信區間？

回答

相關問題