# 第 100 章 不同實驗/研究設計時適用的貝葉斯模型

1. $$Y$$ 是要研究的疾病，$$X$$ 是要研究的暴露因素；
2. 給概率 (probabilities)，和患病率 (prevalences)提供均一先驗概率分佈 (uniform priors)。給作用效果的評價指標例如危險度比 (risk ratios, RR)，危險度差 (risk difference, RD)，或者是比值比 (odds ratios, OR)提供中性的，圍繞着無效爲中心的先驗概率分佈 (“neutral” priors centred on no effect)。

## 100.1 隊列研究設計時的貝葉斯模型

model{
Y0 ~ dbin(r0, X0) # Data model for unexposed
Y1 ~ dbin(r1, X1) # Data model for exposed

# Priors for risks of exposed and unexposed
r0 ~ dbeta(a.unexp, b.unexp)
r1 ~ dbeta(a.exp, b.exp)

# Computation of comparison statistics
rd <-  r1 - r0   # for Risk Difference
rr <- r1 / r0    # for Risk Ratio
}

list(X0 = 2000,
X1 = 1000,
Y0 = 100,
Y1 = 150,
a.unexp = 1, b.unexp = 1,
a.exp   = 1, b.exp   = 1)

model{
Y0 ~ dbin(r0, X0) # Data model for unexposed
Y1 ~ dbin(r1, X1) # Data model for exposed

# define r1 as the unexposed risk plus an effect "rd"
r1 <- r0 + rd

# Priors for risks of exposed and unexposed
r0 ~ dbeta(a.unexp, b.unexp)
rd ~ dunif(min.rd, max.rd)

# Computation of comparison statistics
rr <- r1 / r0    # for Risk Ratio
}

list(X0 = 2000,
X1 = 1000,
Y0 = 100,
Y1 = 150,
a.unexp = 1, b.unexp = 1,
min.rd = -30, max.rd = 30)

model{
Y0 ~ dbin(r0, X0)   # Data model for unexposed
logit(r0) <- lr0    # define the logit of r0
Y1 ~ dbin(r1, X1)   # Data model for exposed
logit(r1) <- lr1    # define the logit of r1

# define lr1 as the unexposed logit(risk) plus log(OR):
lr1 <- lr0 + lor

# priors for the logit of risk of unexposed and of log(OR):
lr0 ~ dnorm(mean.lr0, pr.lr0)
lor ~ dnorm(mean.lor, pr.lor)

# computations of comparison statistics:
rd <- r1 - r0      # For risk difference
rr <- r1 / r0      # For risk ratio
}

list(X0 = 2000,
X1 = 1000,
Y0 = 100,
Y1 = 150,
mean.lr0 = 0, pr.lr0 = 0.3,
mean.lor = 0, pr.lor = 0.33)

## 100.2 病例對照研究設計時的貝葉斯模型

model{
X0 ~ dbin(p0, Y0) # data model for CONTROLS
logit(p0) <- lp0  # define the logit of p0
X1 ~ dbin(p1, Y1) # data model for CASES
logit(p1) <- lp1  # define the logit of p1

# define lp1 as the unexposed logit(exposure) plus log(OR)
lp1 <- lp0 + lor

# Priors for logit of exposure and of log(OR)
lp0 ~ dnorm(mean.lp0, pr.lp0)
lor ~ dnorm(mean.lor, pr.lor)

# Computation of comparison statistics
or <- exp(lor)    # For Odds Ratio
}

list(Y0 = 2000,
Y1 = 1000,
X0 = 800,
X1 = 600,
mean.lp0 = 0, pr.lp0 = 0.3,
mean.lor = 0, pr.lor = 0.33)

## 100.3 橫斷面研究設計時的貝葉斯模型

• 25名吸煙者患有癌症；
• 15名非吸煙者患有癌症；
• 150名吸煙者無癌症；
• 剩餘210名非吸煙者無癌症。

$N \sim \text{Multinomial}([\theta_{11}, \theta_{10}, \theta_{01}, \theta_{00}], 400)$

$\mathbf{\theta} \sim \text{Dirichlet}(\alpha_1, \alpha_2, \dots. \alpha_k) \text{ where } \alpha_i >0, \text{ and we denote } \sum_{i = 1}^k \alpha_i = \alpha$

1. 狄利克雷分布中的任意一個分類的(category)邊際先驗概率分布(marginal prior distribution)是一個Beta分布：$$\theta_i \sim \text{Beta}(\alpha_i, \alpha - \alpha_i)$$
2. 狄利克雷分布中的任意一個分類的參數 $$\theta_i$$ 的均值（期望）是：$$E(\theta_i) = \frac{\alpha_i}{\alpha}$$
3. 狄利克雷分布中的任意一個分類的參數 $$\theta_i$$ 的方差是：$$V(\theta_i) = \frac{E(\theta_i)(1-E(\theta_i))}{\alpha + 1}$$

