Peter Crowther
2014-Feb-25 14:54 UTC
[R] mcr mcreg: Why are BCa and quantile se values calculated, stored, then the stored values set to NA?
In mcreg, there's the following snippet for bootstrap CI method of BCa (and a very similar one for quantile): else if (method.bootstrap.ci == "BCa") { bootB0 <- mc.calc.bca(Xboot = B0, Xjack = B0jack, xhat = glob.coef[1], alpha) bootB1 <- mc.calc.bca(Xboot = B1, Xjack = B1jack, xhat = glob.coef[2], alpha) bootB0$se <- glob.sigma[1] bootB1$se <- glob.sigma[2] bootB0$se <- NA bootB1$se <- NA } [...] mc.res <- mc.make.CIframe(b0 = glob.coef[1], b1 = glob.coef[2], se.b0 = bootB0$se, se.b1 = bootB1$se, CI.b0 = bootB0$CI, CI.b1 = bootB1$CI) I'd quite like to use in my code the SE values that were calculated in mcreg. Can any kind soul shed light on why they're calculated, stored, then the stored values set to NA? Thanks in advance for any insight! - Peter -- Peter Crowther, Director, Melandra Limited
Fabian Model
2014-Mar-17 15:31 UTC
[R] mcr mcreg: Why are BCa and quantile se values calculated, stored, then the stored values set to NA?
This piece of code is indeed confusing. Generally quantile and BCa bootstrap do not estimate a global SE, so the glob.sigma slot you want to access is not really meaningful for those methods. The result object of a BCa bootstrap as calculated by the mc.bootstrap function will contain a glob.sigma slot that is set to the analytical SE estimates of the complete data set for regression methods where analytical estimates of SE are available (LinReg, WLinReg, Deming). For the other regression methods (WDeming, PaBa) glob.sigma will already be NA. Since the analytical SE estimates are not relevant when you perform quantile or BCa bootstrap we decided at some point to always set them NA for these bootstrap methods - independent of regression method. The code copying the glob.sigma results is pointless and we'll remove it in the next release. Now, if you want to calculate the mean or the sd of the individual bootstrap estimates, you can directly access the MCResultResampling and MCResultBCa objects. They store the estimates of regression coefficients for each bootstrap sample in slots B0 (intercept) and B1 (slope). Example: library(mcr) data(creatinine) m <- mcreg(as.matrix(creatinine), method.bootstrap.ci = "quantile", na.rm = TRUE) mean.slope.est<- mean(m at B1) sd.slope.est<- sd(m at B1) mean.intercept.est <- mean(m at B0) sd.intercept.est <- sd(m at B1) However, note that these estimates have no direct relationship to the BCa confidence intervals. Best regards, Fabian On Tuesday, February 25, 2014 3:54:11 PM UTC+1, Peter Crowther wrote:> In mcreg, there's the following snippet for bootstrap CI method of BCa > (and a very similar one for quantile): > > else if (method.bootstrap.ci == "BCa") { > bootB0 <- mc.calc.bca(Xboot = B0, Xjack = B0jack, > xhat = glob.coef[1], alpha) > bootB1 <- mc.calc.bca(Xboot = B1, Xjack = B1jack, > xhat = glob.coef[2], alpha) > bootB0$se <- glob.sigma[1] > bootB1$se <- glob.sigma[2] > bootB0$se <- NA > bootB1$se <- NA > } > [...] > mc.res <- mc.make.CIframe(b0 = glob.coef[1], b1 = glob.coef[2], > se.b0 = bootB0$se, se.b1 = bootB1$se, CI.b0 = bootB0$CI, > CI.b1 = bootB1$CI) > > I'd quite like to use in my code the SE values that were calculated in > mcreg. Can any kind soul shed light on why they're calculated, > stored, then the stored values set to NA? > > Thanks in advance for any insight! > > - Peter > -- > Peter Crowther, Director, Melandra Limited > > ______________________________________________ > R-h... at r-project.org <javascript:> mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide > http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. >