【Andrew Gelman多元统计】（基于R）

这段时间导师开虐了，要10月份内搞完统计书。在这里做一下笔记。。直接先从中间开始了，前面的也在写。强行放到一页上Ctrl+F就方便很多。

---

宗旨：并非旨在全面概述线性模型，而是对各种回归模型有所了解可以适合我们在应用程序中看到的不同数据结构。

一、广义线性模型（Generalized linear model）

广义线性模型一般以$\hat{y}=g^{-1}(X\beta )$为形式，模型的关键在于式子中的变换 $g$ 和数据分布$p(y|\hat{y})$。
在线性回归中， $g(u)=u$， $p(y|\hat{y})$服从正态分布，可以从样本中估计出 $\sigma$。
在逻辑回归中， $g^{-1}(u)=logi{t}^{-1}(u)$ ，数据分布由二进制数据的概率定义：$Pr(y=1)=\hat{y}$。

这一章会有一些其他的分布模型：

6.2 泊松分布Poisson
泊松模型（第6.2节）用于计数数据；也就是说，每个数据点 $y_i$ 可以等于0、1、2…。这里常用的变换 $g$ 是对数，因此 $g(u)=exp(u)$ 将连续线性预测变量$X_i\beta$转换为正$y_i$，数据分布为泊松分布。
通常可以向该模型添加参数以捕获过度分散，以判断数据的变化是否超出了仅根据泊松分布所能预测的范围。

6.3 Logistic二项式模型
Logistic二项式模型用于
每个数据点 $y_i$ 表示 $n_i$次尝试中的成功次数（此$n_i$，即数据点 $i$ 的尝试次数，与 $n$ 不同）。 在此模型中，变换 $g$ 为逆对数，数据分布为二项式。

6.4 Probit概率模型
与logistic回归用法相同，但是logit函数由正态累积分布代替，或者等效地由正态分布代替了潜在数据错误中的logistic。

6.5 多项式logit和probit模型
多项式logit和probit模型是对分类数据进行logistic和probit回归。这里分类数据具有两个以上的选择，例如问卷中选择强烈同意、同意、不同意、非常不同意。这些模型使用logit或概率变换和多项式分布，并且需要附加参数来对数据的多种可能性进行建模。
多项式模型进一步分为有序模型（例如，强烈同意，…，强烈不同意）或无序模型（例如，香草，巧克力，草莓，其他）。

6.6 Robust回归模型
Robust回归模型用其他可能偶尔出现极值的分布1（通常是所谓的Student-t系列模型）代替了通常的逻辑模型。

线性建模函数——lm()
glm()——logistic二项式，概率和泊松回归，适当的情况下可以纠正过度分散。
polr()——有序logit和概率回归
mnp包——拟合无序的概率模型
hett包——t模型。（详见Andrew Gelman多元统计附录C。）
另外，各种错误都可以适用于Bug。

（一）Poisson regression, exposure, and overdispersion 泊松回归，暴露和过度分散

在泊松模型中，每个单位 $i$ 对应于观察到 $y_i$ 事件的设置（通常是空间位置或时间间隔）。
与线性和逻辑回归一样，可以使用线性预测变量 $X$ 来解释 $y$ 的变化。在交通事故示例中，这些预测变量可以包括：常数项，对十字路口附近平均交通速度的度量以及指示交叉路口是否有交通信号，基本的泊松回归模型为 $y_i～Poisson({\theta }_i)$.（表示第 $i$ 个路口的交通事故数量 $y_i$ 服从参数为 $\theta$ 的泊松分布）
参数 ${\theta }_i$ 必须为正，因此需要在对数刻度上拟合线性回归。${\theta }_i=exp(X_i\beta )$，即对自变量的方程作对数变换。

可以对系数β求幂，并将其视为乘积效应，例如，假设交通事故模型为
$y_i\sim Poisson(exp(2.8+0.012X_{i1}-0.20X_{i2}))$,

其中 $X_{i1}$ 是附近街道的平均速度（以英里/小时或mph为单位），如果交叉路口有交通信号，则 $X_{i2}$ = 1，否则为0，然后我们可以将每个系数解释如下：

• 常数项给出了回归的截距，即 $X_{i1}$ 和$X_{i2}=0$ 时的预测。由于这是不可能的（没有街道的平均速度为0），我们将不尝试解释常数项。
• $X_{i1}$ 的系数是每增加1 mph交通速度时y的预期差值（对数刻度），因此预期乘数增加值为 $e^{0.012}=1.012$，意思是每单位$X_{i1}$ 的增加将导致该单位下交通事故发生率提高1.2％。若将 $X_{i1}$ 的单位定义为以十英里/小时，它的系数为0.12，那么作指数转化得$e^{0.12}=1.127$，得每十英里每小时的事故发生率提高12.7％。
• $X_{i2}$ 的系数告诉我们，可以发现交通信号的预测差异是事故率乘以 $e^{-0.20}=0.82$，即设置交通信号将使交通事故的概率减少18％。

1 具有$exposure$的泊松回归Poisson regression with an exposure input

在大多数Poisson回归应用中，可以相对于某些基线或“$exposure$”来解释计数，如通过交叉路口的车辆数。量在一般的Poisson回归模型中，我们将 $y_i$ 视为案例数在具有速率 ${\theta }_i$ 和曝光 $u_i$ 的过程中。

$y_i\sim Poisson(u_i{\theta }_i)$

像以前一样，${\theta }_i=exp(X_i\beta )$。$exposure$的对数 $log(u_i)$ 在广义线性模型术语中称为偏移。回归系数 $\beta$ 概括了预测变量与 ${\theta }_i$ 之间的关联（在我们的示例中，为每辆车的交通事故发生率）。
在 $Poisson$ 回归中包括对数化的 $exposure$ 作为预测变量，如模型$y_i\sim Poisson(u_i{\theta }_i)$所示，将 $exposure$ 的对数作为偏移量放入模型中，等同于将其作为回归预测变量包括在内，但系数固定取值1。另一种选择是将其包括为预测变量，并根据数据估算其系数；在某些设置中，这样做可以使数据更好地拟合。在其他设置中，则更简单，只是将其保留为偏移量，以便估算率 $\theta$ 可以更直接地解释。

泊松模型类似于计数数据的二项式模型（参见第6.3节），但适用于略有不同的情况：

• 如果每个数据点 $y_i$ 都可以解释为ni次试验中“成功”的次数，则标准使用伯努利模型或逻辑模型（如6.3节中所述）或其overdispersed generalization。
• 如果每个数据点 $y_i$ 没有自然极限-它不是基于许多独立试验的-那么使用泊松/对数回归模型（如此处所述）或其overdispersed generalization。

🌰：警察按族裔群体停车
对于警察的分析：

• $i$ 的单位是辖区和族裔群体$（i=1，...，n=3\times 75）$。
• 结果 $y_i$ 是该区域中该族裔成员的站点数。
• 暴露 $u_i$ 是刑事司法服务部（DCJS）记录的该地区上一个族裔在上一年被捕的人数。
• 预测变量是地区和种族指数。有74个区指标（例如，区2–75，以区1为基准）和2个种族指标（例如，西班牙裔和白人，以黑人为基准）。

我们分三步说明模型拟合：首先，我们仅使用偏移量和常数项拟合模型：

glm(formula = stops ~ 1, family=poisson, offset=log(arrests)) 

结果：

           coef.est coef.se
(Intercept) -3.4 0.0
n = 225, k = 1
residual deviance = 44877, null deviance = 44877 (difference = 0)

接下来，我们添加种族指标：

glm(formula = stops ~ factor(eth), family=poisson, offset=log(arrests))

结果：

          coef.est coef.se
(Intercept)  -3.30 0.00
factor(eth2)  0.06 0.01
factor(eth3) -0.18 0.01
n = 225, k = 3
residual deviance = 44133, null deviance = 44877 (difference = 744.1)

这两个种族系数在统计上均具有很高的统计学意义，并且偏差减少了744，远远超过了模型中没有种族解释力的情况下预期的2。与基线DCJS类别1（黑人）相比，我们看到了eth2（西班牙裔）的停留次数增加6％，eth3（白人）的停留次数减少18％。

现在，我们添加了75个区：

glm(formula = stops ~ factor(eth) + factor(precinct), family=poisson, offset=log(arrests))