• 狄利克雷分布$$\text{Dirichlet}(1,2,1)$$，其實是樣本量爲$$\alpha = 4$$，均值是 $$\text{Means}[0.25, 0.5, 0.25]$$的分布；
• 狄利克雷分布$$\text{Dirichlet}(1,2,1)$$，其實是樣本量爲$$\alpha = 40$$，均值也是$$\text{Means}[0.25, 0.5, 0.25]$$的分布。

model{
N[1:4] ~ dmulti(p[], S)  # data model for sample
p[1:4] ~ ddirch(alpha[]) # Dirichlet prior for vector of probabilities

# Computation of comparison statistics:
px0 <- p[2] + p[4]       # proportion of non-exposed
px1 <- 1 - px0           # proportion of exposed
r0 <- p[2] / px0         # risk in the non-exposed
r1 <- p[1] / px1         # risk in the exposed
rr <- r1 / r0            # risk ratio, RR
rd <- r1 - r0            # risk difference, RD
or <- (p[1]*p[4]) / (p[2]*p[3])   # odds ratio, OR
p.crit <- step(or - 1)   # =1 if or >= 1, 0 otherwise

# calculate the total sample size
S <- sum(N[])
}

list(N = c(25, 15, 150, 210),
alpha = c(1, 1, 1, 1)) # uniform prior for p_ij

Iterations = 1001:26000
Thinning interval = 1
Number of chains = 2
Sample size per chain = 25000

1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:

Mean       SD  Naive SE Time-series SE
deviance 18.81623 2.453005 1.097e-02      1.105e-02
or        2.43862 0.861082 3.851e-03      3.886e-03
p[1]      0.06436 0.012207 5.459e-05      5.459e-05
p[2]      0.03959 0.009694 4.335e-05      4.335e-05
p[3]      0.37378 0.024005 1.074e-04      1.078e-04
p[4]      0.52227 0.024829 1.110e-04      1.110e-04
p.crit    0.99372 0.078998 3.533e-04      3.630e-04
r0        0.07046 0.016980 7.594e-05      7.593e-05
r1        0.14689 0.026570 1.188e-04      1.188e-04
rd        0.07643 0.031454 1.407e-04      1.407e-04
rr        2.21362 0.703480 3.146e-03      3.172e-03

2. Quantiles for each variable:

2.5%      25%      50%      75%    97.5%
deviance 16.04000 17.02000 18.19000 19.92000 25.11000
or        1.18700  1.83300  2.29700  2.88600  4.52103
p[1]      0.04272  0.05571  0.06368  0.07215  0.09039
p[2]      0.02285  0.03268  0.03880  0.04563  0.06050
p[3]      0.32720  0.35750  0.37360  0.38980  0.42150
p[4]      0.47350  0.50560  0.52240  0.53890  0.57060
p.crit    1.00000  1.00000  1.00000  1.00000  1.00000
r0        0.04094  0.05839  0.06914  0.08114  0.10730
r1        0.09882  0.12820  0.14570  0.16410  0.20270
rd        0.01569  0.05514  0.07586  0.09718  0.13970
rr        1.16700  1.71900  2.10500  2.58700  3.90200

Iterations = 1001:26000
Thinning interval = 1
Number of chains = 2
Sample size per chain = 25000

1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:

Mean       SD Naive SE Time-series SE
deviance 86.75869 10.43878 4.67e-02       4.67e-02
or        1.41233  0.48014 2.15e-03       2.15e-03
p[1]      0.01854  0.00361 1.61e-05       1.61e-05
p[2]      0.01141  0.00283 1.27e-05       1.28e-05
p[3]      0.53491  0.01330 5.95e-05       5.92e-05
p[4]      0.43514  0.01325 5.93e-05       5.96e-05
p.crit    0.81588  0.38759 1.73e-03       1.74e-03
r0        0.02554  0.00630 2.82e-05       2.86e-05
r1        0.03350  0.00646 2.89e-05       2.89e-05
rd        0.00796  0.00901 4.03e-05       4.03e-05
rr        1.39672  0.46134 2.06e-03       2.06e-03

2. Quantiles for each variable:

