第 12 章　平均因果效应的双重稳健估计量，或增广逆倾向得分加权估计量

在可忽略性 \(Z\perp\!\!\!\perp\{Y(1),Y(0)\}\mid X\) 与重叠 \(0<e(X)<1\) 下，第 11 章给出了平均因果效应 \(\tau=\mathrm{E}\{Y(1)-Y(0)\}\) 的两套识别公式。第一套是结果回归：

\[ \tau=\mathrm{E}\{\mu_1(X)\}-\mathrm{E}\{\mu_0(X)\}, \tag{12.1} \]

其中

\[ \mu_1(X)=\mathrm{E}\{Y(1)\mid X\}=\mathrm{E}(Y\mid Z=1,X), \]

\[ \mu_0(X)=\mathrm{E}\{Y(0)\mid X\}=\mathrm{E}(Y\mid Z=0,X) \]

分别是处理下与对照下、给定协变量后的结果条件均值。第二套是 IPW：

\[ \tau = \mathrm{E}\left\{\frac{ZY}{e(X)}\right\} - \mathrm{E}\left\{\frac{(1-Z)Y}{1-e(X)}\right\}, \tag{12.2} \]

其中 \(e(X)=\operatorname{pr}(Z=1\mid X)\) 就是第 11 章的倾向得分。

结果回归估计量要为结果对处理与协变量建模；结果模型设定正确，它就相合。IPW 估计量要为处理对协变量建模；倾向得分模型设定正确，它就相合。

数学上，把 (12.1) 与 (12.2) 拼来拼去，能生出许多不同的识别公式。下面这一份特别拼法，理论性质很讨人喜欢：倾向得分模型与结果模型，只要对上其中一个，估计就还是相合的。它叫双重稳健估计量（doubly robust estimator），由 James Robins 力推（Scharfstein et al., 1999; Bang and Robins, 2005）。

12.1　双重稳健估计量

12.1.1　总体版本

为结果的条件均值设定模型 \(\mu_1(X,\beta_1)\) 与 \(\mu_0(X,\beta_0)\)，参数是 \(\beta_1\) 与 \(\beta_0\)。例如，工作模型里条件均值是线性或逻辑的，参数就是回归系数。若结果模型设定正确，则 \(\mu_1(X,\beta_1)=\mu_1(X)\)、\(\mu_0(X,\beta_0)=\mu_0(X)\)。再为倾向得分设定工作模型 \(e(X,\alpha)\)，参数是 \(\alpha\)。例如工作模型是逻辑的，\(\alpha\) 就是回归系数。若倾向得分模型设定正确，则 \(e(X,\alpha)=e(X)\)。实践里，两个模型都可能设错。正因为可能设错，它们有时被叫做工作模型（working model）。

定义

\[ \tilde\mu_1^{\mathrm{dr}} = \mathrm{E}\left[ \frac{Z\bigl\{Y-\mu_1(X,\beta_1)\bigr\}}{e(X,\alpha)} + \mu_1(X,\beta_1) \right], \tag{12.3} \]

\[ \tilde\mu_0^{\mathrm{dr}} = \mathrm{E}\left[ \frac{(1-Z)\bigl\{Y-\mu_0(X,\beta_0)\bigr\}}{1-e(X,\alpha)} + \mu_0(X,\beta_0) \right], \tag{12.4} \]

也可以写成

\[ \tilde\mu_1^{\mathrm{dr}} = \mathrm{E}\left[ \frac{ZY}{e(X,\alpha)} - \frac{Z-e(X,\alpha)}{e(X,\alpha)}\,\mu_1(X,\beta_1) \right], \tag{12.5} \]

\[ \tilde\mu_0^{\mathrm{dr}} = \mathrm{E}\left[ \frac{(1-Z)Y}{1-e(X,\alpha)} - \frac{e(X,\alpha)-Z}{1-e(X,\alpha)}\,\mu_0(X,\beta_0) \right]. \tag{12.6} \]

(12.3) 与 (12.4) 是在结果回归估计上，用残差的逆倾向得分加权去增广；(12.5) 与 (12.6) 是在 IPW 估计上，用填补出来的结果去增广。因此，双重稳健估计量也叫增广逆倾向得分加权（augmented inverse propensity score weighting, AIPW）估计量。

这份增广，把理论性质加强成下面这样。

