diff --git a/08-Intro-to-Bayes.Rmd b/08-Intro-to-Bayes.Rmd index 3e2fa7a3d..78736521a 100644 --- a/08-Intro-to-Bayes.Rmd +++ b/08-Intro-to-Bayes.Rmd @@ -1,4 +1,4 @@ -# (PART) 貝葉斯統計 {-} +# (PART) 貝葉斯統計 Introduction to Bayesian Statistics {-} # 貝葉斯統計入門 {#intro-Bayes} @@ -507,6 +507,12 @@ pi <- Vectorize(function(theta) dbeta(theta, 2,8)) curve(pi, xlab=~ theta, ylab="Density", main="Beta prior: a=2, b=8",frame=FALSE, lwd=2) ``` + +你會發現，beta分佈的圖形特徵由它的兩個超參數 (hyperparameter) `a, b` 決定。當相對地 `a` 比較大的時候，beta分佈的概率多傾向於較靠近橫軸右邊，也就是1的位置（較高概率），相對地 `b` 比較大的時候，beta分佈的曲線多傾向於靠近橫軸左邊，也就是0的位置（較低概率）。如果 `a` 和 `b`同時都在變大的話，beta分佈的曲線就變得比較“瘦”。這樣決定概率分佈形狀的參數，又被叫做**形狀參數 (shape parameters)** + + + + 我們定義 $a>0$ 時[伽馬方程](https://zh.wikipedia.org/wiki/%CE%93%E5%87%BD%E6%95%B0)爲 $$\Gamma(a)=\int_0^\infty x^{a-1}e^{-ax}\text{ d}x$$ diff --git a/12-Bayesian-stats.Rmd b/12-Bayesian-stats.Rmd index 23a5ef494..8efedd4b0 100644 --- a/12-Bayesian-stats.Rmd +++ b/12-Bayesian-stats.Rmd @@ -2685,13 +2685,13 @@ Dat <- list( # fit use rjags package -inits <- function (){ - list (beta0=rnorm(1), beta1=runif(1) ) -} +# inits <- function (){ +# list (beta0=rnorm(1), beta1=runif(1) ) +# } post.jags <- jags( data = Dat, - inits = inits, + # inits = inits, parameters.to.save = c("beta0", "beta1", "theta[6]"), n.iter = 5000, model.file = paste(bugpath, @@ -2866,7 +2866,7 @@ jagsModel <- jags.model( sep = ""), data = Dat, n.chains = 2, - inits = inits, + # inits = inits, quiet = TRUE) # Step 2 update 10000 iterations @@ -2904,7 +2904,7 @@ Dat <- list( ) post.jags <- jags( data = Dat, - inits = inits, + # inits = inits, parameters.to.save = c("beta0", "beta1", "theta[6]"), n.iter = 20000, model.file = paste(bugpath, @@ -2960,7 +2960,8 @@ model{ ``` -```{r OpenBUGSPractical0424, message=TRUE, warning=FALSE, echo=TRUE} +```{r OpenBUGSPractical0424, message=TRUE, warning=FALSE, echo=TRUE, eval=TRUE, cache=TRUE} + # R2JAGS codes: Dat <- list( y = c(1, 3, 6, 8, 11, 15, 17, 19), @@ -2968,9 +2969,10 @@ Dat <- list( x = c(30, 32, 34, 36, 38, 40, 42, 44), N = 8 ) + post.jags <- jags( data = Dat, - inits = inits, + # inits = inits, parameters.to.save = c("ED95", "OR", "P35", "beta0", "beta1"), n.iter = 20000, model.file = paste(bugpath, @@ -2986,6 +2988,8 @@ post.jags <- jags( print(post.jags) ``` + + 所以，藥物每增加劑量1mg，有療效的比值比是OR = 1.46 (95%CrI: 1.31, 1.63)。能夠達到95%患者都有療效的劑量是 45.02 mg (95% CrI: 42.7, 48.2 mg)。如果給予患者藥物劑量爲 35 mg，患者的疼痛能夠得到緩解(有療效)的概率是 32.3% (95% CrI: 22.7%, 42.3%)。跟着看到這裏的你是不是也覺得貝葉斯的結果和過程能夠更加豐富地回答我們想知道的問題呢？ ```{r R-OpenBUGSPractical0425, cache=FALSE, fig.width=3.5, fig.height=5, fig.cap='Posterior density plots of ED95, OR, and P35.', fig.align='center', out.width='80%', message=TRUE, warning=FALSE} @@ -3216,7 +3220,31 @@ list(alpha = 10, beta = 0, gamma = -5, logsigma2 = 5) 一開始我們先採集1000個事後概率分布樣本。模型中的四個未知參數$\alpha, \beta, \gamma, \sigma^2$的1000次事後MCMC採樣的歷史痕跡圖繪制如圖 \@ref(fig:BayesianChapter0505)。我們可以看到模型收斂的速度很快。刨除前50次採樣的圖 (Fig. \@ref(fig:BayesianChapter0506))，可以對採樣過程看得更加清楚。像圖 \@ref(fig:BayesianChapter0506) 這樣粗粗的有點像毛毛蟲一樣的歷史痕跡圖通常象徵已經達到理想的收斂。 -```{r BayesianChapter0504, cache=TRUE, message=TRUE, warning=FALSE, echo=TRUE} + + + + + + + + + + + + + + + + + + + + + + + + +```{r BayesianChapter0505, cache=FALSE, fig.width=7, fig.height=8, fig.cap='History plots for iterations 1-1000 for the Gambia example.', fig.align='center', out.width='80%', message=FALSE, warning=FALSE, echo=FALSE} library(BRugs) # Step 1 check model modelCheck(paste(bugpath, "/backupfiles/gambia-model.txt", sep = "")) @@ -3237,10 +3265,6 @@ samplesSet(parameters) # Generate 1000 iterations modelUpdate(1000) -``` - - -```{r BayesianChapter0505, cache=FALSE, fig.width=7, fig.height=8, fig.cap='History plots for iterations 1-1000 for the Gambia example.', fig.align='center', out.width='80%', message=FALSE, warning=FALSE, echo=FALSE} #### PUT THE SAMPLED VALUES IN ONE R DATA FRAME: chain <- data.frame(alpha = samplesSample("alpha"), beta = samplesSample("beta"), diff --git a/backupfiles/logistic-reg-model-centred-stat.txt b/backupfiles/logistic-reg-model-centred-stat.txt index 7a6c2fa92..d63ebdf85 100644 --- a/backupfiles/logistic-reg-model-centred-stat.txt +++ b/backupfiles/logistic-reg-model-centred-stat.txt @@ -7,8 +7,10 @@ model{ } # priors beta0 ~ dunif(-100, 100) - beta1 ~ dunif(-100, 100) + beta1 ~ dnorm(1, 0.000000001) + + # generated values OR <- exp(beta1) # odds ratio of positive response per 1 mg increase in dose - ED95 <- (logit(0.95) - beta0)/beta1 + mean(x[]) # dose that gives 95% of maximal response + ED95 <- mean(x[]) + ((logit(0.95) - beta0) / (beta1)) # dose that gives 95% of maximal response logit(P35) <- beta0 + beta1 * (35 - mean(x[])) -} \ No newline at end of file +} diff --git a/docs/ANOVA.html b/docs/ANOVA.html index da51124bb..004bdcb66 100644 --- a/docs/ANOVA.html +++ b/docs/ANOVA.html @@ -6,7 +6,7 @@ 第 28 章 方差分析 Introduction to Analysis of Variance | 醫學統計學 - + @@ -24,7 +24,7 @@ - + @@ -61,41 +61,66 @@ @@ -575,7 +600,7 @@
  • 38.1 定義
  • 39 The sandwich estimator
    • -
  • VII 貝葉斯統計
    • +
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
  • 40 貝葉斯統計入門
  • 39 The sandwich estimator
    • -
  • VII 貝葉斯統計
    • +
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
  • 40 貝葉斯統計入門
  • 39 The sandwich estimator
    • -
  • VII 貝葉斯統計
    • +
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
  • 40 貝葉斯統計入門 -

    帶方向的連接線 (directed links)被用來表示邏輯依賴關系的方向(natural ordering of the dependence association)。 其實這連接線所帶的方向本身就是回歸模型中表示預測變量和結果變量之間的依賴關系。貝葉斯統計學中常用的有向無環圖 (directed acyclic graphs, DAGs)，是我們常用的輔助建立貝葉斯模型的示意圖。

    +

    帶方向的連接線 (directed links)被用來表示邏輯依賴關系的方向(natural ordering of the dependence association)。 +其實這連接線所帶的方向本身就是回歸模型中表示預測變量和結果變量之間的依賴關系。貝葉斯統計學中常用的有向無環圖 (directed acyclic graphs, DAGs)，是我們常用的輔助建立貝葉斯模型的示意圖。

    例如一個描述吸煙，肥胖和心髒病發病可能的關系的的有向無環圖如下圖 84.2。途中帶方向的連線（箭頭）表示預測變量和結果變量之間的依賴關系。

    DAG for heart disease example @@ -1598,26 +1622,26 @@

    84.7 Practical Bayesian Statistic y ~ dbin(theta, 20)# sampling distribution for 20 observed patients #y <- 15 } -
    # Codes for R2JAGS
-
-
-Dat <- list(
-y = c(15)                    
-)
-
-# fit use R2jags package
-
-post.jags <- jags(
-  data = Dat,
-  parameters.to.save = c("theta"),
-  n.iter = 1100,
-  model.file = paste(bugpath, 
-                     "/backupfiles/MCdrugPractical04.txt", 
-                     sep = ""),
-  n.chains = 1,
-  n.thin = 1,
-  n.burnin = 100,
-  progress.bar = "none")
    + 
    ## Compiling model graph
 ##    Resolving undeclared variables
 ##    Allocating nodes
@@ -1627,7 +1651,7 @@ 
    84.7 Practical Bayesian Statistic
 ##    Total graph size: 4
 ## 
 ## Initializing model
    -
    print(post.jags)
    + 
    ## Inference for Bugs model at "/home/ccwang/Documents/LSHTMlearningnote/backupfiles/MCdrugPractical04.txt", fit using jags,
 ##  1 chains, each with 1100 iterations (first 100 discarded)
 ##  n.sims = 1000 iterations saved
@@ -1658,29 +1682,29 @@ 
    84.7 Practical Bayesian Statistic
 
      
 
    1. 嘗試使用JAGS/OpenBUGS來跑這個模型。記得你需要把數據加載到軟件裏面去。
      2. 
 
    
-
    # Codes for R2JAGS
-
-
-Dat <- list(
-     y = c(1, 3, 6, 8, 11, 15, 17, 19), 
-     n = c(20, 20, 20, 20, 20, 20, 20, 20), 
-     x = c(30, 32, 34, 36, 38, 40, 42, 44), 
-     N = 8 
-)
-
-# fit use R2jags package
-
-post.jags <- jags(
-  data = Dat,
-  parameters.to.save = c("beta0", "beta1", "theta[6]"),
-  n.iter = 1100,
-  model.file = paste(bugpath, 
-                     "/backupfiles/logistic-reg-model.txt", 
-                     sep = ""),
-  n.chains = 1,
-  n.thin = 1,
-  n.burnin = 100,
-  progress.bar = "none")
    
+
 
    ## Compiling model graph
 ##    Resolving undeclared variables
 ##    Allocating nodes
@@ -1690,7 +1714,7 @@ 
    84.7 Practical Bayesian Statistic
 ##    Total graph size: 53
 ## 
 ## Initializing model
    
-
    print(post.jags)
    
+
 
    ## Inference for Bugs model at "/home/ccwang/Documents/LSHTMlearningnote/backupfiles/logistic-reg-model.txt", fit using jags,
 ##  1 chains, each with 1100 iterations (first 100 discarded)
 ##  n.sims = 1000 iterations saved
@@ -1706,8 +1730,8 @@ 
    84.7 Practical Bayesian Statistic
 
      
 
    1. 對未知參數，也就是邏輯回歸的回歸系數 beta0, beta1，和第六次試驗的治療成功百分比 theta[6] 設定監測追蹤其採樣軌跡。試着執行1000次MCMC採樣，把這三個跟蹤的未知參數的採樣歷史軌跡繪制下來。
      2. 
 
    
-
    #### PLOT THE MCMC CHAINS:
-mcmcplots::traplot(post.jags, c("beta0", "beta1", "theta[6]"))
    
+
 
    
 History plot showing 1000 samples of beta0
 
    
@@ -1718,19 +1742,19 @@ 
    84.7 Practical Bayesian Statistic
 
      
 
    1. 接下來，我們給上面同一個邏輯回歸模型增加另一條MCMC採樣鏈，重復跑相同的模型1000次，繪制同樣是這三個未知參數的事後分布MCMC採樣歷史痕跡圖。
      2. 
 
    
-
    # fit use R2jags package
-
-post.jags <- jags(
-  data = Dat,
-  parameters.to.save = c("beta0", "beta1", "theta[6]"),
-  n.iter = 1000,
-  model.file = paste(bugpath, 
-                     "/backupfiles/logistic-reg-model.txt", 
-                     sep = ""),
-  n.chains = 2,
-  n.thin = 1,
-  n.burnin = 0,
-  progress.bar = "none")
    
+
 
    ## Compiling model graph
 ##    Resolving undeclared variables
 ##    Allocating nodes
@@ -1740,7 +1764,7 @@ 
    84.7 Practical Bayesian Statistic
 ##    Total graph size: 53
 ## 
 ## Initializing model
    
-
    mcmcplots::traplot(post.jags, c("beta0", "beta1", "theta[6]"))
    
+
 
    
 History plot showing 1000 samples of beta0, beta1, and theta[6]
 
    
@@ -1748,13 +1772,13 @@ 
    84.7 Practical Bayesian Statistic
 
    
 
    
 
    使用兩個不同的起始值作爲MCMC的採樣起點時，每個未知參數分別獲得兩條不同顏色的歷史痕跡圖。和之前只有一條MCMC採樣鏈時一樣，採樣的起始部分(大約100次左右)都有一些不穩定的值。等到100次採樣過後，我們發現每個參數的採樣結果都趨向於比較穩定，也就是方差，變化小了很多。但是藍色鏈，紅色鏈一直到200-300次採樣之後才逐漸互相交叉重疊。下面的圖把起初的500次採樣刨除了之後重新對每個未知參數的MCMC採樣繪制歷史痕跡。
    
-
    # samplesHistory("*", mfrow = c(3,1), beg = 501, ask = FALSE)
-Simulated <- coda::as.mcmc(post.jags)
-ggSample <- ggs(Simulated)
-
-ggSample %>% 
-  filter(Iteration >= 500 & Parameter %in% c("beta0", "beta1", "theta[6]")) %>% 
-  ggs_traceplot()
    
+
 
    
 History plot showing 1000 samples of beta0, beta1, and theta[6], iteration between 501-1000
 
    
@@ -1763,9 +1787,9 @@ 
    84.7 Practical Bayesian Statistic
 
    
 
    500-1000次之間的隨機採樣被放大了觀察之後，我們發現 beta0, beta1 的兩條MCMC鏈條的重疊程度並不理想，不像theta[6]那樣呈現令人滿意地重疊，兩條採樣鏈上下扭動，且在一些地方差異較大。
    
 
    繪制每個參數的MCMC鏈的自相關圖 (autocorrelation)，下面的圖中只繪制了500-1000次範圍內採樣的自相關圖。
    
-
    ggSample %>% 
-  filter(Iteration >= 500 & Parameter %in% c("beta0", "beta1", "theta[6]")) %>% 
-  ggs_autocorrelation()
    
+
 
    
 Autocorrelation plot for beta0
 
    
@@ -1776,42 +1800,42 @@ 
    84.7 Practical Bayesian Statistic
 
      
 
    1. 重新在OpenBUGS裏跑這個邏輯回歸模型，這一次，嘗試自己給 beta0, beta1 設置起始值。然後更新模型，採集MCMC鏈樣本10000次。這次嘗試同時使用繪制歷史痕跡圖（視覺檢查），和計算Brooks-Gelman-Rubin診斷參數及其圖示來判斷事後概率分布採樣是否達到收斂。你能判斷這個模型需要拋出掉多少一開始採集的樣本嗎(burn-in)？
      2. 
 
    
-
    Dat <- list(
-     y = c(1, 3, 6, 8, 11, 15, 17, 19), 
-     n = c(20, 20, 20, 20, 20, 20, 20, 20), 
-     x = c(30, 32, 34, 36, 38, 40, 42, 44), 
-     N = 8 
-)
-
-# fit use rjags package
-
-inits <- function (){
-  list (beta0=rnorm(1), beta1=runif(1) )
-}
-
-post.jags <- jags(
-  data = Dat,
-  inits = inits,
-  parameters.to.save = c("beta0", "beta1", "theta[6]"),
-  n.iter = 5000,
-  model.file = paste(bugpath,
-                     "/backupfiles/logistic-reg-model.txt",
-                     sep = ""),
-  n.chains = 2,
-  n.thin = 1,
-  n.burnin = 0,
-  progress.bar = "none")
-
-mcmcplots::traplot(post.jags, c("beta0", "beta1", "theta[6]"))
    
+
 
    
 History plots based on first 10000 iterations.
 
    
 圖 84.9: History plots based on first 10000 iterations.
 
    
 
    