2.5%      25%      50%     75%    97.5%
deviance 67.64000 79.50000 86.26000 93.5700 108.5000
or        0.70919  1.07600  1.33600  1.6600   2.5680
p[1]      0.01216  0.01598  0.01833  0.0208   0.0262
p[2]      0.00654  0.00940  0.01118  0.0132   0.0176
p[3]      0.50880  0.52600  0.53480  0.5439   0.5611
p[4]      0.40930  0.42620  0.43510  0.4440   0.4613
p.crit    0.00000  1.00000  1.00000  1.0000   1.0000
r0        0.01473  0.02106  0.02504  0.0294   0.0393
r1        0.02205  0.02894  0.03313  0.0376   0.0472
rd       -0.01009  0.00205  0.00804  0.0140   0.0255
rr        0.71649  1.07400  1.32400  1.6360   2.5040

## 100.4 把不同實驗設計的數據用貝葉斯模型連接起來

### 100.4.1 Linking sub-models throug common parameters

model{
# Data model for unexposed group
Y0c ~ dbin(r0, X0c)
logit(r0) <- lr0

# Data model for exposed group
Y1c ~ dbin(r1, X1c)
logit(r1) <- lr1

# lor is log(OR)
lr1 <- lr0 + lor

# Prior for logit of unexposed risk
lr0 ~ dnorm(0, 0.3)
# Prior for log(OR)
lor ~ dnorm(0, 0.33)
}

model{
# Data model for controls
X0cc ~ dbin(p0, Y0cc)
logit(p0) <- lp0

# Data model for cases
X1cc ~ dbin(p1, Y1cc)
logit(p1) <- lp1

# lor is log(OR)
lp1 <- lp0 + lor

# Prior for logit of probability of exposure for controls
lp0 ~ dnorm(0, 0.3)
# Prior for log(OR)
lor ~ dnorm(0, 0.33)
}

model{
# Cohort sub-model
# Data model for unexposed group
Y0c ~ dbin(r0, X0c)
logit(r0) <- lr0

# Data model for exposed group
Y1c ~ dbin(r1, X1c)
logit(r1) <- lr1

# lor is log(OR)
lr1 <- lr0 + lor

# Prior for logit of unexposed risk
lr0 ~ dnorm(0, 0.3)
# Case-control sub-model
# Data model for controls
X0cc ~ dbin(p0, Y0cc)
logit(p0) <- lp0

# Data model for cases
X1cc ~ dbin(p1, Y1cc)
logit(p1) <- lp1

# lor is log(OR)
lp1 <- lp0 + lor

# Prior for logit of probability of exposure for controls
lp0 ~ dnorm(0, 0.3)

# Prior for common log(OR)
lor ~ dnorm(0, 0.33)
}

## 100.5 Practical Bayesian Statistics 06

### 100.5.1 The GREAT Trial

• 治療組：當患者被家庭醫生 (General Practitioners) 發現心肌梗塞時(myocardial infarctin, MI)，在家中立刻給予阿尼普酶藥物，等到患者被救護車帶到醫院以後給患者服用安慰劑。
• 對照組：同樣情況下，家庭醫生在家中發現患者有心肌梗塞時先不給予藥物治療，而是先讓患者服用安慰劑，等救護車帶患者抵達醫院之後再讓患者服用真正的阿尼普酶藥物。

• 治療組163人，13例死亡；
• 對照組148人，23例死亡。

model {
for (i in 1:2) {
deaths[i] ~ dbin(p[i],n[i])
logit(p[i]) <- alpha + beta*treat[i]
}
alpha ~ dunif(-100,100)
beta ~ dunif(-100,100)
OR <- exp(beta)
}

treat[i] 是一個指示变量 (indicator variable)，i = 1 時表示對照組，i = 0 時表示治療組；
death[i] 是第 i 組對象的死亡病例數；
n[i]     是第 i 組對象的總人數；
p[i]     是第 i 組試驗組患者死亡的概率。
• 解釋上面模型語句中 alpha, beta 的含義是什麼。

alpha 是治療組的對數比值 (log odds for the treatment arm)
beta 是對照組和治療組相比較的死亡比值比 (log odds ratio of death in the control arm compared to the treatment arm)，當 beta 大於1時，說明對照組死亡比值高，結果對治療組有利。

• 用兩組起始值文件來跑這個貝葉斯模型。
# Data file for the model
Dat <- list(
n = c(163, 148),
deaths = c(13, 23),
treat = c(0, 1)
)

# initial values for the model
# the choice is arbitrary

inits <- list(
list(alpha = 1, beta = -1),
list(alpha = -1, beta = 1)
)