定理 12.1　假定可忽略性 \(Z\perp\!\!\!\perp\{Y(1),Y(0)\}\mid X\) 与重叠 \(0<e(X)<1\)。

  1. 若 \(e(X,\alpha)=e(X)\) 或 \(\mu_1(X,\beta_1)=\mu_1(X)\)，则 \(\tilde\mu_1^{\mathrm{dr}}=\mathrm{E}\{Y(1)\}\)。
  2. 若 \(e(X,\alpha)=e(X)\) 或 \(\mu_0(X,\beta_0)=\mu_0(X)\)，则 \(\tilde\mu_0^{\mathrm{dr}}=\mathrm{E}\{Y(0)\}\)。
  3. 若 \(e(X,\alpha)=e(X)\)，或 \(\{\mu_1(X,\beta_1)=\mu_1(X),\,\mu_0(X,\beta_0)=\mu_0(X)\}\)，则 \(\tilde\mu_1^{\mathrm{dr}}-\tilde\mu_0^{\mathrm{dr}}=\tau\)。

由定理 12.1，只要倾向得分模型或结果模型对上其中一个，\(\tilde\mu_1^{\mathrm{dr}}-\tilde\mu_0^{\mathrm{dr}}\) 就等于 \(\tau\)。这就是「双重稳健」这个名字的来历。

定理 12.1 的证明　我只证 \(\mu_1=\mathrm{E}\{Y(1)\}\) 这一侧；\(\mu_0=\mathrm{E}\{Y(0)\}\) 的证明类似。

有分解

\[ \begin{aligned} &\tilde\mu_1^{\mathrm{dr}}-\mathrm{E}\{Y(1)\}\\ &= \mathrm{E}\left[ \frac{Z\bigl\{Y(1)-\mu_1(X,\beta_1)\bigr\}}{e(X,\alpha)} - \bigl\{Y(1)-\mu_1(X,\beta_1)\bigr\} \right] \quad\text{（按定义）}\\ &= \mathrm{E}\left[ \frac{Z-e(X,\alpha)}{e(X,\alpha)} \bigl\{Y(1)-\mu_1(X,\beta_1)\bigr\} \right] \quad\text{（合并同类项）}\\ &= \mathrm{E}\left( \mathrm{E}\left[ \frac{Z-e(X,\alpha)}{e(X,\alpha)} \bigl\{Y(1)-\mu_1(X,\beta_1)\bigr\} \Bigm| X \right] \right) \quad\text{（塔性质）}\\ &= \mathrm{E}\left[ \mathrm{E}\left\{\frac{Z-e(X,\alpha)}{e(X,\alpha)}\Bigm| X\right\} \times \mathrm{E}\bigl\{Y(1)-\mu_1(X,\beta_1)\mid X\bigr\} \right] \quad\text{（可忽略性）}\\ &= \mathrm{E}\left[ \frac{e(X)-e(X,\alpha)}{e(X,\alpha)} \times \bigl\{\mu_1(X)-\mu_1(X,\beta_1)\bigr\} \right]. \end{aligned} \]

因此，只要 \(e(X,\alpha)=e(X)\) 或 \(\mu_1(X,\beta_1)=\mu_1(X)\)，就有 \(\tilde\mu_1^{\mathrm{dr}}-\mathrm{E}\{Y(1)\}=0\)。\(\square\)

12.1.2　样本版本

有了 \(\tilde\mu_1^{\mathrm{dr}}\) 与 \(\tilde\mu_0^{\mathrm{dr}}\) 的总体版本，就可以用样本类似物去构造 \(\tau\) 的双重稳健估计。

定义 12.1（平均因果效应的双重稳健估计量）　基于数据 \((X_i,Z_i,Y_i)_{i=1}^{n}\)，可按下面几步得到 \(\tau\) 的双重稳健估计：1

  1. 得到倾向得分的拟合值：\(e(X_i,\hat\alpha)\)；
  2. 得到结果均值的拟合值：\(\mu_1(X_i,\hat\beta_1)\) 与 \(\mu_0(X_i,\hat\beta_0)\)；
  3. 构造双重稳健估计量：\(\hat\tau^{\mathrm{dr}}=\hat\mu_1^{\mathrm{dr}}-\hat\mu_0^{\mathrm{dr}}\)，其中

\[ \hat\mu_1^{\mathrm{dr}} = \frac{1}{n}\sum_{i=1}^{n} \left[ \frac{Z_i\bigl\{Y_i-\mu_1(X_i,\hat\beta_1)\bigr\}}{e(X_i,\hat\alpha)} + \mu_1(X_i,\hat\beta_1) \right], \]

