Rでメタ・アナリシス / Reproduction of Tango (2016)
Sho Kuroda
July 2017
丹後『新版 メタ・アナリシス入門』（朝倉書店，2016）の再現．丹後先生のサイトでR/S+の分析プログラムが公開されている．
Standard approach for OR (Tango 2016, 3.1)
まずはオッズ比 (OR = odds ratio)．前提とするクロス表のイメージと変数の定義（丹後 2016 表3.1）．
cross_tab <- rbind(c("ai", "bi", "n1i"), c("ci", "di", "n0i"), c("m1i", "m0i", "ni"))
rownames(cross_tab) <- c("treatment", "control", "sum")
colnames(cross_tab) <- c("event", "no-event", "sample size")
print(cross_tab, quote = F)## event no-event sample size
## treatment ai bi n1i
## control ci di n0i
## sum m1i m0i niData: beta blockade
Data source: Yusuf et al. (1985, PCD) Table 10 … 心筋梗塞の2次予防を目的としたβブロッカー投与のRCT
- ai = number of event in treatment
- n1i = sample size of treatment (= ai + bi)
- ci = number of event in control
- n2i = sample size of control (= ci + di) ※ cross-tabの「n0i」に相当
library(metafor)## Warning: package 'metafor' was built under R version 3.3.3## Loading required package: Matrix## Loading 'metafor' package (version 2.0-0). For an overview
## and introduction to the package please type: help(metafor).data(dat.yusuf1985, package = "metafor")
# str(dat.yusuf1985)
yusuf1985t10 <- dat.yusuf1985[dat.yusuf1985$table == 10, ]
# str(yusuf1985t10)ちなみに1つのstudyに対するORと信頼区間は，たとえば以下のように計算される（後述する variance-based method に相当する）．
yusuf1985t10[1, ] # first study/case## table id trial ai n1i ci n2i
## 51 10 5.1 Reynolds (1 y) 3 38 3 39OR1 <- ((3/38) / (1 - 3/38)) / ((3/39) / (1 - 3/39)) # odds ratio
SE1 <- sqrt(1/3 + 1/35 + 1/3 + 1/36) # standard error
OR1; exp(c(log(OR1) - 1.96 * SE1, log(OR1) + 1.96 * SE1)) # 95% CI## [1] 1.028571## [1] 0.194286 5.445369丹後（2016）表3.2に倣ってYusufの最初の15caseに，さらに2case追加．
tango_t32 <- rbind(yusuf1985t10[1:15, ], c(0, 0, 0, 65, 1195, 62, 1200), c(0, 0, 0, 17, 298, 34, 309))
tango_t32$trial[16:17] <- c("LIT", "Boissel") # trial
tail(tango_t32, 4)## table id trial ai n1i ci n2i
## 64 10 5.14 EIS (1 y) 57 858 45 883
## 65 10 5.15 Rehnqvist (3 y) 25 154 31 147
## 16 0 0 LIT 65 1195 62 1200
## 17 0 0 Boissel 17 298 34 309Peto’s method for OR
頁68fのアルゴリズム3.1.
\[ \ln \hat{OR}_i = \frac{O_i - E_i}{V_i}, \quad w_i = V_i, \quad \hat{OR} = \exp(\frac{\sum w_i \ln \hat{OR}_i}{\sum w_i}) = \exp(\frac{\sum (O_i - E_i)}{\sum V_i}) \]
algo31 <- function(data) {
ai <- data$ai # event
n1i <- data$n1i # sample size of treatment
ci <- data$ci # event
n0i <- data$n2i # sample size of control
m1i <- data$ai + data$ci # total # of event
ni <- data$n1i + data$n2i # N
m0i <- ni - m1i
# step 1: approximation of log OR
Ei <- n1i * m1i / ni
Vi <- n1i * n0i * m1i * m0i / (ni^2 * (ni - 1))
logORi <- (ai - Ei) / Vi
# step 2: standard error
# SEi <- 1 / sqrt(Vi) # not used
# step 3: weight
# wi <- Vi # not used
# step 4: combined OR
logOR <- sum(ai - Ei) / sum(Vi)
OR <- exp(logOR)
# step 5: 95% CI
disp_CI95 <- 1.96 * sqrt(1 / sum(Vi))
CI95 <- exp(c(logOR - disp_CI95, logOR + disp_CI95))
# step 6: homogeneity test
Q1 <- sum((ai - Ei)^2 / Vi) - (sum(ai - Ei))^2 / sum(Vi)
Q1df <- nrow(data) - 1
Q1_prob <- 1 - pchisq(q = Q1, df = Q1df)
# step 7: Mantel-Haenszel test for signif.
Q2 <- (abs(sum(ai - Ei)) - .5)^2 / sum(Vi)
Q2_prob <- 1 - pchisq(q = Q2, df = 1)
# return(list = c(OR = OR, CI95 = CI95, Q1 = Q1, Q1df = Q1df, Q1p = Q1_prob, Q2 = Q2, Q2p = Q2_prob))
return(list(OR = OR, CI = CI95, Q1 = c(Q1, Q1df, Q1_prob), Q2 = c(Q2, Q2_prob)))
}
algo31(data = tango_t32)## $OR
## [1] 0.7819482
##
## $CI
## [1] 0.7063996 0.8655766
##
## $Q1
## [1] 21.3816140 16.0000000 0.1643038
##
## $Q2
## [1] 2.226670e+01 2.372857e-06metafor::rma.petoで答え合わせ．
# rma.peto(ai = tango_t32$ai, bi = tango_t32$n1i - tango_t32$ai, ci = tango_t32$ci, di = tango_t32$n2i - tango_t32$ci, n1i = tango_t32$n1i, n2i = tango_t32$n2i)
tango_t32_peto <- rma.peto(ai = ai, bi = n1i - ai, ci = ci, di = n2i - ci, n1i = n1i, n2i = n2i, data = tango_t32)
summary(tango_t32_peto) # estimate は log(OR) であることに注意##
## Fixed-Effects Model (k = 17)
##
## logLik deviance AIC BIC AICc
## -4.2770 21.3816 10.5540 11.3873 10.8207
##
## Test for Heterogeneity:
## Q(df = 16) = 21.3816, p-val = 0.1643
##
## Model Results (log scale):
##
## estimate se zval pval ci.lb ci.ub
## -0.2460 0.0518 -4.7447 <.0001 -0.3476 -0.1444
##
## Model Results (OR scale):
##
## estimate ci.lb ci.ub
## 0.7819 0.7064 0.8656exp(tango_t32_peto$beta) # OR## intrcpt
## 0.7819482forest(tango_t32_peto, slab = tango_t32$trial, atransf = exp) # forest plot
funnel(tango_t32_peto) # funnel plot
Variance-based method for OR
頁73fの漸近分散法．
\[ \ln \hat{OR}_i = \ln \frac{a_i d_i}{b_i c_i}, \quad w_i = (\frac{1}{a_i} + \frac{1}{b_i} + \frac{1}{c_i} + \frac{1}{d_i})^{-1}, \quad \hat{OR} = \exp(\frac{\sum w_i \ln \hat{OR}_i}{\sum w_i}) \]
algo32 <- function(data) {
ai <- data$ai # event
n1i <- data$n1i # sample size of treatment
bi <- n1i - ai
ci <- data$ci # event
n0i <- data$n2i # sample size of control
di <- n0i - ci
logORi <- log(ai * di / (bi * ci))
SEi <- sqrt(1/ai + 1/bi + 1/ci + 1/di)
wi <- 1 / SEi^2
OR <- exp(sum(wi * logORi) / sum(wi))
disp_CI95 <- 1.96 * sqrt(1 / sum(wi))
CI95 <- exp(c(log(OR) - disp_CI95, log(OR) + disp_CI95))
Q1 <- sum(wi * (logORi - log(OR))^2)
Q1_prob <- 1 - pchisq(q = Q1, df = nrow(data) - 1)
Q2 <- (log(OR))^2 * sum(wi)
Q2_prob <- 1 - pchisq(q = Q2, df = 1)
list(OR = OR, CI = CI95, Q1 = c(Q1, Q1_prob), Q2 = c(Q2, Q2_prob))
}
algo32(data = tango_t32)## $OR
## [1] 0.783085
##
## $CI
## [1] 0.7067045 0.8677207
##
## $Q1
## [1] 21.479776 0.160797
##
## $Q2
## [1] 2.180634e+01 3.016019e-06答え合わせ（FE = fixed effect を指定）．
tango_t32_es <- escalc(measure = "OR", ai = ai, bi = n1i - ai, ci = ci, di = n2i - ci, data = tango_t32) # add effect size (yi) and sampling variance (vi)
tango_t32_fe <- rma.uni(tango_t32_es, method = "FE")
summary(tango_t32_fe)##
## Fixed-Effects Model (k = 17)
##
## logLik deviance AIC BIC AICc
## -4.6208 21.4798 11.2416 12.0748 11.5082
##
## Test for Heterogeneity:
## Q(df = 16) = 21.4798, p-val = 0.1608
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.2445 0.0524 -4.6697 <.0001 -0.3471 -0.1419 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1exp(tango_t32_fe$beta)## [,1]
## intrcpt 0.783085forest(tango_t32_fe, slab = tango_t32$trial, atransf = exp)
DerSimonian-Laird method for OR
変量効果モデルの代表らしい．アルゴリズムは頁78fより．
\[ w_i^* = \frac{1}{SE_i^2 + \hat{\tau}^2} , \quad \hat{OR} = \exp(\frac{\sum w_i^* \ln \hat{OR}_i}{\sum w_i^*}) \]
algo34 <- function(data) {
ai <- data$ai # event
n1i <- data$n1i # sample size of treatment
bi <- n1i - ai
ci <- data$ci # event
n0i <- data$n2i # sample size of control
di <- n0i - ci
# step 1 and 2
logORi <- log(ai * di / (bi * ci))
SEi <- sqrt(1/ai + 1/bi + 1/ci + 1/di)
wi <- 1 / SEi^2
OR <- exp(sum(wi * logORi) / sum(wi))
Q1 <- sum(wi * (logORi - log(OR))^2)
## -- ここまでは漸近分散法と同じ --
# step 3
tau2 <- (Q1 - (nrow(data) - 1)) / (sum(wi) - (sum(wi^2)) / sum(wi))
tau2 <- max(0, tau2)
# step 4
I2 <- (Q1 - (nrow(data) - 1)) / Q1
# step 5
wi_star <- 1 / (SEi^2 + tau2)
# step 6
OR <- exp(sum(wi_star * logORi) / sum(wi_star))
# step 7
dist_IC95 <- 1.96 * sqrt(1 / sum(wi_star))
CI95 <- exp(c(log(OR) - dist_IC95, log(OR) + dist_IC95))
# step 8
Q2 <- (log(OR))^2 * sum(wi_star)
Q2_prob <- 1 - pchisq(q = Q2, df = 1)
list(tau2 = tau2, I2 = I2, OR = OR, CI = CI95, Q2 = c(Q2, Q2_prob))
}
algo34(data = tango_t32)## $tau2
## [1] 0.01686299
##
## $I2
## [1] 0.2551133
##
## $OR
## [1] 0.7907509
##
## $CI
## [1] 0.6948917 0.8998338
##
## $Q2
## [1] 1.26794e+01 3.69707e-04metafor::rma.uniで答え合わせ（method = “DL” に対応）．
tango_t32_es <- escalc(measure = "OR", ai = ai, bi = n1i - ai, ci = ci, di = n2i - ci, data = tango_t32) # 再掲
tango_t32_re <- rma.uni(tango_t32_es, method = "DL")
summary(tango_t32_re)##
## Random-Effects Model (k = 17; tau^2 estimator: DL)
##
## logLik deviance AIC BIC AICc
## -3.9845 20.2071 11.9689 13.6353 12.8260
##
## tau^2 (estimated amount of total heterogeneity): 0.0169 (SE = 0.0240)
## tau (square root of estimated tau^2 value): 0.1299
## I^2 (total heterogeneity / total variability): 25.51%
## H^2 (total variability / sampling variability): 1.34
##
## Test for Heterogeneity:
## Q(df = 16) = 21.4798, p-val = 0.1608
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.2348 0.0659 -3.5608 0.0004 -0.3640 -0.1055 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1exp(tango_t32_re$beta)## [,1]
## intrcpt 0.7907509forest(tango_t32_re, slab = tango_t32$trial, atransf = exp)
Cumulative meta-analysis
metafor::cumul関数を利用（cumul関数を利用するとforest関数にslabを引数として渡せないようなので，rma関数で統合する際にslabを指定しておく）．
tango_t32_re_c <- rma.uni(tango_t32_es, method = "DL", slab = tango_t32$trial)
forest(cumul(tango_t32_re_c), atransf = exp)
Standard approach for mean diff (Tango 2016, 3.2)
- treatment … n1i=例数; X1i=平均値; s1i^2=不偏分散
- control … n0i=例数; X0i=平均値; s0i^2=不偏分散
- AD = absolute difference; SMD = standardized mean difference
\[ AD = \bar{X}_{1i} - \bar{X}_{0i} \] \[ SMD = \frac{\bar{X}_{1i} - \bar{X}_{0i}}{s_i}, \quad s_i^2 = \frac{(n_{1i} - 1) s_{1i}^2 + (n_{0i} - 1) s_{0i}^2}{n_{1i} + n_{0i} - 2} \]
Data: length of hospital stay
脳卒中患者の入院日数（特別なケアの有無別）．丹後（2016）表3.4より．出所については Normand (1999, Stat Med) Table II を参照（Cochran Database of Systematic Reviews 1995 よりデータを作成している模様）．
tango_t34 <- data.frame(
t_n = c(155, 31, 75, 18, 8, 57, 34, 110, 60), # number of treatment group
t_m = c(55, 27, 64, 66, 14, 19, 52, 21, 30), # mean
t_sd = c(47, 7, 17, 20, 8, 7, 45, 16, 27), # SD
c_n = c(156, 32, 71, 18, 13, 52, 33, 183, 52), # number (sample size) of control group
c_m = c(75, 29, 119, 137, 18, 18, 41, 31, 23), # mean
c_sd = c(64, 4, 29, 48, 11, 4, 34, 27, 20) # SD
)
rownames(tango_t34) <- c("Ed", "OMi", "OMo", "OS", "MH", "MT", "N", "Um", "Up")
head(tango_t34, 2)## t_n t_m t_sd c_n c_m c_sd
## Ed 155 55 47 156 75 64
## OMi 31 27 7 32 29 4FE model for mean diff
頁90f．
\[ \hat{AD}_i = \bar{X}_{1i} - \bar{X}_{0i}, \quad w_i = \frac{1}{SE_i^2} = ((\frac{1}{n_{1i}} + \frac{1}{n_{0i}}) s_i^2)^{-1}, \quad \hat{AD} = \frac{\sum w_i \hat{AD}_i}{\sum w_i} \]
algo39 <- function(data) {
n1i <- data$t_n # number of group
n0i <- data$c_n
x1i <- data$t_m # mean
x0i <- data$c_m
s1i <- data$t_sd # standard dev.
s0i <- data$c_sd
si2 <- ((n1i - 1) * s1i^2 + (n0i - 1) * s0i^2) / (n1i + n0i - 2) # si^2
ADi <- x1i - x0i
SEi <- sqrt((1/n1i + 1/n0i) * si2)
wi <- 1 / SEi^2
AD <- sum(wi * ADi) / sum(wi)
disp_CI95 <- 1.96 * sqrt(1 / sum(wi))
CI95 <- c(AD - disp_CI95, AD + disp_CI95)
Q1 <- sum(wi * (ADi - AD)^2)
Q1_prob <- 1 - pchisq(q = Q1, df = nrow(data) - 1)
Q2 <- AD^2 * sum(wi)
Q2_prob <- 1 - pchisq(q = Q2, df = 1)
list(AD = AD, CI = CI95, Q1 = c(Q1, Q1_prob), Q2 = c(Q2, Q2_prob))
}
algo39(data = tango_t34)## $AD
## [1] -3.493861
##
## $CI
## [1] -5.026498 -1.961225
##
## $Q1
## [1] 241.059 0.000
##
## $Q2
## [1] 1.996390e+01 7.891828e-06metafor::rmaで確認．MD (mean difference) を指定して効果量を計算．重み付けの定義か何かが異なるようで，上記の丹後本とは結果が少し異なる．
tango_t34_es <- escalc(measure = "MD", n1i = t_n, n2i = c_n, m1i = t_m, m2i = c_m, sd1i = t_sd, sd2i = c_sd, data = tango_t34) # effect size
tango_t34_md <- rma.uni(tango_t34_es, measure = "MD", method = "FE")
summary(tango_t34_md)##
## Fixed-Effects Model (k = 9)
##
## logLik deviance AIC BIC AICc
## -140.0210 238.9158 282.0421 282.2393 282.6135
##
## Test for Heterogeneity:
## Q(df = 8) = 238.9158, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -3.4636 0.7648 -4.5286 <.0001 -4.9626 -1.9646 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ここまで計算してから，データがパッケージ内で用意されていることに気付いたので，そちらも使ってみる．
data(dat.normand1999, package = "metafor")
head(dat.normand1999, 2)## study source n1i m1i sd1i n2i m2i sd2i
## 1 1 Edinburgh 155 55 47 156 75 64
## 2 2 Orpington-Mild 31 27 7 32 29 4normand99_es <- escalc(measure = "MD", m1i = m1i, sd1i = sd1i, n1i = n1i, m2i = m2i, sd2i = sd2i, n2i = n2i, data = dat.normand1999)
normand99_fe <- rma(normand99_es, measure = "MD", method = "FE")
# summary(normand99_fe)
normand99_fe$beta## [,1]
## intrcpt -3.463613forest(normand99_fe, slab = dat.normand1999$source)
DerSimonian-Laird method for mean diff
頁93f．
\[ w_i^* = \frac{1}{SE_i^2 + \hat{\tau}^2}, \quad \hat{AD} = \frac{\sum w_i^* \hat{AD}_i}{\sum w_i^*} \]
normand99_re <- rma(normand99_es, measure = "MD", method = "DL", slab = dat.normand1999$source)
summary(normand99_re)##
## Random-Effects Model (k = 9; tau^2 estimator: DL)
##
## logLik deviance AIC BIC AICc
## -44.5501 47.9739 93.1001 93.4946 95.1001
##
## tau^2 (estimated amount of total heterogeneity): 205.4094 (SE = 179.7183)
## tau (square root of estimated tau^2 value): 14.3321
## I^2 (total heterogeneity / total variability): 96.65%
## H^2 (total variability / sampling variability): 29.86
##
## Test for Heterogeneity:
## Q(df = 8) = 238.9158, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -13.9817 5.1267 -2.7272 0.0064 -24.0299 -3.9336 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1forest(normand99_re)
Meta-regression analysis (Tango 2016, 5.4)
Data: 心疾患リスクのRCT
丹後（2016）表5.6よりデータ入力，出所はThompson (1993, SMMR)．
- trial = 試験 1, 2, …, 25
- int = 介入 intervention … 1=drug, 2=dietary, 3=surgery
- pat = 患者 patient … 1=primary, 2=secondary, 3=mixed
- chol = serum cholesterol reduction
- dur = 追跡期間 duration (in years)
- trt = 治療群 treatment … 0=control, 1=treated
- ihd = 虚血性心疾患 ischaemic heart disease
- subj = 患者数 subject
thompson93 <- data.frame(
trial = as.factor(c(1, rep(1:25, each = 2))),
int = as.factor(c(rep(1, 7), rep(2, 4), 3, 3, rep(1, 6), rep(2, 4), rep(1, 6), rep(2, 4), 1, 1, 2, 2, 1, 1, rep(2, 4), rep(1, 8))),
pat = as.factor(c(rep(2, 3), rep(1, 6), rep(2, 4), 1, 1, rep(2, 4), 1, 1, rep(2, 6), 3, 3, rep(2, 4), 3, 3, rep(2, 8), 3, 3, rep(2, 6))),
chol = c(0, .68, .42, 0, .55, 0, .65, 0, .7, 0, .26, 0, 1.43, 0, .69, 0, .84, 0, .85, 0, .87, 0, 1.13, 0, .68, 0, .59, 0, .49, 0, .95, 0, .57, 0, 1.08, 0, .31, 0, .85, 0, .61, 0, 1.06, 0, .68, 0, 1.35, 0, 1.48, 0, .56),
dur = c(6.2, rep(c(6.2, 5.3, 7.4, 1, 2, 9.7, 5, 5, 3.6, 3.5, 5, 4.5, 5, 1.9, 3.8, 3, 1, 2.5, 5, 1.6, 3.3, 1, 2, 2.7, 1.6), each = 2)),
trt = as.factor(c(0, 1, 1, rep(0:1, 24))),
ihd = c(936, 322, 354, 210, 173, 193, 157, 121, 131, 144, 132, 125, 82, 84, 56, 101, 73, 75, 54, 65, 52, 81, 61, 85, 54, 69, 42, 42, 36, 51, 46, 50, 47, 20, 62, 24, 37, 11, 6, 11, 8, 5, 3, 2, 2, 5, 1, 0, 2, 0, 1),
subj = c(2789, 1119, 1103, 5296, 5331, 1900, 1906, 4516, 4541, 1015, 1018, 417, 421, 2030, 2051, 276, 279, 367, 350, 422, 424, 229, 229, 253, 244, 284, 145, 1129, 1149, 194, 199, 134, 130, 1663, 6582, 237, 221, 72, 71, 52, 28, 30, 60, 30, 88, 94, 94, 52, 94, 29, 23)
)
thompson93$pij <- thompson93$ihd / thompson93$subj # 疾患発生率
head(thompson93)## trial int pat chol dur trt ihd subj pij
## 1 1 1 2 0.00 6.2 0 936 2789 0.33560416
## 2 1 1 2 0.68 6.2 1 322 1119 0.28775693
## 3 1 1 2 0.42 6.2 1 354 1103 0.32094288
## 4 2 1 1 0.00 5.3 0 210 5296 0.03965257
## 5 2 1 1 0.55 5.3 1 173 5331 0.03245170
## 6 3 1 1 0.00 7.4 0 193 1900 0.10157895FE Peto/variance-based method
tango_t55 <- data.frame(
trial = 1:25,
ai = c(676, thompson93$ihd[thompson93$trt == 1][3:26]),
bi = c(2222-676, (thompson93$subj - thompson93$ihd)[thompson93$trt == 1][3:26]),
ci = thompson93$ihd[thompson93$trt == 0],
di = (thompson93$subj - thompson93$ihd)[thompson93$trt == 0],
n1i = c(2222, thompson93$subj[thompson93$trt == 1][3:26]),
n2i = thompson93$subj[thompson93$trt == 0]
)Peto
thompson93_peto <- rma.peto(ai = ai, bi = bi, ci = ci, di = di, n1i = n1i, n2i = n2i, data = tango_t55)
# summary(thompson93_peto)
exp(thompson93_peto$beta) # OR## intrcpt
## 0.8214615forest(thompson93_peto, atransf = exp)
Variance-based
tango_t55_es <- escalc(measure = "OR", ai = ai, bi = bi, ci = ci, di = di, data = tango_t55)
tango_t55_fe <- rma.uni(tango_t55_es, method = "FE")
# summary(tango_t55_fe)
exp(tango_t55_fe$beta) # OR## [,1]
## intrcpt 0.821934# forest(tango_t55_fe, atransf = exp)Logit
\[ \ln \frac{p_{ij}}{1 - p_{ij}} = \mu + TRIAL_i + TRT_j, \quad j = {\rm{control, treated}} \] \[ \hat{OR} = \exp(\hat{TRT}_{\rm{treated}}) \]
以下の結果を見ると，丹後本では（説明されていないものの）サンプルサイズで重み付けした結果が用いられていると思われる．
# Trial固定効果なし
thompson93_logit0 <- glm(pij ~ trt, data = thompson93, family = quasibinomial(link = "logit"))
# summary(thompson93_logit0)
exp(thompson93_logit0$coefficients["trt1"])## trt1
## 0.8610383# Trial固定効果あり
thompson93_logit1 <- glm(pij ~ trial + trt, data = thompson93, family = quasibinomial(link = "logit"))
# summary(thompson93_logit1)
exp(thompson93_logit1$coefficients["trt1"]) # OR## trt1
## 0.7993361# 総患者数によって重み付け
thompson93_logit2 <- glm(pij ~ trial + trt, data = thompson93, weights = subj, family = quasibinomial(link = "logit"))
# summary(thompson93_logit2)
summary(thompson93_logit2)$coef["trt1", ] # このSEを利用するとCIも計算可## Estimate Std. Error t value Pr(>|t|)
## -0.1973087815 0.0459266487 -4.2961719800 0.0002308719exp(thompson93_logit2$coefficients["trt1"]) # OR## trt1
## 0.8209371CHOL変数（コレステロールの減少量）を調整（丹後本ではこのモデルを支持していないが…）．
thompson93_logit3 <- glm(pij ~ trial + trt + chol, data = thompson93, weights = subj, family = quasibinomial(link = "logit"))
# summary(thompson93_logit3)
exp(thompson93_logit3$coefficients["trt1"]) # OR## trt1
## 1.143054Sub-group analysis using logit model
\[ \ln \frac{p_{ijk}}{1 - p_{ijk}} = \mu + TRIAL_i + (TRT * INT)_{jk}, \quad k = 1, 2, 3, \quad \hat{OR} = \exp (\hat{(TRT * INT)}_{\rm{treated}, k}) \]
コレステロール低減を目的としたRCTのサブ・グループ解析（model = pij ~ trial + trt:int をイメージすればよい．ただしcontrolをbaselineにするために I(trt != 1) という面倒なことに…）．
介入の種類別
thompson93_logit4 <- glm(pij ~ trial + I(trt != 1):int, data = thompson93, weights = subj, family = quasibinomial(link = "logit"))
summary(thompson93_logit4)##
## Call:
## glm(formula = pij ~ trial + I(trt != 1):int, family = quasibinomial(link = "logit"),
## data = thompson93, weights = subj)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6692 -0.6732 0.0000 0.6886 1.8666
##
## Coefficients: (3 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.65050 0.04393 -14.807 3.00e-13 ***
## trial2 -2.53105 0.07600 -33.302 < 2e-16 ***
## trial3 -1.53409 0.08049 -19.061 1.38e-15 ***
## trial4 -2.86878 0.10058 -28.524 < 2e-16 ***
## trial5 -1.16601 0.10156 -11.481 5.31e-11 ***
## trial6 -0.19794 0.14185 -1.395 0.176211
## trial7 -2.58200 0.11508 -22.436 < 2e-16 ***
## trial8 -0.02479 0.12179 -0.204 0.840478
## trial9 -0.76279 0.12861 -5.931 4.79e-06 ***
## trial10 -1.14451 0.13942 -8.209 2.75e-08 ***
## trial11 -0.11480 0.14118 -0.813 0.424475
## trial12 -0.19011 0.13193 -1.441 0.163075
## trial13 -0.33022 0.14451 -2.285 0.031852 *
## trial14 -2.58356 0.15036 -17.182 1.29e-14 ***
## trial15 -0.43016 0.15970 -2.694 0.012966 *
## trial16 0.14138 0.17203 0.822 0.419615
## trial17 -3.77832 0.14620 -25.844 < 2e-16 ***
## trial18 -1.18938 0.18339 -6.486 1.28e-06 ***
## trial19 -1.24784 0.32848 -3.799 0.000926 ***
## trial20 -0.49183 0.33570 -1.465 0.156434
## trial21 -1.63076 0.47272 -3.450 0.002179 **
## trial22 -2.53958 0.64330 -3.948 0.000641 ***
## trial23 -2.65730 0.52506 -5.061 4.01e-05 ***
## trial24 -3.48988 0.89937 -3.880 0.000757 ***
## trial25 -3.18961 1.27488 -2.502 0.019911 *
## I(trt != 1)FALSE:int1 -0.22024 0.05116 -4.305 0.000263 ***
## I(trt != 1)TRUE:int1 NA NA NA NA
## I(trt != 1)FALSE:int2 -0.06970 0.08467 -0.823 0.418838
## I(trt != 1)TRUE:int2 NA NA NA NA
## I(trt != 1)FALSE:int3 -0.57084 0.20568 -2.775 0.010763 *
## I(trt != 1)TRUE:int3 NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasibinomial family taken to be 1.592281)
##
## Null deviance: 5809.032 on 50 degrees of freedom
## Residual deviance: 37.537 on 23 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 5exp(c(-0.22024, -0.06970, -0.57084)) # OR## [1] 0.8023262 0.9326736 0.5650506患者のタイプ別
thompson93_logit5 <- glm(pij ~ trial + I(trt != 1):pat, data = thompson93, weights = subj, family = quasibinomial(link = "logit"))
# summary(thompson93_logit5)
exp(c(-0.18448, -0.20193, -0.23241)) # OR## [1] 0.8315366 0.8171521 0.7926211追跡期間別（2年以下，2-5年，5年以上で分割）．
thompson93$dur_f <- cut(thompson93$dur, breaks = c(0, 2, 5, 10))
thompson93_logit6 <- glm(pij ~ trial + I(trt != 1):dur_f, data = thompson93, weights = subj, family = quasibinomial(link = "logit"))
# summary(thompson93_logit6)
exp(c(-0.059419, -0.275191, -0.205175)) # OR## [1] 0.9423119 0.7594271 0.8145048Last updated on July 10, 2017.