# fit the model in jags
jagsModel <- jags(data = Dat,
model.file = paste(bugpath,
"/backupfiles/great-model1.txt",
sep = ""),
parameters.to.save = c("OR", "alpha", "beta", "p"),
n.iter = 1000,
n.chains = 2,
inits = inits,
n.burnin = 0,
n.thin = 1,
progress.bar = "none")
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 2
##    Unobserved stochastic nodes: 2
##    Total graph size: 17
##
## Initializing model
# Show the trace plot
Simulated <- coda::as.mcmc(jagsModel)
ggSample <- ggs(Simulated)
ggSample %>%
filter(Parameter %in% c("OR", "alpha", "beta", "p[1]", "p[2]")) %>%
ggs_traceplot()
gelman.diag(Simulated)
## Potential scale reduction factors:
##
##          Point est. Upper C.I.
## OR             1.05       1.13
## alpha          1.02       1.06
## beta           1.02       1.06
## deviance       1.02       1.07
## p[1]           1.01       1.04
## p[2]           1.00       1.01
##
## Multivariate psrf
##
## 1.02
gelman.plot(Simulated)

• 記錄這時候我們獲得的比值比 OR，迴歸係數，及兩組實驗組的事後死亡概率。
jagsModel <- jags(data = Dat,
model.file = paste(bugpath,
"/backupfiles/great-model1.txt",
sep = ""),
parameters.to.save = c("OR", "alpha", "beta", "p"),
n.iter = 26000,
n.chains = 2,
n.burnin = 1000,
n.thin = 1,
inits = inits,
progress.bar = "none")
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 2
##    Unobserved stochastic nodes: 2
##    Total graph size: 17
##
## Initializing model
print(jagsModel)
## Inference for Bugs model at "/Users/chaoimacmini/Downloads/LSHTMlearningnote/backupfiles/great-model1.txt", fit using jags,
##  2 chains, each with 26000 iterations (first 1000 discarded)
##  n.sims = 50000 iterations saved
##          mu.vect sd.vect   2.5%    25%    50%    75%  97.5%  Rhat n.eff
## OR         2.315   0.915  1.058  1.674  2.144  2.760  4.606 1.001 35000
## alpha     -2.480   0.295 -3.101 -2.669 -2.467 -2.277 -1.938 1.001  6000
## beta       0.769   0.373  0.057  0.515  0.763  1.015  1.527 1.001 35000
## p[1]       0.080   0.021  0.043  0.065  0.078  0.093  0.126 1.001  6000
## p[2]       0.155   0.030  0.102  0.135  0.154  0.174  0.218 1.001 16000
## deviance  11.164   2.024  9.197  9.725 10.540 11.938 16.601 1.001 50000
##
## 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).
##
## DIC info (using the rule, pD = var(deviance)/2)
## pD = 2.0 and DIC = 13.2
## DIC is an estimate of expected predictive error (lower deviance is better).
• 現在我們把先驗概率分佈修改成：

\begin{aligned} \alpha & \sim \text{Logistic}(0,1) \\ \beta & \sim \text{Normal}(0, 0.33)\; 0.33 \text{ is the precision}\\ \end{aligned}

jagsModel <- jags(data = Dat,
model.file = paste(bugpath,
"/backupfiles/great-model1_alt.txt",
sep = ""),
parameters.to.save = c("OR", "alpha", "beta", "p"),
n.iter = 1000,
n.chains = 2,
n.burnin = 0,
n.thin = 1,
inits = inits,
progress.bar = "none")
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 2
##    Unobserved stochastic nodes: 2
##    Total graph size: 18
##
## Initializing model
# Show the trace plot
Simulated <- coda::as.mcmc(jagsModel)
ggSample <- ggs(Simulated)
ggSample %>%
filter(Parameter %in% c("OR", "alpha", "beta", "p[1]", "p[2]")) %>%
ggs_traceplot()

## Potential scale reduction factors:
##
##          Point est. Upper C.I.
## OR            1.020      1.078
## alpha         1.025      1.109
## beta          1.014      1.069
## deviance      1.000      1.000
## p[1]          1.021      1.099
## p[2]          0.999      0.999
##
## Multivariate psrf
##
## 1.02
jagsModel <- jags(data = Dat,
model.file = paste(bugpath,
"/backupfiles/great-model1_alt.txt",
sep = ""),
parameters.to.save = c("OR", "alpha", "beta", "p"),
n.iter = 26000,
n.chains = 2,
n.burnin = 1000,
n.thin = 1,
inits = inits,
progress.bar = "none")
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 2
##    Unobserved stochastic nodes: 2
##    Total graph size: 18
##
## Initializing model
print(jagsModel)
## Inference for Bugs model at "/Users/chaoimacmini/Downloads/LSHTMlearningnote/backupfiles/great-model1_alt.txt", fit using jags,
##  2 chains, each with 26000 iterations (first 1000 discarded)
##  n.sims = 50000 iterations saved
##          mu.vect sd.vect   2.5%    25%    50%    75%  97.5%  Rhat n.eff
## OR         2.086   0.779  0.977  1.536  1.950  2.478  3.993 1.001  4700
## alpha     -2.393   0.279 -2.979 -2.575 -2.381 -2.199 -1.879 1.001  5100
## beta       0.671   0.357 -0.023  0.429  0.668  0.907  1.385 1.001  4700
## p[1]       0.086   0.022  0.048  0.071  0.085  0.100  0.133 1.001  5200
## p[2]       0.154   0.029  0.101  0.133  0.153  0.173  0.215 1.001 30000
## deviance  11.129   1.976  9.192  9.718 10.526 11.893 16.403 1.001 38000
##
## 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).
##
## DIC info (using the rule, pD = var(deviance)/2)
## pD = 2.0 and DIC = 13.1
## DIC is an estimate of expected predictive error (lower deviance is better).
• 對改變了先驗概率前後的貝葉斯分析結果，你有什麼看法？