\[ \hat\mu_0^{\mathrm{dr}} = \frac{1}{n}\sum_{i=1}^{n} \left[ \frac{(1-Z_i)\bigl\{Y_i-\mu_0(X_i,\hat\beta_0)\bigr\}}{1-e(X_i,\hat\alpha)} + \mu_0(X_i,\hat\beta_0) \right]. \]

由定义 12.1，也可以把双重稳健估计量写成

\[ \hat\tau^{\mathrm{dr}} = \hat\tau^{\mathrm{reg}} + \frac{1}{n}\sum_{i=1}^{n} \frac{Z_i\bigl\{Y_i-\mu_1(X_i,\hat\beta_1)\bigr\}}{e(X_i,\hat\alpha)} - \frac{1}{n}\sum_{i=1}^{n} \frac{(1-Z_i)\bigl\{Y_i-\mu_0(X_i,\hat\beta_0)\bigr\}}{1-e(X_i,\hat\alpha)}. \]

与 (12.5)、(12.6) 平行，还可以改写成

\[ \hat\tau^{\mathrm{dr}} = \hat\tau^{\mathrm{ipw}} - \frac{1}{n}\sum_{i=1}^{n} \frac{Z_i-e(X_i,\hat\alpha)}{e(X_i,\hat\alpha)}\,\mu_1(X_i,\hat\beta_1) + \frac{1}{n}\sum_{i=1}^{n} \frac{e(X_i,\hat\alpha)-Z_i}{1-e(X_i,\hat\alpha)}\,\mu_0(X_i,\hat\beta_0). \]

Funk et al. (2011) 建议：从 \((Z_i,X_i,Y_i)_{i=1}^{n}\) 做非参数自助，去近似 \(\hat\tau^{\mathrm{dr}}\) 的方差。

12.2　双重稳健估计量的更多直觉与理论

本章开头说：从结果回归与 IPW 的基本识别公式出发，立刻能生出无穷多种别的识别公式。可 (12.3) 与 (12.4) 这种特别的双重稳健形式，并不是一眼就能想到的。它们最初的动机相当理论，依赖高等数理统计里的半参数效率理论（semiparametric efficiency theory; Bickel et al., 1993）。那超出了本书的程度。下面我从两个更直觉的角度，把 (12.3) 与 (12.4) 重新造一遍。12.2.1 与 12.2.2 两节都只谈 \(\mathrm{E}\{Y(1)\}\) 的估计；\(\mathrm{E}\{Y(0)\}\) 由对称性平行可得。

12.2.1　减小 IPW 估计量的方差

基于

\[ \mu_1 = \mathrm{E}\left\{\frac{ZY}{e(X)}\right\} \]

的 IPW 估计，完全不看 \(Y\) 的结果模型。好处是：不必假定任何结果模型，它仍然相合。可若协变量能预测结果，即便工作结果模型是错的，残差的方差通常也比结果本身更小。对一个可能设错的结果模型 \(\mu_1(X,\beta_1)\)，有一个平凡分解：

\[ \mu_1=\mathrm{E}\{Y(1)\}=\mathrm{E}\bigl\{Y(1)-\mu_1(X,\beta_1)\bigr\}+\mathrm{E}\bigl\{\mu_1(X,\beta_1)\bigr\}. \]

把 \(Y(1)-\mu_1(X,\beta_1)\) 看成处理下的「伪潜在结果」，对第一项再用 IPW 公式，上面就变成

\[ \begin{align} \mu_1 &= \mathrm{E}\left\{\frac{Z\bigl\{Y-\mu_1(X,\beta_1)\bigr\}}{e(X)}\right\} + \mathrm{E}\bigl\{\mu_1(X,\beta_1)\bigr\} \tag{12.7}\\ &= \mathrm{E}\left\{\frac{Z\bigl\{Y-\mu_1(X,\beta_1)\bigr\}}{e(X)}+\mu_1(X,\beta_1)\right\}, \tag{12.8} \end{align} \]

只要倾向得分模型正确，它就成立——不必假定结果模型正确。用工作模型去换效率，是抽样调查里的老想法。Little and An (2004) 与 Lumley et al. (2011) 指出过它与双重稳健估计量的联系。

12.2.2　减小结果回归估计量的偏倚

