详论生物统计学概论

我都给统辅课笔记起了些什么破标题

选课的时候我还侥幸以为生统概论可以放松心情，来随便听听简单了解一下那种。可能是前置没学完的原因，倒也没那么水，ANOVA 还是花了不少时间看的。

wty 老师的英语口语真不错，可惜我听力一般，一走神就寄了。还是要感谢 V1ncent19 的统辅笔记救我，不过他看起来是前置都学完了来上的这个，写的有点抽象，我想写个废话多一点具体点的。

Analysis of Variance

示例——Treatment of Lead-Exposed Children Trial

100 名被试被分成两组分别使用安慰剂和药剂，在实验开始前、实验第 1, 4, 6 周分别检验了血液铅含量水平。第 i 名被试的实验开始前水平为 \(y_{i1}\)，实验结束后水平为 \(y_{i2}\)，变化水平为 \(d_i = y_{i2}-y_{i1}\)，希望利用这些数据观察药剂的有效程度。

假设线性回归模型为 \(d_i = \beta_0 + \beta_1 x_i +e_i\)，注意此处的 \(x_i\) 是一个分类统计量，也就是说被试使用安慰剂时取 \(x_i = 0\)，被试使用试验药剂时 \(x_i =1\)。做一些假设检验：

• \(H_0: \beta_0 =0\) 检验了安慰剂组的平均变化是否为 \(0\)
• \(H_0:\beta_1 = 0\) 检验了安慰剂组和药剂组之间的平均差异是否为 \(0\)，也即药剂是否对治疗有作用

这一模型也可以写作 \(y_{i2} = y_{i1} + \beta_0 +\beta_1 x_i + e_i\)，但实际上 \(y_{i1}\) 对于 \(y_{i2}\) 的影响也有待商榷，未必在 \(y_{i1}\) 变化一个单位时，\(y_{i2}\) 也随之变化一个单位，所以模型可以修正为 \(y_{i2} = \beta_0 + \beta_1 x_i +\beta_2 y_{i1} +e_i\)。这样又多出来一个假设检验：

• \(H_0 :\beta _2 = 0\) 检验了初始状态对末态是否有影响

不同的检验模型会得到一些不同的结果。类似的场景是我们的出发点。

One Factor ANOVA

首先给出一个希望做检验的场景：\(n_T\) 个实验对象被分成 \(r\) 组，每组有 \(n_i\) 个实验对象，有 \(n_T = \Sigma_{i=1} ^r n_i\)。由此我们得到 \(n_T\) 个数据 \(Y_{ij}\)，\(i\) 表示组别，\(1\leq i \leq r\)，\(j\) 表示在某一组内的编号，\(1 \leq j \leq n_i\)。

Cell means model

模型假设是 \(Y_{ij} = \mu_i + \varepsilon_{ij}\)。

其中，\(\mu_i\) 是第 \(i\) 组的理论均值，\(\varepsilon _{ij}\) i.i.d. \(\sim N(0,\sigma^2)\)。注意到在这一模型假设中有 \(r+1\) 个参数，分别是 \(\mu_1,\mu_2,...,\mu_r,\sigma^2\)，我们需要用得到的数据来对这些未知参数进行估计。考虑一些统计量作为参数的估计量：

\[\bar{Y}_{i.} = \frac{1}{n_i} \Sigma_{j=1 }^{n_i} Y_{ij} = \hat \mu _i\]

\[\bar{Y}_{..} = \frac{1}{n_T} \Sigma_{i=1}^r \Sigma_{j=1} ^{n_i} Y_{ij} = \frac{1}{n_T} \Sigma_{i=1}^r n_i \bar{Y}_{i.}\]

\[s_i ^2 = \Sigma(Y_{ij} - \bar Y_{i.})^2 / (n_i -1) \]

\[s^2 = \frac{1}{n_T- r} \Sigma_{i=1}^r (n_i -1) s_i ^2 = \hat \sigma^2\]

在这一模型中，我们关注的假设检验是 \(i\) 组实验之间是否存在差异，假设检验表示为 \(H_0 : \mu_1= \mu_2 = ... = \mu_r = \mu\)，对应的备择假设即为 \(\{ \mu_i \}_{i=1} ^r\) 中存在不同的项。检验最经典的方法即为 ANOVA，analysis of variance。核心是以下的分解：

\[\begin{aligned} SSTO = \Sigma_i \Sigma_j(Y_{ij} - \bar Y_{..})^2 &= \Sigma_i \Sigma_j (Y_{ij} - \bar Y_{i.} +\bar Y_{i.} - \bar Y_{..})^2 \\& = \Sigma_i n_i (\bar Y_{i.} - \bar Y_{..})^2 + \Sigma_i \Sigma_j (Y_{ij} - \bar Y_{i.})^2 \\ &=SSTR + SSE \end{aligned}\]

可以观察到，\(SSTR\) 是组间差距，体现了不同组别之间的差别，\(SSE\) 是组内差距，体现了同一组内各数据的偏差。注意 \(SSTR\) 的自由度是 \(r-1\)，\(SSE\) 的自由度是 \(n_T - r\)，\(SSTO\) 的自由度是 \(n_T - 1\)。由此可以求出 \(MSE= \frac{SSE}{n_T -r}\)，\(MSTR = \frac{SSTR}{r-1}\)。作业里证明了两个统计量的期望，\(\mathbb E(MSE) = \sigma^2\)，\(\mathbb E(MSTR) = \sigma^2 +\frac{\Sigma_i n_i (\mu _i -\mu_.)^2}{r-1}\)，其中 \(\mu_. = \frac{\Sigma_i n_i \mu_i}{n_T}\)。

在 \(H_0\) 成立时，\(\frac{SSE}{\sigma^2} \sim \chi^2 _{n_T - r}\)，\(\frac{SSTR}{\sigma^2} \sim \chi^2_{r-1}\)。因此 \(F=\frac{MSTR}{MSE} \sim F_{r-1,n_T -r}\) 作为最终的检验统计量。

当 \(F^* > F(1-\alpha , r-1 , n_T -r)\) 时拒绝原假设，否则接受；\(Power = P(F^* > F(1-\alpha , r-1,n_T -r)| \delta)\)，其中 \(\delta\) 是一个非中心偏移量，\(\delta = \frac{1}{\sigma} \sqrt{\frac 1 r \Sigma_i n_i (\mu_i - \mu_.)^2}\)。

Factor Effects Model

Factor Effects Model 是 Cell Means Model 的一个重新参数化的结果。模型假设是 \(Y_{ij} = \mu + \tau_i + \varepsilon_{ij}\)。

其中，\(\mu_i\) 是整体的理论均值，\(\varepsilon _{ij} i.i.d. \sim N(0,\sigma^2)\)。它的参数比 cell mean model 多一个，分别是 \(\mu , \tau_1 , ..., \tau_r , \sigma^2\)，但是自由度是相同的，因为 \(\{ \tau _i \}_{i=1}^r\) 存在一个约束 \(\Sigma_{i=1}^r \tau_i = 0\)，如果没有这个约束会导致存在多组解。考虑一些统计量作为参数的估计量：

\[\bar{Y}_{i.} = \frac{1}{n_i} \Sigma_{j=1 }^{n_i} Y_{ij} = \hat \tau _i + \hat \mu\]