• 爲了分析改變先驗概率分佈前後到底哪種給予了模型更多的信息，我們可以把兩個模型改寫成不含數據，只有預測模型的語句，從而可以看到只有先驗概率分佈時的結果是怎樣的。這時候我們把數據加載這一步省略掉，然後對 alpha, beta, p 進行軌跡檢測，繪製他們的預測分佈密度曲線。
library(runjags)
library(rjags)
library(R2jags)

# Step 1 check model
jagsModel <- jags.model(file = paste(bugpath, "/backupfiles/great-model1_forward.txt", sep = ""),
data = list(treat = c(0, 1)),
n.chains = 2)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 0
##    Unobserved stochastic nodes: 2
##    Total graph size: 13
##
## Initializing model
# Step 2 update 1000 iterations
update(jagsModel, n.iter = 25000, progress.bar="none")

# Step 3 set monitor variables

codaSamples <- coda.samples(
jagsModel,
variable.names = c("alpha", "beta", "p[1]", "p[2]"),
n.iter = 25000, progress.bar="none")

# Show density plot

mcmcplots::denplot(codaSamples)
# ggSample <- ggs(codaSamples)
# ggSample %>%
#   filter(Parameter %in% c("alpha", "beta", "p[1]", "p[2]")) %>%
#   ggs_density() 

# Step 1 check model
jagsModel <- jags.model(file = paste(bugpath, "/backupfiles/great-model1_altforward.txt", sep = ""),
data = list(treat = c(0, 1)),
n.chains = 2)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 0
##    Unobserved stochastic nodes: 2
##    Total graph size: 14
##
## Initializing model
# Step 2 update 1000 iterations
update(jagsModel, n.iter = 25000, progress.bar="none")

# Step 3 set monitor variables

codaSamples <- coda.samples(
jagsModel,
variable.names = c("alpha", "beta", "p[1]", "p[2]"),
n.iter = 25000, progress.bar="none")

# Show density plot

mcmcplots::denplot(codaSamples)

1. 作爲替代方案，我們可以把上述模型重新用對數比值比(Log odds ratio, LOR)作爲未知參數來重新建模：

$LOR = \log(\frac{p[1]/(1 - p[1])}{p[2]/(1 - p[2])}) = \text{logit}(p[1]) - \text{logit}(p[2])$

model {
for (i in 1:2) {
deaths[i] ~ dbin(p[i],n[i])
}
logit(p[1]) <- logit(p[2]) + LOR
p[2] ~ dbeta(1,1)
LOR ~ dnorm(0,0.33)
OR <- exp(LOR)
}
• 再看一遍上面寫好的模型，確定你能夠理解其含義，請確認先驗概率分布的實際意義。