结果：

                   coef.est coef.se 
(Intercept)           -4.03 0.05 
factor(eth)2           0.00 0.01 
factor(eth)3          -0.42 0.01 
factor(precinct)2     -0.06 0.07 
factor(precinct)3      0.54 0.06
...
factor(precinct)75     1.41 0.08

n = 225, k = 77
residual deviance = 2828.6, null deviance = 44877 (difference = 42048.4) overdispersion parameter = 18.2

偏差的减少幅度从44,000减小到2800，这是巨大的-远大于如果区域因素是随机噪声时所预期的74的减少。在控制了区域之后，种族系数有所改变——黑人和西班牙裔（ 第1类和第2类）的停止率大致相同，而白人（第3类）的机会比少数族裔低大约42％（更准确地说，白人的指数系数为 $e^{-0.42}=0.66$ ，因此他们被阻止的机会实际上要低34％-近似 $e^{-\beta }\approx 1-\beta$ 仅在β接近0时才是准确的。）。

因此，控制选区实际上增加了白人和少数族裔在停车率上的差异，我们将在15.1节中进一步探讨这个问题。

我们还可以查看回归中的区域系数，例如，控制种族后的停车率（每DCJS逮捕）在区域2中大约降低6％，在区域3中，$e^{0.54}=1.72$ 倍， …，与基准区域1相比，区域75的 $e^{1.41}=4.09$ 倍。

$exposure$ 的输入

在此示例中，将警察逮捕次数与上一年的逮捕次数进行了比较，因此，如果某个种族的人，如“西班牙”或“白人”，比上一年黑人的逮捕数量要高，那么该族的系数将大于1。 同样，2-75区中，与上一年第1区的逮捕率相比更频繁发生逮捕的那些区，指标的系数将超过1。在第15.1节中，我们将考虑另一种可能的分析方法，该方法使用人口而不是前一年的逮捕作为 $exposure$。

（a）残差与预测值之间的关系 （b）标准化残差与预测值之间的关系

Poisson回归模型中的过度分散测试：
正如模型所预期的那样，残差的方差随预测值的增加而增加。标准残差应具有： 平均值为0，标准偏差为1的性质，（±2处的线表示大约95％的误差范围）。标准残差的方差远大于1，表明存在大量的过度分散。

$Overdispersion$过度离势

泊松回归不提供独立的方差参数 $\sigma$，结果可能会过度分散。在泊松分布下，方差等于均值，标准偏差等于平均值​​的平方根。在模型 $y_i\sim Poisson(u_i{\theta }_i)$ 中，$E(y_i)=u_i{\theta }_i$ ，而 $sd(y_i)=\sqrt{u_i{\theta }_i}$
标准化残差定义如下：

$z_i=\frac{y_i-{\hat{y}}_i}{sd({\hat{y}}_i)}=\frac{y_i-{\hat{u}}_i{\theta }_i}{sd({\hat{u}}_i{\theta }_i)}$

其中 $\hat{\theta }=e^{X_i\hat{\theta }}$。如果泊松模型为真，则zi应该是近似独立的（不是完全独立的，因为在计算所有这些时都使用相同的估计系数 $\hat{\beta }$），均值为0，标准差为1。 如果存在过度分散， $z_i$的绝对值就会更大，反映出超出Poisson模型所预测的额外变化。

我们可以通过计算n个标准化残差的平方和 ${\sum }_{i=1}^n{z}_i^2$ 并将其与 ${\chi }_{n-k}^2$分布进行比较来测试经典Poisson回归中的超分散，这是我们在模型下的预期（使用n-k） 而不是用n个自由度来解释k个回归系数的估计）。${\chi }_{n-k}^2$分布具有平均值 $n-k$ ，因此该比率为
$estimated-overdispersion=\frac{1}{n-k}{\sum }_{i=1}^n{z}_i^2$

例如，警察局的经典泊松回归具有 $n=225$ 个数据点和 $k=77$ 个线性预测因子，下图绘制了残差 $y_i-\hat{y}i$ 和标准化残差 $z_i=(y_i-{\hat{y}}_i)/sd({\hat{y}}_i)$ ，如 Poisson回归模型的预测值的函数。正如Poisson模型所期望的那样，残差的方差随着预测值的增加而增加，标准化残差的方差近似恒定。但是标准化残差的方差远远大于1，表明严重的过度分散。

（a）残差与预测值之间的关系 （b）标准化残差与预测值之间的关系

在R中编程过度分散测试：

yhat <- predict (glm.police, type="response")
z <- (stops-yhat)/sqrt(yhat)
cat ("overdispersion ratio is ", sum(z^2)/(n-k), "\n")
cat ("p-value of overdispersion test is ", pchisq (sum(z^2), n-k), "\n")

标准化残差的平方和为 $n_i=1,z_i^2=2700$，相比之下，预期的 $n-k=148$ 。估计的超分散因子为 $\frac{2700}{148}=18.2$，p值为1，表明 来自 ${\chi }_{148}^2$ 分布的随机变量最大可能为2700的概率基本上为零。总而言之，警察站点的数据过度分散了18倍，这是巨大的，甚至过分散系数2也被认为是大的，并且也具有统计意义。

调整过度分散的参考方法

在这个例子中，对超分散的基本修正是将所有回归标准误差乘以$\sqrt{18.2}=4.3$，幸运的是，事实证明我们的主要推论并未受到严重影响，主要关注的参数是${\alpha }_3$ 白人相对于黑人的逮捕率的对数-之前估计为-0.42±0.01，现在变为-0.42±0.04。

转换回原始比例后，白人相对于黑人被逮捕的概率为66％ ，大约50％的可能在
$e^{-0.42\pm (2/3)0.04}$ $=[0.64,0.67]$之间，大约95％的可能在 $e^{-0.42\pm 2\cdot 0.04}$$=[0.61，0.71]$之间。

拟合过度分散的泊松或负二项式模型

更简单地讲，我们可以使用拟泊松族拟合过度分散的模型：

glm(formula = stops ~ factor(eth) + factor(precinct), family=quasipoisson, offset=log(arrests))

结果：

                   coef.est coef.se 
(Intercept)           -4.03 0.21 
factor(eth)2           0.00 0.03 
factor(eth)3          -0.42 0.04 
factor(precinct)2     -0.06 0.30 
factor(precinct)3      0.54 0.24
...
factor(precinct)75     1.41 0.33

n = 225, k = 77
residual deviance = 2828.6, null deviance = 44877 (difference = 42048.4) overdispersion parameter = 18.2

这个模型可以写作：

$y_i\sim overdispersedPoisson(u_{i}exp(X_i\beta ),\omega )$

其中 $\omega$ 是超分散参数（在这种情况下估计为18.2）严格来说，“超分散泊松”不是单个模型，而是描述了任何数据方差为 $\omega$ 乘以均值的计数数据模型。如果ω= 1，则overdispersed Poisson模型就是Poisson模型。在这种情况下，通常使用的特定模型是负二项式分布：

$y_i\sim Negative-binomial(mean=u_{i}exp(X_i\beta ),overdispersion=\omega ).$

不幸的是，负二项式分布的传递不是基于其均值和超分散，而是根据参数a和b表示的，其中分布的均值为a / b，超分散为1 +1 / b。 在拟合此类模型时进行参数化，并且可以通过模拟来自拟合模型的数据集并检查它们看起来像实际数据来进行双重检查（详见第8.3节）。

我们返回警察停靠站示例，在15.1节中使用多级模型纠正过度分散。

6.3 逻辑二项式模型

第5章讨论了二进制（是/否或0/1）数据的逻辑回归，该逻辑模型也可以用于计数数据，使用二项分布来模拟指定范围内可能的成功次数，成功的可能性适合逻辑回归。

计数数据的二项式模型，应用于死刑判决。