12.2.1 节从 IPW 出发，用工作结果模型去换效率。也可以反过来：从结果回归

\[ \tilde\mu_1=\mathrm{E}\bigl\{\mu_1(X,\beta_1)\bigr\} \]

出发；结果模型可能是错的，它也就未必等于 \(\mu_1\)。这份估计的偏倚是 \(\mathrm{E}\{\mu_1(X,\beta_1)-Y(1)\}\)，若倾向得分模型正确，可用 IPW 去估：

\[ B = \mathrm{E}\left\{\frac{Z\bigl\{\mu_1(X,\beta_1)-Y\bigr\}}{e(X)}\right\}. \]

于是去偏之后的估计是 \(\tilde\mu_1-B\)，与 (12.8) 是同一个东西。

12.3　例子

12.3.1　\(\tau\) 的几种典范估计量小结

下面的 R 代码实现 \(\tau\) 的结果回归、HT、Hajek 与双重稳健估计。它们都可以方便地从 glm 的拟合值算出来。倾向得分模型默认用逻辑模型；结果模型默认用线性模型，out.family = gaussian。2 若结果是二值的，也可以指定 out.family = binomial，去拟合逻辑模型。

OS_est = function(z, y, x, out.family = gaussian,
                  truncps = c(0, 1))
{
  ## fitted propensity score
  pscore = glm(z ~ x, family = binomial)$fitted.values
  pscore = pmax(truncps[1], pmin(truncps[2], pscore))
  ## fitted potential outcomes
  outcome1 = glm(y ~ x, weights = z,
                 family = out.family)$fitted.values
  outcome0 = glm(y ~ x, weights = (1 - z),
                 family = out.family)$fitted.values
  ## outcome regression estimator
  ace.reg = mean(outcome1 - outcome0)
  ## IPW estimators
  y.treat     = mean(z*y/pscore)
  y.control   = mean((1 - z)*y/(1 - pscore))
  one.treat   = mean(z/pscore)
  one.control = mean((1 - z)/(1 - pscore))
  ace.ipw0    = y.treat - y.control
  ace.ipw     = y.treat/one.treat - y.control/one.control
  ## doubly robust estimator
  res1      = y - outcome1
  res0      = y - outcome0
  r.treat   = mean(z*res1/pscore)
  r.control = mean((1 - z)*res0/(1 - pscore))
  ace.dr    = ace.reg + r.treat - r.control
  return(c(ace.reg, ace.ipw0, ace.ipw, ace.dr))
}

上面这些估计量的解析方差公式写起来很烦。自助法提供了方便的近似：从 \(\{Z_i,X_i,Y_i\}_{i=1}^{n}\) 再抽样即可。在 OS_est 之上，下面的函数同时返回点估计与自助法标准误。

OS_ATE = function(z, y, x, n.boot = 2*10^2,
                  out.family = gaussian, truncps = c(0, 1))
{
  point.est = OS_est(z, y, x, out.family, truncps)
  ## nonparametric bootstrap
  n = length(z)
  x = as.matrix(x)
  boot.est = replicate(n.boot, {
    id.boot = sample(1:n, n, replace = TRUE)
    OS_est(z[id.boot], y[id.boot], x[id.boot, ],
           out.family, truncps)
  })
  boot.se = apply(boot.est, 1, sd)
  res = rbind(point.est, boot.se)
  rownames(res) = c("est", "se")
  colnames(res) = c("reg", "HT", "Hajek", "DR")
  return(res)
}

12.3.2　模拟

我用模拟看这些估计量的有限样本性质，四种情形：

  1. 倾向得分模型与结果模型都正确；
  2. 倾向得分模型错、结果模型对；
  3. 倾向得分模型对、结果模型错；
  4. 两个模型都错。

报告的是：多次模拟下估计量的平均偏倚、真实标准误，以及平均估计标准误。

第 1 种情形的数据生成过程是

x       = matrix(rnorm(n*2), n, 2)
x1      = cbind(1, x)
beta.z  = c(0, 1, 1)
pscore  = 1/(1 + exp(-as.vector(x1 %*% beta.z)))
z       = rbinom(n, 1, pscore)
beta.y1 = c(1, 2, 1)
beta.y0 = c(1, 2, 1)
y1      = rnorm(n, x1 %*% beta.y1)
y0      = rnorm(n, x1 %*% beta.y0)
y       = z*y1 + (1 - z)*y0

第 2 种情形里，把倾向得分改成非线性：