jagsModel <- jags(data = list(n = c(148, 163), deaths = c(23, 13)),
model.file = paste(bugpath,
"/backupfiles/great-model2.txt",
sep = ""),
parameters.to.save = c("OR", "p"),
n.iter = 26000,
n.chains = 2,
n.burnin = 1000,
n.thin = 1,
inits = list(list(LOR = 0.5, p = c(NA, 0.2)),
list(LOR = 5, p = c(NA, 0.8))),
progress.bar = "none")
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 2
##    Unobserved stochastic nodes: 2
##    Total graph size: 13
##
## Initializing model
print(jagsModel)
## Inference for Bugs model at "/Users/chaoimacmini/Downloads/LSHTMlearningnote/backupfiles/great-model2.txt", fit using jags,
##  2 chains, each with 26000 iterations (first 1000 discarded)
##  n.sims = 50000 iterations saved
##          mu.vect sd.vect  2.5%   25%    50%    75%  97.5%  Rhat n.eff
## OR         2.085   0.777 0.988 1.541  1.949  2.476  3.975 1.001  7700
## p[1]       0.154   0.029 0.101 0.133  0.152  0.173  0.215 1.001 21000
## p[2]       0.086   0.021 0.049 0.071  0.084  0.099  0.132 1.001 12000
## deviance  11.119   1.978 9.193 9.716 10.511 11.864 16.425 1.001  5500
##
## 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).
##
## DIC info (using the rule, pD = var(deviance)/2)
## pD = 2.0 and DIC = 13.1
## DIC is an estimate of expected predictive error (lower deviance is better).
• 修改模型的代碼使之用於向前採集樣本用於預測：
model {
#   for (i in 1:2) {
#       deaths[i] ~ dbin(p[i],n[i])
#   }
logit(p[1]) <- logit(p[2]) + LOR
p[2] ~ dbeta(1,1)
LOR ~ dnorm(0,0.33)
OR <- exp(LOR)
}
# Step 1 check model
jagsModel <- jags.model(file = paste(bugpath, "/backupfiles/great-model2_forward.txt", sep = ""),
n.chains = 2)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 0
##    Unobserved stochastic nodes: 2
##    Total graph size: 9
##
## Initializing model
# Step 2 update 1000 iterations
update(jagsModel, n.iter = 25000, progress.bar="none")

# Step 3 set monitor variables

codaSamples <- coda.samples(
jagsModel,
variable.names = c("p[1]", "p[2]"),
n.iter = 25000, progress.bar="none")

# Show density plot

mcmcplots::denplot(codaSamples)

### 100.5.2 吸煙與癌症

#### 100.5.2.1 隊列研究設計

model{
# Data model for non-smokers
Y0c ~ dbin(r0, X0c)
logit(r0) <- lr0

# Data model for smokers
Y1c ~ dbin(r1, X1c)
logit(r1) <- lr1

lr1 <- lr0 + lor   # lor is log(OR)
OR <- exp(lor)  # comparison statistic

# Priors
lr0 ~ dnorm(0, 0.3)  # priors for logit of non-smokers
lor ~ dnorm(0, 0.33)  # prior for log(OR)
}

# Y0c number of non-smokers developed cancer
# X0c number of nonpsmokers
# Y1c number of smokers developed cancer
# X1c number of smokers 

list(X0c = 2000,
Y0c = 100,
X1c = 1000,
Y1c = 150) 

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

# Data file for the model
Dat <- list(
X0c = 2000,
Y0c = 100,
X1c = 1000,
Y1c = 150
)

# initial values for the model
# the choice is arbitrary

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

# fit the model in jags
jagsModel <- jags(data = Dat,
model.file = paste(bugpath,
"/backupfiles/cohort-model.txt",
sep = ""),
parameters.to.save = c("OR", "lor"),
n.iter = 1000,
n.chains = 2,
inits = inits,
n.burnin = 0,
n.thin = 1,
progress.bar = "none")
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 2
##    Unobserved stochastic nodes: 2
##    Total graph size: 13
##
## Initializing model
# Show the trace plot
Simulated <- coda::as.mcmc(jagsModel)
ggSample <- ggs(Simulated)
ggSample %>%
filter(Parameter %in% c("OR", "lor")) %>%
ggs_traceplot()
## Potential scale reduction factors:
##
##          Point est. Upper C.I.
## OR             1.02       1.02
## deviance       1.10       1.30
## lor            1.01       1.01
##
## Multivariate psrf
##
## 1.04
jagsModel <- jags(data = Dat,
model.file = paste(bugpath,
"/backupfiles/cohort-model.txt",
sep = ""),
parameters.to.save = c("OR", "lor"),
n.iter = 26000,
n.chains = 2,
inits = inits,
n.burnin = 1000,
n.thin = 1,
progress.bar = "none")
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 2
##    Unobserved stochastic nodes: 2
##    Total graph size: 13
##
## Initializing model
# Summary Statistics
print(jagsModel)
## Inference for Bugs model at "/Users/chaoimacmini/Downloads/LSHTMlearningnote/backupfiles/cohort-model.txt", fit using jags,
##  2 chains, each with 26000 iterations (first 1000 discarded)
##  n.sims = 50000 iterations saved
##          mu.vect sd.vect   2.5%    25%    50%    75%  97.5%  Rhat n.eff
## OR         3.354   0.870  2.554  3.023  3.310  3.625  4.306 1.008  9800
## lor        1.198   0.143  0.938  1.106  1.197  1.288  1.460 1.001 19000
## deviance  15.252   7.075 13.130 13.666 14.482 15.836 20.431 1.007  3900
##
## 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).
##
## DIC info (using the rule, pD = var(deviance)/2)
## pD = 25.0 and DIC = 40.3
## DIC is an estimate of expected predictive error (lower deviance is better).