\[\bar{Y}_{..} = \frac{1}{n_T} \Sigma_{i=1}^r \Sigma_{j=1} ^{n_i} Y_{ij} = \frac{1}{n_T} \Sigma_{i=1}^r n_i \bar{Y}_{i.} = \hat{\mu}\]

\[s_i ^2 = \Sigma(Y_{ij} - \bar Y_{i.})^2 / (n_i -1) \]

\[s^2 = \frac{1}{n_T- r} \Sigma_{i=1}^r (n_i -1) s_i ^2 = \hat \sigma^2\]

在这一模型中，我们关注的假设检验仍然是 \(i\) 组实验之间是否存在差异，假设检验表示为 \(H_0 : \tau_1= \tau_2 = ... = \tau_r = 0\)，对应的备择假设即为 \(\{ \mu_i \}_{i=1} ^r\) 中存在不同的项。factor effects model 在参数的含义上比 cell mean model 更清晰。

Two Factors ANOVA

首先给出一个希望做检验的场景：\(nab\) 个实验对象被分成 \(a\times b\) 组，每组有 \(n\) 个实验对象。第 \(ij\) 组的实验条件是 \(A\) 因素的等级为 \(i\)，\(B\) 因素的等级为 \(j\)，其中有 \(1 \leq i \leq a, 1\leq j \leq b\)。由此我们得到 \(nab\) 个数据 \(Y_{ijk}\)，\(i\) 表示以 \(A\) 因素分类的组别，\(j\) 表示以 \(B\) 因素分类的组别，\(k\) 表示在某一组内的编号，\(1 \leq k \leq n\)。

每一组都是 \(n\) 个人，这是一个 balanced design。

Factor Effects Model

模型假设是 \(Y_{ijk} = \mu +\alpha_i + \beta_j +(\alpha \beta)_{ij} + \varepsilon_{ijk}\)。

其中，\(\mu\) 是整体的均值，\(\alpha_i\) 代表只和等级为 \(i\) 的因素 \(A\) 有关的变化，\(\beta_j\) 表示只和等级为 \(j\) 的因素 \(B\) 有关的变化，\((\alpha \beta)_{ij}\) 表示和两个因素同时相关的变化，相当于一个交叉项。

这里面有 \((a+1)(b+1)\) 个参数，对其也有一定的约束：\(\Sigma_i \alpha_i = \Sigma_j \beta_j = \Sigma_i (\alpha \beta)_{ij} = \Sigma_j (\alpha \beta)_{ij} =0\)。做一些其他的参数假设，让参数含义更清晰：

\[\mu_{i.} = \alpha_i +\mu_{..}\]

\[\mu_{.j} = \beta _j + \mu_{..}\]

\[\mu_{ij} = \mu_{..} +\alpha_i + \beta_j +(\alpha \beta)_{ij}\]

如果 \((\alpha \beta)_{ij}=0\) 则说明因素 \(A,B\) 之间不存在相互作用的关系，这样的模型称为 additive model。

考虑一些统计量作为参数的估计量：

\[\bar{Y}_{ij.} = \frac{1}{n} \Sigma_{k=1 }^{n} Y_{ijk} = \hat \mu_{ij} = \hat \mu + \hat \alpha _i + \hat \beta_j + \hat{(\alpha \beta)}_{ij}\]

\[\bar{Y}_{i..} = \frac{1}{bn} \Sigma_{j=1}^b \Sigma_{k=1} ^{n} Y_{ijk} = \hat{\mu} +\hat \alpha_i = \hat \mu_{i.}\]

\[\bar{Y}_{.j.} = \frac{1}{an} \Sigma_{i=1}^a \Sigma_{k=1} ^{n} Y_{ijk} = \hat{\mu} +\hat \beta_j = \hat \mu_{.j}\]

\[\bar{Y}_{...} = \frac{1}{abn}\Sigma_{i=1}^a \Sigma_{j=1}^b \Sigma_{k=1} ^{n} Y_{ijk} = \hat{\mu} \]

Two factors ANOVA 的方差分解更复杂一些：

\[\begin{aligned} SSTO &= \Sigma_i \Sigma_j \Sigma_k (Y_{ijk} - \bar Y_{...})^2 \\ &= \Sigma_i \Sigma_j \Sigma_k ((\bar Y_{i..} - \bar Y_{...})+(\bar{Y}_{.j.} - \bar Y_{...}) +(\bar{Y}_{ij.} - \bar Y_{i..} -\bar Y_{.j.} +\bar Y_{...}) +(Y_{ijk} - \bar Y_{ij.}))^2 \\&=bn\Sigma_i (\bar Y_{i..} - \bar Y_{...})^2+an\Sigma_j (\bar{Y}_{.j.} - \bar Y_{...})^2+n\Sigma_i \Sigma_j (\bar{Y}_{ij.} - \bar Y_{i..} -\bar Y_{.j.} +\bar Y_{...})^2+ \Sigma_{i}\Sigma_j \Sigma_k (Y_{ijk} - \bar Y_{ij.})^2 \\&=SSA +SSB+SSAB+SSE \end{aligned}\]

其中，\(SSA\) 的自由度是 \(a-1\)，\(SSB\) 的自由度是 \(b-1\)，\(SSAB\) 的自由度是 \((a-1)(b-1)\)，\(SSE\) 的自由度是 \(ab(n-1)\)。在这一个复杂问题中我们关心不同的问题，可以做出三种不同的假设检验，有对应的检验统计量。

• \(A\) 因素是否会导致差异？\(H_0: \alpha_1 =\alpha_2=...=\alpha_a\)，在 \(H_0\) 下有 \(F^* = \frac{SSA/(a-1)}{SSE/(ab(n-1))} \sim F_{a-1,ab(n-1)}\)
• \(B\) 因素是否会导致差异？\(H_0: \beta_1 =\beta_2=...=\beta_b\)，在 \(H_0\) 下有 \(F^* = \frac{SSB/(b-1)}{SSE/(ab(n-1))} \sim F_{b-1,ab(n-1)}\)
• \(A,B\) 是否联合作用？\(H_0:(\alpha \beta)_{ij}=0, \forall 1 \leq i \leq a,1\leq j \leq b\)，在 \(H_0\) 下有 \(F^*=\frac{SSAB/(a-1)(b-1)}{SSE/(ab(n-1))} \sim F_{(a-1)(b-1),ab(n-1)}\)

Random Effect Test

Why random effects?

事实上，有的时候对于一个 cell means model 或者 factor effects model，我们并不是很关心每个 \(\mu_i\) 的检验，而是关心整体的 \(\mu_.\) 的情况；这种情况下一般 level of factor 也是随机从一个总体里取出的。所以需要改变一些模型假设，进行 random effect 的检验。

Cell means model for random effect

回顾一下 cell means model：

模型假设是 \(Y_{ij} = \mu_i + \varepsilon_{ij}\)。

其中，\(\mu_i\) 是第 \(i\) 组的理论均值，\(\varepsilon _{ij}\) i.i.d. \(\sim N(0,\sigma^2)\)。注意到在这一模型假设中有 \(r+1\) 个参数，分别是 \(\mu_1,\mu_2,...,\mu_r,\sigma^2\)，我们需要用得到的数据来对这些未知参数进行估计。

对于 cell means model for random effect，需要把 \(\mu_i\) 视作一组随机变量而非参数：

