首页 文章

R:固定设计回归的Bootstrap BCa置信区间

提问于
浏览
0

我有一个固定的设计回归问题,我试图获得引导BCa置信区间,使用R.这是一个例子(使用lmRob),但这仅用于说明:

require(robust)
data(stack.dat)
stack.rob <- lmRob(Loss ~ ., data = stack.dat)

summary(stack.rob)

Call:
lmRob(formula = Loss ~ ., data = stack.dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-8.6299 -0.6713  0.3594  1.1507  8.1740 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -37.65246    5.00256  -7.527 8.29e-07 ***
Air.Flow      0.79769    0.07129  11.189 2.91e-09 ***
Water.Temp    0.57734    0.17546   3.291  0.00432 ** 
Acid.Conc.   -0.06706    0.06512  -1.030  0.31757    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.837 on 17 degrees of freedom
Multiple R-Squared: 0.6205 

Test for Bias:
            statistic p-value
M-estimate      2.751  0.6004
LS-estimate     2.640  0.6197

R中有引导和引导程序包(并且代码为here,但它们都导出了非参数引导程序BCa置信区间 . 但这是一个固定设计的回归设置 . 因此我想知道是否有R软件可用于bootstrap BCa置信区间固定设计回归 . 使用lm的R软件包或类似软件的例子也可以 .

谢谢!

1 回答

  • 0

    我相信我有一个答案,改编自2002年R Journal的例子 .

    require(boot)
    require(robust)
    data(stackloss)
    
    stack.rob <- lmRob(stack.loss ~ ., data = stackloss)
    
    lmRob.coef <- function(data, y, pred) {
        mod <- lmRob(formula = as.formula(eval(paste(y,"~", paste(pred,collapse="+")))) , data = data,  control = lmRob.control(mxr = 1000, mxf = 1000, mxs = 1000))
        coef(mod)
    }
    
    lmRob.results <- function(data, y, pred) {
        mod <- lmRob(formula = as.formula(eval(paste(y,"~", paste(pred,collapse="+")))) , data = data)
        data.frame(fitted = fitted(mod), residuals = resid(mod))
    }
    
    fit.dat <- lmRob.results(data = stackloss, y = "stack.loss", pred = c("Air.Flow", "Water.Temp", "Acid.Conc."))
    
    model.fun <- function(data, i, y, pred, fitted.results) {
        dat <- cbind(data, fitted.results)
        dat[, y] <- dat$fitted + dat$residuals[i]
        lmRob.coef(data = dat, y = y, pred = pred)
    }
    
    lmRob.boot <- boot(data = stackloss, statistic = model.fun, R = 999, y = "stack.loss", pred = c("Air.Flow", "Water.Temp", "Acid.Conc."), fitted.results = fit.dat)
    
    boot.ci(boot.out=lmRob.boot, type="bca", index=4)
    
    BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
    Based on 999 bootstrap replicates
    
    CALL : 
    boot.ci(boot.out = lmRob.boot, type = "bca", index = 4)
    
    Intervals : 
    Level       BCa          
    95%   (-0.2644,  0.2963 )  
    Calculations and Intervals on Original Scale
    

    下一个目标是包括模型中不同情况的权重 .

相关问题