我们在1973-1995年这23年间的34个州中的每个州推翻了死刑判决的比例的研究中，说明了二项式逻辑回归。该分析的单位是34×23 = 784个州， 年（实际上，我们的分析中只有n = 450个州年，因为不同的州自1973年以来在不同的时间重新开始了死刑）。对于每个州年 $i$，我们将 $n_i$ 标记为死刑年数。 指出当年的收入，以 $y_i$ 为单位，这些裁决随后被高等法院推翻。
$y_i\sim Binomial(n_i,p_i)$
$p_i=logi{t}^{-1}(X_i\beta )$ （6.6）
其中X是预测变量的矩阵。首先，我们使用
常数项
33个州指标
时间趋势（即变量1等于1973年、2等于1974年、3等于1975年，依此类推）。
该模型也可以写成
$y_{st}\sim Binomial(n_{st},p_{st})$
$p_{st}=logi{t}^{-1}(\mu +{\alpha }_s+\beta t)$
在状态的下标 $i$ 和时间的下标 $t$ （即时间-1972）之间，我们更喜欢前者（6.6），因为它具有更大的通用性，理解两个公式是有意义的。

过度离势

当采用逻辑回归对数据进行计数时，事实上，通常情况下，数据的变化可能会超出模型所解释的范围，因此出现了过度分散的问题，因为逻辑回归模型没有方差参数σ。更具体地说，如果数据y服从参数为n的和p的二项分布，那么y的平均值为np，y的标准偏差为np(1- p)。
与模型（6.4）一样，我们将每个数据点i的标准化残差定义为：

$z_i=\frac{y_i-{\hat{y}}_i}{sd({\hat{y}}_i)}=\frac{y_i-n_i{\hat{p}}_i}{\sqrt{n_i{\hat{p}}_i(1-{\hat{p}}_i)}}$

其中 $p_i=logi{t}^{-1}(X_i\hat{\beta })$。如果二项式模型为真，则 $z_i$ 应该近似独立，均值均值为0，标准差为1。与泊松模型一样，我们可以计算出估计过度离势 $\frac{1}{n-k}{\sum }_{i=1}^n{z}_i^2$，并通过比较 ${\sum }_{i=1}^n{z}_i^2$ 与 ${\chi }_{n-k}^2$分布来对过度离势进行规范的检验（此处的n表示数据点的数量，与模型（6.6）和（6.7）中的 $n_i$ 表示形式无关，后者指的是州 $i$ 年的病例数）。
实际上，过度分散在逻辑回归中非常普遍（或6.2节中讨论的应用于计数数据的Poisson回归）。
因此也可以采用更一般的过度分散模型，使偏差为 $\omega np(1-p)$，其中ω> 1被称为超分散参数。当ω= 1时，超分散模型就变成了二项式对数回归。

调整过度分散的参考方法

与Poisson回归一样，对过度离势的一种简单校正是将所有系数估算的标准误差乘以估算超分散的平方根（6.5）。如果不进行此调整，则置信区间将太窄，并且推断将过度离势的二项式回归可以使用glm()函数和quasibinomial(link="logit")族在R中拟合，相应的分布是beta-binominal。

二进制数据模型作为计数数据模型的特例

如第5章所述，二进制数据的逻辑回归是二项式（6.6）的特例，其中所有 $i$ 均为 $n_i\equiv 1$.在二进制模型中不会发生单个数据点级别的过度分散，这就是为什么我们没有 在第5章中介绍过度分散的模型。

计数数据模型是二进制数据模型的特例

相反地，二项式模型（6.6）也可以以二进制数据形式（5.1）表示。可以通过将 $n_i$ 个情况中的每一个视为单独的数据点。该扩展回归的样本量为 ${\sum }_i{n}_i$，数据点分别为0和1：每个状态 $i$ 对应于 $y_i$ 和 $n_i-y_i$。最后，将X矩阵展开为具有 $n_i$ 行，其中原始X矩阵的第 $i$ 行成为展开矩阵中的 $n_i$ 个相同行。这时，可以在多级模型中包括过度离势，通过为原始测量创建索引变量（在死刑示例中，采用值1，…，450）并在该级别包含变化的系数或误差项。

平均值为0且标准偏差为1.6的正态密度函数图

平均值为0且标准偏差为1.6的正态密度函数，在大多数实际情况下，与逻辑密度是相同的，因此我们可以将Probit模型中的系数解释为逻辑回归系数除以1.6。

6.4 概率回归：正态分布的潜在数据

Probit模型与logit相同，只是它用正态分布代替了logistic，我们可以将模型直接写为
$Pr(y_i=1)=\Phi (X_i\beta )$
$\Phi$是正态累积分布函数。
$y_i=\{\begin{array}{l}1if{z}_i>0\\ 0if{z}_i<0\end{array}$
$z_i=X_i\beta +ϵ$
${ϵ}_i～N(0,1)$ （6.8）
即潜在误差的正态分布，平均值为0，标准差为1。
更一般而言，模型可以具有误差方差，因此（6.8）的最后一行被替换为${ϵ}_i～N(0,{\sigma }^2)$.
但是 $\sigma$ 是不确定的，因为如果将 $\sigma$ 乘以某个常数 $c$ 然后再将向量 $\beta$ 乘以 $c$ 则模型不变，因此我们需要对参数进行一些限制，标准方法是将 $\sigma$ 固定为1 （6.8）。

概率还是logit

Probit模型接近对数，其残差标准偏差设置为1而不是1.6，因此，Probit回归中的系数通常接近 Logistic回归系数除以1.6。例如，这是的Logistic回归模型的Probit版本：

glm(formula = switch ~ dist100, family=binomial(link="probit"))

结果：

            coef.est coef.se
(Intercept)     0.38 0.04 
dist100        -0.39 0.06
n = 3020, k = 2
residual deviance = 4076.3, null deviance = 4118.1 (difference = 41.8)

对于我们所看到的示例，选择logit或probit模型是一种品味或便利性的问题，例如，在解释probit模型的潜在正态误差时，当我们看到probit回归系数时，只需将它们乘以1.6即可。 获得等效的逻辑系数，例如，我们刚刚拟合的模型 $Pr(y=1)=\Phi (0.38-0.39x)$，基本上等于逻辑模型$Pr(y=1)=logi{t}^{-1}(1.6(0.38-0.39x))=logi{t}^{-1}(0.61-0.62x))$。

6.5 多项式回归

有序和无序分类结果
逻辑和概率回归可以扩展到多个类别，可以有序排列或无序排列，有序分类结果的示例包括民主党，独立，共和党；是的，也许是，否；总是，经常，有时，很少，从不。 无序的分类结果包括自由主义，劳动，保守主义；足球，篮球，棒球，曲棍球；火车，公共汽车，汽车，步行；白人，黑人，西班牙裔，亚洲人，其他；我们首先讨论有序类别，包括扩展示例 ，然后简要讨论无序分类变量的回归模型。

有序多项式logit模型

考虑一个可以取值1,2，…，K的分类结果y。有序逻辑模型可以用两种等效的方式编写。首先我们将其表示为一系列逻辑回归：
$Pr(y>1)=logi{t}^{-1}(X\beta )$
$Pr(y>2)=logi{t}^{-1}(X\beta -c_2)$
$Pr(y>3)=logi{t}^{-1}(X\beta -c_3)$
…
$Pr(y>K-1)=logi{t}^{-1}(X\beta -c_{k-1})$ （6.9）
参数$c_k$（由于我们稍后将解释的原因称为阈值或切点）是递增的：$0=c_1<c_2<\cdot \cdot \cdot <c_{K-1}$，因为假定（6.9）中的概率严格减小，而所有K个结果都具有非零发生概率。由于 $c_1$ 定义为0，因此除 $\beta$ 之外，具有 $K$ 个类别的模型还具有 $K-2$ 个自由参数 $c_k$ 。这对于通常的logistic回归而言是合理的，因为K = 2 仅需估计 $\beta$。
可以使用系数 $\beta$ 同时使用最大似然估计切点 $c_2，...，c_{K-1}$，但是对于某些数据集，参数可能无法确定，就像对二进制数据进行逻辑回归一样。
可以减去（6.9）中的表达式以获得单个结果的概率：

$Pr(y=k)=Pr(y>k-1)-Pr(y>k)=logi{t}^{-1}(X\beta -c_{k-1})-logi{t}^{-1}(X\beta -c_k).$

有序分类逻辑模型中的切点的说明。在此示例中，有K = 4个类别，切点为$c_1=0，c_2=0.8，c_3=1.8$。这三个图说明了对应于三个的潜在结果 $z$ 的分布，线性预测变量 $X\beta$ 的值各不相同。每个切点都显示结果y等于1、2、3或4。

带有切点的潜在变量解释

