第 39 章 混杂的调整，交互作用，和模型的可压缩性
临床医学，流行病学研究的许多问题，需要我们通过数据来评估某些结果变量 (outcome) 和某些预测变量 (predictors/exposures) 之间的关系 (甚至是因果关系)。这些问题的最佳解决方法应该说是随机临床试验 (ramdomized clinical trial, RCT)。但是有更多的时候 (由于违反医学伦理，或者现状所困，甚至是知识有限) 我们无法设计 RCT 来解决这些问题，就只能借助于观察性研究 (observational study)。观察性研究最大的局限性在于无法像 RCT 那样从实验设计阶段把混杂因素排除或者降到最低，所以观察数据在分析的时候，混杂 (confounding) 是必须要加以考虑的一大要因。在简单线性回归章节 (Section 26.5)，详细讨论过混杂因素的定义及条件：
对于一个预测变量是否够格被叫做混杂因子，它必须满足下面的条件：
- 与关心的预测变量相关 (i.e. \(\delta_1 \neq 0\))；
- 与因变量相关 (当关心的预测变量不变时，\(\beta_2\neq0\) )；
- 不在预测变量和因变量的因果关系 (如果有的话) 中作媒介。 Not be on the causal pathway between the predictor of interest and the dependent variable.
下面的统计数据来自一个比较手术和超声碎石术对于肾结石治疗结果的评价。已知大多数医生都公认，肾结石的直径小于 2 公分时治疗成功的概率较高。
| Group | Surgery | Lithotripsy | Surgery | Lithotripsy |
|---|---|---|---|---|
| Success | 81.00 | 234 | 192.00 | 55 |
| Failure | 6.00 | 36 | 71.00 | 25 |
| Total | 87.00 | 270 | 263.00 | 80 |
| Odds Ratios | 2.08 | 1.23 |
可以看到，在上面的分组表格中，左右两边的四格表分别统计了肾结石尺寸小于2 cm 和大于2 cm 时，手术摘除肾结石的成功/失败次数，以及超声碎石术的成功/失败次数。这个表格告诉我们，无论肾结石的尺寸是大于还是小于 2 cm，手术的成功的比值比都大于超声碎石术。但是如果我们没有把数据按照肾结石尺寸区分时，数据就被压缩 (collapsed) 成了下面表格总结的样子：
| Outcome | Surgery | Lithotripsy |
|---|---|---|
| Success | 273 (78%) | 289 (83%) |
| Failure | 77 | 61 |
| Total | 350 | 350 |
| Odds ratio | 0.75 |
当肾结石尺寸被忽略时，数据却显示超声碎石术的成功比值比要高于手术，和之前的结果是矛盾的，你会信哪个？
不要慌，数据不会撒谎，撒谎(有意掩盖事实)的永远是人类。我们来把手术次数，超声碎石术次数，以及肾结石尺寸之间的关系再列个表格：
| Size of the Stone | Surgery | Lithotripsy |
|---|---|---|
| \(< 2\) cm | 87 (33%) | 270 (77%) |
| \(\geqslant 2\) cm | 263 | 80 |
| Total | 350 | 350 |
可见医生也都不是傻子，明明肾结石尺寸很大还要用超声碎石术的人很少，有 77% 的肾结石尺寸小的患者选择了超声碎石术。然后我们再列一个表格来看看肾结石的尺寸大小和超声碎石术成功与否有没有关系：
| Outcome | \(< 2\) cm | \(\geqslant 2\) cm |
|---|---|---|
| Success | 234 (87%) | 55 (69%) |
| Failure | 36 | 25 |
| Total | 370 | 80 |
可见结石尺寸较大时，超声碎石术的成功率 (69%)，明显低于尺寸小的时候的成功率 (87%)。在选择做外科手术的患者中，大多数人的结石尺寸大于 2 cm，成功率仍然达到了 73%。
| Outcome | \(< 2\) cm | \(\geqslant 2\) cm |
|---|---|---|
| Success | 81 (93%) | 192 (73%) |
| Failure | 6 | 71 |
| Total | 87 | 263 |
看到这里，你是否也发现了，肾结石尺寸大小，同时和手术类型的选择有关(尺寸小的倾向于选择超声碎石术)，尺寸大小，同样也和手术结果的成功与否，高度相关(结石小的人成功率高)。所以，分析手术类型和结石手术的成功率的关系的时候，患者体内结石尺寸的大小，对这个关系是产生了混杂效应的。他们三者之间的关系，可以用下面的图 39.1 来理解：
图 39.1: Confounding in kidney stones example
当数据被压缩(忽略了肾结石尺寸时)，比较两种手术类型的成功率的时候，接受外科手术患者的尺寸大多数都较大的事实被“人为的掩盖住了” ，但是当数据按照结石大小分层以后，你会看见外科手术不论是对付大的结石，还是小的结石，成功率都高于超声碎石术。这个例子是混杂效应直接逆转真实相关关系的极佳的实例。同时也提示我们，总体的比值比 (overall odds ratio) 不能通过简单地把分层表格直接压缩 (collapsed table) 获得的数字来计算。
39.1 混杂因素的调整
在前面的肾结石手术的例子中，我们利用已有的背景知识(小尺寸结石的成功率高)，和初步的统计分析(变量之间两两列表分析其内在关系) 发现了混杂因素(结石尺寸)，并且其结果也让我们认定了要暴露因素和结果变量之间的关系，混杂因素必须被调整(adjusted)。
如肾结石数据这样简单的情境下(被认为是混杂因素的变量和预测变量/暴露变量都只是一个二进制变量)，我们可以把变量两两捉对列表分析其相互关系，确定了混杂效应以后把暴露变量和结果变量按照混杂因素的有无列表之后，就可以求得混杂因素被调整后的分层的比值比。这些分层比值比，在暴露变量与结果变量的关系保持混杂因素的层之间保持不变的前提下，可以被“平均化”(简单地说)以后求得总体/合并的比值比(overall /pooled odds ratio)。这就是 Mantel-Haenszel 法或者 Woolf 法的合并比值比的思想出发点。我们来回顾一下 Woolf 法的全部计算过程：
现在假设我们关心的是某种疾病的患病与否(是/否)，和某个暴露变量(是/否) 之间的关系，但是同时想要调整另一个具有\(k\) 个分层的混杂因素变量。
| \(X=0\) | \(X=1\) | ||
|---|---|---|---|
| \(D=0\) | \(X=0\) | \(n_{00}\) | \(n_{10}\) |
| \(X=1\) | \(n_{01}\) | \(n_{10}\) |
39.1.1 Woolf 法估算合并比值比
对想要调整的一个\(k\) 组的混杂因素，把数据按照它的分组分层整理以后，可以得到上面的\(k\) 个四格表(每个分层四格表都是暴露变量和结果变量结合的四个观察值\(a_i, b_i, c_i, d_i, (i=1,\cdots, k)\))。每个分层四格表的观测比值比为 \(\hat\Psi_i = \frac{a_id_i}{c_ib_i}\)，且可以证明方差为
\[ \text{Var}(\text{log}\hat\Psi_i) \approx \frac{1}{a_i} + \frac{1}{b_i} + \frac{1}{c_i} + \frac{1} {d_i} = \frac{1}{w_i} \]
合并比值比的对数 \(\text{log}\hat\Psi_w\) 的 Woolf 的计算方法就是
\[ \text{log}\hat\Psi_w = \frac{\sum w_i\text{log}\hat\Psi_i}{\sum w_i} \]
所以，每个分层的对数比值比乘以自己的方差的倒数 (权重weights) 之后求和，再除以所有权重之和，获得的是合并后的对数比值比，然后再逆运算回来就是Woolf 法计算合并比值比的原理是。
这个合并比值比的对数的方差是
\[ \text{Var}(\text{log}\hat\Psi_w) = \frac{1}{\sum w_i} \]
有了加权后的合并比值比，和方差，就可以求加权后的合并比值比的置信区间了。值得注意的是，分层之后每个分层四格表中的四个数字 (四个样本量) 都不能太小。肾结石手术数据满足这些条件，那么不妨跟我一起手算一下结石尺寸调整后的手术与超声碎石术成功与否的比值比：
\[ \begin{aligned} \hat\Psi_1 = 2.08 ;&\; \hat\Psi_2 = 1.23 \\ \text{Var}(\text{log}\hat\Psi_1) = \frac{1}{81} & + \frac{1}{234} + \frac{1}{6} + \frac{1} {36} = 0.2111 \\ \text{Var}(\text{log}\hat\Psi_2) = \frac{1}{192} & + \frac{1}{55} + \frac{1}{71} + \frac{1} {25} = 0.0775 \\ w_1 = \frac{1}{\text{Var}(\text{log}\hat\Psi_1)} = 4.74 ; \;& w_2 = \frac{1}{\text{Var}(\text{log} \hat\Psi_2)} = 12.91 \\ \text{log}\hat\Psi_w = & \frac{4.74\times\text{log(2.08)} + 12.91\times\text{log(1.23)}}{4.74 + 12.91} \\ = & 0.3481 \\ \Rightarrow \hat\Psi_w =& e^{0.3481} = 1.42\\ \text{Var}(\hat\Psi_w) =& \frac{1}{4.74+12.91} = 0.0567 \\ \Rightarrow 95\% \text{ CI} = & e^{0.3481 \pm 1.96\times\sqrt{0.0567}} \\ = & (0.89, 2.26) \end{aligned} \]
Woolf 的计算调整后的合并比值比的方法是在线性回归和广义线性回归被发现之前诞生的，但是其想法之精妙，确实令人感叹。可惜其最大的缺陷是无法用这样的方法进行连续型变量的调整，也很难同时进行多个变量的调整，所以现在这一算法已经逐渐被淘汰。现在我们有了广义线性回归模型这一更强大的工具，只要把变量加入广义线性模型进行调整就可以计算曾经难以计算和扩展的调整后的合并比值比。从下面的代码计算获得的调整后比值比 \(1.43 (0.91, 2.34)\) 也可以看出，Woolf 方法的计算结果也是足够令人满意的。
size <- c("< 2cm", "< 2cm", ">= 2cm", ">= 2cm")
treatment <- c("Surgery","Lithotripsy","Surgery","Lithotripsy")
n <- c(87, 270, 263, 80)
Success <- c(81, 234, 192, 55)
Stone <- data.frame(size, treatment, n, Success)
ModelStone <- glm(cbind(Success, n - Success) ~ treatment + size,
family = binomial(link = logit), data = Stone)
summary(ModelStone)
Call:
glm(formula = cbind(Success, n - Success) ~ treatment + size,
family = binomial(link = logit), data = Stone)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.9366 0.1704 11.361 < 2e-16 ***
treatmentSurgery 0.3572 0.2291 1.559 0.119
size>= 2cm -1.2606 0.2390 -5.274 1.33e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 33.1239 on 3 degrees of freedom
Residual deviance: 1.0082 on 1 degrees of freedom
AIC: 26.355
Number of Fisher Scoring iterations: 3| term | estimate | std.error | statistic | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|
| (Intercept) | 1.94 | 0.17 | 11.36 | 0.00 | 1.61 | 2.28 |
| treatmentSurgery | 0.36 | 0.23 | 1.56 | 0.12 | -0.09 | 0.81 |
| size>= 2cm | -1.26 | 0.24 | -5.27 | 0.00 | -1.74 | -0.80 |
所以，此次分析的结论是，在 5% 的显著性水平下，数据无法提供有效证据证明，当调整了结石尺寸之后，外科手术和超声碎石术治疗肾结石有差别。
39.2 交互作用
我们决定调整一个混杂因素的时候，其实同时包含了“在不同混杂因素的程度下，暴露变量和结果变量之间的关系不变” 的假设。但是，一个混杂因素，它可能同时还扮演了改变暴露变量和结果变量之间关系的角色 (effect modifier/交互作用效应)。另外还有的情况下，某变量可能改变了暴露变量和结果变量之间的关系，却不一定是一个混杂因素。此时我们说这个起到改变关系的变量和暴露变量之间发生了交互作用。
如果暴露变量在某个分组变量的不同层 (strata) 之间是不同质的 (hetergeneous, not constant)，我们建议要分析且报告不同层各自的比值比。惟一的例外是 RCT 临床试验的时候，我们更加关心调整后合并比值比。因为一项治疗药物对不同的 “个体” 有不同的疗效是必然的，但是，RCT 的目的是要评价的其实是这个药物或者新疗法在整个人群中的疗效是怎样的。
我们给肾结石数据加上治疗方案和结石尺寸大小的交互作用项，结果如下：
ModelStone2 <- glm(cbind(Success, n - Success) ~ treatment*size,
family = binomial(link = logit), data = Stone)
summary(ModelStone2)
Call:
glm(formula = cbind(Success, n - Success) ~ treatment * size,
family = binomial(link = logit), data = Stone)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.8718 0.1790 10.455 < 2e-16 ***
treatmentSurgery 0.7309 0.4594 1.591 0.11163
size>= 2cm -1.0833 0.3004 -3.606 0.00031 ***
treatmentSurgery:size>= 2cm -0.5245 0.5372 -0.976 0.32882
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 3.3124e+01 on 3 degrees of freedom
Residual deviance: 1.6653e-14 on 0 degrees of freedom
AIC: 27.347
Number of Fisher Scoring iterations: 3| term | estimate | std.error | statistic | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|
| (Intercept) | 1.87 | 0.18 | 10.46 | 0.00 | 1.53 | 2.24 |
| treatmentSurgery | 0.73 | 0.46 | 1.59 | 0.11 | -0.10 | 1.73 |
| size>= 2cm | -1.08 | 0.30 | -3.61 | 0.00 | -1.67 | -0.49 |
| treatmentSurgery:size>= 2cm | -0.52 | 0.54 | -0.98 | 0.33 | -1.65 | 0.48 |
交互作用项的回归系数是否为零的检验结果是\(p = 0.329\)，提示数据无法提供足够的证据证明结石尺寸对治疗方案和手术成功与否之间的关系造成有意义的交互作用(There is no evidence of an interaction between stone size and surgery)。所以此次的数据分析我们可以报告结石尺寸调整后的手术成功比值比就可以了。其中，交互作用项的回归系数\(-0.5245 = \text{log}(0.59)\)，的含义是当结石尺寸\(\geqslant 2 \text{cm}\) 时，外科手术和超声碎石术成功比值比的对数的差。我们可以看到一开始我们计算的分层比值比的比值 \(\frac{1.23}{2.08} = 0.59\)。还注意到这已经是一个饱和模型 (模型偏差为零)，模型的拟合值和我们的观测值是完全相同的。
另外一点不得不提的是，交互作用项的回归系数的点估计\(0.59\) 其实低于零假设时的\(1\) 挺多的，它的\(95\%\) 置信区间也相当的宽$(0.21,1.70) $，所以其实这里我们没有获得足够的证据证明替代假设(有交互作用)，很难说不是因为样本量不足导致的统计效能较低，所以我们没有办法从这个数据完全排除结石尺寸对治疗方案的选择和手术成功的关系之间的交互作用。 (We really cannot be sure that there is no interaction in truth - the data are consistent with quite large interactions between size and surgery effect.) 因此，有些统计学家可能会倾向于报告分层的比值比，因为我们没有办法从这个样本数据排除有较强交互作用存在的可能性。
39.3 可压缩性 collapsibility
模型可压缩性的概念可以这样来理解：
当我们把一个我们认为很重要的混杂因子 加到模型中去时，自然而然我们会期待其对结果变量的效果估计(effect estimate) (回归系数)在调整前后发生变化。如果是反过来，当我们把一个我们认为不重要的非混杂因子 加到模型中去时，我们也会自然而然地期待其对结果变量的效果估计(effect estimate)在调整前后不会发生改变才是。不幸的是，我们这种理想中的想当然，仅仅在某些情境下成立，其中之一是简单线性回归 (Section 23)。
39.3.1 线性回归的可压缩性
证明
令 \(Y\) 标记结果变量，\(X\) 标记暴露变量，\(Z\) 则标记我们想要调整的某个混杂因子：
\[ Y = \alpha + \beta_X X + \beta_Z Z + \varepsilon, \text{ where } \varepsilon \sim N(0, \sigma^2) \]
然后把等式两边同时取以暴露变量 \(X\) 为条件的期待值：
\[ E(Y | X) = \alpha + \beta_X X + \beta_Z E(Z|X) \]
如果\(Z\) 和\(X\) 是相互独立的(即，不是\(X, Y\) 之间关系的混淆因子)，那么\(E(Z|X) = E(Z) = \mu_Z\)，因为$X $ 不能提供任何\(Z\) 的有效信息。所以，等式就变为：
\[ E(Y|X) = \alpha + \beta_Z \mu_Z + \beta_X X \]
所以，当我们用简单线性回归来拟合\(X,Y\) 之间的关系时，如果\(Z,X\) 是相互独立的，模型中增加了\(Z\)，不会影响\(X\) 的回归系数。线性回归的这个性质被定义为模型的可压缩性 (linear regression models are collapsible)。
39.3.2 逻辑链接方程时的不可压缩性
使用对数链接方程 (\(\text{log link}\)) 的回归模型，同样具有和线性回归类似的可压缩性。但是，逻辑链接方程 (\(\text{logit link}\)) 的回归模型则不具有可压缩性。下面举例的分层表格和压缩表格，证明了逻辑链接方程不具有可压缩性：
| Strata 1 | Strata 2 | |||
|---|---|---|---|---|
| Outcome | Exposure \(+\) | Exposure \(-\) | Exposure \(+\) | Exposure \(-\) |
| Success | 90 | 50 | 50 | 10 |
| Failure | 10 | 50 | 50 | 90 |
| Total | 100 | 100 | 100 | 100 |
| Odds Ratios | 9 | 9 | ||
从表格的数据计算可知，要被调整的分组变量的两层数据中，暴露变量和结果变量的关系相同，比值比都是\(9\)，没有交互作用，也没有混杂效应(每一层中暴露与非暴露均占相同的\(50\%\))。但是，你如果把这个观测数据合并：
| Outcome | Exposure \(+\) | Exposure \(-\) |
|---|---|---|
| Success | 140 | 60 |
| Failure | 60 | 140 |
| Total | 200 | 200 |
| Odds ratio | 5.4 |
既然已知我们拿来分层的变量对暴露和结果的关系既没有交互作用，也没有混杂效应，那么压缩数据以后的合并比值比应该和分层比值比一样才说的通，可惜我们获得了截然不同的合并比值比(非调整的)。所以在应用逻辑链接方程建立广义线性回归模型的时候，一定不能忘记其不可压缩性的特征。所以，即使加入一个非混淆因子，暴露变量的逻辑回归的效应估计 (系数) 也会发生改变。调整了 \(Z\) 以后的比值比被叫做条件 (或直接使用分层) 比值比。如同表格中实例所示，条件比值比会比边缘比值比 (非调整) 看起来要大一些。
逻辑回归的不可压缩性给我们的启示是，加入某个变量进入逻辑回归模型前后，其他预测变量的回归系数发生的变化可能仅仅是由于模型的不可压缩性导致的变化，而非由于新加入的变量对原先变量与结果之间的关系起到了交互作用或者混杂效应。所以，拟合回归模型的时候，如果你要考虑对某个因素进行调整，必须做的一件事是，先分析它和其他模型中已有的预测变量之间的关系，从而有助于分析并判断把要调整的变量放进模型前后的回归系数变化是否是真的来自于交互作用或者混杂效应。
39.4 交互作用对尺度的依赖性
GLM 模型中的交互作用检验，对选用的尺度 (比值比 OR，还是危险度比 RR) 依赖性极高。用模型可压缩性的数据例子也可以说明交互作用对尺度的依赖性。上文书说到，两个分层中的比值比都是 9，该分层变量既没有交互作用，也不是混淆因子 (当使用比值比的时候)。如果我们改用危险度比(risk ratio, RR)，在分类变量的第一层(Strata 1) 中，暴露的危险度比是\(\frac{90/100}{50/100} = 1.8\)；分类变量的第二层(Strata 2) 中，暴露的危险度比是\(\frac{50/100}{10/100} = 5\)。所以使用危险度比作为评价指标的时候，被调整的分类变量就突然摇身一变变成了有交互作用的因子。这里，我们用数据，证明了交互作用的存在与否，对尺度的选用依赖性极高。 这就导致我们在描述一个变量是否对我们关心的暴露和结果之间的关系有交互作用时，必须明确指出所使用的是比值比，还是危险度比进行的交互作用评价。
39.5 GLM例7
在本次练习中，我们用一个计算机模拟的有5000名高血压患者的RCT实验数据。该数据只有三个随机产生的二进制变量：
treat表示患者是接收了新疗法 (1 = new)，或者继续维持现有的疗法 (0 = current)；basecontrol表示患者在刚进入实验时 (baseline) 血压原本控制的状态 (0 = bad; 1 = good)；fupcontrol表示患者在实验过程的随访结果中 (followup) 血压原本控制的状态 (0 = bad; 1 = good)；
39.5.1 使用你熟悉的统计学软件拟合一个由 fupcontrol 作为结果变量，treat 作为唯一预测变量的广义线性回归模型。根据报告的结果，写一段适用于医学/流行病学文献杂志的报告。
highbp <- haven::read_dta("../Datas/highbp.dta")
m0 <- glm(fupcontrol ~ treat, data = highbp,
family = binomial(link = logit))
summary(m0); jtools::summ(m0, exp = TRUE, confint = TRUE, digits = 7)
Call:
glm(formula = fupcontrol ~ treat, family = binomial(link = logit),
data = highbp)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.81512 0.04337 -18.80 <2e-16 ***
treat 0.82632 0.05900 14.01 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 6749.1 on 4999 degrees of freedom
Residual deviance: 6548.2 on 4998 degrees of freedom
AIC: 6552.2
Number of Fisher Scoring iterations: 4| Observations | 5000 |
| Dependent variable | fupcontrol |
| Type | Generalized linear model |
| Family | binomial |
| Link | logit |
| χ²(1) | 200.8606211 |
| Pseudo-R² (Cragg-Uhler) | 0.0531595 |
| Pseudo-R² (McFadden) | 0.0297611 |
| AIC | 6552.2390052 |
| BIC | 6565.2733916 |
| exp(Est.) | 2.5% | 97.5% | z val. | p | |
|---|---|---|---|---|---|
| (Intercept) | 0.4425851 | 0.4065197 | 0.4818502 | -18.7953325 | 0.0000000 |
| treat | 2.2849008 | 2.0353891 | 2.5649993 | 14.0057422 | 0.0000000 |
| Standard errors: MLE |
There was strong evidence (p < 0.001) that treatment is associated with the probability of having BP controlled at followup, with patients randomised to the treatment having odds of BP controlled of 2.28 (95%CI 2.04 to 2.56) higher than those randomised to the current treatment.
39.5.2 分析 treat 和 basecontrol 之间的关系，结果是否如你的预期那样？
m1 <- glm(basecontrol ~ treat, data = highbp,
family = binomial(link = logit))
summary(m1); jtools::summ(m1, exp = TRUE, confint = TRUE, digits = 7)
Call:
glm(formula = basecontrol ~ treat, family = binomial(link = logit),
data = highbp)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.01760 0.04000 0.440 0.660
treat 0.02881 0.05658 0.509 0.611
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 6930.2 on 4999 degrees of freedom
Residual deviance: 6929.9 on 4998 degrees of freedom
AIC: 6933.9
Number of Fisher Scoring iterations: 3| Observations | 5000 |
| Dependent variable | basecontrol |
| Type | Generalized linear model |
| Family | binomial |
| Link | logit |
| χ²(1) | 0.2592686 |
| Pseudo-R² (Cragg-Uhler) | 0.0000691 |
| Pseudo-R² (McFadden) | 0.0000374 |
| AIC | 6933.9324824 |
| BIC | 6946.9668687 |
| exp(Est.) | 2.5% | 97.5% | z val. | p | |
|---|---|---|---|---|---|
| (Intercept) | 1.0177563 | 0.9410103 | 1.1007613 | 0.4399943 | 0.6599412 |
| treat | 1.0292268 | 0.9211969 | 1.1499256 | 0.5091777 | 0.6106277 |
| Standard errors: MLE |
这个分析结果提示我们无证据表明该数据中 basecontrol 与治疗方案有关。模型比值比十分接近 1。这与完全符合预期，因为如何选择治疗方案在这个实验中是完全随机分配的，它不应与基线时的血压控制情况有任何关联。随机过程把基线血压控制的情况在两个治疗组之间平衡了。
39.5.3 已知模型中如果增加调整基线变量可能对 fupcontrol 有一定的预测效果。在你的模型中增加基线血压控制情况的变量。与 m0 的结果 (治疗效果 treatment effect；参数标准误 standard error；和 p 值)。重新修改之前用于发表在医学杂志上关于这个分析结果的报告描述。
m2 <- glm(fupcontrol ~ treat + basecontrol, data = highbp,
family = binomial(link = logit))
summary(m2); jtools::summ(m2, exp = TRUE, confint = TRUE, digits = 7)
Call:
glm(formula = fupcontrol ~ treat + basecontrol, family = binomial(link = logit),
data = highbp)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.99454 0.06748 -29.56 <2e-16 ***
treat 1.00376 0.06655 15.08 <2e-16 ***
basecontrol 1.95677 0.06798 28.78 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 6749.1 on 4999 degrees of freedom
Residual deviance: 5588.2 on 4997 degrees of freedom
AIC: 5594.2
Number of Fisher Scoring iterations: 4| Observations | 5000 |
| Dependent variable | fupcontrol |
| Type | Generalized linear model |
| Family | binomial |
| Link | logit |
| χ²(2) | 1160.9142180 |
| Pseudo-R² (Cragg-Uhler) | 0.2797289 |
| Pseudo-R² (McFadden) | 0.1720102 |
| AIC | 5594.1854084 |
| BIC | 5613.7369879 |
| exp(Est.) | 2.5% | 97.5% | z val. | p | |
|---|---|---|---|---|---|
| (Intercept) | 0.1360766 | 0.1192192 | 0.1553176 | -29.5584135 | 0.0000000 |
| treat | 2.7285122 | 2.3948517 | 3.1086597 | 15.0827896 | 0.0000000 |
| basecontrol | 7.0764349 | 6.1936470 | 8.0850477 | 28.7828138 | 0.0000000 |
| Standard errors: MLE |
After adjusting for baseline BP control, the estimated log OR for the new vs current treatment is 1.00 (95%CI 0.87 to 1.13). This is quite a bit larger than the corresponding estimate from part 1, which was 0.83 (0.71, 0.94) . The standard error has, perhaps contrary to expectation, increased from 0.059 to 0.067. The p-values are highly significant from both analyses, but the z-statistics is larger in the baseline adjusted analysis, which indicates this result is more statistically significant.
The increase in the log OR estimate is due to the fact that the baseline adjusted analysis is estimating a different parameter. This is because, although there is no confounding, odds ratios are not collapsible - a conditional odds ratio has a different interpretation from a marginal one. The baseline adjusted analyses estimates the odds ratio for two patients who have the same value of baseline BP control, and this differs from the unconditional OR.
39.5.4 你更推荐使用哪个模型作为最终主要结果的汇报？
都可以。取决于你想拿这结果用于怎样的解释。因为他们回答的是不同的问题。条件比值比 (conditional OR) 的一个缺点是，它的解释需要考虑这个比值比究竟调整了哪些变量（有哪些条件）。这就带来一个缺点，如果相同的治疗方法选取了不同的条件，也就是调整了不同的变量的话，他们就会完全不同。因为条件比值比回答的问题更加关心的是病人个体水平的层面(individual level characteristics) 的问题，也就是说他/她在基线时是否是血压控制的良好的，或者增加其他的变量(性别或者年龄等)。相反的，无条件，或者粗比值比回答的问题是更加广泛的人群水平的治疗效果而不拘泥与究竟调整了哪些变量。后者对于需要采纳针对整体人群政策的决策者来说更加有参考价值。前者对于有（或没有）某些具体特征的病人（或非病人）来说则更加有意义。
39.5.5 实验研究者更想知道新的治疗方案是否由于基线时患者的血压控制情况而有不同。为了回答这个问题，请拟合对应的广义线性回归模型。根据结果回答这个问题。
m3 <- glm(fupcontrol ~ treat + basecontrol + treat*basecontrol,
data = highbp, family = binomial(link = logit))
summary(m3); jtools::summ(m3, exp = TRUE, confint = TRUE, digits = 7)
Call:
glm(formula = fupcontrol ~ treat + basecontrol + treat * basecontrol,
family = binomial(link = logit), data = highbp)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.02870 0.08864 -22.886 <2e-16 ***
treat 1.05611 0.10941 9.653 <2e-16 ***
basecontrol 2.00490 0.10502 19.090 <2e-16 ***
treat:basecontrol -0.08351 0.13795 -0.605 0.545
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 6749.1 on 4999 degrees of freedom
Residual deviance: 5587.8 on 4996 degrees of freedom
AIC: 5595.8
Number of Fisher Scoring iterations: 4| Observations | 5000 |
| Dependent variable | fupcontrol |
| Type | Generalized linear model |
| Family | binomial |
| Link | logit |
| χ²(3) | 1161.2815922 |
| Pseudo-R² (Cragg-Uhler) | 0.2798075 |
| Pseudo-R² (McFadden) | 0.1720647 |
| AIC | 5595.8180342 |
| BIC | 5621.8868070 |
| exp(Est.) | 2.5% | 97.5% | z val. | p | |
|---|---|---|---|---|---|
| (Intercept) | 0.1315068 | 0.1105340 | 0.1564591 | -22.8862857 | 0.0000000 |
| treat | 2.8751646 | 2.3202252 | 3.5628315 | 9.6525024 | 0.0000000 |
| basecontrol | 7.4253853 | 6.0439727 | 9.1225341 | 19.0899856 | 0.0000000 |
| treat:basecontrol | 0.9198831 | 0.7019588 | 1.2054622 | -0.6053666 | 0.5449354 |
| Standard errors: MLE |
Likelihood ratio test
Model 1: fupcontrol ~ treat + basecontrol + treat * basecontrol
Model 2: fupcontrol ~ treat + basecontrol
#Df LogLik Df Chisq Pr(>Chisq)
1 4 -2793.9
2 3 -2794.1 -1 0.3674 0.5444从模型本身的交互作用项结果(p = 0.54)和两个模型的似然比检验结果(p = 0.54) 来看，均无证据表示基线时的血压控制情况和疗效之间存在有意义的交互作用。
39.5.6 换一个模型，先不考虑 basecontrol，使用危险度比 (risk ratio) 来评价不同治疗方案之间的疗效。
m4 <- glm(fupcontrol ~ treat,
data = highbp, family = binomial(link = log))
summary(m4); jtools::summ(m4, exp = TRUE, confint = TRUE, digits = 7)
Call:
glm(formula = fupcontrol ~ treat, family = binomial(link = log),
data = highbp)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.18156 0.03006 -39.31 <2e-16 ***
treat 0.49400 0.03604 13.71 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 6749.1 on 4999 degrees of freedom
Residual deviance: 6548.2 on 4998 degrees of freedom
AIC: 6552.2
Number of Fisher Scoring iterations: 5| Observations | 5000 |
| Dependent variable | fupcontrol |
| Type | Generalized linear model |
| Family | binomial |
| Link | log |
| χ²(1) | 200.8606211 |
| Pseudo-R² (Cragg-Uhler) | 0.0531595 |
| Pseudo-R² (McFadden) | 0.0297611 |
| AIC | 6552.2390052 |
| BIC | 6565.2733916 |
| exp(Est.) | 2.5% | 97.5% | z val. | p | |
|---|---|---|---|---|---|
| (Intercept) | 0.3068000 | 0.2892478 | 0.3254173 | -39.3095666 | 0.0000000 |
| treat | 1.6388526 | 1.5270777 | 1.7588090 | 13.7062367 | 0.0000000 |
| Standard errors: MLE |
链接方程换为对数而非逻辑方程的广义线性回归模型时，新治疗方案和现行治疗方案的血压控制 RR 是 1.64 (95%CI: 1.52, 1.76)
39.5.7 在前一模型m4中加入 basecontrol，与未加入该变量时模型的输出结果相比，有什么不同？
m5 <- glm(fupcontrol ~ treat + basecontrol,
data = highbp, family = binomial(link = log))
summary(m5); jtools::summ(m5, exp = TRUE, confint = TRUE, digits = 7)
Call:
glm(formula = fupcontrol ~ treat + basecontrol, family = binomial(link = log),
data = highbp)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.87003 0.04514 -41.43 <2e-16 ***
treat 0.44211 0.03175 13.92 <2e-16 ***
basecontrol 1.11661 0.04324 25.83 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 6749.1 on 4999 degrees of freedom
Residual deviance: 5613.9 on 4997 degrees of freedom
AIC: 5619.9
Number of Fisher Scoring iterations: 6| Observations | 5000 |
| Dependent variable | fupcontrol |
| Type | Generalized linear model |
| Family | binomial |
| Link | log |
| χ²(2) | 1135.2087599 |
| Pseudo-R² (Cragg-Uhler) | 0.2742121 |
| Pseudo-R² (McFadden) | 0.1682015 |
| AIC | 5619.8908665 |
| BIC | 5639.4424461 |
| exp(Est.) | 2.5% | 97.5% | z val. | p | |
|---|---|---|---|---|---|
| (Intercept) | 0.1541189 | 0.1410694 | 0.1683756 | -41.4273054 | 0.0000000 |
| treat | 1.5559807 | 1.4620984 | 1.6558911 | 13.9236039 | 0.0000000 |
| basecontrol | 3.0544926 | 2.8063141 | 3.3246190 | 25.8258341 | 0.0000000 |
| Standard errors: MLE |
由于对数链接方程和逻辑链接方程不同，是具有可压缩性的 (collapisble)。所以如果增加的预测变量和结果变量是无关的，那么我们应该认为评价治疗效果的危险度比在加入无关变量前后不会有太大的变化。但实际结果我们看见危险度比略微减少了一些。这主要是因为这里的模型本身是有问题的。
39.5.8 给上述模型增加交互作用项。对于危险度比作为指标时的交互作用分析结果，和使用比值比时相比，你有怎样的思考和结论？
m6 <- glm(fupcontrol ~ treat + basecontrol + treat*basecontrol,
data = highbp, family = binomial(link = log))
summary(m6); jtools::summ(m6, exp = TRUE, confint = TRUE, digits = 7)
Call:
glm(formula = fupcontrol ~ treat + basecontrol + treat * basecontrol,
family = binomial(link = log), data = highbp)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.15225 0.07834 -27.473 < 2e-16 ***
treat 0.85895 0.09112 9.426 < 2e-16 ***
basecontrol 1.44713 0.08336 17.359 < 2e-16 ***
treat:basecontrol -0.48113 0.09705 -4.958 7.14e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 6749.1 on 4999 degrees of freedom
Residual deviance: 5587.8 on 4996 degrees of freedom
AIC: 5595.8
Number of Fisher Scoring iterations: 6| Observations | 5000 |
| Dependent variable | fupcontrol |
| Type | Generalized linear model |
| Family | binomial |
| Link | log |
| χ²(3) | 1161.2815922 |
| Pseudo-R² (Cragg-Uhler) | 0.2798075 |
| Pseudo-R² (McFadden) | 0.1720647 |
| AIC | 5595.8180342 |
| BIC | 5621.8868070 |
| exp(Est.) | 2.5% | 97.5% | z val. | p | |
|---|---|---|---|---|---|
| (Intercept) | 0.1162228 | 0.0996798 | 0.1355112 | -27.4728008 | 0.0000000 |
| treat | 2.3606845 | 1.9745766 | 2.8222919 | 9.4262970 | 0.0000000 |
| basecontrol | 4.2509087 | 3.6101303 | 5.0054218 | 17.3593671 | 0.0000000 |
| treat:basecontrol | 0.6180868 | 0.5110253 | 0.7475780 | -4.9576328 | 0.0000007 |
| Standard errors: MLE |
Likelihood ratio test
Model 1: fupcontrol ~ treat + basecontrol + treat * basecontrol
Model 2: fupcontrol ~ treat + basecontrol
#Df LogLik Df Chisq Pr(>Chisq)
1 4 -2793.9
2 3 -2806.9 -1 26.073 3.288e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1我们发现当我们把评价指标换成危险度比之后，治疗方案和基线血压控制情况之间的交互作用变得高度有统计学意义 (p < 0.0001)。这和之前使用逻辑链接函数的模型大相径庭。这个结果提示我们交互作用存在与否是取决于你使用的评价指标的 (scale dependent)。这里我们发现本数据中 log risk ratio 的尺度上交互作用存在且有意义， log odds ratio 的尺度上则并不存在有意义的交互作用。所以，采用危险度比作为评价指标时，我们发现之前我们作出的结论应该需要被修正。
39.5.9 如果说不考虑一个RCT的统计分析不能在收集完数据之后修改这一事实，你认为危险度比模型和比值比模型更应该使用哪一个来总结本数据的结果呢？
Since there is no interaction on the log odds scale, using a logisitic regression is probably preferable, as it gives a simpler model which correctly models the data. Whether the unadjusted or baseline adjusted results are presented is a question which is still hotly debated. The latter has increased power to detect a treatment effect, but as we have seen it estimates a different parameter to the marginal unadjusted OR.
39.5.10 证明危险度比模型是可以压缩的。 Prove that the log-link models are collapsible.
令 \(Y\) 为结果变量，\(X\) 为主要暴露变量 (main exposure)，\(Z\) 为与 \(X\) 相互独立的变量，也就是 \(Z\) 本身并不是一个 \(X\) 的混杂因子。那么为了计算 \(E(Y|X)\)，我们可以使用期望的数学算法来证明：
\[ \begin{aligned} \text{Assume that } E(Y | X, Z) & = \exp(\beta_0 + \beta_X X + \beta_Z Z) \\ \text{then } E(Y|X) & = E(E(Y | X,Z) | X) \text{ by the law of total expectation}\\ & = E(e^{\beta_0}e^{\beta_XX}e^{\beta_ZZ} | X) \\ & = e^{\beta_0}e^{\beta_XX}E(e^{\beta_ZZ} |X) \\ \because X \perp \!\!\! \perp Z & \therefore E(e^{\beta_ZZ} |X) = E(e^{\beta_ZZ}) \text{ is a constant}\\ & = e^{\beta_0}e^{\beta_XX} \times k \\ \Rightarrow \log(E(Y|X)) & = \beta_0 + \beta_XX + \log(k) \end{aligned} \]
由此，我们证明了在危险度比的对数链接方程的广义线性回归模型中，增加非混杂因子变量并不会改变主要暴露变量的回归系数的参数估计。 log risk ratio without adjusting for \(Z\) is \(\beta_X\), the same value as in the adjusted value. Therefore, a GLM with a log-link do possess the collapsibility property.
如果你觉得上面的证明你看不懂，我下面自己写的这个证明可能相对容易一些：
假定 \(Z\) 也是二进制变量。
\[ \begin{aligned} E(Y | X) & = \text{Pr}(Y = 1 | X) \\ & = \text{Pr}(Y = 1 | X, Z = 1) \text{Pr}(Z = 1| X) + \text{Pr}(Y = 1 | X, Z = 0)\text{ Pr}(Z = 0 | X) \\ & = \sum_Z \text{Pr}(Y = 1 | X,Z) \text{Pr}(Z | X) \\ & = \exp(\beta_0 + \beta_XX + \beta_Z) \text{Pr}(Z = 1 | X) + \exp(\beta_0 + \beta_XX ) \text{Pr}(Z = 0 | X) \\ & = \exp(\beta_0 + \beta_XX) \color{red}{\{\exp(\beta_Z)\text{Pr}(Z = 1 | X) +\text{Pr}(Z = 0 | X) \}} \\ & \text{the part in red color is a constant denote as }k \\ & = \exp(\beta_0 + \beta_XX + \log(k)) \\ \Rightarrow & \beta_X \text{ was not affected with or without adjusting for }Z \end{aligned} \]