#### 100.5.2.2 病例對照研究設計下的模型

• 病例組中600人過去20年中有過吸煙史；
• 對照組中800人過去20年中有過吸煙史。

model{
# Data model for non-cancer controls
X0cc ~ dbin(p0, Y0cc)
logit(p0) <- lp0

# Data model for cancer cases
X1cc ~ dbin(p1, Y1cc)
logit(p1) <- lp1

lp1 <- lp0 + lor # lor is log(OR)
OR <- exp(lor)  # comparison statistic

# Priors
lp0 ~ dnorm(0, 0.3)  # prior for logit of probability of exposure for controls
lor ~ dnorm(0, 0.33) # prior for log(OR)
}

# X0cc indicates number of smokers among controls (without cancer)
# Y0cc indicates number of controls
# X1cc is the number of smokers among cancer cases
# Y1cc is the number of cancer cases 

• X0cc 表示對照組中有吸煙史的人數；
• Y0cc 表示對照組的總人數；
• X1cc 表示癌症患者組中有吸煙史的人數；
• Y1cc 表示癌症患者組的總人數。

list(X0cc = 800,
Y0cc = 2000,
X1cc = 600,
Y1cc = 1000)

# Data file for the model
Dat <- list(X0cc = 800,
Y0cc = 2000,
X1cc = 600,
Y1cc = 1000)

# initial values for the model
# the choice is arbitrary

inits <- list(
list(lp0 = -2, lor = 0),
list(lp0 = 2, lor = 5)
)

# fit the model in jags
jagsModel <- jags(data = Dat,
model.file = paste(bugpath,
"/backupfiles/case-control-model.txt",
sep = ""),
parameters.to.save = c("OR", "lor"),
n.iter = 1000,
n.chains = 2,
inits = inits,
n.burnin = 0,
n.thin = 1,
progress.bar = "none")
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 2
##    Unobserved stochastic nodes: 2
##    Total graph size: 13
##
## Initializing model
# Show the trace plot
Simulated <- coda::as.mcmc(jagsModel)
ggSample <- ggs(Simulated)
ggSample %>%
filter(Parameter %in% c("OR", "lor")) %>%
ggs_traceplot()
## Potential scale reduction factors:
##
##          Point est. Upper C.I.
## OR             1.01       1.06
## deviance       1.03       1.09
## lor            1.02       1.07
##
## Multivariate psrf
##
## 1.03
jagsModel <- jags(data = Dat,
model.file = paste(bugpath,
"/backupfiles/case-control-model.txt",
sep = ""),
parameters.to.save = c("OR", "lor"),
n.iter = 26000,
n.chains = 2,
inits = inits,
n.burnin = 1000,
n.thin = 1,
progress.bar = "none")
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 2
##    Unobserved stochastic nodes: 2
##    Total graph size: 13
##
## Initializing model
print(jagsModel)
## Inference for Bugs model at "/Users/chaoimacmini/Downloads/LSHTMlearningnote/backupfiles/case-control-model.txt", fit using jags,
##  2 chains, each with 26000 iterations (first 1000 discarded)
##  n.sims = 50000 iterations saved
##          mu.vect sd.vect   2.5%    25%    50%    75%  97.5%  Rhat n.eff
## OR         2.280   1.628  1.925  2.131  2.250  2.373  2.635 1.054 28000
## lor        0.812   0.099  0.655  0.756  0.811  0.864  0.969 1.002 50000
## deviance  18.199  52.611 15.381 15.909 16.742 18.138 22.836 1.039 26000
##
## 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).
##
## DIC info (using the rule, pD = var(deviance)/2)
## pD = 1383.8 and DIC = 1402.0
## DIC is an estimate of expected predictive error (lower deviance is better).

#### 100.5.2.3 聯合模型 joint model

# Joint model

model{
# cohort sub-model
Y0c ~ dbin(r0, X0c) # data model for non-smokers
logit(r0) <- lr0
Y1c ~ dbin(r1, X1c) # data model for smokers
logit(r1) <- lr1
lr1 <- lr0 + lor # lor is log(OR)
# prior for cohort sub-model
lr0 ~ dnorm(0, 0.3) # prior for logOdds of nonsmokers

# case-control sub-model
X0cc ~ dbin(p0, Y0cc) # data model for non-cancer controls
logit(p0) <- lp0
X1cc ~ dbin(p1, Y1cc) # data model for cancer cases
logit(p1) <- lp1
lp1 <- lp0 + lor # lor is log(OR)
# prior for case-control sub-model
lp0 ~ dnorm(0, 0.3) # prior for logOdds of exposure for controls

# Common code
lor ~ dnorm(0, 0.33) # prior for common log(OR)
OR <- exp(lor) # comparison statistic
}

