library(MKinfer)
library(resample)
library(tidyverse)
QK Box 2.3
This box continues with the Low et al. anesthetic example from Box 2.2
Preliminaries
Packages: MKinfer, resample
Use low data:
<- read_csv("data/lowco2.csv") low
Get jackknife SE for two groups
<- subset(low,anesth=="iso")
low1 jackknife(low1$co2,mean)
## Call:
## jackknife(data = low1$co2, statistic = mean)
## Replications: 12
##
## Summary Statistics:
## Observed SE Mean Bias
## mean 50 3.2891 50 0
<- subset(low,anesth=="ac")
low2 jackknife(low2$co2,mean)
## Call:
## jackknife(data = low2$co2, statistic = mean)
## Replications: 11
##
## Summary Statistics:
## Observed SE Mean Bias
## mean 70.90909 6.090909 70.90909 1.421085e-13
Get bootstrap SE and 95%CI
<- bootstrap(low1$co2,mean,R=9999)
low1boot
low1boot## Call:
## bootstrap(data = low1$co2, statistic = mean, R = 9999)
## Replications: 9999
##
## Summary Statistics:
## Observed SE Mean Bias
## mean 50 3.146432 49.96038 -0.03962063
CI.percentile(low1boot, probs=c(0.025,0.975))
## 2.5% 97.5%
## mean 43.66667 58.11949
CI.bca(low1boot, probs=c(0.025,0.975))
## 2.5% 97.5%
## mean 44.58333 60.16667
<- bootstrap(low2$co2,mean,R=9999)
low2boot
low2boot## Call:
## bootstrap(data = low2$co2, statistic = mean, R = 9999)
## Replications: 9999
##
## Summary Statistics:
## Observed SE Mean Bias
## mean 70.90909 5.850436 70.92512 0.01602888
CI.percentile(low2boot, probs=c(0.025,0.975))
## 2.5% 97.5%
## mean 59 85.63636
CI.bca(low2boot, probs=c(0.025,0.975))
## 2.5% 97.5%
## mean 60.27273 88.09091
Get bootstrap SE and CI on difference
<- bootstrap2(low$co2,mean,treatment=low$anesth,R=9999,ratio=FALSE)
lowboot
lowboot## Call:
## bootstrap2(data = low$co2, statistic = mean, treatment = low$anesth,
## R = 9999, ratio = FALSE)
## Replications: 9999
## Two samples, sample sizes are 11 12
##
## Summary Statistics for the difference between samples 1 and 2:
## Observed SE Mean Bias
## mean: ac-iso 20.90909 6.666543 20.86282 -0.04626826
CI.percentile(lowboot, probs=c(0.025,0.975))
## 2.5% 97.5%
## mean: ac-iso 6.016849 36.92163
Randomization test
perm.t.test(co2~anesth, data=low, R=9999, paired= FALSE)
##
## Permutation Welch Two Sample t-test
##
## data: co2 by anesth
## (Monte-Carlo) permutation p-value = 0.005001
## permutation difference of means (SE) = 20.92096 (7.964221)
## 95 percent (Monte-Carlo) permutation percentile confidence interval:
## 5.75000 36.06818
##
## Results without permutation:
## t = 3.0206, df = 15.485, p-value = 0.008362
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 6.194866 35.623316
## sample estimates:
## mean in group ac mean in group iso
## 70.90909 50.00000