-
    Simulated <- coda::as.mcmc(post.jags)
-# postsamples <- buildMCMC("*")
-gelman.diag(Simulated)
    
+
 
    ## Potential scale reduction factors:
 ## 
 ##          Point est. Upper C.I.
@@ -1823,7 +1847,7 @@ 
    84.7 Practical Bayesian Statistic
 ## Multivariate psrf
 ## 
 ## 2.08
    
-
    gelman.plot(Simulated)
    
+
 
    
 Brooks-Gelman-Rubin diagnostic graph
 
    
@@ -1834,29 +1858,29 @@ 
    84.7 Practical Bayesian Statistic
 
      
 
    1. 確認了你想要刨除的初始樣本量(burn-in)，繼續再進行MCMC採樣直到獲得1,000,000個事後概率分布樣本。此時你對獲得的事後概率分布樣本量提供的參數估計精確度滿意否？嘗試報告此時獲得的參數們的事後均值，及其95%可信區間。
      2. 
 
    
-
    Dat <- list(
-     y = c(1,   3,  6,  8, 11, 15, 17, 19), 
-     n = c(20, 20, 20, 20, 20, 20, 20, 20), 
-     x = c(30, 32, 34, 36, 38, 40, 42, 44), 
-     N = 8 
-)
-
-# inits <- function (){
-#   list (beta0=rnorm(1), beta1=runif(1) )
-# }
-
-post.jags <- jags(
-  data = Dat,
-  # inits = inits,
-  parameters.to.save = c("beta0", "beta1", "theta[6]"),
-  n.iter = 502000,
-  model.file = paste(bugpath,
-                     "/backupfiles/logistic-reg-model.txt",
-                     sep = ""),
-  n.chains = 2,
-  n.thin = 1,
-  n.burnin = 2000,
-  progress.bar = "none")
    
+
 
    ## Compiling model graph
 ##    Resolving undeclared variables
 ##    Allocating nodes
@@ -1934,36 +1958,36 @@ 
    84.7 Practical Bayesian Statistic
   beta0 ~ dunif(-100, 100)
   beta1 ~ dunif(-100, 100)
 }
    
-
    # JAGS codes:
-# 
-Dat <- list(
-     y = c(1, 3, 6, 8, 11, 15, 17, 19), 
-     n = c(20, 20, 20, 20, 20, 20, 20, 20), 
-     x = c(30, 32, 34, 36, 38, 40, 42, 44), 
-     N = 8 
-)
-
-# step 1 check model
-jagsModel <- jags.model(
-                     file = paste(bugpath, 
-                                  "/backupfiles/logistic-reg-model-centred.txt",
-                                  sep = ""),
-                     data = Dat,
-                     n.chains = 2, 
-                     inits = inits,
-                     quiet = TRUE)
-# Step 2 update 10000 iterations
-
-update(jagsModel, n.iter = 1, progress.bar = "none")
-
-# Step 3 set monitor variables
-
-codaSamples <- coda.samples(
-  jagsModel, variable.names = c("beta0", "beta1", "theta[6]"),
-  n.iter = 10000, progress.bar = "none"
-)
-
-mcmcplots::traplot(codaSamples, c("beta0", "beta1", "theta[6]"))
    
+
 
    
 History plots based on first 10000 iterations with centred covariates.
 
    
@@ -1971,7 +1995,7 @@ 
    84.7 Practical Bayesian Statistic
 
    
 
    
 
    用中心化之後的模型我們發現需要10000或更多的起始樣本來達到事後概率分布的收斂。BGR診斷圖 84.14 也提示我們大概在10000次之後的重復採樣才能達到收斂。所以我們決定要刨除起始10000次採集的樣本。(burn-in = 10000)
    
-
    gelman.diag(codaSamples)
    
+
 
    ## Potential scale reduction factors:
 ## 
 ##          Point est. Upper C.I.
@@ -1982,7 +2006,7 @@ 
    84.7 Practical Bayesian Statistic
 ## Multivariate psrf
 ## 
 ## 1
    
-
    gelman.plot(codaSamples)
    
+
 
    
 Brooks-Gelman-Rubin diagnostic graph
 
    
@@ -2004,10 +2028,10 @@ 
    84.7 Practical Bayesian Statistic
 ##  2 chains, each with 20000 iterations (first 10000 discarded)
 ##  n.sims = 20000 iterations saved
 ##          mu.vect sd.vect   2.5%    25%    50%    75%  97.5%  Rhat n.eff
-## beta0      0.002   0.206 -0.398 -0.139  0.001  0.140  0.412 1.001 20000
-## beta1      0.377   0.057  0.269  0.337  0.375  0.413  0.494 1.001 19000
-## theta[6]   0.753   0.049  0.654  0.721  0.755  0.787  0.844 1.001 13000
-## deviance  25.782   2.046 23.759 24.302 25.143 26.601 31.325 1.001 20000
+## beta0     -0.002   0.208 -0.404 -0.141 -0.002  0.136  0.408 1.002  2300
+## beta1      0.376   0.057  0.272  0.336  0.374  0.413  0.493 1.001 20000
+## theta[6]   0.752   0.049  0.651  0.719  0.754  0.786  0.842 1.001  4100
+## deviance  25.769   2.068 23.758 24.304 25.140 26.538 31.408 1.001 20000
 ## 
 ## For each parameter, n.eff is a crude measure of effective sample size,
 ## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
@@ -2042,46 +2066,47 @@ 
    84.7 Practical Bayesian Statistic
   ED95 <- (logit(0.95) - beta0)/beta1 + mean(x[]) # dose that gives 95% of maximal response
   logit(P35) <- beta0 + beta1 * (35 - mean(x[])) # probability of positive reaction if given 35mg drugs
 }
    
-
    # R2JAGS codes:
-Dat <- list(
-     y = c(1, 3, 6, 8, 11, 15, 17, 19), 
-     n = c(20, 20, 20, 20, 20, 20, 20, 20), 
-     x = c(30, 32, 34, 36, 38, 40, 42, 44), 
-     N = 8 
-)
-post.jags <- jags(
-  data = Dat,
-  inits = inits,
-  parameters.to.save = c("ED95", "OR", "P35", "beta0", "beta1"),
-  n.iter = 20000,
-  model.file = paste(bugpath,
-                     "/backupfiles/logistic-reg-model-centred-stat.txt",
-                     sep = ""),
-  n.chains = 2,
-  n.thin = 1,
-  n.burnin = 10000,
-  progress.bar = "none")
    
+
 
    ## Compiling model graph
 ##    Resolving undeclared variables
 ##    Allocating nodes
 ## Graph information:
 ##    Observed stochastic nodes: 8
 ##    Unobserved stochastic nodes: 2
-##    Total graph size: 74
+##    Total graph size: 76
 ## 
 ## Initializing model
    
-
    # mcmcplots::traplot(post.jags, c("beta0", "beta1", "theta[6]"))
-print(post.jags)
    
+
 
    ## Inference for Bugs model at "/home/ccwang/Documents/LSHTMlearningnote/backupfiles/logistic-reg-model-centred-stat.txt", fit using jags,
 ##  2 chains, each with 20000 iterations (first 10000 discarded)
 ##  n.sims = 20000 iterations saved
 ##          mu.vect sd.vect   2.5%    25%    50%    75%  97.5%  Rhat n.eff
-## ED95      45.037   1.376 42.754 44.081 44.895 45.833 48.134 1.001 14000
-## OR         1.458   0.083  1.310  1.400  1.453  1.509  1.636 1.001 20000
-## P35        0.322   0.050  0.227  0.288  0.322  0.356  0.422 1.001  3800
-## beta0     -0.002   0.202 -0.402 -0.138 -0.003  0.136  0.393 1.002  2500
-## beta1      0.375   0.057  0.270  0.337  0.373  0.412  0.492 1.001 20000
-## deviance  25.723   2.037 23.760 24.288 25.103 26.451 31.239 1.001 16000
+## ED95      45.021   1.374 42.700 44.073 44.883 45.800 48.139 1.001 20000
+## OR         1.459   0.083  1.309  1.402  1.453  1.510  1.636 1.001 12000
+## P35        0.322   0.050  0.227  0.288  0.321  0.355  0.422 1.001 12000
+## beta0     -0.002   0.203 -0.400 -0.137 -0.002  0.133  0.395 1.001 20000
+## beta1      0.376   0.057  0.269  0.338  0.374  0.412  0.492 1.001 12000
+## deviance  25.733   2.053 23.762 24.272 25.085 26.530 31.161 1.001 20000
 ## 
 ## For each parameter, n.eff is a crude measure of effective sample size,
 ## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
@@ -2090,24 +2115,24 @@ 
    84.7 Practical Bayesian Statistic
 ## pD = 2.1 and DIC = 27.8
 ## DIC is an estimate of expected predictive error (lower deviance is better).
    
 
    所以，藥物每增加劑量1mg，有療效的比值比是OR = 1.46 (95%CrI: 1.31, 1.63)。能夠達到95%患者都有療效的劑量是 45.02 mg (95% CrI: 42.7, 48.2 mg)。如果給予患者藥物劑量爲 35 mg，患者的疼痛能夠得到緩解(有療效)的概率是 32.3% (95% CrI: 22.7%, 42.3%)。跟着看到這裏的你是不是也覺得貝葉斯的結果和過程能夠更加豐富地回答我們想知道的問題呢？
    
-
    #### PLOT THE DENSITY and HISTOGRAMS OF THE SAMPLED VALUES
-par(mfrow=c(3,1))
-Simulated <- coda::as.mcmc(post.jags)
-ggSample <- ggs(Simulated)
-
-ED95 <- ggSample %>% 
-  filter(Parameter == "ED95")
-OR <- ggSample %>% 
-  filter(Parameter == "OR")
-P35 <- ggSample %>% 
-  filter(Parameter == "P35")
-
-plot(density(ED95$value), main = "ED95 sample", 
-     ylab = "P(ED95)", xlab = "Distribution of ED95", col = "red")
-plot(density(OR$value), main = "OR sample", 
-     ylab = "P(OR)", xlab = "Distribution of OR", col = "red")
-plot(density(P35$value), main = "P35 sample", 
-     ylab = "P(P35)", xlab = "Distribution of P35", col = "red")
    
+
 
    
 Posterior density plots of ED95, OR, and P35.
 
    