通过将潜在变量公式（5.4）归纳为K个类别，最容易理解有序分类模型：
$y_i=\{\begin{array}{l}1if{z}_i<0\\ 2if{z}_i\in (0,c_2)\\ 3ifzi\in (c_2,c_3)\\ ...\\ K-1ifzi\in (c_{K-2},c_{K-3})\\ Kif{z}_i>c_{K-1}\end{array}$ （6.10）

$z_i=X_i\beta +{ϵ}_i$
具有（5.4）中的逻辑分布的独立错误i。
潜在变量模型显示，两个相邻切点 $c_{k-1}$，$c_k$之间的距离将影响$y=k$的概率。我们还可以看到，如果线性预测变量$X\beta$足够高，y几乎可以肯定取尽可能高的值，如果$X\beta$足够低，则y几乎可以肯定等于最低可能值。

示例：可存储的投票

我们以实验性经济学为研究对象，对有序分类数据分析进行了说明，该主题涉及“可投票数”。这个例子有些复杂，说明了有序逻辑模型的使用和潜在局限性。 总共有两个游戏。在每个游戏中，一组k个玩家在两个问题上都可以投票，每个玩家总共获得4票。 他们可以选择投1，2或3票，其余票投第二期，每期的获胜方由多数票决定，届时获胜方的每个参与者都将获得正数收益，这些收益来自于在间隔[1，100]上均匀分布。

为了增加预期的收益，玩家应该遵循一种策略，针对潜在收益较高的问题投更多票，在进行此实验的方式中，告知参与者可能的收益分配，并告知他们每项收益的潜在收益。因此，在选择要在第一期投票中投多少票时，每个玩家仅知道自己在该投票中的潜在收益，然后，告知玩家他们在第二期投票中的潜在收益， 但由于此时他们将自动花费所有剩余的选票，因此不涉及任何选择，因此，根据潜在的收益x，玩家的策略可以概括为他们选择的初始选票y = 1、2或3。

可存储选票研究中的一些示例数据：垂直线显示了估计的临界点，曲线显示了使用有序逻辑回归估计的预期响应；左两个图显示了与模型相当吻合的数据；其他两个图则在某种程度上拟合了模型 但并不完美。上图绘制了实验中大约100位学生中有6位学生的回答，并选择了这6位来代表几种不同的数据模式，我们并不惊讶地看到这些回答通常是单调的，也就是说，当他们的潜在收益更高时，学生倾向于花费更多的选票，但有趣的是看到了所选择的各种近似单调策略。
大多数个人的行为可以归纳为三个参数：1和2之间的临界点，2和3之间的临界点以及这些划分的模糊性。这两个临界点表征了所选的单调策略， 划分的敏锐度表明了该策略的一致性。
有序逻辑模型的三个参数化，使用有序logit建模响应很方便，其参数化与模型（6.10）略有不同，以适应我们对单调策略的理解。

$y_i=\{\begin{array}{l}1if{z}_i<c_{1.5}\\ 2if{z}_i\in (c_{1.5},c_{2.5})\\ 3if{z}_i>c_{2.5}\end{array}$                  ~~~~~~~~~~~~~~~~                 （6.11）

$z_i～logistic(x_i,{\sigma }^2)$
在此模型中，切点 $c_{1.5}$ 和 $c_{2.5}$ 处于数据x的1-100范围内，而误差的比例σ对应于切点的模糊性。
该模型具有与常规参数化（6.10）相同的参数数量-消失了两个回归系数，同时增加了一个附加的自由割点和误差方差，这里是模型（6.10），其中K = 3个类别，一个 预测变量x。

$y_i=\{\begin{array}{l}1if{z}_i<0\\ 2if{z}_i\in (0,c_2)\\ 3if{z}_i>c_2\end{array}$                  ~~~~~~~~~~~~~~~~                 （6.12）
$z_i=a+\beta x+{ϵ}_i$
其中${ϵ}_i～logistic(0,1)$

另外一种形式的模型保留了两个不同的分割点，去掉了常数项。

$y_i=\{\begin{array}{l}1if{z}_i<c_{1/2}\\ 2if{z}_i\in (0,c_{2/3})\\ 3if{z}_i>c_{2/3}\end{array}$                  ~~~~~~~~~~~~~~~~                 （6.13）
$z_i=\beta x+{ϵ}_i$

其中${ϵ}_i～logistic(0,1)$
这三个模型实际上是等价的，在（6.13）中为$z_i/\beta$，在（6.12）中为$(z_i-\alpha )/\beta$
对应于（6.11）中的 $z_i$ 。匹配的参数如下：

我们更喜欢参数化（6.11），因为我们可以直接将c1.5和c2.5解释为输入x的标度上的阈值，而σ对应于从1到2以及从2到3的过渡的渐进性。 但是，使用标准参数化（6.12）和（6.13）拟合模型非常方便，因此能够在模型之间来回切换会很有帮助。
在R中拟合模型。我们可以使用polr（“比例比对数逻辑回归”）函数拟合有序logit（或概率）模型，该函数是R中MASS软件包的一部分。我们用来自可储存票研究数据进行说明 ：

polr (factor(y) ~ x)

结果：

Coefficients: x
0.07911799
Intercepts:
1|2      2|3
1.956285 4.049963

从输出中我们可以看到它拟合了形式为（6.13）的模型，其估计值$\hat{\beta }=0.079$，${\hat{c}}_{1|2}=1.96$和${\hat{c}}_{2|3}=4.05$。使用这三个模型的表转换为模型（6.11）， 我们得到${\hat{c}}_{1.5}=1.96/0.079=24.8$，${\hat{c}}_{2.5}=4.03/0.079=51.3$，$\hat{\sigma }=1/0.079=12.7$。

显示拟合的模型图6.4显示了切点c1.5，c2.5和预期投票数 E(y) 作为 x 的函数，这是根据多个数据的估计得出的。
根据模型（6.11），预期投票可以写成

这样就能计算出顺序分类下y所获得的预期选票了。其中$logi{t}^{-1}(x)=\frac{e^x}{(1+e^x)}$。表达式（6.14）看起来很复杂，但很容易在R中作为函数进行编程：

expected <- function (x, c1.5, c2.5, sigma){ 
p1.5 <- invlogit ((x-c1.5)/sigma)
p2.5 <- invlogit ((x-c2.5)/sigma)
return ((1*(1-p1.5) + 2*(p1.5-p2.5) + 3*p2.5))
}

然后，可以将图6.4中的数据，切点和曲线绘制如下：

plot (x, y, xlim=c(0,100), ylim=c(1,3), xlab="Value", ylab="Vote") 
lines (rep (c1.5, 2), c(1,2))
lines (rep (c2.5, 2), c(2,3))
curve (expected (x, c1.5, c2.5, sigma), add=TRUE)

在显示了这些估计值之后，下一步是研究参数在总体中的分布，以了解学生采用的策略的范围，在这种情况下，数据具有多层次的结构-每个观测值30次观察数位学生的学习方式，我们将在多级广义线性模型这一章的15.2节中进一步研究该示例。

建模有序分类数据的替代方法

排序的分类数据可以通过几种方式建模，包括：

• 如前所述，具有K-1切点参数的有序logit模型。
• 概率模型中的相同模型。
• 简单的线性回归（可能先对结果值进行简单的转换），如果类别的数量很大并且可以认为它们之间的间距相等，则这是一个好主意，前提是实际上类别的合理范围是例如，如果等级为1到10，但实际上总是等于9或10，则线性模型可能无法很好地工作。
• 单独的逻辑回归-y = 1或者y = 2，…，K的逻辑回归模型；若y≥2，则y = 2活着y = 3，…，K，以此类推，如果y≥K -1且y = K -1或y = K。或可以使用Probit模型进行设置。单独的逻辑（或Probit）回归具有灵活地拟合数据的优势，但缺点是我们无法简单地描述切点模型的潜在变量解释。
• 最后，我们在第6.6节中讨论的robit回归是logistic回归的竞争者，后者考虑了偶然的异常数据，如离群值。
  
  图6.5 使用逻辑回归拟合的假设数据：（a）具有“异常值”的数据集（左上角附近的意外y = 1值。；（b）从逻辑回归模型模拟的数据，没有异常值。 在每个图中，虚线和实线分别显示拟合的logit和robit回归。在每种情况下，尤其是对于受污染的数据，robit line都更陡峭，因为它有效地减轻了看起来不适合模型的点的影响。

无序分类回归

正如在6.5节开始时讨论的那样，有时将离散结果建模为无序是适当的，我们研究中出现的一个例子就是换井问题，如5.4节所述，拥有不安全井的家庭可以选择改用安全的井，但实际的替代方案更为复杂，可以概括为：（0）什么都不做，（1）改用现有的私人井，（2）改用现有的社区井，（3）自己安装新井。如果将它们编码为0、1、2、3，则可以对Pr（y≥1），Pr（y≥2 | y≥1），Pr（y = 3 | y≥2）进行建模。选项可以被认为是以某种方式排序的，应用排序的多项式logit或probit模型是没有意义的，因为不同的因素可能会影响三个不同的决策，相反，拟合单独的logit（或probit）更有意义模型要由以下三个组成部分决定：（a）是切换还是不执行任何操作？（B）如果切换，则切换到n现有的井或自己建造新的井（C）如果您改用现有的井，是私人井还是社区井？有关此重要模型类别的更多信息，请参见本章末尾的参考资料。

6.6 使用t模型的稳健回归

t分布而不是正态分布

当回归模型偶尔可能会有非常大的误差时，对模型对误差项的回归通常更适合使用Student-t而不是正态分布。回归的基本形式不变：$y=X\beta +ϵ$，但对于 $ϵ$ 的分布不同，因此估计 $\beta$ 的方法略有不同（请参见第18章中关于最大似然估计的讨论），并且预测的分布也有所不同。使用 $t$ 模型估计的回归被认为是稳健的，因为系数估计值受单个外部数据点的影响较小，可以使用 $R$ 中 $hett$ 包中的tlm()函数来拟合t误差的回归。

Robit代替logit或probit

Logistic回归（和本质上等效的Probit回归）可以灵活方便地对二进制数据进行建模，但是它们可能会遇到离群值问题，离群值通常被认为是极端观察。图6.5a展示了从logistic回归模拟的数据，其中极点从0切换到1。在logistic模型的上下文中，对于某些x的值，y = 1的观察将是极不可能，但在实际数据中肯定会发生这种“错误分类”，因此，即使该模型不完全合适，该图也代表了我们可能适合进行逻辑回归的数据类型。
有关数据点异常的逻辑回归的另一个说明，请参见图6.4的右上图，这是一个具有三个结果的示例；为简单起见，我们将注意力集中在二进制结果上。
通过归纳潜在数据公式（5.4），可以方便地“加强”逻辑回归：

$y_i=\{\begin{array}{l}1if{z}_i>0\\ 0if{z}_i<0\end{array}$

$z_i=X_i\beta +{ϵ}_i$

给出潜在误差 $ϵ$ 的 $t$ 分布：

${ϵ}_i～t_v(0,\frac{v-2}{v})$

根据数据估算自由度参数 $\nu >2$ 并缩放 $t$ 分布，使其标准偏差等于1。
$i$ 的 $t$ 模型允许偶尔进行意外的预测——如对于线性预测变量 $X\beta$ 的高度负值，$z$ 为正值，或者与之相反的极端情况。图6.5a通过模拟的“污染”数据集进行了说明：实线显示 $Pr(y=1)$ 是拟合robit回归的x的函数，比拟合logistic模型要陡得多。$t$ 分布有效地降低了不一致的数据点的权重，因此该模型更好地拟合了数据的主要部分。
图6.5b显示了来自逻辑模型的数据：这里，robit模型接近logit，这是有道理的，因为它没有发现很多需要省略的点。
在数学上，可以将robit模型视为probit的泛化和logit的近似泛化，probit对应于自由度 $\nu =\infty$，而 logit 与 ν=7 的 robit 模型非常接近。

6.7 建立更复杂的广义线性模型

到目前为止，我们已经考虑的模型可以在实践中处理许多回归问题，对于连续数据，我们从具有正误差的线性回归开始，考虑第4章中讨论的适当转换和交互作用，并为偶尔有大误差的数据切换到错误模型。 对于二进制数据，我们使用logit，probit或robit，再次转换输入变量并考虑残差图，如第5章中所述。对于计数数据，起点是过度分散的二项式和泊松分布，对于离散结果则为我们可以使用不止两个类别来拟合有序或无序多项式logit或概率回归，在这里我们简要描述一些可以考虑其他模型的情况。

混合离散/连续数据

收入是一个具有离散和连续方面的结果变量的示例，在第4章的收入和身高回归中，我们通过删除所有收益为零的受访者来对数据进行预处理，但是通常可以对变量建模。例如收益，分两个步骤：首先对 $Pr(y>0)$ 进行逻辑回归，然后 $log(y)$ 进行线性回归，并以 $y>0$ 为条件。然后，必须分两步对此类模型进行预测，这最方便使用模拟（请参见第7章）。
在分几个步骤对结果进行建模时，有时需要进行编程工作才能将推论转换为原始数据规模，例如，在用于预测身高和性别的收入的两步模型中，我们首先使用逻辑回归来预测收益是否为正：

earn.pos <- ifelse (earnings>0, 1, 0)
fit.1a <- glm (earn.pos ~ height + male, family=binomial(link="logit"))

结果：

            coef.est coef.se
(Intercept)    -3.85 2.07
height          0.08 0.03
male            1.70 0.32
n = 1374,residual deviance = 988.3,k = 3
null deviance = 1093.2 (difference = 104.9)

然后，我们将线性回归拟合为正收益的对数：

log.earn <- log(earnings)
fit.1b <- lm (log.earn ~ height + male, subset = earnings>0)

结果：

              coef.est coef.se
(Intercept)       8.12 0.60
height            0.02 0.01
male              0.42 0.07
n = 1187, k =3
residual sd = 0.88 R-Squared = 0.09

因此，举例来说，身高66英寸的女性获得正收入的概率为 $logi{t}^{-1}(-3.85+0.08\cdot 66+1.70\cdot 0)=0.81$，即81％的机会。如果她的收入为正， 它们的预测值是 $e^{8.12+0.02\cdot 66+0.42\cdot 0}=12600$。
潜在数据模型。混合数据建模的另一种方法是通过潜在数据，例如，假设“基本”收入水平 $z_i$（即我被雇用后将拥有的收入），只有 $y_i>0$ 时才能观察到。 Tobit回归就是这样一种模型，在计量经济学中很流行。

蟑螂和zero-inflated Poisson 模型

二项式和泊松模型及其过度分散的概括都可以用潜在的连续概率或事件发生率来表示，但是有时，潜在率本身具有离散的方面。
例如，在对城市公寓中蟑螂侵扰的研究中，我为每个公寓设置了几天的陷阱，我们将 $u_i$ 标记为陷阱天数，将 $y_i$ 标记为被捕获的蟑螂数，目标是在给定预测变量X（包括公寓居民的收入和种族，邻里指标以及公寓质量的衡量指标）的情况下预测蟑螂出没，我们将从模型开始：
$y_i\sim overdispersedPoisson(u_i{e}^{X_i\beta },\omega )$             ~~~~~~~~~~~             (6.16)

但是，实际数据可能有比模型预测结果有更多的零（即，公寓 $i$ 的蟑螂数为 $y_i=0$ ）。一个自然的解释是，有些公寓的确有为零（或非常接近零） 的蟑螂数量，还有一些公寓蟑螂的数量只是因为离散计数而写作零，零膨胀模型zero-inflated Poisson 将（6.16）放入混合模型中：

$y\{\begin{array}{l}=0if{S}_i=0\\ ～overdispersedPoisson(u_i{e}^{X_i\beta },\omega )if{S}_i=1\end{array}$

在这里，Si可以指示我的公寓中是否有蟑螂，并且可以使用逻辑回归建模：

$Pr(S_i=1)=logi{t}^{-1}(X_i\gamma )$

其中 $\gamma$ 是模型这一部分的一组新的回归系数。估计此两阶段模型并不简单-没有观测到 $S_i$ ，因此无法直接估计 $\gamma$ ；并且我们不知道哪个零观测值对应于 $S_i=0$ 且与泊松分布的结果相对应，因此我们无法直接估计 $\beta$ 。已经编写了一些 $R$ 函数来拟合此类模型，并且也可以使用 $Bugs$ 对其进行拟合。

在我们的实际示例中，过度分散的Poisson模型在预测零的数量方面做得很合理；请参见第161页。但是在其他类似的数据集中，零膨胀模型既可以理解数据，又可以很好地拟合数据，因此在此进行介绍。

其他模型

线性，逻辑和泊松模型的基本选择，以及这些模型的取消以及它们的过度分散，鲁棒性和多项式概括，可以解决许多回归问题，但是，其他分布形式已用于特定类型的数据；这些包括用于等待时间数据的指数模型，伽玛模型和威布尔模型，以及用于生存数据的危害模型更一般而言，非参数模型（包括广义加性模型，神经网络和许多其他模型）已被开发出来，以超越广义线性建模框架。 允许输入和数据之间的数据拟合非线性关系。

6.8 建设性选择模型

到目前为止，我们已经将回归建模作为一种描述性工具，用于研究如何在给定输入变量的情况下预测结果，这是一种完全不同的方法，有时适用于选择数据，例如前文所述的关于逻辑回归和广义线性的示例模型，就是将决策建模为目标或效用的平衡。
我们以孟加拉国的油井转换为例（参见第5.4节）证明了这一想法，我们如何理解距离，砷含量与转移决定之间的关系？有意义的是，砷含量较高的人更有可能转移，但是我们应该期望什么系数值？该关系应该是对数还是线性尺度？实际健康风险被认为是砷浓度呈线性；这是否意味着对数模型不合适？可以使用怎样的方法解决这些问题？

要建立选择模型，我们必须指定一个值函数，该值函数表示一个决策相对于另一个决策的偏好强度，在这种情况下，表示相对于不切换的切换偏好。 在上例中，正值对应于转换的偏好，负值导致不转换，因此该模型类似于逻辑回归的潜伏数据解释；实际上，该模型是一种特殊情况，例如我们将在这里看到。

逻辑或概率回归作为一维选择模型

有简单的一维选择模型，可以用一个预测变量将其简化为概率回归或对数回归，正如我们通过将给定距离切换到最近的井的模型所说明的那样，对数回归为

glm(formula = switch ~ dist100, family=binomial(link="logit")) 
           coef.est coef.se
(Intercept)    0.61 0.06 
dist100       -0.62 0.10
n = 3020, k = 2
residual deviance = 4076.2, null deviance = 4118.1 (difference = 41.9)

现在让我们从首要原则考虑它作为决策问题。

$a_i=从不安全的井切换到安全的井的好处$

$b_i+c_i{x}_i=转换到距离xi远的新井的成本$

我们假设一种效用理论，其中收益（降低疾病风险）可以与成本（无法再使用自己的油井带来的不便，加上距离所成比例的额外努力）所能表达的规模相同。
Logit模型，在本实用新型下，如果 $a_i>b_i+c_i{x}_i$，则家庭 $i$ 会切换，但是，我们无法直接测量 $a_i$，$b_i$ 和 $c_i$ 的值，我们可以从数据中学到的就是切换为 $a$ 的概率。 $x_i$ 的函数；这样就将 $a_i$，$b_i$，$c_i$视为随机变量，其分布由总体中这些参数的（未知）值确定。

$Pr(switch)=Pr(y_i=1)=Pr(a_i>b_i+c_i{x}_i)$      ~~~~     $(6.17)$

表达式（6.17）可以写成

$Pr(y_i=1)=Pr(\frac{a_i-b_i}{c_i}>x_i)$

重新表达很有用，因为它会将所有随机变量放在同一位置，并揭示 $y$ 和 $x$ 之间的总体关系取决于 $\frac{a-b}{c}$ 在总体中的分布。
为了方便起见，标签 $d_i=\frac{a_i-b_i}{c_i}$：切换到相邻井的净收益，除以每口到新井的距离成本，如果 $d_i$ 在总体中具有逻辑分布，并且 $d$ 与 $x$ 无关，则 $Pr(y=1)$ 将具有对 $x$ 的逻辑回归的形式，正如我们将在此处显示的那样。
如果 $d_i$ 具有中心 $\mu$ 和比例为 $\sigma$ 的逻辑分布，则 $d_i=\mu +{\sigma }_i$，其中 $i$ 为单位逻辑密度；

这只是一个系数为 $\frac{\mu }{\sigma }$ 和 $\frac{-1}{\sigma }$ 的逻辑回归，然后我们可以拟合该逻辑回归并求解 $\mu$ 和 $\sigma$，例如，良好转换模型 $Pr(y=1)=logi{t}^{-1}(0.61-0.62x)$，对应于$\frac{\mu }{\sigma }=0.61$和$\frac{-1}{\sigma }=-0.62$；因此 $\sigma =1/0.62=1.6$ 和 $\mu =0.61/0.62=0.98$。图6.6显示了 $d$ 的分布以及 $Pr(d>x)$ 的曲线是 $x$ 的函数。

图6.6（a）总体中 $d_i=\frac{a_i-b_i}{c_i}$的假设逻辑分布，以及（b）切换给定距离的概率的对应logistic回归曲线，它们都对应于模型 $Pr(y_i=1)=Pr(d_i>x_i)=logi{t}^{-1}(0.61-0.62x)$。（b）中曲线的暗部对应于井转换数据中x的范围（以100米为单位的距离）；

Probit模型：从效用参数的正态分布开始，得到类似的模型：$d〜N(\mu ,{\sigma }^2)$。

这只是一个Probit回归，模型 $Pr(y=1)=\Phi (0.38-0.39x)$对应 $\mu /\sigma =0.38$和 $-1/\sigma =-0.39$；因此 $\sigma =1/0.39=2.6$ ,$\mu =0.38/0.39=0.98$。图6.7显示了该模型，该模型与图6.6中所示的逻辑模型几乎相同。

图6.7（a）假想的正态分布 $d_i=(a_i-b_i)/c_i$，均值为0.98，标准差为2.6，并且（b）给出了切换给定距离的概率的相应概率回归曲线，它们都对应于模型$Pr(y_i=1)=Pr(d_i>x_i)=\Phi (0.38-0.39x)$与图6.6进行比较。

选择模型，离散数据回归和潜在数据Choice models, discrete data regressions, and latent data

通常将逻辑回归和广义线性模型作为估计不同结果y的方法（给定预测变量x）。拟合模型代表整个总体，模型中的“错误”通过概率而不是0或 1（因此，数据点与拟合曲线之间的间隙，例如下图）。

拟合逻辑回归的图形表示，
$Pr(switchingwells)=logi{t}^{-1}(0.61-0.62\cdot dist100)$
相比之下，选择模型是在个人的层次上定义的，正如我们在很好的转换示例中看到的那样，其中每个家庭 $i$ 都有自己的数据 $X_i$，$y_i$ 和自己的参数 $a_i$，$b_i$ ，$c_i$来确定它的效用函数，从而决定是否切换。

Logit或Probit回归作为多维选择模型

通过将井水中的砷含量纳入模型中，我们可以将 well-switching model 扩展到多个维度。
$a_i\times (As{)}_i=$从砷含量不安全的井切换到安全井的好处
（将收益与当前砷含量成正比是有意义的，因为据信风险基本上与累积接触砷。）
$b_i+c_i{x}_i=$转换到距离xi远的新井的成本。
如果$a_i\times (As{)}_i>b_i+c{x}_i$，则家庭 $i$ 应进行转移。因此，决策取决于家庭的砷水平$(As{)}_i$，其到最近井的距离 $x_i$ 及其效用参数 $a_i$，$b_i$，$c_i$。
图6.8显示了一个家庭的决定空间，取决于其砷含量和到最近安全井的距离，给定 $a_i$，$b_i$，$c_i$该模型下的决定是确定性的，但是 $a_i$，$b_i$，$c_i$不能直接观察到-我们所看到的只是针对家庭的决策（$y_i=0$或$1$），因为他们的砷水平为$(As{)}_i$，到最近的安全井的距离 $x_i$。
例如，如果 $a_i$ 和 $c_i$ 是常数，而 $b_i/a_i$ 具有独立于$(As{)}_i$和 $x_i$ 的逻辑分布，则总体中（a，b，c）的某些分布会减少到拟合的logistic回归。 如果因素相加，系数在总体中不变，并且存在固定成本（在此示例中为 $b_i$），且在总体中具有逻辑分布，则选择模型可简化为逻辑回归。
（a，b，c）的其他分布也是可能的。可以拟合相应的模型，并将这些效用参数视为潜在数据。在R中没有使用glm（）拟合此类模型的简便方法（除非特殊情况是减少到logit和probit），但它们可以适合Bug（请参阅练习17.7）。

决策模型的见解

选择模型即使没有正式拟合也可以提供一些见解，例如，在拟合逻辑回归中，我们发现距离可以很好地用作线性预测变量，砷水平在对数尺度上更合适。暗示这两个因素应该线性地出现，砷的转化表明人们（错误地）以对数尺度感知风险-例如，看到4到8之间的差异不比1和2之间的差异更糟。 2.（此外，我们的残差图表明，人们似乎低估了砷水平非常接近0.5的风险。在此之后的简化假设是，砷水平低于0.5的所有井都是“安全的”。）
我们也可以使用本实用新型模型来解释模型中的教育系数，即受过良好教育的人更容易转换，表明他们的转换成本较低，或者他们从减少砷暴露中获得的感知收益较高。人口中潜在的效用参数。
还可以对模型进行详细考虑，以考虑各种选择方案，包括不进行任何操作，改用现有的私人井，改用现有的社区井或挖掘新的私人井，具体取决于成本。步行，对健康风险的感知，财务资源和未来计划。

PartB 处理回归推断

现在，我们讨论如何超越简单地查看回归系数，首先通过使用仿真总结和传播推断不确定性，然后考虑如何将回归用于因果推断。

三、概率模型和统计推断的仿真

每当我们使用点估计和标准误差来表示参数的推论时，我们都会执行数据约简，如果估计值是正态分布的，则此摘要不会丢弃任何信息，因为正态分布完全由均值和方差定义。 但是在其他情况下，通过一组随机模拟来表示参数估计中的不确定性可能是有用的，这些随机模拟表示参数向量的可能值（在模拟中更可能出现值）。 那么，我们的意思是通过随机数而不是通过点估计和标准误差来概括推理。

7.1 概率模型的仿真

在本节中，我们将介绍两个简单概率模型的仿真，本章的其余部分讨论如何使用仿真来总结和理解回归并推广线性模型，下一章将仿真应用于模型检查和验证。 还准备了多级模型，我们使用基于仿真的推理进行拟合，如第2B部分所述。

离散预测仿真的简单示例

400个婴儿中有多少个女孩？婴儿是男孩还是女孩的概率分别为48.8％或51.2％。假设在某年在医院出生的婴儿为400个，那么女孩将是多少？
我们可以使用二项分布来模拟400次出生：

n.girls <- rbinom (1, 400, .488) 
print (n.girls)

为了了解可能发生的情况的分布，我们对过程进行了1000次模拟（首先创建矢量n.girls以存储模拟结果之后）：这向我们展示了400例婴儿中可能发生的情况。

n.sims <- 1000
n.girls <- rep (NA, n.sims)  '建立数组’'
for (s in 1:n.sims){    '为数组赋值，每个值服从预设参数的二项分布’'
n.girls[s] <- rbinom (1, 400, .488)} 
hist (n.girls)

这样就得到了图7.1所示的直方图，该直方图代表了女婴出生数量的概率分布，经过1000次模拟可以捕获不确定性。

从模拟结果看到女孩出生数量的分布主要集中在190-210。

如果我们要一个医院中女孩的个数，需要考虑会影响女孩个数的几种情况。这一胎是否是双胞胎，是异卵双胞胎还是同卵双胞胎。

birth.type <- sample (c("fraternal twin","identical twin","single birth"), size=400, replace=TRUE, prob=c(1/125, 1/300, 1 - 1/125 - 1/300))
girls <- rep (NA, 400) 
for (i in 1:400){
  if (birth.type[i]=="single birth"){ girls[i] <- rbinom (1, 1, .488)}
  else if (birth.type[i]=="identical twin"){ girls[i] <- 2*rbinom (1, 1, .495)}
  else if (birth.type[i]=="fraternal twin"){ girls[i] <- rbinom (1, 2, .495)}
}
n.girls <- sum (girls)

将这样的模拟进行1000次

n.sims <- 1000
n.girls <- rep (NA, n.sims) 
for (s in 1:n.sims)
{
  birth.type <-sample(c("fraternal twin","identical twin","single birth"), 
                      size=400, replace=TRUE, prob=c(1/125, 1/300, 1 - 1/125 - 1/300))
  girls <- rep (NA, 400) 
  for (i in 1:400){
    if (birth.type[i]=="single birth"){ girls[i] <- rbinom (1, 1, .488)}
    else if (birth.type[i]=="identical twin"){ girls[i] <- 2*rbinom (1, 1, .495)}
    else if (birth.type[i]=="fraternal twin"){ girls[i] <- rbinom (1, 2, .495)}
  }
  n.girls[s] <- sum (girls)
}
hist(n.girls)

形状和上面的均匀分布大同小异。

连续随机变量的模拟

n.sims <- 1000
avg.height <- rep (NA, n.sims) 
for (s in 1:n.sims){
  sex <- rbinom (10, 1, .52)
  height <- ifelse (sex==0, rnorm (10, 69.1, 2.9), rnorm (10, 64.5, 2.7)) 
  avg.height[s] <- mean (height)
}
hist (avg.height, main="Average height of 10 adults")

那10个人的最大身高呢？要确定这个高度，只需在循环中添加以下行：

max.height[s] <- max (height)

当然在此之前要初始化max.height数组

max.height <- rep (NA, n.sims)

总的来一遍：

sex <- rbinom (10, 1, .52)
height <- ifelse (sex==0, rnorm (10, 69.1, 2.9), rnorm (10, 64.5, 2.7)) 
avg.height <- mean (height)
print (avg.height)

n.sims <- 1000
avg.height <- rep (NA, n.sims) 
max.height <- rep (NA, n.sims)

for (s in 1:n.sims){
  sex <- rbinom (10, 1, .52)
  height <- ifelse (sex==0, rnorm (10, 69.1, 2.9), rnorm (10, 64.5, 2.7)) 
  avg.height[s] <- mean (height)
  max.height[s] <- max(height)
}
hist (avg.height, main="Average height of 10 adults")
hist(max.height, main="Max height of 10 adults")

使用定制功能在R中进行仿真

如果我们在R中将单个模拟的步骤表示为函数，则模拟的编码会变得更加清晰。我们以平均高度的模拟为例进行说明。

Height.sim <- function (n.adults){
  sex <- rbinom (n.adults, 1, .52)
  height <- ifelse (sex==0, rnorm (10, 69.5, 2.9), rnorm (10, 64.5, 2.7)) 
  return (mean(height))
}

（为简单起见，我们对女性的比例以及男性和女性身高的均值和标准差进行了“硬编码”，但更笼统地说，这些可以作为函数的参数提供。）
然后，我们使用replicate()函数调用Height.sim()1000次：

有关在伪数据模拟中执行功率计算的R中函数使用的更详细示例，请参见20.5节。

7.2 使用模拟总结线性回归：朴素贝叶斯方法

在回归设置中，我们可以使用模拟来捕获预测不确定性（回归模型中的误差项）和推论不确定性（系数的标准误差和有关残差的不确定性），我们首先讨论模拟预测的最简单情况错误，然后考虑推论不确定性以及两种变异来源的组合。

仿真代表预测不确定性

我们使用第63页的模型（4.2）来说明预测不确定性以及预测身高68英寸的男人的收入的问题。
自动获取点和区间预测，使用R中的回归软件可以轻松访问预测的估计和置信区间：

x.new <- data.frame (height=68, male=1)
pred.interval <- predict (earn.logmodel.3, x.new, interval="prediction",
level=.95)

在R中，我们可以为新数据创建一个data frame，然后使用predict()函数。其中earn.logmodel.3为前面建立的模型：

earn.logmodel.3 <- lm (log.earn ~ height + male + height:male)
               coef.est  coef.se
(Intercept)       8.388  0.844
height            0.017  0.013 
male             -0.079  1.258 
height:male       0.007  0.019
n = 1192, k = 4 residual sd = 0.88, R-Squared = 0.09

然后求幂以得出原始（未记录）收入规模的预测：

exp (pred.interval)

使用模拟构建预测间隔现在，我们讨论如何使用从拟合回归模型得到的模拟来“手动”获得预测间隔。 但是，仿真是一种通用工具，我们将能够将其应用于更复杂的预测设置中，如本章后面和下一章所述。

对数收益的点估计为$8.4+0.017\cdot 68-0.079\cdot 1+0.007\cdot 68\cdot 1=9.95$，标准差为0.88。将这些值放在原始（未记录）标度上，我们求幂以得出几何平均值的 $e^{9.95}=21000$ ，几何标准偏差 $e^{0.88}=2.4$。
然后，68％的预测间隔为 $[21000/2.4，21000\cdot 2.4]=[8800，50000]$，而95％的预测间隔为 $[21000/2.42，21000\cdot 2.42]=[3600，121000]$。

模拟预测是一组随机数，其对数的平均值为9.95，标准偏差为0.88，例如，在R中，我们可以使用以下命令来汇总预测分布：

pred <- exp (rnorm (1000, 9.95, .88))

它告诉R从正态分布中取1000个随机数，均值9.95，方差0.88，然后对这些值取幂。
我们可以将模拟显示为直方图（参见图7.2），还可以计算各种数值汇总，例如，
•均值：mean(pred)
•中位数: quantile(pred.5)
•50％间隔： quantile(pred，c(.25，.75))
•95％间隔： quantile(pred，c(025，.975))
（这些计算忽略了回归参数中的不确定性，因此只是近似的；我们将在本节后面介绍更完整的计算过程。）

mean(pred)
quantile(pred,.5)
quantile(pred,c(.25,.75))
quantile(pred,c(.025,.975))

为什么我们需要模拟来进行预测性推论？

在许多情况下，点估计，标准误差以及从R中的predict()函数获得的间隔就足够了，因为中心极限定理可确保对于除最小样本量外的所有样本以及合理地表现良好的误差分布，系数估计为考虑到标准误差估计中的不确定性，具有 $n-k$ 个自由度的 $t$ 分布（其中 $k$ 是模型中的预测变量数）是适当估计的可靠近似解析。程序也可用于获得参数和预测的线性组合的不确定性（预测的线性组合的一个示例是使用第3章中的一种模型来预测某项的平均测试得分）。一组100名孩子的母亲的教育程度和智商很高。

但是，对于更笼统的预测，计算不确定性的最简单，最可靠的方法是通过仿真，例如，假设我们有一个身高68英寸的女人和一个身高68英寸的男人，我们想使用模型（4.2）来预测他们的收益差异作为点估计，我们可以使用点预测的差异：
$e^{8.4+0.017\cdot 68-0.079\cdot 1+0.007\cdot 68\cdot 1}-e^{8.4+0.017\cdot 68-0.079\cdot 0+0.007\cdot 68\cdot 0}=6900$。
为此预测获取标准误差或不确定性区间的最简单方法是使用仿真：

pred.man <- exp (rnorm (1000, 8.4 + .017*68 - .079*1 + .007*68*1, .88)) 
pred.woman <- exp (rnorm (1000, 8.4 + .017*68 - .079*0 + .007*68*0, .88)) 
pred.diff <- pred.man - pred.woman
pred.ratio <- pred.man/pred.woman

模拟表示回归系数的不确定性

拟合回归的通常总结给出了每个系数的标准误以及估计值，这些都给人一种估计上的不确定性的感觉（请参阅第40页的图3.7）。但是，当超出单个系数的推论时， 通过仿真总结推论是有帮助的，这使我们在传播有关参数和预测组合的不确定性方面具有完全的灵活性。
对于经典线性回归和广义线性模型，我们使用R中的sim()函数实现这些模拟。例如，如果执行以下代码

n.sims <- 1000
fit.1 <- lm (log.earn ~ height + male + height:male) 
sim.1 <- sim (fit.1, n.sims)

那么sim.1$beta是一个具有1000行4列的矩阵（表示向量的1000个独立模拟 $(\beta 0,\beta 1,\beta 2,\beta 3))$，而 sim.1$sigma 是一个长度为1000的向量（表示估计） 残留标准偏差参数 $\sigma$ 中的不确定性。
我们可以检查这些模拟是否等同于回归计算，例如通过R中的以下命令，这些命令将在拟合模型中打印出系数的平均值，标准偏差和95％的区间：

height.coef <- sim.1$beta[,2]
mean (height.coef)
sd (height.coef)
quantile (height.coef, c(.025,.975))

举一个更有趣的示例，请考虑以下问题：在这个交互模型中，关于男人之间的身高系数可以说什么？我们不能使用回归输出直接回答这个问题：男人的斜率是身高和身高的总和 ：男性系数，根据回归表中的信息，没有简单的方法来计算其标准误差，最直接的方法是直接从推论模拟中计算出95％的区间：

height.for.men.coef <- sim.1$beta[,2] + sim.1$beta[,4] 
quantile (height.for.men.coef, c(.025,.975))

结果为 $[-0.003，0.049]$ ，即 $[-0.3％，4.9％]$ ，统计意义不是分析的目标-如果例如以95％的间隔为单位，则我们的结论不应有太大变化 $[0.1％，5.3％]$ 但是要知道估计的不确定性很重要，并且使用推论模拟可以方便地做到这一点，将推论模拟的能力与预测结合起来可以更有效地展示。

模拟程序的细节

为了获得n.sims仿真绘图（例如，通常1000个绰绰有余；请参阅第17章），我们基于贝叶斯推论应用了以下过程（请参阅第18章）。

1.使用 $k$ 个预测变量上的 $n$ 个数据点的经典回归来计算估计参数的矢量 $\beta$,、未缩放的估计协方差矩阵 $V\beta$ 和残差 ${\hat{\sigma }}^2$。

2.创建系数向量 $\beta$ 和残差标准偏差 $\sigma$ 的 n.sims随机模拟。
（a）模拟$\sigma =\hat{\sigma }\sqrt{(n-k)/X}$，其中X是来自具有 $n-k$ 自由度的 ${\chi }^2$ 分布的随机抽取。
（b）给定 $\sigma$ 的随机抽取，从均值 $\hat{\beta }$ 和方差矩阵 ${\sigma }^{2}V\beta$ 的多元正态分布模拟 $\beta$。
这些模拟以估计值 $\hat{\beta }$ 和 $\hat{\sigma }$ 为中心，并在参数中表示估计的不确定性（例如，${\beta }_1$的大约68％的模拟将在 ${\hat{\beta }}_1$ 的±1标准误差范围内，大约95％将在误差在±2个标准误差之内，依此类推。）

这些步骤由我们的R函数 sim()自动执行，该函数从拟合的线性模型中提取 $n，k，\beta ,V\beta ，\sigma$，然后对n.sims模拟执行循环：

for (s in 1:n.sims){
  sigma[s] <- sigma.hat*sqrt((n-k)/rchisq(1,n-k)) 
  beta[s,] <- mvrnorm (1, beta.hat, V.beta*sigma[s]^2)
}

sim()函数最后返回 $\sigma$ 模拟的向量和 $\beta$ 模拟的 $n_{sims}\times k$ 矩阵:

return (list (beta=beta, sigma=sigma))

为列表项指定了名称，以便可以从仿真对象中使用这些名称访问这些功能，该功能对广义线性模型（如logistic和Poisson回归）的工作原理类似，可以通过使用系数估计的标准误差来调整任何过度分散。如果模型中包含超分散，则对其进行缩放。

朴素贝叶斯推理

贝叶斯推理是指将未知参数（以及缺失和潜在数据）建模为随机变量的统计程序，如第18.3节中所述，贝叶斯推理从对未知参数的先验分布开始，并用概率的先验分布对其进行更新。数据，产生后验分布，用于推理和预测。
本书的第2部分讨论了贝叶斯推理如何适用于多级建模，在这种模型中，很自然地将概率分布拟合到一批参数中。对于第1部分中考虑的经典模型，贝叶斯推理较为简单，通常以“ “非信息性”或未知参数上的均匀先验分布。除了注意到此处提供的模拟与这些非信息性先验分布相对应之外，我们在此不再进一步探讨技术问题。将这些模拟视为代表配置也可能会有帮助。与观察到的数据兼容的参数和预测的定义-在同一个意义上，经典的置信区间包含与数据不矛盾的一系列参数值。