list(X0c  = 2000,
Y0c  = 100,
X1c  = 1000,
Y1c  = 150,
X0cc = 800,
Y0cc = 2000,
X1cc = 600,
Y1cc = 1000
)

# Data file for the model
Dat <- list(X0c  = 2000,
Y0c  = 100,
X1c  = 1000,
Y1c  = 150,
X0cc = 800,
Y0cc = 2000,
X1cc = 600,
Y1cc = 1000)

# initial values for the model
# the choice is arbitrary

inits <- list(
list(lr0 = -1, lp0 = -2, lor = 0),
list(lr0 = -5, lp0 = 2, lor = 5)
)

# fit the model in jags
jagsModel <- jags(data = Dat,
model.file = paste(bugpath,
"/backupfiles/jointmodel.txt",
sep = ""),
parameters.to.save = c("OR", "lor"),
n.iter = 1000,
n.chains = 2,
inits = inits,
n.burnin = 0,
n.thin = 1,
progress.bar = "none")
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 4
##    Unobserved stochastic nodes: 3
##    Total graph size: 21
##
## Initializing model
# Show the trace plot
Simulated <- coda::as.mcmc(jagsModel)
ggSample <- ggs(Simulated)
ggSample %>%
filter(Parameter %in% c("OR", "lor")) %>%
ggs_traceplot()
## Potential scale reduction factors:
##
##          Point est. Upper C.I.
## OR             1.01       1.03
## deviance       1.02       1.03
## lor            1.01       1.04
##
## Multivariate psrf
##
## 1.01
jagsModel <- jags(data = Dat,
model.file = paste(bugpath,
"/backupfiles/jointmodel.txt",
sep = ""),
parameters.to.save = c("OR", "lor"),
n.iter = 26000,
n.chains = 2,
inits = inits,
n.burnin = 1000,
n.thin = 1,
progress.bar = "none")
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 4
##    Unobserved stochastic nodes: 3
##    Total graph size: 21
##
## Initializing model
print(jagsModel)
## Inference for Bugs model at "/Users/chaoimacmini/Downloads/LSHTMlearningnote/backupfiles/jointmodel.txt", fit using jags,
##  2 chains, each with 26000 iterations (first 1000 discarded)
##  n.sims = 50000 iterations saved
##          mu.vect sd.vect   2.5%    25%    50%    75%  97.5%  Rhat n.eff
## OR         2.509   1.099  2.178  2.376  2.487  2.603  2.842 1.065 15000
## lor        0.912   0.086  0.779  0.865  0.911  0.957  1.045 1.006 50000
## deviance  38.958  57.785 35.151 36.168 37.317 39.067 44.113 1.082  5000
##
## 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).
##
## DIC info (using the rule, pD = var(deviance)/2)
## pD = 1669.3 and DIC = 1708.2
## DIC is an estimate of expected predictive error (lower deviance is better).

#### 100.5.2.4 你是否能證明兩個研究的比值比是相同的？

Cancer
Y N
Smoking Y S
N NS
C NC

$\text{OR}_{\text{cohort}} = \frac{\text{Pr}(C|S)\times\text{Pr}(NC|NS)}{\text{Pr}(C|NS)\times\text{Pr}(NC|S)}$

Smoking
Y N
Cancer Y C
N NC
S NS

$\text{OR}_{\text{case-control}} = \frac{\text{Pr}(S|C)\times\text{Pr}(NS|NC)}{\text{Pr}(S|NC)\times\text{Pr}(NS|C)}$

\begin{aligned} \text{OR}_{\text{cohort}} & = \frac{\text{Pr}(C|S)\times\text{Pr}(NC|NS)}{\text{Pr}(C|NS)\times\text{Pr}(NC|S)} \\ & = \frac{\frac{\text{Pr}(S|C)\text{Pr}(C)}{\text{Pr}(S)}\times\frac{\text{Pr}(NS|NC)\text{Pr}(NC)}{\text{Pr}(NS)}}{\frac{\text{Pr}(NS|C)\text{Pr}(C)}{\text{Pr}(NS)}\times\frac{\text{Pr}(S|NC)\text{Pr}(NC)}{\text{Pr}(S)}} \\ & = \frac{\text{Pr}(S|C)\times\text{Pr}(NS|NC)}{\text{Pr}(S|NC)\times\text{Pr}(NS|C)} \\ & = \text{OR}_{\text{case-control}} \end{aligned}