##讀取檔案
data(anorexia, package="MASS")

##anorexia存入dta
dta <- anorexia
head(dta)
##   Treat Prewt Postwt
## 1  Cont  80.7   80.2
## 2  Cont  89.4   80.1
## 3  Cont  91.8   86.4
## 4  Cont  74.0   86.3
## 5  Cont  78.1   76.1
## 6  Cont  88.3   78.1
##資料概覽
summary(dta)
##   Treat        Prewt           Postwt      
##  CBT :29   Min.   :70.00   Min.   : 71.30  
##  Cont:26   1st Qu.:79.60   1st Qu.: 79.33  
##  FT  :17   Median :82.30   Median : 84.05  
##            Mean   :82.41   Mean   : 85.17  
##            3rd Qu.:86.00   3rd Qu.: 91.55  
##            Max.   :94.90   Max.   :103.60
##資料維度
dim(dta)
## [1] 72  3
##列名
row.names(dta)
##  [1] "1"  "2"  "3"  "4"  "5"  "6"  "7"  "8"  "9"  "10" "11" "12" "13" "14" "15"
## [16] "16" "17" "18" "19" "20" "21" "22" "23" "24" "25" "26" "27" "28" "29" "30"
## [31] "31" "32" "33" "34" "35" "36" "37" "38" "39" "40" "41" "42" "43" "44" "45"
## [46] "46" "47" "48" "49" "50" "51" "52" "53" "54" "55" "56" "57" "58" "59" "60"
## [61] "61" "62" "63" "64" "65" "66" "67" "68" "69" "70" "71" "72"
##回歸
anrx_lm <- lm(Postwt ~ Prewt + Treat, data=dta)

##顯示回歸結果
summary(anrx_lm)
## 
## Call:
## lm(formula = Postwt ~ Prewt + Treat, data = dta)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.1083  -4.2773  -0.5484   5.4838  15.2922 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  49.7711    13.3910   3.717  0.00041 ***
## Prewt         0.4345     0.1612   2.695  0.00885 ** 
## TreatCont    -4.0971     1.8935  -2.164  0.03400 *  
## TreatFT       4.5631     2.1333   2.139  0.03604 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.978 on 68 degrees of freedom
## Multiple R-squared:  0.2777, Adjusted R-squared:  0.2458 
## F-statistic: 8.713 on 3 and 68 DF,  p-value: 5.719e-05
##相關
coef(anrx_lm)
## (Intercept)       Prewt   TreatCont     TreatFT 
##  49.7711090   0.4344612  -4.0970655   4.5630627

設定拔靴法公式，改為while-loop

bootstrap_function  <- function(string_formula, nc, model_data, ndraws) {
      coeff_mtx <- matrix(0, nrow = ndraws,  ncol = nc) #coeff_mtx為矩陣
      i <- 1
      while (i < ndraws) { #抽取<ndraws次
        i <- i+1
        bootstrap_ids   <- sample(seq(nrow(model_data)), nrow(model_data), replace = TRUE)
        #從model_data中抽取，抽取後放入
        bootstrap_data <- model_data[bootstrap_ids,]
        bootstrap_model <- lm(as.formula(string_formula), bootstrap_data)
        coeff_mtx[i,]   <- bootstrap_model$coefficients
        if (i == 1) {
          print(bootstrap_model$coefficients)
        }
      }     
      return(coeff_mtx)
    }
    
##按照設定的公式輸入資料與數字
bootstrap_estimates <- bootstrap_function("Postwt ~ Prewt + Treat", nc=4, dta, 500) 

## 列出前後2.5%
bootstrap_CIs <- matrix(0, nrow = 4, ncol = 2)
    for (i in 1:4) {
      bootstrap_CIs[i,1]  <- quantile(bootstrap_estimates[,i], 0.025)
      bootstrap_CIs[i,2]  <- quantile(bootstrap_estimates[,i], 0.975)
    }
## 產出bootstrap_CIs
print(bootstrap_CIs)
##            [,1]       [,2]
## [1,] 14.5694440 75.4969944
## [2,]  0.1196479  0.8410811
## [3,] -7.7418322 -0.3263648
## [4,]  0.2109885  8.6719603
##信賴區間
confint(anrx_lm)
##                  2.5 %     97.5 %
## (Intercept) 23.0498681 76.4923499
## Prewt        0.1128268  0.7560955
## TreatCont   -7.8754712 -0.3186599
## TreatFT      0.3060571  8.8200682

將繪圖指令改為for-loop

pict <- function(j){
  pic <- par(mfrow=c(2,2))
  for (j in 1:4)
  hist(bootstrap_estimates[,j])
  return(pic)
}
pict(1)

## $mfrow
## [1] 1 1
##