x1     = cbind(1, x, exp(x))
beta.z = c(-1, 0, 0, 1, -1)
pscore = 1/(1 + exp(-as.vector(x1 %*% beta.z)))

第 3 种情形里，把结果模型改成非线性：

beta.y1 = c(1, 0, 0,  0.2, -0.1)
beta.y0 = c(1, 0, 0, -0.2,  0.1)
y1      = rnorm(n, x1 %*% beta.y1)
y0      = rnorm(n, x1 %*% beta.y0)

第 4 种情形里，倾向得分与结果模型都改。

样本量取 \(n=500\)，按上面的数据生成过程独立生成 500 套数据。第 1 种情形：

reg HT Hajek DR
ave.bias \(0.00\) \(0.02\) \(0.03\) \(0.01\)
true.se \(0.11\) \(0.28\) \(0.26\) \(0.13\)
est.se \(0.10\) \(0.25\) \(0.23\) \(0.12\)

所有估计量几乎无偏。两个加权估计量的方差更大。第 2 种情形：

reg HT Hajek DR
ave.bias \(0.00\) \(-0.76\) \(-0.75\) \(-0.01\)
true.se \(0.12\) \(0.59\) \(0.47\) \(0.18\)
est.se \(0.13\) \(0.50\) \(0.38\) \(0.18\)

两个加权估计量偏倚很重，因为倾向得分模型设错了。结果回归与双重稳健几乎无偏。

第 3 种情形：

reg HT Hajek DR
ave.bias \(-0.05\) \(0.00\) \(-0.01\) \(0.00\)
true.se \(0.11\) \(0.15\) \(0.14\) \(0.14\)
est.se \(0.11\) \(0.14\) \(0.13\) \(0.14\)

结果回归的偏倚比另外三个大，因为结果模型设错了。加权与双重稳健几乎无偏。

第 4 种情形：

reg HT Hajek DR
ave.bias \(-0.08\) \(0.11\) \(-0.07\) \(0.16\)
true.se \(0.13\) \(0.32\) \(0.20\) \(0.41\)
est.se \(0.13\) \(0.25\) \(0.16\) \(0.26\)

两个模型都错，所有估计量都有偏。HT 与双重稳健的偏倚最大。两边都错时，双重稳健看起来像双重脆弱。

上面所有情形里，只要估计量对真平均因果效应几乎无偏，自助法标准误就与真实标准误接近。

12.3.3　应用

回到例 10.3，得到如下估计与自助法标准误：

reg HT Hajek DR
est \(-0.017\) \(-1.516\) \(-0.156\) \(-0.019\)
se \(0.230\) \(0.492\) \(0.246\) \(0.233\)

两个加权估计量比另外两个大得多。把估计倾向得分截断在 \([0.1,0.9]\)，得到：

reg HT Hajek DR
est \(-0.017\) \(-0.713\) \(-0.054\) \(-0.043\)
se \(0.223\) \(0.422\) \(0.235\) \(0.231\)

Hajek 估计量变得更靠近结果回归与双重稳健，Horvitz–Thompson 仍是离群的那个。

12.4　再多说几句

回想定理 12.1 的证明：双重稳健的关键，是乘积结构

\[ \tilde\mu_1^{\mathrm{dr}}-\mathrm{E}\{Y(1)\} = \mathrm{E}\left[ \frac{e(X)-e(X,\alpha)}{e(X,\alpha)} \times \bigl\{\mu_1(X)-\mu_1(X,\beta_1)\bigr\} \right], \]

它保证：只要 \(e(X)=e(X,\alpha)\) 或 \(\mu_1(X)=\mu_1(X,\beta_1)\)，估计误差就是零。这份精致的结构也让双重稳健在两个模型都设错时，可能变成双重脆弱。两个误差相乘，有可能生出大得多的误差。12.3.2 节的模拟，印证了这一点。

Kang and Schafer (2007) 用模拟批评过双重稳健估计量。他们发现：有限样本里，它的表现甚至可能比简单的结果回归与 IPW 更野。尽管有这样的批评，自 Scharfstein et al. (1999) 那篇开创性工作以来，双重稳健仍是因果推断里的标准策略。近年它在理论统计与计量文献里又活了过来，还换了个更时髦的名字——双重机器学习（double machine learning; Chernozhukov et al., 2018）。基本想法是：把倾向得分与结果的工作模型，换成机器学习工具；相对传统参数模型，它们可以看成更灵活的模型。