@@ -2183,7 +2208,7 @@ 
    References
    
   (function () {
     var script = document.createElement("script");
     script.type = "text/javascript";
-    var src = "";
+    var src = "true";
     if (src === "" || src === "true") src = "https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-MML-AM_CHTML";
     if (location.protocol !== "file:")
       if (/^https?:/.test(src))
diff --git a/docs/author.html b/docs/author.html
index 291b67325..be4492b41 100644
--- a/docs/author.html
+++ b/docs/author.html
@@ -6,7 +6,7 @@
   
   我是誰 | 醫學統計學
   
-  
+  
 
   
   
@@ -24,7 +24,7 @@
 
 
 
-
+
 
   
   
@@ -61,41 +61,66 @@
 
 
 
 
 
@@ -575,7 +600,7 @@
 
  • 38.1 定義
    •

  • 39 The sandwich estimator
    • -
  • VII 貝葉斯統計
    • +
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
  • 40 貝葉斯統計入門
  • 39 The sandwich estimator
    • -
  • VII 貝葉斯統計
    • +
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
  • 40 貝葉斯統計入門
  • 39 The sandwich estimator
    • -
  • VII 貝葉斯統計
    • +
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
  • 40 貝葉斯統計入門
  • 39 The sandwich estimator
    • -
  • VII 貝葉斯統計
    • +
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
  • 40 貝葉斯統計入門 -

    從 \(59\) 個號碼中隨機取出任意 \(6\) 個號碼的方法有 \(^{59}C_6\) 種。 \[^{59}C_6=\frac{59!}{6!(59-6)!}=45,057,474\]

    +

    從 \(59\) 個號碼中隨機取出任意 \(6\) 個號碼的方法有 \(^{59}C_6\) 種。 +\[^{59}C_6=\frac{59!}{6!(59-6)!}=45,057,474\]

    每次選取六個號碼做爲一組的可能性相同，所以，你買了一組樂透號碼，能中獎的概率就是 \(1/45,057,474 = 0.00000002219\)。你還會再去買彩票麼？

    5.3.1 如果我只想中其中的 \(3\) 個號碼，概率有多大？

    @@ -1433,7 +1457,7 @@

    5.3.1 如果我只想中其中的 (function () { var script = document.createElement("script"); script.type = "text/javascript"; - var src = ""; + var src = "true"; if (src === "" || src === "true") src = "https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-MML-AM_CHTML"; if (location.protocol !== "file:") if (/^https?:/.test(src)) diff --git a/docs/bookdown_files/figure-html/06-RobustStatistic-22-1.png b/docs/bookdown_files/figure-html/06-RobustStatistic-22-1.png index eb818482c..3a5d62085 100644 Binary files a/docs/bookdown_files/figure-html/06-RobustStatistic-22-1.png and b/docs/bookdown_files/figure-html/06-RobustStatistic-22-1.png differ diff --git a/docs/bookdown_files/figure-html/06-RobustStatistic-24-1.png b/docs/bookdown_files/figure-html/06-RobustStatistic-24-1.png index 1d7680988..05888a81c 100644 Binary files a/docs/bookdown_files/figure-html/06-RobustStatistic-24-1.png and b/docs/bookdown_files/figure-html/06-RobustStatistic-24-1.png differ diff --git a/docs/bookdown_files/figure-html/12-Bayesian-stats-6-1.png b/docs/bookdown_files/figure-html/12-Bayesian-stats-6-1.png index 4a552e931..67266531a 100644 Binary files a/docs/bookdown_files/figure-html/12-Bayesian-stats-6-1.png and b/docs/bookdown_files/figure-html/12-Bayesian-stats-6-1.png differ diff --git a/docs/bookdown_files/figure-html/12-Bayesian-stats-7-1.png b/docs/bookdown_files/figure-html/12-Bayesian-stats-7-1.png index 1156d23d2..bf4b3d730 100644 Binary files a/docs/bookdown_files/figure-html/12-Bayesian-stats-7-1.png and b/docs/bookdown_files/figure-html/12-Bayesian-stats-7-1.png differ diff --git a/docs/bookdown_files/figure-html/12-JAGS-stats-4-1.png b/docs/bookdown_files/figure-html/12-JAGS-stats-4-1.png index 475bf244a..279dc0586 100644 Binary files a/docs/bookdown_files/figure-html/12-JAGS-stats-4-1.png and b/docs/bookdown_files/figure-html/12-JAGS-stats-4-1.png differ diff --git a/docs/bookdown_files/figure-html/JAGS00-1.png b/docs/bookdown_files/figure-html/JAGS00-1.png index 1a5a12e8a..7d594f793 100644 Binary files a/docs/bookdown_files/figure-html/JAGS00-1.png and b/docs/bookdown_files/figure-html/JAGS00-1.png differ diff --git a/docs/bookdown_files/figure-html/JAGS01-1.png b/docs/bookdown_files/figure-html/JAGS01-1.png index 4427f9203..62d9829f5 100644 Binary files a/docs/bookdown_files/figure-html/JAGS01-1.png and b/docs/bookdown_files/figure-html/JAGS01-1.png differ diff --git a/docs/bookdown_files/figure-html/JAGS02-1.png b/docs/bookdown_files/figure-html/JAGS02-1.png index b38d4c655..72ee3f09d 100644 Binary files a/docs/bookdown_files/figure-html/JAGS02-1.png and b/docs/bookdown_files/figure-html/JAGS02-1.png differ diff --git a/docs/bookdown_files/figure-html/OpenBUGSPractical0407-1.png b/docs/bookdown_files/figure-html/OpenBUGSPractical0407-1.png index ee72f0ee9..a336c8b83 100644 Binary files a/docs/bookdown_files/figure-html/OpenBUGSPractical0407-1.png and b/docs/bookdown_files/figure-html/OpenBUGSPractical0407-1.png differ diff --git a/docs/bookdown_files/figure-html/OpenBUGSPractical040808-1.png b/docs/bookdown_files/figure-html/OpenBUGSPractical040808-1.png index ddbec52e3..043bc8e63 100644 Binary files a/docs/bookdown_files/figure-html/OpenBUGSPractical040808-1.png and b/docs/bookdown_files/figure-html/OpenBUGSPractical040808-1.png differ diff --git a/docs/bookdown_files/figure-html/OpenBUGSPractical0415-1.png b/docs/bookdown_files/figure-html/OpenBUGSPractical0415-1.png index 3575414a4..c06104d91 100644 Binary files a/docs/bookdown_files/figure-html/OpenBUGSPractical0415-1.png and b/docs/bookdown_files/figure-html/OpenBUGSPractical0415-1.png differ diff --git a/docs/bookdown_files/figure-html/OpenBUGSPractical0418-1.png b/docs/bookdown_files/figure-html/OpenBUGSPractical0418-1.png index 970dc3a46..639ed85cb 100644 Binary files a/docs/bookdown_files/figure-html/OpenBUGSPractical0418-1.png and b/docs/bookdown_files/figure-html/OpenBUGSPractical0418-1.png differ diff --git a/docs/bookdown_files/figure-html/OpenBUGSPractical0419-1.png b/docs/bookdown_files/figure-html/OpenBUGSPractical0419-1.png index a787abe5b..b0c500d6a 100644 Binary files a/docs/bookdown_files/figure-html/OpenBUGSPractical0419-1.png and b/docs/bookdown_files/figure-html/OpenBUGSPractical0419-1.png differ diff --git a/docs/bookdown_files/figure-html/OpenBUGSPractical0420-1.png b/docs/bookdown_files/figure-html/OpenBUGSPractical0420-1.png index c10b5c5a1..86a9acdd3 100644 Binary files a/docs/bookdown_files/figure-html/OpenBUGSPractical0420-1.png and b/docs/bookdown_files/figure-html/OpenBUGSPractical0420-1.png differ diff --git a/docs/bookdown_files/figure-html/OpenBUGSPractical0420-2.png b/docs/bookdown_files/figure-html/OpenBUGSPractical0420-2.png index e713058bf..bdb0e9b8a 100644 Binary files a/docs/bookdown_files/figure-html/OpenBUGSPractical0420-2.png and b/docs/bookdown_files/figure-html/OpenBUGSPractical0420-2.png differ diff --git a/docs/bookdown_files/figure-html/R-OpenBUGSPractical0425-1.png b/docs/bookdown_files/figure-html/R-OpenBUGSPractical0425-1.png index 1fffb803a..2a7afc38a 100644 Binary files a/docs/bookdown_files/figure-html/R-OpenBUGSPractical0425-1.png and b/docs/bookdown_files/figure-html/R-OpenBUGSPractical0425-1.png differ diff --git a/docs/bookdown_files/figure-html/propscore00-1.png b/docs/bookdown_files/figure-html/propscore00-1.png index a8d28dbe8..71b05f5e5 100644 Binary files a/docs/bookdown_files/figure-html/propscore00-1.png and b/docs/bookdown_files/figure-html/propscore00-1.png differ diff --git "a/docs/causal-languages-\345\233\240\346\236\234\346\216\250\346\226\267\347\232\204\350\252\236\346\263\225.html" "b/docs/causal-languages-\345\233\240\346\236\234\346\216\250\346\226\267\347\232\204\350\252\236\346\263\225.html" index bd45c16ba..77a222e25 100644 --- "a/docs/causal-languages-\345\233\240\346\236\234\346\216\250\346\226\267\347\232\204\350\252\236\346\263\225.html" +++ "b/docs/causal-languages-\345\233\240\346\236\234\346\216\250\346\226\267\347\232\204\350\252\236\346\263\225.html" @@ -6,7 +6,7 @@ 第 91 章 Causal Languages 因果推斷的語法 | 醫學統計學 - + @@ -24,7 +24,7 @@ - + @@ -61,41 +61,66 @@ @@ -575,7 +600,7 @@
  • 38.1 定義

  • 39 The sandwich estimator
    • -
  • VII 貝葉斯統計
    • +
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
  • 40 貝葉斯統計入門
  • 39 The sandwich estimator
    • -
  • VII 貝葉斯統計
    • +
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
  • 40 貝葉斯統計入門
  • 39 The sandwich estimator
    • -
  • VII 貝葉斯統計
    • +
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
  • 40 貝葉斯統計入門
  • 39 The sandwich estimator
    • -
  • VII 貝葉斯統計
    • +
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
  • 40 貝葉斯統計入門
  • 39 The sandwich estimator
    • -
  • VII 貝葉斯統計
    • +
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
  • 40 貝葉斯統計入門
  • 39 The sandwich estimator
    • -
  • VII 貝葉斯統計
    • +
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
  • 40 貝葉斯統計入門
  • 39 The sandwich estimator
    • -
  • VII 貝葉斯統計
    • +
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
  • 40 貝葉斯統計入門
  • 39 The sandwich estimator
    • -
  • VII 貝葉斯統計
    • +
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
  • 40 貝葉斯統計入門
  • 39 The sandwich estimator
    • -
  • VII 貝葉斯統計
    • +
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
  • 40 貝葉斯統計入門
  • 39 The sandwich estimator
    • -
  • VII 貝葉斯統計
    • +
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
  • 40 貝葉斯統計入門
  • 39 The sandwich estimator
    • -
  • VII 貝葉斯統計
    • +
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
  • 40 貝葉斯統計入門
  • 39 The sandwich estimator
    • -
  • VII 貝葉斯統計
    • +
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
  • 40 貝葉斯統計入門

    推薦書目：

      -
    1. “Principles of Statistical Inference” by D.R. Cox (D. Cox 2006)
      2. +
    2. “Principles of Statistical Inference” by D.R. Cox (Cox 2006)
    3. “Bayesian Data Analysis” by Gelman, Carlin, Stern, Dunson, Vehtari, and Rubin (Gelman et al. 2013), website for the book
    4. “Bayesian Biostatistics” by Vehtari and Rubin (Lesaffre and Lawson 2012)

    貝葉斯統計推斷，提供了不同於概率論推斷的另一種考察和解決問題的思路。所有的思考，都源於貝葉斯定理 Bayes’ Theorem (Section 2)。起源於英國統計學家托馬斯貝葉斯 (Thomas Bayes) 死後被好友 Richard Price 整理發表的論文: “An essay towards solving a problem in the doctrine of chances.”

    概率論推斷與貝葉斯推斷的中心都圍繞似然 likelihood (Section 12) 的概念。然而二者對似然提供的信息之理解和解釋完全不同。即在對於觀察數據提供的信息的理解，和如何應用已有信息來影響未來決策（或提供預測）的問題上常常被認爲是統計學中形成鮮明對比的兩種哲學理念。過去幾個世紀二者之間孰優孰劣的爭論相當激烈。但是，從實際應用的角度來看，我們目前更關心哪種思維能更加實用地描述和模擬真實世界。幸運地是，多數情況下，二者的差距不大。所以無法簡單地從一個實驗或者一次爭論中得出誰更出色的結論。現在的統計學家們通常不再如同信仰之爭那樣的互相水火不容，而是從實用性角度來判斷一些實際情況下，採用哪種思想能使計算過程更加簡便或者計算結果更加接近真實情況。

    -

    請思考如下的問題： 什麼是概率？ What is probability?

    +

    請思考如下的問題： +什麼是概率？ What is probability?

    1. 概率論思想下的定義：某事件在多次重複觀察實驗結果中發生次數所佔的比例。
      The probability of an event is the limit of its relative frequency in a large number of trials."
      2. -
    2. 貝葉斯思想下的定義：概率是你相信某事件會發生的可能性。
      Probability is a measure of the degree of belief about an event.
      3. +
    3. 貝葉斯思想下的定義：概率是你相信某事件會發生的可能性。
      Probability is a measure of the degree of belief +about an event.

    40.1 概率論推斷的複習

    @@ -1400,7 +1425,7 @@

    40.2.1 演繹推理 deductive rea

    40.2.2 如何給可能性定量 Quantifying plausibility

    -

    進行可能性定量之前，R.T. Cox 制定了如下的規則(R. T. Cox 1946)：

    +

    進行可能性定量之前，R.T. Cox 制定了如下的規則(Cox 1946)：

    1. \(\text{plausibility}(A)\) 是一個有邊界的實數；
    2. 傳遞性，transitivity：如果 @@ -1442,7 +1467,8 @@

      40.2.2 如何給可能性定量 Q

    40.3 貝葉斯推理的統計學實現

    -

    在經典概率論中，概率分佈的標記 \(f_X(x;\theta)\) 的涵義爲： 對於一個隨機變量 \(X\)，它在我們假設的某種固定的真實（上帝才知道是多少的）參數 \(\theta\) 的分佈框架下，不斷重複相同的實驗之後獲得的概率分佈。

    +

    在經典概率論中，概率分佈的標記 \(f_X(x;\theta)\) 的涵義爲： +對於一個隨機變量 \(X\)，它在我們假設的某種固定的真實（上帝才知道是多少的）參數 \(\theta\) 的分佈框架下，不斷重複相同的實驗之後獲得的概率分佈。

    在貝葉斯統計推理中，一切都被看作是一個服從概率分佈的隨機變量。利用貝葉斯定理，我們將先驗隨機概率分佈 (prior probability distribution)，和觀察數據作條件概率 (condition on the observed data)，從而獲得事後概率分佈 (posterior probability distribution)。

    40.3.1 醫學診斷測試 diagnostic testing

    @@ -1507,7 +1533,7 @@

    40.3.3 說點小歷史

    圖 40.1: Sir Ronald Fisher

    -

    Ronald Aylmer Fisher (1890-1962) 推動了統計學在20世紀前半頁的重大發展。他鞏固了概率論統計學堅實的基礎，並且積極提倡這一套理論(R. A. Fisher 1922)。但是 Fisher 本人對於統計學的“統計學意義, level of significance” 的認識卻是隨着時間和他年齡的變化而變化的：

    +

    Ronald Aylmer Fisher (1890-1962) 推動了統計學在20世紀前半頁的重大發展。他鞏固了概率論統計學堅實的基礎，並且積極提倡這一套理論(Fisher 1922)。但是 Fisher 本人對於統計學的“統計學意義, level of significance” 的認識卻是隨着時間和他年齡的變化而變化的：

    @@ -1578,7 +1607,8 @@

    40.4 練習題

    \end{aligned} \]

    其中分子部分 \(\frac{f(\theta,x)}{\pi(\theta)}\) 就是條件概率 \(f(x|\theta)\)。

    -

    分母的 \(f(x)\) 部分 \[ +

    分母的 \(f(x)\) 部分 +\[ \begin{aligned} f(x) &= \int f(x,\theta) \text{d}\theta \\ &= \int \frac{f(x,\theta)}{\pi(\theta)} \cdot \pi(\theta) \text{d}\theta \\ @@ -1705,25 +1735,25 @@

    40.4 練習題

    References

    -

    Cox, D.R. 2006. Principles of Statistical Inference. Cambridge University Press. https://books.google.co.jp/books?id=nRgtGZXi2KkC.

    +

    Cox, D.R. 2006. Principles of Statistical Inference. Cambridge University Press. https://books.google.co.jp/books?id=nRgtGZXi2KkC.

    -

    Cox, R. T. 1946. “Probability, Frequency and Reasonable Expectation.” American Journal of Physics 14 (1): 1–13. doi:10.1119/1.1990764.

    +

    Cox, R. T. 1946. “Probability, Frequency and Reasonable Expectation.” American Journal of Physics 14 (1): 1–13. https://doi.org/10.1119/1.1990764.

    -

    Fisher, R. A. 1922. “On the Mathematical Foundations of Theoretical Statistics.” Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character 222. The Royal Society: 309–68. http://www.jstor.org/stable/91208.

    +

    Fisher, R. A. 1922. “On the Mathematical Foundations of Theoretical Statistics.” Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character 222. The Royal Society: 309–68. http://www.jstor.org/stable/91208.

    -

    Gelman, A., J.B. Carlin, H.S. Stern, D.B. Dunson, A. Vehtari, and D.B. Rubin. 2013. Bayesian Data Analysis, Third Edition. Chapman & Hall/Crc Texts in Statistical Science. Taylor & Francis. https://books.google.co.uk/books?id=ZXL6AQAAQBAJ.

    +

    Gelman, A., J.B. Carlin, H.S. Stern, D.B. Dunson, A. Vehtari, and D.B. Rubin. 2013. Bayesian Data Analysis, Third Edition. Chapman & Hall/Crc Texts in Statistical Science. Taylor & Francis. https://books.google.co.uk/books?id=ZXL6AQAAQBAJ.

    -

    Lesaffre, E., and A.B. Lawson. 2012. Bayesian Biostatistics. Statistics in Practice. Wiley. https://books.google.co.uk/books?id=WV7KVjEQnJMC.

    +

    Lesaffre, E., and A.B. Lawson. 2012. Bayesian Biostatistics. Statistics in Practice. Wiley. https://books.google.co.uk/books?id=WV7KVjEQnJMC.

      -

    1. http://www.senns.demon.co.uk/wdict.html↩

      2. +

    2. http://www.senns.demon.co.uk/wdict.html↩

    @@ -1787,7 +1817,7 @@

    References

    (function () { var script = document.createElement("script"); script.type = "text/javascript"; - var src = ""; + var src = "true"; if (src === "" || src === "true") src = "https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-MML-AM_CHTML"; if (location.protocol !== "file:") if (/^https?:/.test(src)) diff --git a/docs/intro.html b/docs/intro.html index 1885698dd..45662544d 100644 --- a/docs/intro.html +++ b/docs/intro.html @@ -6,7 +6,7 @@ 第 1 章 概率論入門：定義與公理 | 醫學統計學 - + @@ -24,7 +24,7 @@ - + @@ -61,41 +61,66 @@ @@ -575,7 +600,7 @@
  • 38.1 定義
  • 39 The sandwich estimator
    • -
  • VII 貝葉斯統計
    • +
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
  • 40 貝葉斯統計入門
  • 39 The sandwich estimator
    • -
  • VII 貝葉斯統計
    • +
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
  • 40 貝葉斯統計入門
  • 39 The sandwich estimator
    • -
  • VII 貝葉斯統計
    • +
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
  • 40 貝葉斯統計入門
  • 39 The sandwich estimator
    • -
  • VII 貝葉斯統計
    • +
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
  • 40 貝葉斯統計入門
    • 表 40.1: Fisher’s interpretation of ‘level of significance’ and the Neyman-Pearson interpretation @@ -1528,13 +1554,16 @@

    40.3.3 說點小歷史
    -統計學有意義的水平（傳統上使用 \(\alpha=5\%\)），必須在實施統計檢驗之前就被決定。因此，統計學意義的水平是相應統計學檢驗本身的性質之一。
    Thus, the level of significance is a property of the test. +統計學有意義的水平（傳統上使用 \(\alpha=5\%\)），必須在實施統計檢驗之前就被決定。因此，統計學意義的水平是相應統計學檢驗本身的性質之一。
    +Thus, the level of significance is a property of the test.     		-統計學意義的水平，應該被精確計算並且在報告中明確 \(p\) 值的大小，故統計學意義的水平本身是在實施了統計檢驗之後計算的。它應該是屬於觀察數據的固有性質。
    Here the level of significance is a property of the data. +統計學意義的水平，應該被精確計算並且在報告中明確 \(p\) 值的大小，故統計學意義的水平本身是在實施了統計檢驗之後計算的。它應該是屬於觀察數據的固有性質。
    +Here the level of significance is a property of the data.     		-\(\alpha\) 和 \(\beta\) 作爲統計檢驗的第一類錯誤和第二類錯誤指標，應該在實施統計檢驗之前被決定。所以 \(\alpha, \beta\) 是屬於統計檢驗的性質。
    Yet, to determine \(\alpha, \beta\) no convention is required, but rather a cost-benefit estimation of the severity of the two kinds of error. +\(\alpha\) 和 \(\beta\) 作爲統計檢驗的第一類錯誤和第二類錯誤指標，應該在實施統計檢驗之前被決定。所以 \(\alpha, \beta\) 是屬於統計檢驗的性質。
    +Yet, to determine \(\alpha, \beta\) no convention is required, but rather a cost-benefit estimation of the severity of the two kinds of error.
    @@ -2333,7 +2358,7 @@

    13.5.3 Q3

    (function () { var script = document.createElement("script"); script.type = "text/javascript"; - var src = ""; + var src = "true"; if (src === "" || src === "true") src = "https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-MML-AM_CHTML"; if (location.protocol !== "file:") if (/^https?:/.test(src)) diff --git a/docs/lm.html b/docs/lm.html index 882b1b5c5..82f70de78 100644 --- a/docs/lm.html +++ b/docs/lm.html @@ -6,7 +6,7 @@ 第 26 章 簡單線性迴歸 Simple Linear Regression | 醫學統計學 - + @@ -24,7 +24,7 @@ - + @@ -61,41 +61,66 @@ @@ -575,7 +600,7 @@
  • 38.1 定義
  • 39 The sandwich estimator
    • -
  • VII 貝葉斯統計
    • +
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
  • 40 貝葉斯統計入門
    • @@ -1793,9 +1816,9 @@

    26.7 R 演示 例 2： 表

    擬合簡單線性迴歸也是小菜一碟：

    -
    wk_age <- lm(Age ~ Group, data=Walk)
-
-summary(wk_age)
    +     
    ## 
 ## Call:
 ## lm(formula = Age ~ Group, data = Walk)
@@ -1815,7 +1838,7 @@ 
    26.7 R 演示 例 2： 表
    -
    anova(wk_age)
    +     
    ## Analysis of Variance Table
 ## 
 ## Response: Age
@@ -1827,7 +1850,7 @@ 
    26.7 R 演示 例 2： 表這裏的 \(\hat\alpha=10.125\)，意爲參照組 (此處，“exercise” 被默認設定爲參照組，而 “control” 被默認拿來和參照組相比較) 的兒童也就是，積極練習走路的小朋友這組能夠獨立行走的平均年齡是 \(10.125\) 個月。
    
 
    \(\hat\beta=2.225\)，意爲和參照組 (積極練習組) 相比，對照組兒童能夠自己行走的年齡平均要晚 \(2.225\) 個月。所以對照組兒童能夠直立行走的平均年齡就是 \(10.125+2.225=12.35\) 個月。
    
 
    上述結果，你如果拿來和下面的兩樣本 \(t\) 檢驗的結果相比就知道，是完全一致的。其中統計量 \(t^2=2.9285^2=F_{1,9}=8.58\)。
    
-
    t.test(Age~Group, data=Walk, var.equal=TRUE)
    
+    
 
    ## 
 ##  Two Sample t-test
 ## 
@@ -1843,9 +1866,9 @@ 
    26.7 R 演示 例 2： 表
 
    26.8 練習
    
 
    使用的數據內容爲：兩次調查同一樣本，99 名健康男性的血清膽固醇水平，間隔一年。
    
-
    # 數據讀入
-Chol <- read_dta("backupfiles/chol.dta")
-summary(Chol)
    
+    
 
    ##        id           chol1         chol2    
 ##  Min.   : 1.0   Min.   :152   Min.   :170  
 ##  1st Qu.:25.5   1st Qu.:235   1st Qu.:240  
@@ -1853,8 +1876,8 @@ 
    26.8 練習
    
 ##  Mean   :50.0   Mean   :265   Mean   :264  
 ##  3rd Qu.:74.5   3rd Qu.:290   3rd Qu.:290  
 ##  Max.   :99.0   Max.   :360   Max.   :355
    
-
    # Alternative Descriptive Statistics using psych package
-describe(Chol)
    
+
 
    ## Chol 
 ## 
 ##  3  Variables      99  Observations
@@ -1883,16 +1906,16 @@ 
    26.8 練習
    
 ## 
 ## lowest : 170 190 195 200 205, highest: 320 330 345 350 355
 ## ----------------------------------------------------------------------------------------------------
    
-
    # 兩次膽固醇水平的直方圖 Distribution of the two measures
-par(mfrow=c(1,2))
-hist(Chol$chol1)
-hist(Chol$chol2)
    
+
 
    
-
    # 對兩次膽固醇水平作散點圖
-ggplot(Chol, aes(x=chol1, y=chol2)) + geom_point(shape=20) +
-  scale_x_continuous(breaks=seq(150, 400, 50),limits = c(150, 355))+
-  scale_y_continuous(breaks=seq(150, 400, 50),limits = c(150, 355)) +
-   theme_stata() +labs(x = "Cholesterol at visit 1 (mg/100ml)", y = "Cholesterol at visit 2 (mg/100ml)")
    
+
 
    
 
    
 
    26.8.1 兩次測量的膽固醇水平分別用 \(C_1, C_2\) 來標記的話，考慮這樣的簡單線性迴歸模型：\(C_2=\alpha+\beta C_2 + \varepsilon\)。我們進行這樣迴歸的前提假設有哪些？
    
@@ -1902,14 +1925,14 @@ 
    26.8.1 兩次測量的膽固醇
 
  • 殘差值，在每一個給定的 \(C_1\) 值處呈現正態分佈，且方差不變。
    • 
 
 
    從散點圖來看這些假設應該都能得到滿足。
    
-
    # 計算兩次膽固醇水平的 均值，方差，以及二者的協方差
-mean(Chol$chol1); mean(Chol$chol2)
    
+
 
    ## [1] 264.6
    
 
    ## [1] 263.5
    
-
    var(Chol$chol1); var(Chol$chol2)
    
+
 
    ## [1] 1661
    
 
    ## [1] 1457
    
-
    cov(Chol$chol1, Chol$chol2)
    
+
 
    ## [1] 961.2
    
 
    
 
    
@@ -1921,15 +1944,15 @@ 
    26.8.2 計算普通最小二乘
 &=\frac{1661.061}{961.224}=0.578
 \end{aligned}
 \]
    
-
    cov(Chol$chol1, Chol$chol2)/var(Chol$chol1)
    
+
 
    ## [1] 0.5787
    
 
    \[\hat\alpha=\bar{y}-\hat\beta\bar{x}=263.54-0.578\times264.59=110.425\]
    
-
    mean(Chol$chol2)-mean(Chol$chol1)*cov(Chol$chol1, Chol$chol2)/var(Chol$chol1)
    
+
 
    ## [1] 110.4
    
 
    
 
    
 
    26.8.3 和迴歸模型計算的結果作比較，解釋這些估計值的含義
    
-
    summary(lm(chol2~chol1, data=Chol))
    
+
 
    ## 
 ## Call:
 ## lm(formula = chol2 ~ chol1, data = Chol)
@@ -1955,23 +1978,24 @@ 
    26.8.3 和迴歸模型計算的
 
    
 
    
 
    26.8.4 加上計算的估計值直線 (即迴歸直線)
    
-
    ggplot(Chol, aes(x=chol1, y=chol2)) + geom_point(shape=20, colour="grey40") +
-  stat_smooth(method = lm, se=FALSE, size=0.5) +
-   scale_x_continuous(breaks=seq(150, 400, 50),limits = c(150, 355))+
-  scale_y_continuous(breaks=seq(150, 400, 50),limits = c(150, 355)) +
-   theme_stata() +labs(x = "Cholesterol at visit 1 (mg/100ml)", y = "Cholesterol at visit 2 (mg/100ml)")
    
+
 
    
-
    可以注意到，第一次訪問時膽固醇水平高的人，第二次被測量時膽固醇值高於平均值，但是卻沒有第一次高出平均值的部分多。 相似的，第一次膽固醇水平低的人，第二次膽固醇水平低於平均值，但是卻沒有第一次低於平均值的部分多。這一現象被叫做 “向均數迴歸-regression to the mean”
    
+
    可以注意到，第一次訪問時膽固醇水平高的人，第二次被測量時膽固醇值高於平均值，但是卻沒有第一次高出平均值的部分多。
+相似的，第一次膽固醇水平低的人，第二次膽固醇水平低於平均值，但是卻沒有第一次低於平均值的部分多。這一現象被叫做 “向均數迴歸-regression to the mean”
    
 
    
 
    
 
    26.8.5 下面的代碼用於模型的假設診斷
    
-
    M <- lm(chol2~chol1, data=Chol)
-par(mfrow = c(2, 2))  # Split the plotting panel into a 2 x 2 grid
-plot(M)
    
+
 
    
 
    好心人在 github 上共享了 Check_assumption.R 的代碼，可以使用 ggplot2 來獲取高逼格的模型診斷圖：
    
-
    source("checkassumptions.R")
-check_assumptions(M)
    
+
 
    
 
    
 
@@ -2037,7 +2061,7 @@ 
    26.8.5 下面的代碼用於模
   (function () {
     var script = document.createElement("script");
     script.type = "text/javascript";
-    var src = "";
+    var src = "true";
     if (src === "" || src === "true") src = "https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-MML-AM_CHTML";
     if (location.protocol !== "file:")
       if (/^https?:/.test(src))
diff --git a/docs/missing-data-1.html b/docs/missing-data-1.html
index a9d9d911c..1bbd3cc3e 100644
--- a/docs/missing-data-1.html
+++ b/docs/missing-data-1.html
@@ -6,7 +6,7 @@
   
   第 67 章 Missing data 1 | 醫學統計學
   
-  
+  
 
   
   
@@ -24,7 +24,7 @@
 
 
 
-
+
 
   
   
@@ -61,41 +61,66 @@
 
 
 
 
 
@@ -575,7 +600,7 @@
 
  • 38.1 定義
    • 
 
 
  • 39 The sandwich estimator
    • 
-
  • VII 貝葉斯統計
    • 
+
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
    • 
 
  • 40 貝葉斯統計入門
    • 
 
  • 39 The sandwich estimator
    • 
-
  • VII 貝葉斯統計
    • 
+
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
    • 
 
  • 40 貝葉斯統計入門
    • 
 
  • 39 The sandwich estimator
    • 
-
  • VII 貝葉斯統計
    • 
+
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
    • 
 
  • 40 貝葉斯統計入門
    • 
 
  • 39 The sandwich estimator
    • 
-
  • VII 貝葉斯統計
    • 
+
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
    • 
 
  • 40 貝葉斯統計入門
    • 
 
  • 39 The sandwich estimator
    • 
-
  • VII 貝葉斯統計
    • 
+
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
    • 
 
  • 40 貝葉斯統計入門
    • 
 
    根據方差協方差矩陣進行的主成分分析結果，我們發現主成分 6 和 7 可以忽略不計。相同的計算結果可以在R裏面通過方便的計算包 FactoMineR 來計算並用 factoextra 來實現其分析圖形的美觀展示：
    
-
    # library(FactoMineR)
-org.pca <- PCA(orange[, 1:7], ncp = 7, graph = FALSE)
-
-# library(factoextra)
-eig.val <- get_eigenvalue(org.pca)
-eig.val # eigenvalue (variances of each principal components)
    
+
 
    ##        eigenvalue variance.percent cumulative.variance.percent
 ## Dim.1 4.743692688      67.76703840                   67.767038
 ## Dim.2 1.333289855      19.04699793                   86.814036
 ## Dim.3 0.819841150      11.71201643                   98.526053
 ## Dim.4 0.084023297       1.20033282                   99.726386
 ## Dim.5 0.019153009       0.27361442                  100.000000
    
-
    # eigen vectors:
-org.pca$svd$V
    
+
 
    ##             [,1]        [,2]         [,3]         [,4]         [,5]
 ## [1,]  0.21100074  0.65340689 -0.517409852  0.028573070  0.030958154
 ## [2,]  0.45241413  0.11618305 -0.064606287  0.266760192  0.295222955
@@ -1764,7 +1789,8 @@ 
    68.4 主成分分析數據實例<
 y_5 & = 0.0310x_1 + 0.2952x_2 - 0.2250x_3 + 0.3456x_4 - 0.4106x_5 + 0.6712x_6 + 0.3503x_7 \\
 \end{aligned}
 \]
    
-
    於是，解釋完了如何從原始數據變量根據計算獲得的特徵值向量轉換成爲新的變量之後，要面對的問題是，我們要保留多少主成分？ 我們通常會使用圖 68.6 那樣的碎石圖 (Scree plot) 來輔助判斷。碎石圖通常縱軸是每個主成分能夠解釋的數據總體方差的百分比，然後橫軸是主成分的個數。所以我們會期待出現一個像手肘一樣的形狀提示應該在第幾個主成分的地方停下。通常在統計分析中，我們默認的準則是，至少保留的主成分個數要能夠解釋總體方差的 70%/80% 以上才較爲理想。Kaiser 準則 建議的是，最好保留下特徵值大於等於1(也就是標準化數據之後獲得的主成分變量方差大於等於1)的主成分變量。在我們的橙汁數據實例中，顯然保留前兩個主成分就已經能夠解釋 86.81% 的總體方差，我們認爲這是理想的主成分個數。
    
+
    於是，解釋完了如何從原始數據變量根據計算獲得的特徵值向量轉換成爲新的變量之後，要面對的問題是，我們要保留多少主成分？
+我們通常會使用圖 68.6 那樣的碎石圖 (Scree plot) 來輔助判斷。碎石圖通常縱軸是每個主成分能夠解釋的數據總體方差的百分比，然後橫軸是主成分的個數。所以我們會期待出現一個像手肘一樣的形狀提示應該在第幾個主成分的地方停下。通常在統計分析中，我們默認的準則是，至少保留的主成分個數要能夠解釋總體方差的 70%/80% 以上才較爲理想。Kaiser 準則 建議的是，最好保留下特徵值大於等於1(也就是標準化數據之後獲得的主成分變量方差大於等於1)的主成分變量。在我們的橙汁數據實例中，顯然保留前兩個主成分就已經能夠解釋 86.81% 的總體方差，我們認爲這是理想的主成分個數。
    
 
    
 Orange data: eigenvalues among all variances (varaince explained) by each dimension (principle component) provided by PCA
 
    
@@ -1772,8 +1798,8 @@ 
    68.4 主成分分析數據實例<
 
    
 
    
 
    另外一種輔助的圖形是叫做分數圖 (score plot)，又名個人圖 (graph of individuals)，如果個體的變量特徵相近，他們會在圖中聚在較爲靠近的地方：
    
-
    fviz_pca_ind(org.pca, pointsize = "cos2", pointshape = 21, 
-             fill = "#E7B800", repel = TRUE, labelsize = 2) 
    
+
 
    
 Score plot/individual plot of the orange data.
 
    
@@ -1782,7 +1808,7 @@ 
    68.4 主成分分析數據實例<
 
    
 
    細心觀察的話，你會發現圖 68.7 中各個橙汁 (個體,individual) 的座標其實是來自於PCA分析結果中第一和第二主成分變量的結果，展示在第一和第二主成分變量構成的平面。該平面解釋了總體數據慣性 (inertia) 的 86.82% (= 67.77% + 19.05%)。其中第一個主成分 Dim.1 把 Tropicana fr. 和 Pampryl amb. 兩種橙汁分別歸類在最右邊和最左邊。這是因爲原始數據中 Tropicana fr. 是 Odour.typicality 得分最高 Bitternes 得分倒數第二低，同時 Pampryl amb. 則是在這兩個項目上得分分別是最低和最高。也就是說這兩種橙汁在這兩個項目上得分分別是左右兩種極端，所以首先在第一主成分中把這兩中橙汁分離開來。接下來，第二主成分變量 Dim.2 則是將第一主成分成功分離開的兩個個體(橙汁)從數據中拿掉以後，剩下的四種橙汁的分類。可以看到第二個主成分軸，把 Pampryl fr. 和 Tropicana amb. 兩種橙汁放在了該軸的兩個極端，這是因爲 Pampryl fr. 在 Intensity.of.taste 項目上得分最高，而 Tropicana amb. 在拿掉了第一主成分分離的兩種橙汁之後，在 Odour.intensity 項目上得分最低。
    
 
    回到 R 幫忙分析的主成分結果報告來：
    
-
    summary(org.pca)
    
+
 
    ## 
 ## Call:
 ## PCA(X = orange[, 1:7], ncp = 7, graph = FALSE) 
@@ -1844,7 +1870,7 @@ 
    68.4 主成分分析數據實例<
 Bitterness         | -0.935 18.420  0.874 |  0.188  2.651  0.035 | -0.285  9.936  0.081 |
 Sweetness          |  0.955 19.220  0.912 | -0.159  1.889  0.025 |  0.187  4.246  0.035 |
    
 
    根據這個結果繪製的變量相關關係圖如下：
    
-
    fviz_pca_var(org.pca, repel = TRUE, labelsize = 2) 
    
+
 
    
 Variable plot of the orange data.
 
    
@@ -1857,8 +1883,8 @@ 
    68.4 主成分分析數據實例<
 
  • 每個變量從原點出發時的箭頭長度越長 cos2，代表它在該主成分軸上代表質量更好 (the quality of representation of the variable on the component)
    • 
 
 
    如果你願意，且數據和變量不至於多到眼花繚亂，我們還可以把個體圖和變量圖結合起來觀察：
    
-
    fviz_pca_biplot(org.pca, repel = TRUE, pointsize = "cos2", pointshape = 21, 
-             labelsize = 2) 
    
+
 
    
 Biplot of the orange data.
 
    
@@ -2101,9 +2127,9 @@ 
    68.5 在PCA圖形中加入補充

    我們可以把這些沒有用於計算主成分分析的變量 (active variables)，和其餘的輔助性變量 (supplementary variables) 同時繪製在變量相關係數圓盤圖中。此時我們只需要在進行PCA運算的時候告訴R這些變量是輔助性的連續/分類變量即可：

    -
    org.pca <- PCA(orange, quanti.sup = 8:14, quali.sup = 15:16,
-               graph = FALSE)
-org.pca$quanti.sup
    + 
    ## $coord
 ##                         Dim.1       Dim.2         Dim.3        Dim.4       Dim.5
 ## Glucose          -0.572454497  0.31123036  0.0263849025 -0.208332016 -0.72892600
@@ -2134,7 +2160,7 @@ 
    68.5 在PCA圖形中加入補充
 ## Citric.acid      0.5466683910 0.014786677 0.038314802810 0.0776568810 0.322573248
 ## Vitamin.C        0.0019870119 0.100477991 0.064778491512 0.8191451868 0.013611319

    然後用下面的代碼繪製包含了輔助性變量的變量相關圓盤圖：

    -
    fviz_pca_var(org.pca, repel = TRUE)
    +
    Orange juice data: representation of the active and supplementary variables (in blue).

    @@ -2146,9 +2172,9 @@

    68.5 在PCA圖形中加入補充

    68.5.1 展示分類輔助性變量和個體的關係

    根據不同的儲存方式，兩類的橙汁區別很清楚。

    -
    p <- fviz_pca_ind(org.pca, habillage = 15, 
-             palette = "jco", repel = TRUE)
-p
    +
    Plane representation of the scatterplot of individuals with a supplementary categorical variable (way of preserving).

    @@ -2156,9 +2182,9 @@

    68.5.1 展示分類輔助性變

    根據橙子的產地區分繪製的個人圖：

    -
    p <- fviz_pca_ind(org.pca, habillage = 16, 
-             palette = "jco", repel = TRUE)
-p
    +
    Plane representation of the scatterplot of individuals with a supplementary categorical variable (origin).

    @@ -2183,9 +2209,9 @@

    68.6.1 使用的數據和簡單
    1. 在R裏讀入你的數據，看看這4種生物標幟物的簡單統計量和分佈，它們用的是相同的測量單位嗎？
    -
    plant <- read_dta("backupfiles/plant.dta")
-plant <- plant[, 1:4]
-head(plant)
    + 
    ## # A tibble: 6 x 4
 ##       bm1     bm2     bm3     bm4
 ##     <dbl>   <dbl>   <dbl>   <dbl>
@@ -2195,7 +2221,7 @@ 
    68.6.1 使用的數據和簡單
 ## 4 106      10      44.6    57.6  
 ## 5 141.    122     115.    123.   
 ## 6   0.5     0.800   0.200   0.5
    -
    summ(plant)
    + 
    ## 
 ## No. of observations = 50
 ## 
@@ -2204,7 +2230,7 @@ 
    68.6.1 使用的數據和簡單
 ## 2 bm2       50   53.21  52.7    45.13  0      143.6 
 ## 3 bm3       50   61.43  55.25   51.47  0.2    147.9 
 ## 4 bm4       50   57.43  56.75   45.45  0.1    146.1
    -
    psych::describe(plant)
    + 
    ##     vars  n  mean    sd median trimmed   mad min   max range skew kurtosis   se
 ## bm1    1 50 56.60 48.05  47.55   53.36 66.05 0.0 143.0 143.0 0.27    -1.38 6.80
 ## bm2    2 50 53.21 45.13  52.70   49.62 59.45 0.0 143.6 143.6 0.43    -1.07 6.38
@@ -2215,7 +2241,7 @@ 
    68.6.1 使用的數據和簡單
 
  • 這些生物標幟物能否單獨提供關於該植物的某部分特徵信息呢？思考我們該如何回答這個問題（提示：計算這些指標直接的相關係數）
    • 
 
 
    我們可以通過計算這四個生物標幟物濃度測量值之間的相關係數，來觀察它們之間是否具有相似性或者是否提供了部分相似的信息。
    
-
    cor(plant)
    
+
 
    ##            bm1        bm2        bm3        bm4
 ## bm1 1.00000000 0.49826220 0.59414820 0.26769269
 ## bm2 0.49826220 1.00000000 0.50574946 0.33347350
@@ -2230,13 +2256,13 @@ 
    68.6.1 使用的數據和簡單
 
  • 再次思考問題1.的答案，思考並選擇合適的測量不同樣本個體之間距離 (distance) 的度量衡。嘗試使用簡單的聚類分析命令對樣本進行分類。
    • 
 
 
    由於每個生物標記物在所有樣本中的數值基本在相似的比例或者刻度，每個標幟物在這個樣本中的標準差/方差數值也較爲相近。我們嘗試用連續變量最常見的均值測量距離指標:
    
-
    # prepare hierarchical cluster
-hc <-  hclust(dist(plant), "ave")
-
-
-plot(hc, cex = 0.8, hang = -1, 
-     main = "", ylab = "L2 dissimilarity measure", 
-     xlab = "No. of specimen")
    
+
 
    
 Dendrogram for L2_avlink cluster analysis
 
    
@@ -2247,11 +2273,11 @@ 
    68.6.1 使用的數據和簡單
 
      
 
    1. 從簡單的歐幾里得距離改成歐幾里得距離平方來測量樣本之間的距離的話，圖形會變成什麼樣？
      2. 
 
    
-
    hc <- hclust(dist(plant)^2)
-
-plot(hc, cex = 0.8, hang = -1, 
-     main = "", ylab = "L2squared dissimilarity measure", 
-     xlab = "No. of specimen", sub = "")
    
+
 
    
 Dendrogram for L2sq_avlink cluster analysis
 
    
@@ -2260,8 +2286,8 @@ 
    68.6.1 使用的數據和簡單
 
    
 
    當使用歐幾里得距離的平方作爲樣本間隔的度量衡時，我們發現聚類的過程其實總體來說和使用歐幾里得距離本身並無本質上的區別。只是在差異性較低的地方聚類加速 (squeeze the dissimilarities at the lower end)，並且在較大的聚類區分之間變得更加明顯，視覺效果上更容易區分。
    
 
    如果說，不用歐幾里得平方，而是使用簡單的曼哈頓距離 (L1 度量衡)，那麼圖形則又會呈現爲:
    
-
    plot(cluster::agnes(plant, metric = "manhattan", stand = F), which.plots = 2, hang = -1, 
-     xlab = "No. of specimen", main = "", ylab = "L1 dissimilarity measure", sub = "", cex = 0.8)
    
+
 
    
 Dendrogram for L1_avlink cluster analysis
 
    
@@ -2336,7 +2362,7 @@ 
    68.6.1 使用的數據和簡單
   (function () {
     var script = document.createElement("script");
     script.type = "text/javascript";
-    var src = "";
+    var src = "true";
     if (src === "" || src === "true") src = "https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-MML-AM_CHTML";
     if (location.protocol !== "file:")
       if (/^https?:/.test(src))
diff --git a/docs/profile-log-likelihood.html b/docs/profile-log-likelihood.html
index 6fab65b68..1734ac648 100644
--- a/docs/profile-log-likelihood.html
+++ b/docs/profile-log-likelihood.html
@@ -6,7 +6,7 @@
   
   第 19 章 多個參數時的統計推斷 – 子集似然函數 profile log-likelihoods | 醫學統計學
   
-  
+  
 
   
   
@@ -24,7 +24,7 @@
 
 
 
-
+
 
   
   
@@ -61,41 +61,66 @@
 
 
 
 
 
@@ -575,7 +600,7 @@
 
  • 38.1 定義
    •

  • 39 The sandwich estimator
    • -
  • VII 貝葉斯統計
    • +
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
  • 40 貝葉斯統計入門
  • 39 The sandwich estimator
    • -
  • VII 貝葉斯統計
    • +
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
  • 40 貝葉斯統計入門
  • 39 The sandwich estimator
    • -
  • VII 貝葉斯統計
    • +
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
  • 40 貝葉斯統計入門
  • 39 The sandwich estimator
    • -
  • VII 貝葉斯統計
    • +
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
  • 40 貝葉斯統計入門
  • 39 The sandwich estimator
    • -
  • VII 貝葉斯統計
    • +
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
  • 40 貝葉斯統計入門
  • 39 The sandwich estimator
    • -
  • VII 貝葉斯統計
    • +
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
  • 40 貝葉斯統計入門
  • 39 The sandwich estimator
    • -
  • VII 貝葉斯統計
    • +
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
  • 40 貝葉斯統計入門
  • 39 The sandwich estimator
    • -
  • VII 貝葉斯統計
    • +
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
  • 40 貝葉斯統計入門
  • 39 The sandwich estimator
    • -
  • VII 貝葉斯統計
    • +
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
  • 40 貝葉斯統計入門

    注意我們用beta 分布來描述 p[2] 的先驗概率分布。

    -
    # Step 1 check model
-modelCheck(paste(bugpath, "/backupfiles/great-model2.txt", sep = ""))
    + 
    ## model is syntactically correct
    -
    # Load the data 
-modelData(paste(bugpath, "/backupfiles/great-data-alt.txt", sep = ""))
    + 
    ## data loaded
    -
    # compile the model with two separate chains
-modelCompile(numChains = 2)
    + 
    ## model compiled
    -
    # generate initial values 
-# the choice is arbitrary
-initlist <- list(LOR = 0.5, p = c(NA, 0.2))
-modelInits(bugsData(initlist))
    + 
    ## Initializing chain 1:
    ## initial values loaded and chain initialized but another chain contain uninitialized variables
    -
    initlist1 <- list(LOR = 5, p = c(NA, 0.8))
-modelInits(bugsData(initlist1))
    + 
    ## Initializing chain 2:
    ## model is initialized
    -
    # Set monitors on nodes of interest#### SPECIFY, WHICH PARAMETERS TO TRACE:
-parameters <- c("OR", "p")
-samplesSet(parameters)
    + 
    ## monitor set for variable 'OR'
    ## monitor set for variable 'p'
    -
    # Generate 10000 iterations
-modelUpdate(26000)
    + 
    ## 26000 updates took 0 s
    -
    # Summary Statistics
-sample.statistics <- samplesStats("*", beg = 1001)
-print(sample.statistics)
    + 
    ##        mean      sd  MC_error val2.5pc  median val97.5pc start sample
 ## OR   2.0790 0.77960 0.0073820  0.97320 1.94300    3.9890  1001  50000
 ## p[1] 0.1539 0.02947 0.0001353  0.10060 0.15230    0.2162  1001  50000
@@ -1948,34 +1972,34 @@ 
    86.5.1 The GREAT Trial
    
     LOR ~ dnorm(0,0.33) 
     OR <- exp(LOR)
 }
    -
    # Step 1 check model
-modelCheck(paste(bugpath, "/backupfiles/great-model2_forward.txt", sep = ""))
    + 
    ## model is syntactically correct
    -
    # Load the data 
-# modelData(paste(bugpath, "/backupfiles/great-forward.txt", sep = ""))     
-# compile the model with two separate chains
-modelCompile(numChains = 2)
    + 
    ## model compiled
    -
    # generate initial values 
-# the choice is arbitrary
-initlist <- list(LOR = 0.5, p = c(NA, 0.2))
-modelInits(bugsData(initlist))
    + 
    ## Initializing chain 1:
    ## initial values loaded and chain initialized but another chain contain uninitialized variables
    -
    initlist1 <- list(LOR = 5, p = c(NA, 0.8))
-modelInits(bugsData(initlist1))
    + 
    ## Initializing chain 2:
    ## model is initialized
    -
    # Set monitors on nodes of interest#### SPECIFY, WHICH PARAMETERS TO TRACE:
-parameters <- c("OR", "p")
-samplesSet(parameters)
    + 
    ## monitor set for variable 'OR'
    ## monitor set for variable 'p'
    -
    # Generate 10000 iterations
-modelUpdate(26000)
    + 
    ## 26000 updates took 0 s
    -
    # Density plot the first 1000 run 
-samplesDensity("p", mfrow = c(1,2), ask = FALSE)
    +
    Density plots for parameters prediction in GREAT trial second prior.

    @@ -2021,30 +2045,30 @@

    86.5.2.1 隊列研究設計

    list(lr0 = -1, lor = 0)
 list(lr0 = -5, lor = 5)

    對這麼模型進行50000次事後樣本採集，獲取比值比 OR 的事後概率分布描述：

    -
    # Step 1 check model
-modelCheck(paste(bugpath, "/backupfiles/cohort-model.txt", sep = "")) 
-# Load the data 
-modelData(paste(bugpath, "/backupfiles/dataforcohort.txt", sep = ""))     
-# compile the model with two separate chains
-modelCompile(numChains = 2) 
-# generate initial values 
-# the choice is arbitrary
-initlist <- list(lr0 = -1, lor = 0)
-modelInits(bugsData(initlist))
-initlist1 <- list(lr0 = -5, lor = 5)
-modelInits(bugsData(initlist1))
-
-# Set monitors on nodes of interest#### SPECIFY, WHICH PARAMETERS TO TRACE:
-parameters <- c("OR", "lor")
-samplesSet(parameters)
-
-
-# Generate 1000 iterations
-modelUpdate(1000)
-
-
-# Check trace history for the first 1000 run 
-samplesHistory("*", mfrow = c(2,1), ask = FALSE)
    +
    Density plots for parameters prediction in smoking and cancer cohort study.

    @@ -2066,12 +2090,12 @@

    86.5.2.1 隊列研究設計

    圖 86.12: Gelman-Rubin convergence statistic for smoking and cancer cohort study.

    -
    # Generate 50000 iterations
-modelUpdate(25000)
    + 
    ## 25000 updates took 0 s
    -
    # Summary Statistics
-sample.statistics <- samplesStats("*", beg = 1001)
-print(sample.statistics)
    + 
    ##      mean     sd MC_error val2.5pc median val97.5pc start sample
 ## OR  3.335 0.4513 0.003901   2.5360  3.304     4.302  1001  50000
 ## lor 1.195 0.1347 0.001163   0.9306  1.195     1.459  1001  50000
    @@ -2119,30 +2143,30 @@

    86.5.2.2 病例對照研究設計 X1cc = 600, Y1cc = 1000)

    模型擬合過程和結果分別羅列如下：

    -
    # Step 1 check model
-modelCheck(paste(bugpath, "/backupfiles/case-control-model.txt", sep = "")) 
-# Load the data 
-modelData(paste(bugpath, "/backupfiles/dataforcasecontrol.txt", sep = ""))     
-# compile the model with two separate chains
-modelCompile(numChains = 2) 
-# generate initial values 
-# the choice is arbitrary
-initlist <- list(lp0 = -2, lor = 0)
-modelInits(bugsData(initlist))
-initlist1 <- list(lp0 = 2, lor = 5)
-modelInits(bugsData(initlist1))
-
-# Set monitors on nodes of interest#### SPECIFY, WHICH PARAMETERS TO TRACE:
-parameters <- c("OR", "lor")
-samplesSet(parameters)
-
-
-# Generate 1000 iterations
-modelUpdate(1000)
-
-
-# Check trace history for the first 1000 run 
-samplesHistory("*", mfrow = c(2,1), ask = FALSE)
    +
    Density plots for parameters prediction in smoking and cancer case-control study.

    @@ -2164,12 +2188,12 @@

    86.5.2.2 病例對照研究設計 圖 86.14: Gelman-Rubin convergence statistic for smoking and cancer cohort study.

    -
    # Generate 50000 iterations
-modelUpdate(25000)
    + 
    ## 25000 updates took 0 s
    -
    # Summary Statistics
-sample.statistics <- samplesStats("*", beg = 1001)
-print(sample.statistics)
    + 
    ##       mean      sd  MC_error val2.5pc median val97.5pc start sample
 ## OR  2.2530 0.17780 0.0011660   1.9260 2.2450    2.6210  1001  50000
 ## lor 0.8092 0.07875 0.0005161   0.6555 0.8088    0.9636  1001  50000
    @@ -2214,30 +2238,30 @@

    86.5.2.3 聯合模型 joint model Y1cc = 1000 )

    接下來是對聯合模型的擬合及對OR值事後樣本的採集過程：

    -
    # Step 1 check model
-modelCheck(paste(bugpath, "/backupfiles/jointmodel.txt", sep = "")) 
-# Load the data 
-modelData(paste(bugpath, "/backupfiles/dataforjoint.txt", sep = ""))     
-# compile the model with two separate chains
-modelCompile(numChains = 2) 
-# generate initial values 
-# the choice is arbitrary
-initlist <- list(lr0 = -1, lp0 = -2, lor = 0)
-modelInits(bugsData(initlist))
-initlist1 <- list(lr0 = -5, lp0 = 2, lor = 5)
-modelInits(bugsData(initlist1))
-
-# Set monitors on nodes of interest#### SPECIFY, WHICH PARAMETERS TO TRACE:
-parameters <- c("OR", "lor")
-samplesSet(parameters)
-
-
-# Generate 1000 iterations
-modelUpdate(1000)
-
-
-# Check trace history for the first 1000 run 
-samplesHistory("*", mfrow = c(2,1), ask = FALSE)
    +
    Density plots for parameters prediction in smoking and cancer joint model.

    @@ -2259,12 +2283,12 @@

    86.5.2.3 聯合模型 joint model 圖 86.16: Gelman-Rubin convergence statistic for smoking and cancer cohort study.

    -
    # Generate 50000 iterations
-modelUpdate(25000)
    + 
    ## 25000 updates took 0 s
    -
    # Summary Statistics
-sample.statistics <- samplesStats("*", beg = 1001)
-print(sample.statistics)
    + 
    ##       mean      sd  MC_error val2.5pc median val97.5pc start sample
 ## OR  2.4930 0.17070 0.0009080   2.1760 2.4870     2.844  1001  50000
 ## lor 0.9111 0.06842 0.0003648   0.7777 0.9111     1.045  1001  50000
    @@ -2438,7 +2462,7 @@

    86.5.2.4 你是否能證明兩個 (function () { var script = document.createElement("script"); script.type = "text/javascript"; - var src = ""; + var src = "true"; if (src === "" || src === "true") src = "https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-MML-AM_CHTML"; if (location.protocol !== "file:") if (/^https?:/.test(src)) diff --git "a/docs/\344\272\214\351\240\205\345\210\206\344\275\210\346\225\270\346\223\232\347\232\204\345\273\243\347\276\251\347\267\232\346\200\247\350\277\264\346\255\270\346\250\241\345\236\213-logistic-regression-model.html" "b/docs/\344\272\214\351\240\205\345\210\206\344\275\210\346\225\270\346\223\232\347\232\204\345\273\243\347\276\251\347\267\232\346\200\247\350\277\264\346\255\270\346\250\241\345\236\213-logistic-regression-model.html" index 1b5cfa0d0..729f58a92 100644 --- "a/docs/\344\272\214\351\240\205\345\210\206\344\275\210\346\225\270\346\223\232\347\232\204\345\273\243\347\276\251\347\267\232\346\200\247\350\277\264\346\255\270\346\250\241\345\236\213-logistic-regression-model.html" +++ "b/docs/\344\272\214\351\240\205\345\210\206\344\275\210\346\225\270\346\223\232\347\232\204\345\273\243\347\276\251\347\267\232\346\200\247\350\277\264\346\255\270\346\250\241\345\236\213-logistic-regression-model.html" @@ -6,7 +6,7 @@ 第 45 章 二項分佈數據的廣義線性迴歸模型 logistic regression model | 醫學統計學 - + @@ -24,7 +24,7 @@ - + @@ -61,41 +61,66 @@ @@ -575,7 +600,7 @@
  • 38.1 定義

  • 39 The sandwich estimator
    • -
  • VII 貝葉斯統計
    • +
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
  • 40 貝葉斯統計入門
  • 39 The sandwich estimator
    • -
  • VII 貝葉斯統計
    • +
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
  • 40 貝葉斯統計入門
  • 39 The sandwich estimator
    • -
  • VII 貝葉斯統計
    • +
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
  • 40 貝葉斯統計入門
  • 39 The sandwich estimator
    • -
  • VII 貝葉斯統計
    • +
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
  • 40 貝葉斯統計入門
  • 39 The sandwich estimator
    • -
  • VII 貝葉斯統計
    • +
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
  • 40 貝葉斯統計入門
  • 39 The sandwich estimator
    • -
  • VII 貝葉斯統計
    • +
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
  • 40 貝葉斯統計入門
  • 39 The sandwich estimator
    • -
  • VII 貝葉斯統計
    • +
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
  • 40 貝葉斯統計入門 -

    你是否還記得貝葉斯定理的公式: 如果\(A, B\)分別標記兩個事件，那麼有

    +

    你是否還記得貝葉斯定理的公式: +如果\(A, B\)分別標記兩個事件，那麼有

    \[ p(A|B) = \frac{p(B|A)p(A)}{p(B)} \]

    @@ -1375,44 +1400,45 @@

    83.2 二項分布(似然)數據 & \propto \theta^y(1-\theta)^{n-y}\theta^{a-1}(1-\theta)^{b-1} \\ & = \theta^{y + a -1}(1-\theta)^{n - y + b -1} \end{aligned} -\] 眼尖的人立刻能看出來，這個事後概率分布本身也是一個Beta分布的概率方程，只是它的超參數和先驗概率相比發生了變化(更新)：

    +\] +眼尖的人立刻能看出來，這個事後概率分布本身也是一個Beta分布的概率方程，只是它的超參數和先驗概率相比發生了變化(更新)：

    \[ p(\theta | y,n) = \text{Beta}(a + y, b + n -y) \]

    像這樣，先驗概率和事後概率兩個概率分布都來自相同分佈家族的情況，先驗概率又被叫做和似然成共軛關系的先驗概率(共軛先驗概率, conjugate prior)。

    -
    par(mfrow=c(3,1))
-# Plot exact prior probability density 
-# values of the hyperparameters
-a <- 9.2 
-b <- 13.8
-
-# prior function
-prior <- Vectorize(function(theta) dbeta(theta, a, b))
-# Plot 
-curve(prior, 0, 1, type = "l", main = "Prior for "~theta, ylim = c(0, 4.5), frame = F,
-      xlab = "Probability of positive response", ylab = "Density", lwd = 2, cex.axis = 1.5, cex.lab = 1.5)
-
-# binomial likelihood function (likelihood)
-
-Likelihood <- Vectorize(function(theta) dbinom(15, 20, theta))
-
-# Plot
-curve(Likelihood, 0, 1, type = "l", main = "Likelihood for the data", 
-      frame = FALSE, xlab = "Probability of positive response", ylab = "Density", 
-      lwd = 2, cex.axis =  1.5, cex.lab = 1.5)
-# n <- 0; r <- 0; a <- 9.2; b <- 13.8; np <- 20
-# plot(0:20, BetaBinom(0:20), type = "b", xlab = "r*", ylab = "P(R = r*)", 
-#      main = "Prior predictive: a = 9.2, b = 13.8")
-
-# Posterior function 
-
-posterior <- Vectorize(function(theta) dbeta(theta, a+15, b+20-15))
-
-# Plot
-
-curve(posterior, 0, 1, type = "l", main = "Posterior for "~theta, 
-      frame = FALSE, xlab = "Probability of positive response", ylab = "Density", 
-      lwd = 2, cex.axis = 1.5, cex.lab = 1.5)
    +
    Prior, likelihood, and posterior for Drug example

    @@ -1461,26 +1487,26 @@

    83.2.1 事後概率分布預測 -
    Dat <- list(a = 9.2,             # parameters of prior distribution 
-            b = 13.8, 
-            y = 15,              # number of successes in completed trial
-            n = 20,              # number of patients in completed trial
-            m = 40,              # number of patients in future trial 
-            mcrit = 25)          # critical value of future successes
-
-# fit use R2jags package
-
-post.jags <- jags(
-  data = Dat,
-  parameters.to.save = c("P.crit", "theta", "y.pred"),
-  n.iter = 50100,
-  model.file = paste(bugpath, 
-                     "/backupfiles/MCdrugP29.txt", 
-                     sep = ""),
-  n.chains = 1,
-  n.thin = 10,
-  n.burnin = 100,
-  progress.bar = "none")
    + 
    ## Compiling model graph
 ##    Resolving undeclared variables
 ##    Allocating nodes
@@ -1490,7 +1516,7 @@ 
    83.2.1 事後概率分布預測
    -
    print(post.jags)
    + 
    ## Inference for Bugs model at "/home/ccwang/Documents/LSHTMlearningnote/backupfiles/MCdrugP29.txt", fit using jags,
 ##  1 chains, each with 50100 iterations (first 100 discarded), n.thin = 10
 ##  n.sims = 5000 iterations saved
@@ -1503,39 +1529,39 @@ 
    83.2.1 事後概率分布預測
    -
    Simulated <- coda::as.mcmc(post.jags)
-#### PUT THE SAMPLED VALUES IN ONE R DATA FRAME:
-ggSample <- ggmcmc::ggs(Simulated)
-
-#### PLOT THE DENSITY and HISTOGRAMS OF THE SAMPLED VALUES
-par(mfrow=c(1,2))
-
-Theta <- ggSample %>% 
-  filter(Parameter == "theta")
-Y <- ggSample %>% 
-  filter(Parameter == "y.pred")
-plot(density(Theta$value), main = "theta sample 50000", 
-     ylab = "P(theta)", xlab = "Probability of response", col = "red")
-hist(Y$value, main = "y.pred sample 50000", ylab = "P(y.pred)", 
-     xlab = "Number of success", col = "red", prob = TRUE,xlim = c(0, 40))
    +
    Posterior and predictive distributions for Drug example

    圖 83.2: Posterior and predictive distributions for Drug example

    -
    # # Or simply use the denplot() function from the mcmcplots package
-# mcmcplots::denplot(post.jags, "theta")
    +

    圖83.2左邊的圖，是前一次試驗結果的事後概率分布，20名患者中觀察到15名患者症狀改善。右邊的圖則是對下一次40人的試驗的結果做的預測，平均22.5名患者可能會有症狀改善，這個均值的標準差是4.3。

    -
    mcmcplots::traplot(post.jags, "theta")
    +
    Plot of the MCMC chain of the parameter, Drug example.

    圖 83.3: Plot of the MCMC chain of the parameter, Drug example.

    -
    # plot(Theta$value, main="", type="l", ylab="theta", 
-#      xlab="iteration", col="red", ylim = c(0,1.2))
    +

    爲了比較，我們可以把精確計算獲得的答案和蒙特卡羅算法給出的預測做個比較：

    1. \(\theta:\)均值爲0.563，標準差是0.075；
      2. @@ -1579,34 +1605,34 @@

      83.3 正態分布(似然)數據

    2. 事後概率分布的方差也是和先前提到的先驗概率樣本量有密切關系。它是觀測數據的方差除以觀測樣本量和先驗概率樣本量之和。

    3. 當然，當你的觀測數據樣本量趨向於無窮大時 \(n \rightarrow \infty\)，事後概率分布本身就接近與觀測數據的似然 \(p(\theta | \mathbf{y}) \rightarrow N(\bar y, \sigma^2/n)\)。也就是說觀測數據的信息量佔絕對主導，先驗概率分布，不再能提供太多有價值的信息，可以忽略不計了。

      4. -
      # prior function
-xseq<-seq(80,180,.01)
-densities<-dnorm(xseq, 120,10)
-
-# Plot 
-plot(xseq, densities, col="darkgreen",xlab="mean THM concentration, ug/l (theta)", 
-     ylab="Density", type="l",lwd=2, cex=2, 
-     # main="PDF of Prior, likelihood, and posterior for THM example.", 
-     cex.axis=0.9, cex.lab = 1, ylim = c(0,0.12))
-
-# normal likelihood function (likelihood)
-
-Likelihood <- dnorm(xseq, 130, 5)
-
-# Plot
-points(xseq, Likelihood, col="darkred", type="l",lwd=2, cex=2)
-
-
-# Posterior function 
-
-posterior <-  dnorm(xseq, 128.9, 3.33)
-
-# Plot
-
-points(xseq, posterior, col="black", type="l",lwd=2, cex=2)
-
-legend("topright", c("Prior", "Likelihood", "Posterior"), bty = "n", lwd = 2,
-       lty = c(1,1,1), col = c("darkgreen", "darkred", "black"))
      +
      PDF of Prior, likelihood and posterior for THM example.

      @@ -1783,50 +1809,50 @@

      83.6 Practical Bayesian Statistic m = 40, # future number of trials ncrit = 25) # critical value of future successes

      把這裏的模型保存稱爲 drug-model.txt文件，把數據保存成 drug-data.txt文件，我們來試着用OpenBUGS/JAGS跑這個模型:

      -
      # Codes for OpenBUGS
-
-# # Step 1 check model
-# modelCheck(paste(bugpath, "/backupfiles/drug-model.txt", sep = "")) 
-# # Load the data 
-# modelData(paste(bugpath, "/backupfiles/drug-data.txt", sep = ""))     
-# # compile the model
-# modelCompile(numChains = 1) 
-# # There is no need to provide initial values as 
-# # they are aleady provided in the model specification
-# modelGenInits() #
-# # Set monitors on nodes of interest#### SPECIFY, WHICH PARAMETERS TO TRACE:
-# parameters <- c("P.crit", "theta", "y.pred")
-# samplesSet(parameters)
-
-# # Generate 10000 iterations
-# modelUpdate(10000)
-# #### SHOW POSTERIOR STATISTICS
-# sample.statistics <- samplesStats("*")
-# print(sample.statistics)
-
-# Codes for R2JAGS
-
-
-Dat <- list(a = 9.2,             # parameters of prior distribution 
-            b = 13.8, 
-            y = 15,              # number of successes in completed trial
-            n = 20,              # number of patients in completed trial
-            m = 40,              # number of patients in future trial 
-            ncrit = 25)          # critical value of future successes
-
-# fit use R2jags package
-
-post.jags <- jags(
-  data = Dat,
-  parameters.to.save = c("P.crit", "theta", "y.pred"),
-  n.iter = 10100,
-  model.file = paste(bugpath, 
-                     "/backupfiles/drug-model.txt", 
-                     sep = ""),
-  n.chains = 1,
-  n.thin = 1,
-  n.burnin = 100,
-  progress.bar = "none")
      + 
      ## Compiling model graph
 ##    Resolving undeclared variables
 ##    Allocating nodes
@@ -1836,7 +1862,7 @@ 
      83.6 Practical Bayesian Statistic
 ##    Total graph size: 12
 ## 
 ## Initializing model
      -
      print(post.jags)
      + 
      ## Inference for Bugs model at "/home/ccwang/Documents/LSHTMlearningnote/backupfiles/drug-model.txt", fit using jags,
 ##  1 chains, each with 10100 iterations (first 100 discarded)
 ##  n.sims = 10000 iterations saved
@@ -1851,16 +1877,16 @@ 
      83.6 Practical Bayesian Statistic
 ## DIC is an estimate of expected predictive error (lower deviance is better).

      此時我們獲得我們關心的各個事後概率分佈描述，其中陽性反應率的點估計是0.563，其95%可信區間是(0.412, 0.707)。40名患者中25人或者以上出現有療效反應的概率是32.4%。

      請繪製theta的事後概率分佈的密度曲線，以及y.pred的預測概率分佈:

      -
      #### PLOT THE DENSITY and HISTOGRAMS OF THE SAMPLED VALUES
-par(mfrow=c(1,2))
-Theta <- ggSample %>% 
-  filter(Parameter == "theta")
-Y <- ggSample %>% 
-  filter(Parameter == "y.pred")
-plot(density(Theta$value), main = "theta sample 10000", 
-     ylab = "P(theta)", xlab = "Probability of response", col = "red")
-hist(Y$value, main = "y.pred sample 10000", ylab = "P(y.pred)", 
-     xlab = "Number of success", col = "red", prob = TRUE,xlim = c(0, 40))
      +
      Posterior and predictive distributions for Drug example

      @@ -1908,51 +1934,51 @@

      83.6 Practical Bayesian Statistic m = 40, # future number of trials ncrit = 25 # critical value of future successes ) -
      # Codes for OpenBUGS
-
-# # Step 1 check model
-# modelCheck(paste(bugpath, "/backupfiles/drug-modeluniform.txt", sep = "")) 
-# # Load the data 
-# modelData(paste(bugpath, "/backupfiles/drug-datauniform.txt", sep = ""))     
-# # compile the model
-# modelCompile(numChains = 1) 
-# # There is no need to provide initial values as 
-# # they are aleady provided in the model specification
-# modelGenInits() #
-# # Set monitors on nodes of interest#### SPECIFY, WHICH PARAMETERS TO TRACE:
-# parameters <- c("P.crit", "theta", "y.pred")
-# samplesSet(parameters)
-# 
-# # Generate 10000 iterations
-# modelUpdate(10000)
-# #### SHOW POSTERIOR STATISTICS
-# sample.statistics <- samplesStats("*")
-# print(sample.statistics)
-
-
-# Codes for R2JAGS
-
-
-Dat <- list(
-y = 15,                                # number of successes
-n = 20,                                # number of trials 
-m = 40,                                # future number of trials 
-ncrit = 25                             # critical value of future successes
-)
-
-# fit use R2jags package
-
-post.jags <- jags(
-  data = Dat,
-  parameters.to.save = c("P.crit", "theta", "y.pred"),
-  n.iter = 10100,
-  model.file = paste(bugpath, 
-                     "/backupfiles/drug-modeluniform.txt", 
-                     sep = ""),
-  n.chains = 1,
-  n.thin = 1,
-  n.burnin = 100,
-  progress.bar = "none")
      + 
      ## Compiling model graph
 ##    Resolving undeclared variables
 ##    Allocating nodes
@@ -1962,7 +1988,7 @@ 
      83.6 Practical Bayesian Statistic
 ##    Total graph size: 12
 ## 
 ## Initializing model
      -
      print(post.jags)
      + 
      ## Inference for Bugs model at "/home/ccwang/Documents/LSHTMlearningnote/backupfiles/drug-modeluniform.txt", fit using jags,
 ##  1 chains, each with 10100 iterations (first 100 discarded)
 ##  n.sims = 10000 iterations saved
@@ -2006,27 +2032,27 @@ 
      83.6 Practical Bayesian Statistic
    # theta ~ dunif(-10000, 10000) 
 }

      在BUGS語言中，正態分佈用 dnorm(theta, tau)，其中 theta 爲均值，tau 是精確度(precision = 1/variance)，它是方差的倒數。

      -
      # Codes for R2JAGS
-
-
-Dat <- list(
-y = c(128, 132),                       # observed 
-tau = 1/(5^2)
-)
-
-# fit use R2jags package
-
-post.jags <- jags(
-  data = Dat,
-  parameters.to.save = c("theta"),
-  n.iter = 10100,
-  model.file = paste(bugpath, 
-                     "/backupfiles/thm-model.txt", 
-                     sep = ""),
-  n.chains = 1,
-  n.thin = 1,
-  n.burnin = 100,
-  progress.bar = "none")
      + 
      ## Compiling model graph
 ##    Resolving undeclared variables
 ##    Allocating nodes
@@ -2036,7 +2062,7 @@ 
      83.6 Practical Bayesian Statistic
 ##    Total graph size: 8
 ## 
 ## Initializing model
      -
      print(post.jags)
      + 
      ## Inference for Bugs model at "/home/ccwang/Documents/LSHTMlearningnote/backupfiles/thm-model.txt", fit using jags,
 ##  1 chains, each with 10100 iterations (first 100 discarded)
 ##  n.sims = 10000 iterations saved
@@ -2047,47 +2073,47 @@ 
      83.6 Practical Bayesian Statistic
 ## DIC info (using the rule, pD = var(deviance)/2)
 ## pD = 1.0 and DIC = 12.4
 ## DIC is an estimate of expected predictive error (lower deviance is better).
      -
      # OpenBUGS code:
-# # Step 1 check model
-# modelCheck(paste(bugpath, "/backupfiles/thm-model.txt", sep = "")) 
-#   
-# # compile the model
-# modelCompile(numChains = 1) 
-# # There is no need to provide initial values as 
-# # they are aleady provided in the model specification
-# modelGenInits() #
-# # Set monitors on nodes of interest#### SPECIFY, WHICH PARAMETERS TO TRACE:
-# parameters <- c("theta")
-# samplesSet(parameters)
-
-# # Generate 10000 iterations
-# modelUpdate(10000)
-# #### SHOW POSTERIOR STATISTICS
-# sample.statistics <- samplesStats("*")
-# print(sample.statistics)
      +

      所以這個區域三氯甲烷的事後均值和95%可信區間分別是 128.9 (122.3, 135.4)。和我們之前用精確計算法給出的結果相同/相近。

      如果在這個模型中，我們嘗試沒有信息的先驗概率分佈，結果會怎樣呢？

      -
      # Codes for R2JAGS
-
-
-Dat <- list(
-y = c(128, 132),                       # observed 
-tau = 1/(5^2)
-)
-
-# fit use R2jags package
-
-post.jags <- jags(
-  data = Dat,
-  parameters.to.save = c("theta"),
-  n.iter = 10100,
-  model.file = paste(bugpath, 
-                     "/backupfiles/THM-vaguemodel.txt", 
-                     sep = ""),
-  n.chains = 1,
-  n.thin = 1,
-  n.burnin = 100,
-  progress.bar = "none")
      + 
      ## Compiling model graph
 ##    Resolving undeclared variables
 ##    Allocating nodes
@@ -2097,7 +2123,7 @@ 
      83.6 Practical Bayesian Statistic
 ##    Total graph size: 6
 ## 
 ## Initializing model
      -
      print(post.jags)
      + 
      ## Inference for Bugs model at "/home/ccwang/Documents/LSHTMlearningnote/backupfiles/THM-vaguemodel.txt", fit using jags,
 ##  1 chains, each with 10100 iterations (first 100 discarded)
 ##  n.sims = 10000 iterations saved
@@ -2108,24 +2134,24 @@ 
      83.6 Practical Bayesian Statistic
 ## DIC info (using the rule, pD = var(deviance)/2)
 ## pD = 0.9 and DIC = 12.4
 ## DIC is an estimate of expected predictive error (lower deviance is better).
      -
      # OpenBUGS code:
-# # Step 1 check model
-# modelCheck(paste(bugpath, "/backupfiles/THM-vaguemodel.txt", sep = "")) 
-# 
-# # compile the model
-# modelCompile(numChains = 1) 
-# # There is no need to provide initial values as 
-# # they are aleady provided in the model specification
-# modelGenInits() #
-# # Set monitors on nodes of interest#### SPECIFY, WHICH PARAMETERS TO TRACE:
-# parameters <- c("theta")
-# samplesSet(parameters)
-# 
-# # Generate 10000 iterations
-# modelUpdate(10000)
-# #### SHOW POSTERIOR STATISTICS
-# sample.statistics <- samplesStats("*")
-# print(sample.statistics)
      +

      MC試驗10000次給出的事後均值是129.9，它和樣本均值相同，因爲此時我們給模型一個幾乎不含有效信息的先驗概率分佈時，意味着我們讓試驗獲得的數據做完全主導作用，“data speak for themselves”。

      C. Disease risk in a small area

      下面的BUGS/JAGS模型代碼可以用來進行泊淞－伽馬分布似然的MC計算。在這個例子中，在某個區域我們觀察到５例白血病新病例，已知該區域的年齡性別標準化發病期望病例數是\(E = 2.8\)。注意看我們在代碼中加入了兩種不同類型的先驗概率分布，一個是沒有太多信息的(vague prior distribution)dgamma(0.1, 0.1)，一個則是有一定信息的(informative prior, stribution) dgamma(48, 40)。我們把兩個先驗概率分布同時加在一個模型裏，這是十分便於進行兩種先驗概率對結果的影響的比較的手段。你當然可以把它寫成兩個不同的模型。注意模型代碼是如何表示事後概率分布和計算相對危險度比(relative risk)大於1的概率的。

      @@ -2150,26 +2176,26 @@

      83.6 Practical Bayesian Statistic # y[2] <- 5 # replicate data to fit both models together } -
      # Codes for R2JAGS
-
-
-Dat <- list(
-y = c(5, 5)                      # replicate data to fit both models together
-)
-
-# fit use R2jags package
-
-post.jags <- jags(
-  data = Dat,
-  parameters.to.save = c("P.excess[1]", "P.excess[2]", "lambda[1]", "lambda[2]"),
-  n.iter = 10100,
-  model.file = paste(bugpath, 
-                     "/backupfiles/disease-model.txt", 
-                     sep = ""),
-  n.chains = 1,
-  n.thin = 1,
-  n.burnin = 100,
-  progress.bar = "none")
      + 
      ## Compiling model graph
 ##    Resolving undeclared variables
 ##    Allocating nodes
@@ -2179,7 +2205,7 @@ 
      83.6 Practical Bayesian Statistic
 ##    Total graph size: 15
 ## 
 ## Initializing model
      -
      print(post.jags)
      + 
      ## Inference for Bugs model at "/home/ccwang/Documents/LSHTMlearningnote/backupfiles/disease-model.txt", fit using jags,
 ##  1 chains, each with 10100 iterations (first 100 discarded)
 ##  n.sims = 10000 iterations saved
@@ -2193,59 +2219,59 @@ 
      83.6 Practical Bayesian Statistic
 ## DIC info (using the rule, pD = var(deviance)/2)
 ## pD = 1.1 and DIC = 9.8
 ## DIC is an estimate of expected predictive error (lower deviance is better).
      -
      # # Step 1 check model
-# modelCheck(paste(bugpath, "/backupfiles/disease-model.txt", sep = "")) 
-# 
-# # compile the model
-# modelCompile(numChains = 1) 
-# # There is no need to provide initial values as 
-# # they are aleady provided in the model specification
-# modelGenInits() #
-# # Set monitors on nodes of interest#### SPECIFY, WHICH PARAMETERS TO TRACE:
-# parameters <- c("P.excess[1]", "P.excess[2]", "lambda[1]", "lambda[2]")
-# samplesSet(parameters)
-# 
-# # Generate 10000 iterations
-# modelUpdate(10000)
-# #### SHOW POSTERIOR STATISTICS
-# sample.statistics <- samplesStats("*")
-# print(sample.statistics)
      +

      在沒有信息的先驗概率(dgamma(0.1, 0.1))分布條件下，該地區可能有較高白血病發病率的概率是85%，但是在有參考價值信息的先驗概率分布(dgamma(48, 40))條件下，該地區可能有較高白血病發病率的概率是93%。所以，盡管相對危險度(relative risk)的事後均值在沒有信息的先驗概率分布條件下比較高(lambda[1] = 1.770 > lambda[2] = 1.238)，但是沒有信息的先驗概率分布條件下，這個相對危險度大於1的概率(85%)要比有信息的先驗概率分布條件下相對危險度大於1的概率要低(93%)。這是因爲在沒有信息的先驗概率分布條件下，相對危險度估計本身有太多的不確定性(uncertainty，圖83.7)。

      -
      # 
-# # Step 1 check model
-# modelCheck(paste(bugpath, "/backupfiles/disease-model.txt", sep = "")) 
-# 
-# # compile the model
-# modelCompile(numChains = 1) 
-# # There is no need to provide initial values as 
-# # they are aleady provided in the model specification
-# modelGenInits() #
-# # Set monitors on nodes of interest#### SPECIFY, WHICH PARAMETERS TO TRACE:
-# parameters <- c("P.excess[1]", "P.excess[2]", "lambda[1]", "lambda[2]")
-# samplesSet(parameters)
-# 
-# # Generate 10000 iterations
-# modelUpdate(10000)
-
-
-#### PUT THE SAMPLED VALUES IN ONE R DATA FRAME:
-# chain <- data.frame(P.excess1 = samplesSample("P.excess[1]"), 
-#                     P.excess2 = samplesSample("P.excess[2]"), 
-#                     lambda1 = samplesSample("lambda[1]"), 
-#                     lambda2 = samplesSample("lambda[2]"))
-Simulated <- coda::as.mcmc(post.jags)
-#### PUT THE SAMPLED VALUES IN ONE R DATA FRAME:
-ggSample <- ggmcmc::ggs(Simulated)
-
-Lambda1 <- ggSample %>% 
-  filter(Parameter == "lambda[1]")
-Lambda2 <- ggSample %>% 
-  filter(Parameter == "lambda[2]")
-
-boxplot(Lambda1$value,Lambda2$value, col = "green", ylab = "lambda", 
-        outline=FALSE, main = "boxplot: lambda", ylim = c(0,4), 
-        names = c("1", "2")) 
-abline(h = 1.5)
      +
      Box plots of relative risk (lambda) of leukaemia under different priors (vague = 1, informative = 2).

      @@ -2275,26 +2301,26 @@

      83.6 Practical Bayesian Statistic # y[2] <- 100 # replicate data to fit both models together } -
      # Codes for R2JAGS
-
-
-Dat <- list(
-y = c(100, 100)                    # the observed new cases changed from 5 to 100
-)
-
-# fit use R2jags package
-
-post.jags <- jags(
-  data = Dat,
-  parameters.to.save = c("P.excess[1]", "P.excess[2]", "lambda[1]", "lambda[2]"),
-  n.iter = 10100,
-  model.file = paste(bugpath, 
-                     "/backupfiles/disease-modelupdated.txt", 
-                     sep = ""),
-  n.chains = 1,
-  n.thin = 1,
-  n.burnin = 100,
-  progress.bar = "none")
      + 
      ## Compiling model graph
 ##    Resolving undeclared variables
 ##    Allocating nodes
@@ -2304,7 +2330,7 @@ 
      83.6 Practical Bayesian Statistic
 ##    Total graph size: 15
 ## 
 ## Initializing model
      -
      print(post.jags)
      + 
      ## Inference for Bugs model at "/home/ccwang/Documents/LSHTMlearningnote/backupfiles/disease-modelupdated.txt", fit using jags,
 ##  1 chains, each with 10100 iterations (first 100 discarded)
 ##  n.sims = 10000 iterations saved
@@ -2318,24 +2344,24 @@ 
      83.6 Practical Bayesian Statistic
 ## DIC info (using the rule, pD = var(deviance)/2)
 ## pD = 4.0 and DIC = 20.7
 ## DIC is an estimate of expected predictive error (lower deviance is better).
      -
      # OpenBUGS codes:
-# # Step 1 check model
-# modelCheck(paste(bugpath, "/backupfiles/disease-modelupdated.txt", sep = "")) 
-# 
-# # compile the model
-# modelCompile(numChains = 1) 
-# # There is no need to provide initial values as 
-# # they are aleady provided in the model specification
-# modelGenInits() #
-# # Set monitors on nodes of interest#### SPECIFY, WHICH PARAMETERS TO TRACE:
-# parameters <- c("P.excess[1]", "P.excess[2]", "lambda[1]", "lambda[2]")
-# samplesSet(parameters)
-# 
-# # Generate 10000 iterations
-# modelUpdate(10000)
-# #### SHOW POSTERIOR STATISTICS
-# sample.statistics <- samplesStats("*")
-# print(sample.statistics)
      +

      有了更多的觀察數據的模型，我們看到，相對危險度的事後估計lambda[1], lambda[2]都變得更加精確了，有了更小的標準差。這和我們的預期相同，因爲觀察數據增加自然會提高相對危險度估計的精確度。可是，我們發現在沒有信息的先驗概率分布條件下，相對危險度的事後估計幾乎沒有太大的變化(1.783 v.s. 1.770)。相反地，有信息的先驗概率分布條件下，相對危險度的事後估計比之前增加了不少(1.540 v.s. 1.238)。這主要是因爲，現在我們得到更多的觀察數據，這些數據得到的信息被給予了更多的權重，且觀察數據提示相對危險度應該要比較大(100/56 = 1.786)。

      盡管觀察數據確實給模型提供了較多的有價值的信息，你會發現，我們使用的第二個先驗概率分布（也就是有信息的先驗概率分布）仍然對相對危險度的事後估計起到了相當的影響(the prior variance is 0.03, giving prior sd of 0.17)。這也是因爲這個先驗概率分布給出的信息量，幾乎相當於觀察數據給出的信息量（二者的標準差很接近）。另外一種理解此現象的方法是看事後概率分布推導出的新的伽馬分布的計算式：

      \[ @@ -2347,24 +2373,24 @@

      83.6 Practical Bayesian Statistic \end{aligned} \]

      由於事後伽馬分布的方差（標準差）主要由第二個參數\((b + E)\)決定，從上面的公式推導我們可以看見，事後伽馬分布的第二個參數\((b + E)\)，分別是先驗概率分布的第二個參數\((b = 40)\)，和觀察數據的期望值\((E = 56)\)。由於這兩個數值大小接近，所以我們也可以理解此時先驗概率提供的信息和我們觀察數據提供的信息旗鼓相當。另外，在這個新的觀察數據條件下，我們無論使用哪個先驗概率分布做條件，都獲得了100%的結果也就是這個特定區域的白血病發病率大於期望值的概率是100%。也就是我們此時有100%的把握認爲這個特定區域的白血病發病率較高。

      -
      #### PUT THE SAMPLED VALUES IN ONE R DATA FRAME:
-# chain <- data.frame(P.excess1 = samplesSample("P.excess[1]"), 
-#                     P.excess2 = samplesSample("P.excess[2]"), 
-#                     lambda1 = samplesSample("lambda[1]"), 
-#                     lambda2 = samplesSample("lambda[2]"))
-Simulated <- coda::as.mcmc(post.jags)
-#### PUT THE SAMPLED VALUES IN ONE R DATA FRAME:
-ggSample <- ggmcmc::ggs(Simulated)
-
-Lambda1 <- ggSample %>% 
-  filter(Parameter == "lambda[1]")
-Lambda2 <- ggSample %>% 
-  filter(Parameter == "lambda[2]")
-
-boxplot(Lambda1$value,Lambda2$value, col = "green", ylab = "lambda", 
-        outline=FALSE, main = "boxplot: lambda", ylim = c(0,4), 
-        names = c("1", "2")) 
-abline(h = 1.66)
      +
      Box plots of relative risk (lambda) of leukaemia under different priors (vague = 1, informative = 2) with more observations.

      @@ -2395,26 +2421,26 @@

      83.6 Practical Bayesian Statistic # y <- 13 # observed number of correct taste tests in original experiment } -
      # Codes for R2JAGS
-
-
-Dat <- list(
-y = c(13)                    # Bond had 13 correct taste tests
-)
-
-# fit use R2jags package
-
-post.jags <- jags(
-  data = Dat,
-  parameters.to.save = c("P.ability", "P.Moneyback", "theta", "y.pred"),
-  n.iter = 100100,
-  model.file = paste(bugpath, 
-                     "/backupfiles/bondmodel.txt", 
-                     sep = ""),
-  n.chains = 1,
-  n.thin = 1,
-  n.burnin = 100,
-  progress.bar = "none")
      + 
      ## Compiling model graph
 ##    Resolving undeclared variables
 ##    Allocating nodes
@@ -2424,7 +2450,7 @@ 
      83.6 Practical Bayesian Statistic
 ##    Total graph size: 12
 ## 
 ## Initializing model
      -
      print(post.jags)
      + 
      ## Inference for Bugs model at "/home/ccwang/Documents/LSHTMlearningnote/backupfiles/bondmodel.txt", fit using jags,
 ##  1 chains, each with 100100 iterations (first 100 discarded)
 ##  n.sims = 100000 iterations saved
@@ -2438,23 +2464,23 @@ 
      83.6 Practical Bayesian Statistic
 ## DIC info (using the rule, pD = var(deviance)/2)
 ## pD = 0.9 and DIC = 4.6
 ## DIC is an estimate of expected predictive error (lower deviance is better).
      -
      # # Step 1 check model
-# modelCheck(paste(bugpath, "/backupfiles/bondmodel.txt", sep = "")) 
-# 
-# # compile the model
-# modelCompile(numChains = 1) 
-# # There is no need to provide initial values as 
-# # they are aleady provided in the model specification
-# modelGenInits() #
-# # Set monitors on nodes of interest#### SPECIFY, WHICH PARAMETERS TO TRACE:
-# parameters <- c("P.ability", "P.Moneyback", "theta", "y.pred")
-# samplesSet(parameters)
-# 
-# # Generate 100000 iterations
-# modelUpdate(100000)
-# #### SHOW POSTERIOR STATISTICS
-# sample.statistics <- samplesStats("*")
-# print(sample.statistics)
      +

      1. 第一個問題，可以用 theta 的事後概率分布來回答，十萬次MC計算的結果顯示，邦德能夠準確分辨馬丁尼酒的調制方法的概率是0.78，且這個事後概率分布的95%可信區間是(0.564, 0.933)。

      2. 第二個問題，邦德擁有能準確分辨馬丁尼酒調制方法的能力(不是猜的)的概率是 P.ability = 0.993。如果邦德只是瞎猜，那麼 theta 就只能等於或者十分接近0.5。所以我們相信邦德有這樣一種分辨能力的概率是99.3%。

        3. @@ -2523,7 +2549,7 @@

        83.6 Practical Bayesian Statistic (function () { var script = document.createElement("script"); script.type = "text/javascript"; - var src = ""; + var src = "true"; if (src === "" || src === "true") src = "https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-MML-AM_CHTML"; if (location.protocol !== "file:") if (/^https?:/.test(src)) diff --git "a/docs/\345\206\215\350\250\252-mcmc.html" "b/docs/\345\206\215\350\250\252-mcmc.html" index cb07d8eb0..5c000409b 100644 --- "a/docs/\345\206\215\350\250\252-mcmc.html" +++ "b/docs/\345\206\215\350\250\252-mcmc.html" @@ -6,7 +6,7 @@ 第 89 章 再訪 MCMC | 醫學統計學 - + @@ -24,7 +24,7 @@ - + @@ -61,41 +61,66 @@ @@ -575,7 +600,7 @@
      3. 38.1 定義
      4. 39 The sandwich estimator
        5. -
      5. VII 貝葉斯統計
        6. +
      6. VII 貝葉斯統計 Introduction to Bayesian Statistics
      7. 40 貝葉斯統計入門
      8. 39 The sandwich estimator
        9. -
      9. VII 貝葉斯統計
        10. +
      10. VII 貝葉斯統計 Introduction to Bayesian Statistics
      11. 40 貝葉斯統計入門

      -

      事實上，上面的現象在使用邏輯迴歸的時候基本上都會呈現。在經典論文 (L. D. Robinson and Jewell 1991) 中給出了詳細的論證。所以其實使用邏輯迴歸擬合數據的 RCT 臨牀試驗，我們可以推論，當模型中加入第三個僅和結果變量有關的基線共變量 (baseline covariates)，如果模型估計的對數比值比在調整前後變化不大 (即，不可壓縮性造成的影響很小)，那這樣的調整對於改善分析的統計效能上幾乎也沒有貢獻。(跟使用線性迴歸的 RCT 完全不同！)

      +

      事實上，上面的現象在使用邏輯迴歸的時候基本上都會呈現。在經典論文 (Robinson and Jewell 1991) 中給出了詳細的論證。所以其實使用邏輯迴歸擬合數據的 RCT 臨牀試驗，我們可以推論，當模型中加入第三個僅和結果變量有關的基線共變量 (baseline covariates)，如果模型估計的對數比值比在調整前後變化不大 (即，不可壓縮性造成的影響很小)，那這樣的調整對於改善分析的統計效能上幾乎也沒有貢獻。(跟使用線性迴歸的 RCT 完全不同！)

      由於邏輯迴歸受使用 \(\text{logit}\) 鏈接方程時不可壓縮性的侷限，同時還由於使用 \(\text{log}\) 鏈接方程時獲得的危險度比 (risk ratios) 比比值比 (odds ratios) 更加容易讓人理解，結果變量爲二分類的 RCT 臨牀試驗常常會選用 \(\text{log}\) 鏈接方程的廣義線性迴歸模型 (見 Section 45.3 第 5 條討論)。選用 \(\text{log}\) 鏈接方程的 GLM 最大的問題在於，當模型中加入過多的預測變量時，會導致模型無法收斂 (converge)–無法找到極大似然估計。

      至於使用泊松迴歸模型的時候，預測變量如果放入不合理，那麼很容易違反泊松分佈的前提 (方差和均值相同)。對於違反了泊松分佈前提，模型變得過度離散 (over-dispersed) 的 GLM，加入適當的基線共變量 (baseline covariates) 則有助於減少模型的過度離散，減小參數估計的標準誤 (使之變得更精確些)。和線性迴歸相同的是，泊松迴歸模型不受不可壓縮性 (non-collapsibility) 的影響。

      @@ -1585,13 +1602,13 @@

      51.3.1 不成熟的小策略

      • 第一步，先看看暴露變量和結果變量之間的關係
      -
      lbw <- read_dta("http://www.stata-press.com/data/r12/lbw.dta")
-lbw$race <- factor(lbw$race)
-lbw$smoke <- factor(lbw$smoke)
-lbw$ht <- factor(lbw$ht)
-a <- Epi::stat.table(list("Birthweight <2500g" = low, "History of hypertension"=ht), list(count(),percent(low)), data = lbw, margins = TRUE)
-# We first tabulate the data
-print(a, digits = c(percent = 2))
      + 
      ##  -------------------------------------- 
 ##               -History of hypertension- 
 ##  Birthweight         0       1   Total  
@@ -1610,8 +1627,8 @@ 
      51.3.1 不成熟的小策略
      
 
        
 
      • 第二步，分析母親高血壓病史和嬰兒低出生體重之間的調整前 (粗) 關係。
        • 
 
      
-
      Model0 <- glm(low~ht, data = lbw, family = binomial(link = "logit"))
-summary(Model0); epiDisplay::logistic.display(Model0)
      
+
 
      ## 
 ## Call:
 ## glm(formula = low ~ ht, family = binomial(link = "logit"), data = lbw)
@@ -1647,9 +1664,9 @@ 
      51.3.1 不成熟的小策略
      
 
        
 
      • 接下來，分析潛在的混雜因子是否和主要暴露變量相關：
        • 
 
      
-
      # lwt is the last weight of mothers before pregnancy
-Model1 <- lm(lwt ~ ht, data = lbw)
-summary(Model1); epiDisplay::regress.display(Model1)
      
+
 
      ## 
 ## Call:
 ## lm(formula = lwt ~ ht, data = lbw)
@@ -1675,8 +1692,8 @@ 
      51.3.1 不成熟的小策略
      
 ##                                                    
 ## No. of observations = 189
      
 
      可見，有高血壓病史的母親，孕前體重較高。再看其與結果變量是否有關係：
      
-
      Model2 <- glm(low ~ lwt, data = lbw, family = binomial(link = "logit"))
-summary(Model2); epiDisplay::logistic.display(Model2)
      
+
 
      ## 
 ## Call:
 ## glm(formula = low ~ lwt, family = binomial(link = "logit"), data = lbw)
@@ -1709,8 +1726,8 @@ 
      51.3.1 不成熟的小策略
      
 ## No. of observations = 189
 ## AIC value = 232.7081
      
 
      由此知，母親孕前體重較高的人，有較低的可能剩下低出生體重的嬰兒。這兩個單獨的關係，各自看都具有 5% 的統計學意義，但是這 (或者其他變量分析的結果沒有統計學意義時) 並不是決定模型中是否加入母親孕前體重這一潛在的混雜因子的理由。接下來，我們通過模型中加入母親孕前體重這一變量前後模型的參數估計變化來分析：
      
-
      Model3 <- glm(low ~ ht + lwt, data = lbw, family = binomial(link = "logit"))
-summary(Model3);epiDisplay::logistic.display(Model3)
      
+
 
      ## 
 ## Call:
 ## glm(formula = low ~ ht + lwt, family = binomial(link = "logit"), 
@@ -1754,8 +1771,8 @@ 
      51.3.1 不成熟的小策略
      
 
      上面的分析結果，告訴我們，數據提供了足夠的證據證明母親孕前體重和是否有高血壓病史，在調整了彼此之後，仍然獨立地和嬰兒低出生體重的發生有相關性。這裏，我們可以下結論認爲，模型中加入母親孕前體重作爲混雜因子，是合情合理的。
      
 
      完成了目前爲止的初步分析和混雜因子的判斷以後，下一階段的分析側重於尋找有沒有任何第三方的預測變量，會對主要暴露變量 \(X_1\) (孕期高血壓) 與結果變量 \(Y\) (嬰兒出生體重過低) 之間的關係產生交互作用。如果數據中的預測變量有多個，那可能導致需要分析潛在的交互作用有許多對，通常建議在遇到多個預測變量之間的複雜關係需要討論的時候，建議不要一股腦全部作交互作用的分析，而是限定一個或者幾個最有可能有交互作用的變量就可以了。否則模型過於複雜，反而不利於理解。一般生物醫學的統計分析中考慮的重要交互作用分析，需要有重要的生物學意義，常見的例子是年齡，性別等。
      
 
      本節使用的例子中，令人感興趣的是，母親的孕前體重，會不會對妊娠高血壓的有無與嬰兒出生體重過低之間的關係造成交互作用：
      
-
      Model4 <- glm(low ~ ht*lwt, family = binomial(link = "logit"), data = lbw)
-summary(Model4); epiDisplay::logistic.display(Model4)
      
+
 
      ## 
 ## Call:
 ## glm(formula = low ~ ht * lwt, family = binomial(link = "logit"), 
@@ -1812,11 +1829,11 @@ 
      51.3.2 補充
      
 
      
 
      51.4 分析目的 2 和 3 – 建立預測模型 (predictive models)
      
 
      建立預測模型的過程，其實就是選擇哪個或者那些變量進入模型的過程。方法有很多，可惜的是，沒有哪種是公認完美的。這裏只介紹兩種最常見，也最常被批評的方法 – 前/後 逐步選擇法 (forward stepwise selection/backward elimination)。強調一下，逐步法本身並不是神奇法術，不同的算法選擇的變量自然會有不同，如果你用了逐步選擇法，選出來的模型變量僅僅只能作爲參考，而不能作爲最終結論。
      
-
      vitc <- haven::read_dta("backupfiles/vitC.dta")
-vitc$ctakers <- factor(vitc$ctakers)
-vitc$sex <- factor(vitc$sex)
-
-stats::step(lm(seruvitc~1,data=vitc[complete.cases(vitc),]),direction="forward",scope=~age+height+weight+sex+cigs+ctakers)
      
+
 
      ## Start:  AIC=575.43
 ## seruvitc ~ 1
 ## 
@@ -1866,7 +1883,7 @@ 
      51.4 分析目的 2 和 3 – 建
 ## Coefficients:
 ## (Intercept)     ctakers1         sex1         cigs  
 ##     46.0283      19.8317       9.7642     -12.2550
      
-
      stats::step(lm(seruvitc~.,data=vitc[complete.cases(vitc),]),direction="backward")
      
+
 
      ## Start:  AIC=563.38
 ## seruvitc ~ serial + age + height + cigs + weight + sex + ctakers
 ## 
@@ -1998,7 +2015,7 @@ 
      References
      
   (function () {
     var script = document.createElement("script");
     script.type = "text/javascript";
-    var src = "";
+    var src = "true";
     if (src === "" || src === "true") src = "https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-MML-AM_CHTML";
     if (location.protocol !== "file:")
       if (/^https?:/.test(src))
diff --git "a/docs/\345\210\206\346\236\220\347\255\226\347\225\245\345\222\214\346\250\241\345\236\213\346\252\242\346\237\245-model-checking-survival-analysis.html" "b/docs/\345\210\206\346\236\220\347\255\226\347\225\245\345\222\214\346\250\241\345\236\213\346\252\242\346\237\245-model-checking-survival-analysis.html"
index 8da8a22a5..67951f28a 100644
--- "a/docs/\345\210\206\346\236\220\347\255\226\347\225\245\345\222\214\346\250\241\345\236\213\346\252\242\346\237\245-model-checking-survival-analysis.html"
+++ "b/docs/\345\210\206\346\236\220\347\255\226\347\225\245\345\222\214\346\250\241\345\236\213\346\252\242\346\237\245-model-checking-survival-analysis.html"
@@ -6,7 +6,7 @@
   
   第 75 章 分析策略和模型檢查 Model checking-survival analysis | 醫學統計學
   
-  
+  
 
   
   
@@ -24,7 +24,7 @@
 
 
 
-
+
 
   
   
@@ -61,41 +61,66 @@
 
 
 
 
 
@@ -575,7 +600,7 @@
 
    4. 38.1 定義
      5. 
 
 
    5. 39 The sandwich estimator
      6. 
-
    6. VII 貝葉斯統計
      7. 
+
    7. VII 貝葉斯統計 Introduction to Bayesian Statistics
      8. 
 
    8. 40 貝葉斯統計入門
      9. 
 
    9. 39 The sandwich estimator
      10. 
-
    10. VII 貝葉斯統計
      11. 
+
    11. VII 貝葉斯統計 Introduction to Bayesian Statistics
      12. 
 
    12. 40 貝葉斯統計入門

    -
    set.seed(1234)
-Normal <- rnorm(2500, mean = 120, sd = 8)
-h <- hist(Normal,breaks = 20, col = "lightblue", xlab = "some value" ,
-          ylim = c(0,300))
-xfit<-seq(min(Normal),max(Normal),length=40)
-yfit<-dnorm(xfit,mean=mean(Normal),sd=sd(Normal))
-yfit <- yfit*diff(h$mids[1:2])*length(Normal)
-lines(xfit, yfit, col="blue", lwd=2)
    +
    Appearance of histogram with normal curve

    圖 25.1: Appearance of histogram with normal curve

    -
    qqnorm(Normal,frame=F); qqline(Normal)
    +
    Appearance of normal plot for a normally distributed variable

    @@ -1426,46 +1450,46 @@

    25.4 數學冪轉換 power transf \cdots,x^{-2},x^{-1},x^{-\frac{1}{2}},\text{log}(x),x^{\frac{1}{2}},x^1,x^2,\cdots \]

    上面舉例的數學冪轉換方法，都是常見的手段用於降低原始數據的偏度 (skewness)，相反地，冪轉換卻不一定能夠改變數據的峯度 (kurtosis)。下面的方程，(非常的羅嗦的方程 sorry)，用於實施類似 ladder 在 Stata 中的效果，即對數據進行各種轉換，然後輸出每種冪轉換後的數據是否爲正態分佈的檢驗結果 (使用 shapiro.test())：

    -
    Ladder.x <- function(x){
-    data <- data.frame(x^3,x^2,x,sqrt(x),log(x),1/sqrt(x),1/x,1/(x^2),1/(x^3))
-    names(data) <- c("cubic","square","identity","square root","log","1/(square root)",
-                     "inverse","1/square","1/cubic")
-   # options(scipen=5)
-    test1 <- shapiro.test(data$cubic)
-    test2 <- shapiro.test(data$square)
-    test3 <- shapiro.test(data$identity)
-    test4 <- shapiro.test(data$`square root`)
-    test5 <- shapiro.test(data$log)
-    test6 <- shapiro.test(data$`1/(square root)`)
-    test7 <- shapiro.test(data$inverse)
-    test8 <- shapiro.test(data$`1/square`)
-    test9 <- shapiro.test(data$`1/cubic`)
-    W.statistic <- c(test1$statistic,
-                     test2$statistic,
-                     test3$statistic,
-                     test4$statistic,
-                     test5$statistic,
-                     test6$statistic,
-                     test7$statistic,
-                     test8$statistic,
-                     test9$statistic)
-    p.value <- c(test1$p.value,
-                 test2$p.value,
-                 test3$p.value,
-                 test4$p.value,
-                 test5$p.value,
-                 test6$p.value,
-                 test7$p.value,
-                 test8$p.value,
-                 test9$p.value)
-    Hmisc::format.pval(p.value ,digits=5, eps = 0.00001, scientific = FALSE)
-    Transformation <- c("cubic","square","identity","square root","log","1/(square root)",
-                        "inverse","1/square","1/cubic")
-    Formula <- c("x^3","x^2","x","sqrt(x)","log(x)","1/sqrt(x)","1/x","1/(x^2)","1/(x^3)")
-    (results <- data.frame(Transformation, Formula, W.statistic, p.value))
-  }
    -
    Normal <- rnorm(2500, mean = 120, sd = 8)
-Ladder.x(Normal)
    + + 
    ##    Transformation   Formula W.statistic   p.value
 ## 1           cubic       x^3      0.9901 3.906e-12
 ## 2          square       x^2      0.9968 4.283e-05
@@ -1608,7 +1632,7 @@ 
    25.4.4 百分比的轉換
    
 
    References
    
 
    
 
    
-
    Belle, G. van, L.D. Fisher, P.J. Heagerty, and T. Lumley. 2004. Biostatistics: A Methodology for the Health Sciences. Wiley Series in Probability and Statistics. Wiley. https://books.google.co.uk/books?id=KSh8IOrLPzwC.
    
+
    Belle, G. van, L.D. Fisher, P.J. Heagerty, and T. Lumley. 2004. Biostatistics: A Methodology for the Health Sciences. Wiley Series in Probability and Statistics. Wiley. https://books.google.co.uk/books?id=KSh8IOrLPzwC.
    
 
    
 
    
             
@@ -1672,7 +1696,7 @@ 
    References
    
   (function () {
     var script = document.createElement("script");
     script.type = "text/javascript";
-    var src = "";
+    var src = "true";
     if (src === "" || src === "true") src = "https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-MML-AM_CHTML";
     if (location.protocol !== "file:")
       if (/^https?:/.test(src))
diff --git "a/docs/\345\237\272\346\226\274\347\247\251\346\254\241\347\232\204\351\235\236\345\217\203\346\225\270\346\252\242\351\251\227.html" "b/docs/\345\237\272\346\226\274\347\247\251\346\254\241\347\232\204\351\235\236\345\217\203\346\225\270\346\252\242\351\251\227.html"
index 8e1af826b..25bf674ab 100644
--- "a/docs/\345\237\272\346\226\274\347\247\251\346\254\241\347\232\204\351\235\236\345\217\203\346\225\270\346\252\242\351\251\227.html"
+++ "b/docs/\345\237\272\346\226\274\347\247\251\346\254\241\347\232\204\351\235\236\345\217\203\346\225\270\346\252\242\351\251\227.html"
@@ -6,7 +6,7 @@
   
   第 36 章 基於秩次的非參數檢驗 | 醫學統計學
   
-  
+  
 
   
   
@@ -24,7 +24,7 @@
 
 
 
-
+
 
   
   
@@ -61,41 +61,66 @@
 
 
 
 
 
@@ -575,7 +600,7 @@
 
  • 38.1 定義
    •

  • 39 The sandwich estimator
    • -
  • VII 貝葉斯統計
    • +
  • VII 貝葉斯統計 Introduction to Bayesian Statistics
  • 40 貝葉斯統計入門