模型假设是 \(Y_{ij} = \mu_i + \varepsilon_{ij}\)。其中有 \(i = 1,2,...,r\)，\(j =1,2,...,n\)，暂时设定为一个 balanced test。

其中，\(\mu_i \sim N(\mu_., \sigma_\mu ^2)\) 不要求相互独立，\(\varepsilon_{ij}\) i.i.d. \(\sim N(0,\sigma^2)\)，且 \(\mu_i\) 和 \(\varepsilon_{ij}\) 之间是相互独立的。 因此，这一模型中共有 \(3\) 个参数，分别是 \(\mu_., \sigma_\mu, \sigma\)，而且 \(Y_{ij}\) 不相互独立。

事实上，我们不要求关于特定的 \(\mu_i\) 的结论，而是关心在 population \(\lbrace \mu_i \rbrace_{i=1} ^r\) 上的推断，以推断 \(\mu\) 和 \(\sigma_\mu^2\) 为主。

简单推导可以得到以下结论：

• \(\mathbb E(Y_{ij}) = \mu_.\)
• \(Var(Y_{ij}) = \sigma^2 +\sigma_\mu ^2\)
• \(Cov(Y_{ij} , Y_{ij'}) = \sigma_\mu ^2, j \neq j'\)
• \(Cov(Y_{ij},Y_{i'j'}) = 0, i \neq i', j \neq j'\)
• \(\mathbb E(MSE) = \sigma^2\)
• \(\mathbb E(MSTR) = \sigma^2 + n\sigma_\mu ^2\)

记 \(\sigma_Y^2 = \sigma_\mu^2 + \sigma^2\)，有一个在推断中有意义的参数是 \(\frac{\sigma_\mu ^2}{\sigma_Y ^2} = \frac{\sigma_\mu ^2}{\sigma_\mu ^2 + \sigma^2}\)。这表征了在整体方差中，由 \(Y_{ij}\) 的随机性产生的比例。

对此我们常做的一个检验是 \(H_0 :\sigma_\mu ^2 = 0; H_a: \sigma_\mu ^2 >0\)，原假设的含义在于说明 \(\mu_i\)​ 是相等的，事实上这和 fixed effect 下的 cell means model 所做的检验是相同的。在 \(H_0\) 成立时，\(\frac{SSE}{\sigma^2} \sim \chi^2 _{r(n-1)}\)，\(\frac{SSTR}{\sigma^2} \sim \chi^2_{r-1}\)。因此 \(F=\frac{MSTR}{MSE} \sim F_{r-1,r(n-1)}\) 作为最终的检验统计量。

当 \(F^* > F(1-\alpha , r-1 , r(n-1))\) 时拒绝原假设，否则接受；\(Power = P(F^* > F(1-\alpha , r-1,r(n-1))| \delta)\)，其中 \(\delta\) 是一个非中心偏移量，\(\delta = \frac{1}{\sigma} \sqrt{\frac 1 r \Sigma_i n_i (\mu_i - \mu_.)^2}\)。这里是暂时认为做的是 balanced test。

如果实验并不是一个 balanced test，每一组的数据量为 \(n_1,n_2,...,n_r\) 且有 \(n_T = \Sigma_{i=1} ^r n_i\)，于是有：

\[\mathbb E(MSTR) = \sigma^2 + n' \sigma_\mu ^2, \quad n' = \frac{1}{r-1} [n_T - \frac{\Sigma_{i=1}^r n_i ^2}{n_T}]\]

对这一模型的参数做一些估计，上文中提到 cell means model for random effect 有三个参数 $_., ^2 $ 和 \(\sigma_\mu ^2\)，事实上对于后两个参数我们也关心 \(\frac{\sigma_\mu ^2}{\sigma_\mu ^2 + \sigma^2}\) 的估计，以下还是在 balanced test 的前提下分点来阐述技术细节：

• \(\mu_.\) 的一个无偏估计即为 \(Y_{..} = \frac{\Sigma_{i=1}^r \Sigma_{j=1}^n Y_{ij} }{rn}\)，它的方差是 \(Var(Y_{..}) = \frac{\sigma^2}{rn} + \frac{\sigma_\mu ^2}{r}\)。

  因此 \(s^2(\bar{Y}_{..}) = \frac{MSTR}{rn}\) 是 \(Y_{..}\) 的方差的无偏估计，\(\mu_.\) 的置信区间是 \(\bar{Y}_{..} \pm t(1-\alpha/2;r-1)s(\bar{Y}_{.. } )\)。

• 由于 \(\frac{MSTR/(n\sigma_\mu ^2 + \sigma^2)}{MSE/ \sigma^2} \sim F_{r-1,r(n-1)}\)，因此有：

  \[L = \frac 1n [\frac{MSTR}{MSE} (\frac{1}{F[1-\alpha/2;r-1;r(n-1)]})-1]\]

  \[R = \frac 1n [\frac{MSTR}{MSE} (\frac{1}{F[\alpha/2; r-1,r(n-1)]})-1]\]

  此时 \(\frac{\sigma_\mu ^2}{\sigma_\mu ^2 + \sigma^2}\) 的 \(100(1-\alpha) \%\) 置信区间是 \([\frac{L}{1+L} , \frac{R}{1+R}]\)。

• \(\sigma^2\) 的无偏估计是 \(MSE = \frac{\Sigma_{i=1} ^r \Sigma_{j=1}^n (Y_{ij} - \bar Y_{i.})}{r(n-1)}\)，置信区间的估计可以建立在 \(\frac{r(n-1)MSE}{\sigma^2} \sim \chi^2[r(n-1)]\) 的基础上，于是有 \(\sigma^2\) 的 \(100(1-\alpha) \%\) 置信区间是 \([\frac{r(n-1)MSE}{\chi^2[1-\alpha /2; r(n-1)]}, \frac{r(n-1)MSE}{\chi^2 [\alpha/2;r(n-1)]}]\)

• \(\sigma_\mu^2\) 的无偏估计是 \(s_\mu ^2 = \frac{MSTR - MSE}{n}\)，但有时数据点决定此处点估计值是负的，此时就用 \(0\) 来代替 \(\sigma_\mu ^2\) 的估计值。

Random factor effects model

相应地可以对 cell means model for random effect 做一些改进。回顾一下 factor effects model：

Factor Effects Model 是 Cell Means Model 的一个重新参数化的结果。模型假设是 \(Y_{ij} = \mu + \tau_i + \varepsilon_{ij}\)。

其中，\(\mu_i\) 是整体的理论均值，\(\varepsilon _{ij}\) i.i.d. \(\sim N(0,\sigma^2)\)。它的参数比 cell mean model 多一个，分别是 \(\mu , \tau_1 , ..., \tau_r , \sigma^2\)，但是自由度是相同的，因为 \(\{ \tau _i \}_{i=1}^r\) 存在一个约束 \(\Sigma_{i=1}^r \tau_i = 0\)，如果没有这个约束会导致存在多组解。

Random factor effects model 也是 cell means model for random effects 的一个重新参数化的结果。

模型假设是 \(Y_{ij} = \mu + \tau_i +\varepsilon_{ij}\)，其中 \(i = 1,2,...,a,j=1,2,...,n_i\)。相对于 cell means model for random effect，这里将 \(\tau_i\) 都视作随机变量，满足 \(\tau_i \sim N(0,\sigma_\tau ^2)\) 不要求相互独立，\(\varepsilon_{ij}\) i.i.d. \(\sim N(0,\sigma^2)\)，且 \(\tau_i\) 和 \(\varepsilon_{ij}\) 之间是相互独立的。这一模型中也是三个参数，分别是 \(\mu_., \sigma_\mu, \sigma\)，而且 \(Y_{ij}\) 不相互独立。

由于是随机变量，自然不再要求满足 \(\Sigma_{i=1}^a \tau_i =0\) 的限制。

实际上可以看到，在 random effect 意义下，factor effects model 和 cell means model 的关系更紧密了，只是相差一个常数。因此所有的检验和推断可以类似的进行。

Mixed factors effects model

回顾一下 two factors effects model：

模型假设是 \(Y_{ijk} = \mu +\alpha_i + \beta_j +(\alpha \beta)_{ij} + \varepsilon_{ijk}\)。

其中，\(\mu\) 是整体的均值，\(\alpha_i\) 代表只和等级为 \(i\) 的因素 \(A\) 有关的变化，\(\beta_j\) 表示只和等级为 \(j\) 的因素 \(B\) 有关的变化，\((\alpha \beta)_{ij}\) 表示和两个因素同时相关的变化，相当于一个交叉项。

这里面有 \((a+1)(b+1)\) 个参数，对其也有一定的约束：\(\Sigma_i \alpha_i = \Sigma_j \beta_j = \Sigma_i (\alpha \beta)_{ij} = \Sigma_j (\alpha \beta)_{ij} =0\)。

类似地，mixed factors effects model 也是把模型假设进行了一些修改，将 \(\alpha_i\) 从参数改成了随机变量。注意 \(\alpha_i , \beta_j\) 中只要有一个改为 random effect，其联合作用 \((\alpha \beta)_{ij}\) 也变为随机变量而非参数。

模型假设是 \(Y_{ijk} = \mu +\alpha_i + \beta_j +(\alpha \beta)_{ij} + \varepsilon_{ijk}\)。

其中，\(\mu\) 是整体的均值，\(\alpha_i\) 代表只和等级为 \(i\) 的因素 \(A\) 有关的变化，\(\beta_j\) 表示只和等级为 \(j\) 的因素 \(B\) 有关的变化，\((\alpha \beta)_{ij}\) 表示和两个因素同时相关的变化，相当于一个交叉项。

需要满足一定的约束：

• \(\Sigma_{i=1} ^ a \alpha _i =0\)
• \(\beta_j \sim N(0,\sigma_\beta ^2)\)
• \((\alpha \beta)_{ij} \sim N(0,\sigma_{\alpha \beta}^2)\)

检验方式和统计量见 ANOVA table construction 环节。

ANOVA table construction

e.g: three factors model

\[Y_{ijkl} = \mu + \alpha_i + \beta_j + \gamma _k +(\alpha \beta)_{ij} + (\alpha \gamma)_{ik} + (\beta \gamma)_{jk} + (\alpha \beta \gamma)_{ijk} + \varepsilon_{ijkl}\]

\[i = 1,2,...,a\]

\[j = 1,2,...,b\]

\[k = 1,2,...,c\]

\[l = 1,2,...,n\]

其中 \(a\) 是 fixed effect，而 \(b,c\) 是 random effect。于是模型参数是 \(\theta = \lbrace \mu , \alpha_i^{i=1,2,...,a} , \sigma_\beta^2 , \sigma_\gamma ^2 , \sigma_{\alpha \beta}^2 , \sigma_{\beta \gamma }^2 , \sigma_{\alpha,\gamma}^2, \sigma_{\alpha,\beta,\gamma}^2,\sigma^2 \rbrace\)

ANOVA table construction 的过程详见 Statistics Note by V1ncent19，P212-215，这里仅给出结果：

Random/Fix	F	R	R	R
# level	a	b	c	n
Index	i	j	k	l	\(\mathbb E(MS)\)
\(\alpha_i\)	0	b	c	n	\(\sigma^2 + cn \sigma_{\alpha \beta}^2 + bn \sigma_{\alpha \gamma}^2 + n \sigma_{\alpha \beta \gamma}^2 + bcn \frac{\Sigma \alpha_i ^2}{a-1}\)
\(\beta_j\)	a	1	c	n	\(\sigma^2 + an \sigma_{\beta \gamma}^2 + acn \sigma_{\beta}^2\)
\(\gamma_k\)	a	b	1	n	\(\sigma^2 + an \sigma_{\beta \gamma}^2 + abn \sigma_{\gamma}^2\)
\((\alpha \beta)_{ij}\)	0	1	c	n	\(\sigma^2 + cn \sigma_{\alpha \beta}^2 + n \sigma_{\alpha \beta \gamma}^2\)
\((\alpha \gamma)_{ik}\)	0	b	1	n	\(\sigma ^2 + bn \sigma_{\alpha \gamma}^2 + n \sigma_{\alpha \beta \gamma}^2\)
\((\beta \gamma)_{jk}\)	a	1	1	n	\(\sigma^2 + an\sigma_{\beta \gamma}^2\)
\((\alpha \beta \gamma)_{ijk}\)	0	1	1	n	\(\sigma^2 + n \sigma_{\alpha \beta \gamma}^2\)
\(\varepsilon_{ijkl}\)	1	1	1	1	\(\sigma^2\)

于是可以利用上述期望进行检验：

• \(H_0 : \alpha_1 = \alpha_2 = ...=\alpha_a\)，检验 \(F_{\alpha_i} ^* = \frac{MS_\alpha + MS_{\alpha ]beta \gamma}}{MS_{\alpha \beta} + MS_{\alpha \gamma}} \sim ^{H_0} F_{(a-1)+(a-1)(b-1)(c-1), (a-1)(b-1)+(a-1)(c-1)}\)
• \(H_0 : \sigma^2 =0\)
• \(H_0 : \sigma_{\beta}^2 =0\)
• ......

R code for ANOVA

详见 Mixed ANOVA in R，和课上讲的是同一个东西但是更精细一点（

期末考前又去看我第二次作业的那个 project，绷不住了，我怎么做了个不合理的 mixed effect model 还拿了满分（

实验设计方法

不如说是实验设计小知识（

Designed experiment 的主要目标是研究若干 explanatory factor 和 response variable 之间的因果关系。在实验中变量控制情况（也就是所谓的 treatments）和被试者的匹配之间是完全随机的，例如在试药时安慰剂和特效药是随机分发的，甚至会将实验条件对被试者保密。

Observational & Experimental Study

• Observational study: The investigator use the data observed in the population to make inference on the relationship between the variables

• Experimental study: The investigator intervene in the natural history by actively alternating one of the variables and then making inference on the relationship between the variables based on the outcomes

简单来说，直接研究通过观察收集来的数据就是 observational study，数据只来自于自然的观察，收集者对事物发生不作任何干涉；experimental study 则需要实验者自行控制变量，在所需条件下收集数据，再进行研究。

Crossed and Nested Factors

这是多因素分析中的一个概念。

• Crossed factors: all combinations of the levels of the two factors are included in the study. e.g.

  Factor Combination	Temperature	Solvent Concentration	Outcome
  1	Low	Low	/
  2	Low	High	/
  3	High	Low	/
  4	High	High	/

• Nested factors: the level of one factor are unique to a particular level of another factor. e.g.

  Factor Combination	Temperature	Solvent Concentration	Outcome
  1	Low	Low	/
  2	High	High	/

Completely randomized design

• 对于一个有三个因素的 designed experiment，建立的模型包括三个部分：

  \(Y=\) overall constant + treatment effect + experimental error

  如果采用线性回归模型，以三个因素为例，全体被试者被分为 \(8\) 组，使得

  \[Y_{ij} = \beta_0 + \beta_1 X_{ij1} + \beta_2 X_{ij2} + \beta_3 X_{ij3} + \varepsilon_{ij}\]

  可以轻松地将其和上述的三个模型组成部分对应起来。其中 \(X_1,X_2,X_3\) 类似于示性函数，此处是一种 categorical variable，$X_i = 1 $ if treatment \(i\)，otherwise \(0\)。

  也可以用因子效应模型来考虑这个问题，表示起来会比较容易；另外要注意的是，regression model 关心平均响应的预测，但 factor effect model 更关心平均响应的差别。

• 仍然以上一部分中的 Treatment of Lead-Exposed Children Trial 为例。

  \(100\) 名被试被分成两组分别使用安慰剂和药剂，在实验开始前、实验第 \(1, 4, 6\) 周分别检验了血液铅含量水平。第 i 名被试的实验开始前水平为 \(y_{i1}\)，实验结束后水平为 \(y_{i2}\)，变化水平为 \(d_i = y_{i2}-y_{i1}\)，希望利用这些数据观察药剂的有效程度。

  我们在这一实验中设定三个不同的 explanatory factor，分别是 treatment type (succimer or placebo)，children's ID (from \(1\) to \(100\)，用于展示个体间的差别)，time point (week \(0,1,4,6\))。其中，treatment type 和 time point 可能会协同作用。Outcome，也即 response variable 设定为血铅浓度 \(Y\)。

  模型为：

  \[Y_{ijt} = \mu + \alpha_i + \beta_{j(i)} + \gamma_t + (\alpha \gamma)_{it} + \varepsilon_{ijt}\]

  参数假设为：

  \[\beta_{j(i)} \sim N(0,\sigma_\beta ^2)\]

  \[\varepsilon_{ijk} \sim N(0,\sigma^2)\]

  自由度和方差分解为：

  \[Y_{ijt} - Y_{...} = (\bar Y_{i..} - \bar Y_{...}) +(\bar Y_{..t} - \bar Y_{...}) +(\bar Y_{i.t} - \bar Y_{i..} - \bar Y_{..t} + \bar Y_{...}) + (\bar Y_{ij.} - \bar Y_{i..}) + (Y_{ijt} - \bar Y_{ij.} - \bar Y_{i.t} + \bar Y_{i..})\]

  \[df: Nn-1 = (a-1)+(n-1)+(a-1)(n-1)+(N-a)+(N-a)(n-1)\]

  可以提出并进行检验的问题有：

  • 安慰剂和药剂之间有无明显的区别？在第 \(0,1,4,6\) 周时，两组之间的血铅值是否有明显差异？

    \(H_0: \mu_t(S) = \mu_t(P), \forall t= 1,2,3,4\)，其中 \(\mu_t = \mu + \gamma_t + (\alpha \gamma)_{it}\)，也即取出全体和时间有关的项。

  • 也可以换一种表示方法：安慰剂和药剂是否在降低血铅方面有同等的效果？

    \(H_0 : \mu_t(S) - \mu_1 (S) = \mu_t(P) - \mu_1(P),\forall t =2,3,4\)

    可以推出这等价于 \(H_0: (\alpha \gamma)_{it} = 0,\forall i = 1,2, t=1,2,3,4\)，这是因为参数假设中认为 \(\Sigma (\alpha \gamma)_{it}=0\)。

Clinical Trials

临床试验最准确的定义是 the evaluation of intervention (treatment) on disease in a controlled experimental setting（绷不住了，感觉看原文比翻译过来要靠谱的多），有不同的类型。

Randomized clinical trial

RCT 的目的一般是在不存在因素的互相干扰情况下做对比，目的往往有找出最好的治疗方案，测试药剂相对于安慰剂的效果，等等。因此，实验组必须在各个方面和初始量化指标上都是一致的，在此基础上接受不同的治疗方案，才能控制变量找出方案对于效果的影响。

RCT 是最有力的临床试验方案，因为是对几乎相同的组别随机分配治疗方案，得到的结果是普适性的。

Non-randomized clinical trial

不随机的点有很多，比如被试者可以自行选择想参加实验的组，自由选择使用的是药剂还是安慰剂，或者这一分组是被实验组织者人为设定而不是随机分配的。

Contingency Table in Clinical Trials

Clinical trial 里面也会用到一些 contingency table 和相应的检验。在这里先给出我们关注的问题和一些相应的定义。

此处我们主要关心的情境是，两组被试者分别使用 medicine 和 placebo，可能会出现 adverse reaction 也可能不会出现，统计全体被试者的反应得到以下表格：

	Adverse Reaction (\(D\))	No Adverse Reaction (\(D^C\))	Total
Medicine (\(E\))	\(n_{11}\)	\(n_{12}\)	\(n_{1+}\)
Placebo (\(E^C\))	\(n_{21}\)	\(n_{22}\)	\(n_{2+}\)
Total	\(n_{+1}\)	\(n_{+2}\)	\(n_{++}\)

相应地定义 \(p_1\) 是 medicine 组发生不良反应的概率，\(p_2\) 是 placebo 组发生不良反应的概率。于是有：

\[p_1 = \mathbb P(D | E) , p_2 = \mathbb P(D |E^C)\]

\[\text{Risk Difference} = \Delta = p_1 - p_2\]

\[\text{Relative Risk} = \phi = \frac{p_1}{p_2}\]

\[\text{Odds Ratio} = \theta = \frac{p_1 / (1-p_1)}{p_2 / (1-p_2)}\]

在这样的情境下我们希望做的检验即为 \(H_0 : p_1 = p_2\)，也即观察 medicine 相对 placebo 来说是否发生不良反应的可能性是相同的。事实上由于上述条件概率的定义，也就是希望检验 adverse reaction 和 medicine/placebo 的选择是独立的。这就归结为了 contingency table 的 independence inference，详见下一部分的 Pearson Test 和 Likelihood Ratio Test。

另外，我们希望使用观察到的数据 \(n_{ij}\) 对上述参数做一些估计。但具体的估计与 clinical trial 的类型有关，此处分别考虑两种类型的 clinical trial，分别是 prospective clinical trial 和 retrospective clinical trial。

虽然 contingency table 从理论上来说对于随机变量 \(Y,Z\) 是完全对称的，但是在实验设计中我们需要先确定 \(Y\) 或者 \(Z\)，进行优先级的选择后再收集数据。这可能会导致条件概率不同。

Prospective clinical trial

招募了许多患有某一疾病的被试者进行试药，\(Y = E / E^C\) 是 medicine 或者 placebo 组在实验开始之前就对全体被试者划分好了，于是 \(Z = D / D^C\) 是是否出现不良反应的量度，是我们在实验中需要关注的。

在这一实验中由于是预先招募了被试者，所以 \(n_{1+},n_{2+},n_{++}\) 是确定的数值。

Prospective clinical trial 显然更可控一些，但是也会遇到一些生存分析的问题。另外，如果是关于罕见病的实验，这一方式可能会比较低效。

此时有参数估计：

\[\hat p_1 = \frac{n_{11}}{n_{1+}}, \hat p_2 = \frac{n_{21}}{n_{2+}}\]

\[\hat \phi = \frac{n_{11}/ n_{1+}}{n_{21}/ n_{2+}}\]

\[\hat \theta = \frac{n_{11}n_{22}}{n_{21}n_{12}}\]

\[\hat{Var} ( \log \hat \theta) = \frac{1}{n_{11}} + \frac{1}{n_{12}} + \frac{1}{n_{21}} + \frac{1}{n_{22}}\]

注意一些参数含义：\(\phi\) 和 \(\theta\) 都表示了组间差距，如果这两个参数比较大就说明 \(Y\) 对 \(Z\) 和 \(Z^C\) 的影响不同。

Retrospective clinical trial

招募了许多产生了不良反应和未产生不良反应的患者，\(Z = D/D^C\) 是是否有不良反应的组在实验开始之前就对被试者划分好了。再对他们发放问卷或者进行采访，了解出现不良反应之前的用药历史。

在这一实验中对被试者进行的划分是关于 \(Z\) 的，因此 \(n_{+1},n_{+2},n_{++}\) 是确定的数值。

Retrospective clinical trial 显然是没有那么精准的，我们难以精准控制 \(Y\)，但是数据的收集更快更方便。

此时有参数估计：

\[\hat p_1 = \frac{\rho \frac{n_{11}}{n_{+1}}}{\rho \frac{n_{11}}{n_{+1}}+(1-\rho) \frac{n_{12}}{n_{+2}}},\hat p_2 = \frac{\rho \frac{n_{21}}{n_{+1}}}{\rho \frac{n_{21}}{n_{+1}}+(1-\rho) \frac{n_{22}}{n_{+2}}} \]

其中 \(\rho\) 是人群中的发病率。

\[\hat \phi = \frac{\hat p_1}{\hat p_2}\]

\[\hat \theta = \frac{n_{11}n_{22}}{n_{21}n_{12}}\]

\[\hat{Var} ( \log \hat \theta) = \frac{1}{n_{11}} + \frac{1}{n_{12}} + \frac{1}{n_{21}} + \frac{1}{n_{22}}\]

可以发现 odds ratio 在两种 clinical trial 里是一致的，因此这一参数更佳。

Gehan's Two-Stage Design

如果某一实验看起来完全达不到我们的最低期望，那我们希望能够尽早在临床试验的早期抛弃这一选择。比如说，对于 \(n\) 位患者，我们能够接受的最低治疗成功率是 \(20\%\)，那么在治疗有效的患者数小于 \(m\) 时就可以认为这一治疗方案是无效的，然后尽早放弃继续的探索。Gehan's two-stage design 展示了选取 \(m\) 的方式。

实际上这个过程类似于寻找算法的终止条件。在最初的状态下选择 \(n_0\) 名患者参与治疗，如果这些人里面没有治疗有效的，那么直接停止这一实验过程；否则继续做实验，并且每次加入新的被试者，希望能够达到期望的最低治疗成功率。

比如说，我们认为当有 \(\pi _0 = 20\%\) 的患者症状缓解时，可以认为这一治疗方案是有效的，这也就是可以接受的最低治疗成功率。在 Gehan's two-stage design 里设定 level 是 \(\alpha_0 = 5 \%\)，在此基础下考虑第一次实验的被试数量 \(n_0\) 的取值。

事实上治疗有效果的患者数是一个服从伯努利分布的随机变量 \(X \sim B(\pi)\)，那么 \(P(X_0 = 0) = (1-\pi)^{n_0}\leq \alpha_0\) 对于任意的 \(\pi \geq \pi_0\) 是成立的，那么 \(n_0 \geq \frac{\log (0.05)}{\log (1-0.20)}=14\) 即为第一次实验的最少被试数量。

也可以通过置信区间的方法来确定。我们希望能够有 \(95 \%\) 的信度认为这一实验可以继续下去，允许 \(\pi_0\) 有 \(\pm 15 \%\) 的浮动，则最少的被试参与数满足 \(1.96(\frac{(1-\pi_0)\pi_0}{n})^{\frac 1 2} = 0.15\)，解得 \(n=28\)。

Inference on Contingency Table

简单来说，孟德尔种豌豆那个表格就是一个 contingency table，到底为什么呈现出 \(9:3:3:1\) 的规律也可以用 goodness of fit test 进行检验。Contingency table analysis 本质上是对 categorical data 的检验。

Motivation Example

我们希望知道吸烟和肺癌之间有无确切的联系，采集了一些数据：吸烟且肺癌的有 \(688\) 人，吸烟但没有得肺癌的有 \(650\) 人；不吸烟且肺癌的有 \(21\) 人，不吸烟且没有得肺癌的有 \(59\) 人。样本量共计 \(1418\) 人，列表如下：

	Lung Cancer	Healthy	Total
Smoke	688	650	1338
Not Smoke	21	59	80
Total	709	709	1418

“看起来”吸烟和肺癌是有关系的，因为吸烟者得肺癌的频率看似更高。但是如何从频率得到概率，仍然需要一些检验的逻辑。

Contingency Table

抽象出来的 contingency table 的逻辑是这样的：

	\(Z=1\)	\(Z=2\)	Total
\(Y=1\)	\(n_{11}\)	\(n_{12}\)	\(n_{1+}\)
\(Y=2\)	\(n_{21}\)	\(n_{22}\)	\(n_{2+}\)
Total	\(n_{+1}\)	\(n_{+2}\)	\(n_{++}\)

当然，\(Y\) 和 \(Z\) 都可以有更多的指标，这不影响我们的构造。比如说它可以长成下面这种通用的形式：

	1	...	J	total
1	\(n_{11}\)	...	\(n_{1J}\)	\(n_{1+}\)
...	...	...	...	...
I	\(n_{I1}\)	...	\(n_{IJ}\)	\(n_{I+}\)
total	\(n_{+1}\)	...	\(n_{+J}\)	\(n_{++}\)

与此同时，\(Y\) 和 \(Z\) 的位置是完全对称的，可以交换二者在行和列里的顺序。

在 \(n_{++},n_{1+},n_{2+},n_{+1},n_{+2}\) 都已知的情况下，完全可以把 \(n_{11}\) 视作一个随机变量，当然表格中的 \(n_{ij}\) 都可以被视为随机变量，但四者中有其一被确定后整个表格就完全被确定，因此只需要讨论 \(n_{11}\) 的数值。

事实上，\(n_{11}\) 服从的是超几何分布。\(P(n_{11} = k) = \frac{C_{n_{1+}}^k C_{n_{2+}}^{n_{+1}-k}}{C_{n_{++}}^{n_{+1}}}\) 对于 \(\max(0,n_{1+} + n_{+1}-n) \leq k \leq \min(n_{1+},n_{+1})\) 都成立。

P-value Calculation Example: 女士品茶

对于一个如下所示的表格，希望判断品茶的女士是否真的有品尝先加奶还是先加茶的能力——尽管看起来是有的，\(4\) 次尝试里她只失误了一次，但是我们希望通过一些更严谨的推断方法来证明。

Actually poured first		Lady says poured first
	tea	milk
tea	3	1
milk	1	3

由于 \(P(n_{11} = 3) = \frac{C_4 ^3 C_4 ^1}{C_8 ^4} = \frac{16}{70}\)，\(P(n_{11}=4) = \frac{C_4 ^4 C_4 ^0}{C_8 ^4} = \frac{1}{70}\)，回顾 p-value 的定义是“比所得数据更加极端的情形所出现的概率”，事实上也就是 \(P(n_{11} \geq 3) = \frac{16}{70} + \frac{1}{70} = \frac{17}{70} = 0.243\)，无法给出充足的证据证明这位女士能够品出添加的顺序。

Inference on Independence

Contingency table 的通用形式如下所示：

	1	...	J	total
1	\(n_{11}\)	...	\(n_{1J}\)	\(n_{1+}\)
...	...	...	...	...
I	\(n_{I1}\)	...	\(n_{IJ}\)	\(n_{I+}\)
total	\(n_{+1}\)	...	\(n_{+J}\)	\(n_{++}\)

另一个检验的重点在于 \(Y\) 和 \(Z\) 联合作用于 \(\{ n_{ij} \}\) 的取值时，我们希望了解二者的作用是否独立。

如果二者的作用独立，就有 \(\pi_{ij} = P(Y=i,Z=j) = P(Y=i)P(Z=j) = \pi_{i+} \pi_{+j}\)。因此我们只要检查实际的 contingency table 和理想中符合二者独立时的 contingency table 有何区别，此处用到 Goodness of Fit Test。

Hypothesis

\(H_0:\) The independence model is true, i.e. \(\pi_{ij}= \pi_{i+} \pi_{+j}\) for all pairs of \((i,j)\).

\(H_1 :\) The saturated model is true, i.e. \(\pi_{ij}\neq \pi_{i+} \pi_{+j}\) for at least one pair of \((i,j)\).

Test Statistic

• The standard Pearson Chi-Square statistic:

  \[X^2 = \Sigma_i \Sigma_j \frac{(O_{ij} - E_{ij})^2}{E_{ij}} \sim ^{approx} \chi^2 _{(I-1)(J-1)}\]

• Likelihood Ratio Test statistic:

  $G^2 = -2 () = 2 i j O{ij} {E{ij} } ^{approx} ^2 _{(I-1)(J-1)} $

注意 \(X^2\) 和 \(G^2\) 的分布都是近似分布，二者并不是等价的检验统计量，计算值时常不同，但检验结果往往是类似的。

Survival Analysis

Basic Concepts

首先要明确一些概念，给我的感觉是像黑话一样（

Concept

Event 是指生存分析所研究的事件，大多数情况下是指某种治疗过程中患者死亡这样的事件发生，但其实生存分析的应用很广泛，甚至不只是在生医领域都有这些概念的使用。患者治疗后恢复/犯人被释放后再次犯罪/确诊新冠都可以作为研究所感兴趣的事件，这些事件发生就称为 event。

At risk 是指尚未发生 event 的个体。事实上大多数情况下研究的 event 都是患者去世之类的负面事件，因此尚未发生 event 的个体会被先入为主地定义成 at risk，即为暴露在即将发生 event 的风险下。什么达摩克斯之剑

Function

一些常用的函数：

• survival function：\(S(t) = P(T > t) = 1-P(T \leq t) = 1-F(t)\)
• mean survival time：\(E(T) = \int_{0} ^{\infty} 1-F(t) dt = \int_{0} ^{\infty} S(t ) dt\)
• hazard function：\(\lambda (t) = \lim_{t \downarrow 0} \frac{P(t \leq T < t+\Delta t | T \geq t)}{\Delta t} = \frac{f(t)}{S(t)} = \frac{-S'(t)}{S(t)} = -\frac{d}{dt} \log S(t)\)
• cumulative hazard function：\(\Lambda (t) = \int_0 ^{t} \lambda(s) ds = - \log S(t)\)

事实上都是对 \(T\) 的分布函数 \(F(t)\) 的变形。

Data

生存分析中有两种不同的数据，censored data 和 uncensored data。

Censored data 是指在研究过程中未能记录到 event 发生（未能记录到 failure）的数据，往往是因为直到研究结束 event 都尚未发生，或者研究对象中途退出研究，也可能是没有及时发现跟进从而未能记录准确的时间，等等。Censored data 又分为三种情况，right-censored，left-censored 和 interval-censored 的情况。

• right-censored data 往往是生存期比所记录的时间要长的情况，censor time 在 event time 之前。例如临床试验持续了一段时间，将 event 设定为观察到患者死亡。某一位患者在第一周加入试验，试验持续十周结束，他直到试验结束都保持健康，没有发生 event，则他的生存期是 10+ weeks，这是一个 right-censored data。
• left-censored data 的一个例子是认为 event 是感染新冠，一位患者在试验第七天发现抗原阳性，但在这之前没有做过检测，不清楚是从什么时候开始感染新冠的，只知道 censor time 在 event time 之后。因此这是一个 left-censored data。
• interval-censored data 的一个例子仍然是认为 event 是感染新冠，一位患者在试验第三天抗原阴性，第七天发现抗原阳性，但在这两次检测之间没有做过检测，也还是不清楚从什么时候开始感染新冠，只知道 event time 在两个 censor time 之间。称为一个 interval-censored data。

此处简化情况，只研究 right-censored data，尝试对它进行标记。记 survival time 为 \(T\)，censor time 为 \(C\)，记 right-censored data 为 \((X,\delta)\)，其中 \(X = \min(T,C)\)。\(\delta\) 是一个示性变量，如果 \(X = C\) 则 \(\delta=0\)，也就是 event 未发生；否则 \(\delta=1\)，\(X = T\)，也就是观察到了 event 发生。

每一组数据贡献的 likelihood 是 \(\lambda (x)^\delta S(x)\)，整体的 likelihood function 即为 $ L = _{i=1} ^n (x_i )^{_i} S(x_i)$。

Kaplan & Meier Estimation

Theory

一个估计 survival function 的方法。原理上来说是把 \(S(t) = P(T >t)\) 分解成了一列条件概率的积：

\[S(t) = P(T>t) = P(T>t_1) P(T>t | T>t_1) = P(T>t_1) \{\Pi_{j=2} ^k P(T>t_j | T>t_{j-1}) \} P(T>t|T>t_k)\]

与此同时重新组织数据，计算每个时刻的 items at risk 和 count of events，列表如下：

Distinct Event Times	\(t_0\)	\(t_1\)	...	\(t_k\)
# of events at \(t_j\)	\(0\)	\(d_1\)	...	\(d_k\)
# at risk at \(t_j\)	\(n\)	\(n_1\)	...	\(n_k\)

此时认为条件概率的估计量为 \(\hat P (T> t_j | T> t_{j-1}) = \frac{n_j - d_j}{n_j}\)，因此 survival function 的估计是：

\[\hat S(t) = \Pi_{j=1}^k \frac{n_j - d_j}{n_j} = \Pi_{j=1} ^n (1- \frac{d_j}{n_j})\]

Example

对于一组生存数据：\(4 \quad 5+ \quad 6 \quad 10+ \quad 6 \quad 2 \quad 6\)

重新排序得到：\(2 \quad 4 \quad 5+ \quad 6 \quad 6 \quad 6 \quad 10+\)

按照 event time 重新组织数据，并计算 Kaplan-Meier Estimator：

Distinct Event Times	0	2	4	6
# of events at \(t_j\)	0	1	1	3
# at risk at \(t_j\)	7	7	6	4
\(\hat S(t)\)	1	\(\frac{6}{7}\)	\(\frac{5}{7}\)	\(\frac{5}{28}\)

Cox Proportional Hazard Model

Theory

在生存分析中数据往往是 skewed 的情况，希望分析变量（比如，患者的年龄，所用的药物，实验组别）对于生存期的影响时，一般对 hazard function 采用 \(\lambda(t) = \lambda_0 (t) e^{X^T\beta}\) 模型，类似于一个 transformed linear regression。

其中，\(\lambda_0(t)\) 作为一个调整参数，只和时间有关，和 explanatory variable \(X\) 无关。

以一元的 explanatory variable \(X\) 为例，不同数据之间的 hazard ratio 是：

\[\frac{\lambda_i(t)}{\lambda_j (t)} = \frac{\lambda_0(t)}{\lambda_0(t)} \frac{exp(\beta_0 + \beta_1 X_i)}{exp(\beta_0 + \beta_1 X_j)} = e^{(X_i - X_j)\beta_1}\]

模型假设是 \((T_i ,C_i)\) 在 \(X_i\) 给定的条件下是条件独立的，这是一个针对 censored data 的假设。

另外一个结论是在时间 \(t\) 处，uncensored 对象 \(i\) 发生 event 的概率是 \(P_i(t) = \frac{ exp(x_i \beta)}{\Sigma_j I(T_j \geq t) exp(x_j \beta)}\)。由此可知，用于估计 \(\beta\) 的 likelihood function 是：

\[L(\beta) = \Pi_{j=1}^ k \frac{ exp(x_i \beta)}{\Sigma_j I(T_j \geq t) exp(x_j \beta)} = \Pi_{j=1} ^n (\frac{ exp(x_i \beta)}{\Sigma_j I(T_j \geq t) exp(x_j \beta)})^{\delta _i}\]

仍然是以一元的 explanatory variable \(X\) 为例，hazard ratio 是 \(HR = \frac{\lambda(t | \text{exposure})}{\lambda(t| \text{not exposure})} = \frac{\lambda_0 (t) exp(\beta_1)}{\lambda_0(t)} = exp(\beta_1)\)

Example

分别有五位患者参与两种试药，生存数据如下所示，希望通过 Cox Model 来比较不同种类的药物效果。

Drug A: \(8+ \quad 132 \quad 52 \quad 220 \quad 75+\)

Drug B: \(45 \quad 80 \quad 90+ \quad 20 + \quad 120\)

此处的唯一一个 explanatory variable 即为用药种类，认为 \(X_i = \begin{cases} 1 \quad \text{i use drug A}\\ 0 \quad \text{i use drug B} \end{cases}\)；

对数据重新排序，得到：\(8+ \quad 20+ \quad 45 \quad 52 \quad 75+ \quad 80 \quad 90+ \quad 120 \quad 132 \quad 220\)，对应的 \(X_i\) 是 \(1,0,0,1,1,0,0,0,1,1\)。

Cox model 设定为 \(\lambda_i (t) = \lambda_0 (t) e^{X ^T\beta} = \lambda_0 (t) \exp(\beta_0 +\beta_1 X_i)\)，\(i=1,2,...,10\)。

此时的 partial likelihood function 是 \(PL(\beta) = \Pi_{i=1}^n (\frac{\exp(\beta^T X_i)}{\Sigma_{j: t_j \geq t_i} \exp(\beta^T X_j)})\)，代入数据有：

\[PL(\beta) = \Pi_{i=1}^{10} \frac{\exp(\beta_0 + \beta_1 X_i)}{\Sigma_{j: t_j \geq t_i} \exp(\beta_0 + \beta_1 X_i)} = \frac{exp(2\beta_1)}{360(1+\exp (\beta_1))^3 (5+4\exp(\beta_1)) (3+4\exp(\beta_1)) (2+3\exp(\beta_1)) (1+3\exp(\beta_1))}\]

对其进行 maximum likelihood 处理即可得到 \(\beta_1\) 的估计值。至于 \(\beta_0\) 的估计其实没什么意义，上下直接消掉了。

Why Cox Model

• 体现了 explanatory variable 带来的影响
• baseline hazard \(\lambda_0 (t)\) 事实上并不重要，一般也用不到
• 得到 proportional hazard function 就可以作图比较是否有一个时间点处发生明显的差异
• 单独对不同组别的 event 发生情况使用 p-value （见下一部分）并不够，需要精确到每个个体的生存时间

Comparison of Survival Curves

Log Rank Test

想要比较两组数据的生存期情况，比如说，试药组和安慰剂组之间有无生存周期上的差别，来判断药物是否有效或者有副作用。这里涉及到的问题都属于比较简单的情况，一般用一个 log rank test 就可以了。

	Distinct Event Times	\(t_1\)	...	\(t_i\)	...	\(t_k\)
Pooled sample	# of Events	\(d_1\)	...	\(d_i\)	...	\(d_k\)
	# at Risk at \(t_i\)	\(n_1\)	...	\(n_i\)	...	\(n_k\)
Sample 1	# of Events	\(d_{11}\)	...	\(d_{1i}\)	...	\(d_{1k}\)
	# at Risk at \(t_i\)	\(n_{11}\)	...	\(n_{1i}\)	...	\(n_{1k}\)
Sample 2	# of Events	\(d_{21}\)	...	\(d_{2i}\)	...	\(d_{2k}\)
	# at Risk at \(t_i\)	\(n_{21}\)	...	\(n_{2i}\)	...	\(n_{2k}\)

有以下的近似分布关系，对近似正态分布做检验即可：

\[ U = \Sigma_{i=1} ^k (d_{1i} - \frac{d_i n_{1i}}{n_i})\]

\[\hat{Var} (U) = \Sigma_{i=1} ^k (\frac{d_i n_{1i} (n_i - d_i)(n_i - n_{1i})}{n_i ^2 (n_i -1)})\]

\[\frac{U}{\hat{Var} (U)^{\frac 1 2}} \sim ^{H_0,asy} N(0,1)\]

P-value mode

对于组织成 contingency table 形式的数据，按照逻辑回归部分的 p-value 方法检验即可。一个数据的例子如下所示：

	Drug A	Drug B
Alive	10	16
Dead	8	4

完结撒花

在这之后讲的东西都不考了，我不想做生统，也不是很感兴趣，不想写了。说到底高中的时候连生物都没选考，怎么会想到来尝试生统概论的（。希望大作业没事。

昨天才在路上遇到了 wty 老师和她男朋友/爱人，相比之下（草）她好娇小好可爱，呜呜呜（虽然感觉她也瞟了我好几眼，不会也认识我吧，不会是这个发型太有辨识度了我又每次都苟在最后一排就被记住了吧（

传统艺能一下（不会现在还有人是 FGO 厨吧，不会吧不会吧（

后记：莫名其妙拿了个 A，实际上期末考大寄特寄好几问没做完，我还在想这课能不挂就挺好的了，不知道是老师太善良还是毕业班同学太多帮忙抬上来的（（只不过之前看 V1ncent19 说他这课是 A+，羡慕（x