12.5　习题

12.1　一次理智检查
考虑协变量离散 \(X\in\{1,\ldots,K\}\)、关心参数为 \(\tau\) 的情形。不加任何模型假定时，估计倾向得分 \(\hat e(X)\) 等于 \(\hat e_{[k]}=\widehat{\operatorname{pr}}(Z=1\mid X=k)\)，也就是第 \(k\) 层里接受处理的比例；估计结果均值则是层内样本均值 \(\hat{\bar Y}_{[k]}(1)=\hat{\mathrm{E}}(Y\mid Z=1,X=k)\) 与 \(\hat{\bar Y}_{[k]}(0)=\hat{\mathrm{E}}(Y\mid Z=0,X=k)\)（\(k=1,\ldots,K\)）。证明：分层估计量、结果回归估计量、HT 估计量、Hajek 估计量与双重稳健估计量，在数值上完全相同。

12.2　\(\tau\) 的双重稳健估计量的另一种形式
受 (12.7) 启发，\(\mu_1=\mathrm{E}\{Y(1)\}\) 还有另一种双重稳健形式：

\[ \tilde\mu_1^{\mathrm{dr2}} = \frac{ \mathrm{E}\left[ \dfrac{Z\bigl\{Y-\mu_1(X,\beta_1)\bigr\}}{e(X,\alpha)} \right] }{ \mathrm{E}\left[ \dfrac{Z}{e(X,\alpha)} \right] } + \mathrm{E}\bigl\{\mu_1(X,\beta_1)\bigr\}. \]

证明：若 \(e(X,\alpha)=e(X)\) 或 \(\mu_1(X,\beta_1)=\mu_1(X)\)，则 \(\tilde\mu_1^{\mathrm{dr2}}=\mu_1\)。给出估 \(\mu_0\) 的对偶公式。再给出基于这些公式的、\(\tau\) 的双重稳健估计量的样本版本。

注： 这种形式的双重稳健估计量出现在 Robins et al. (2007)。

12.3　双重稳健估计量偏倚的一个上界
考虑 \(\mathrm{E}\{Y(1)\}\) 的双重稳健估计量的总体版本 \(\tilde\mu_1^{\mathrm{dr}}\)。证明

\[ \bigl|\tilde\mu_1^{\mathrm{dr}}-\mathrm{E}\{Y(1)\}\bigr| \le \sqrt{ \mathrm{E}\left[ \frac{\bigl\{e(X)-e(X,\alpha)\bigr\}^2}{e(X,\alpha)^2} \right] \times \mathrm{E}\bigl[\bigl\{\mu_1(X)-\mu_1(X,\beta_1)\bigr\}^2\bigr] }. \]

再给出 \(\tilde\mu_0^{\mathrm{dr}}\) 对 \(\mathrm{E}\{Y(0)\}\) 的偏倚的对偶上界。

注： 证明时，附录 A.1.4 可能有用。

12.4　例 10.1 的数据分析
用到目前为止讨论过的方法，分析数据集 cps1re74.csv。

12.5　分析卡罗林斯卡学院的一套数据
Rubin (2008) 用数据集 karolinska.txt 说明观察性研究里因果推断的想法。数据包含 1988 至 1995 年间瑞典中部与北部诊断的 158 名贲门癌患者：79 人在高容量医院诊断（该时期治疗贲门癌超过十人的医院），79 人在其余小容量医院诊断。处理 \(z\) 是患者是否在高容量医院诊断的指示。结果 \(y\) 是诊断后是否存活超过一年。协变量 \(x\) 包含年龄、是否来自农村、是否男性。

karolinska = read.table("karolinska.txt", header = TRUE)
z = karolinska$hvdiag
y = 1 - (karolinska$year.survival == 1)
x = as.matrix(karolinska[, c(3, 4, 5)])

用到目前为止讨论过的方法分析这套数据。

12.6　推荐阅读
Lunceford and Davidian (2004) 对第 11、12 章讨论的许多方法做了综述与比较。

  1. 前面不想强调工作模型里的参数时，我用 \(\hat e(X_i)\) 表示 \(e(X_i,\hat\alpha)\)，用 \(\hat\mu_z(X_i)\) 表示 \(\mu_z(X_i,\hat\beta_z)\)。↩︎

  2. glm 比 lm 更一般。当 out.family = gaussian 时，glm 与 lm 相同。↩︎