Tuesday, September 23, 2014

Các hàm R phổ biến 8: Phân tích bootstrap

Mã R dùng cho phân tích bootstrap. Tôi có thêm một ví dụ dùng phương pháp bootstrap cho việc phân tích sự khác biệt giữa 2 tỉ lệ bằng mô hình Bayes. 




Phân tích 
Hàm R
Lấy mẫu ngẫu nhiên

x = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
sample(x, 5, replace=T)
sample(x, 5, replace=F)

Khoảng tin cậy 95% cho trung vị 
x = c(0.05,0.15,0.35,0.25, ...);
n = length(pain);
B=1000;
median = numeric(B);
for (i in 1:B)
{
    bs.sample = sample(x, n, replace=T)
    median[i] = median(bs.sample)
 }
quantile(median, probs=c(0.025, 0.975))
hist(median, col="blue", border="white", main="")

So sánh 2 số trung vị
control = c(52,104,146, ...);
treatment = c(94, 197, 16, ...);
B=1000;
con = rep(NA, B);
trt = rep(NA, B);
for (i in 1:B) {
  con[i] = median(sample(control, replace=T))
  trt[i] = median(sample(treatment, replace=T))
  dif[i] = trt[i]-con[i]
}
quantile(dif, c(0.025, 0.5, 0.975))

So sánh 2 tỉ lệ
n = c(11037, 11034);
disease = c(104, 189);
aspirin = c(rep(1, disease[1]), rep(0, n[1]-disease[1]));
control = c(rep(1, disease[2]), rep(0, n[2]-disease[2]));
B = 10000;
bs1 = rep(NA, B);
bs2 = rep(NA, B);
for (i in 1:B) {
  resample1 = sample(aspirin, n[1], replace=T)
  dis1 = sum(resample1)
  resample2 = sample(control, n[2], replace=T)
  dis2 = sum(resample2)
  bs1[i] = dis1/n[1]
  bs2[i] = dis2/n[2]
}

RR = bs1/bs2;
quantile(RR, c(0.025, 0.50, 0.975)


So sánh 2 tỉ lệ, phương pháp Bayes
x = 11; n1 = 20; alpha1 = 1; beta1 = 1;
y = 5;  n2 = 20; alpha2 = 1; beta2 = 1;
p1 = rbeta(10000, x + alpha1, n1 - x + beta1);
p2 = rbeta(10000, y + alpha2, n2 - y + beta2);
rd = p2 - p1;
plot(density(rd));
quantile(rd, c(.025, 0.5, 0.975))




No comments:

Post a Comment