设随机变量 \(Y\) 的概率密度函数或概率质量函数 \(f(y)\) 可以表达为 f(y; \theta) = \exp\left\{ a(y) b(\theta) + c(\theta) + d(y) \right\}, 且 \(f(\cdot;\theta)\) 的支撑集不依赖于 \(\theta\) ， 则称 \(Y\) 服从 单参数指数族分布 。

单参数指数族分布与广义线性模型有密切联系， 其中 \(b(\cdot)\) 称为经典联系函数， 作为广义线性模型中 \(g(\cdot)\) 的一种重要选择。

对单参数指数族分布随机变量 \(Y\) ， = -\frac{c'(\theta)}{b'(\theta)}, \quad \text{Var}(Y) = \frac{b''(\theta) c'(\theta) - c''(\theta) b'(\theta)}{[b'(\theta)]^3} .

单参数指数族分布包含了许多常见分布， 如方差已知的正态分布、泊松分布、二项分布等。

泊松分布的概率质量函数为 f(y; \lambda) = \frac{\lambda^y}{y!} e^{-\lambda} = \exp\{y \ln \lambda - \lambda - \ln(y!) \}, \ y=0,1,\dots, 是单参数指数族分布。 参数为 \(\lambda = EY\) 。 经典联系函数为 \(g(\theta) = \ln\theta\) ， 一元泊松回归模型可以写成 \[\begin{aligned} Y \sim & \text{Poisson}(\theta), \\ g(\theta) =& \ln(\theta) = a + b x . \end{aligned}\]

二项分布B( \(n,p\) )当 \(n\) 已知， 以 \(p\) 为参数时，是单参数指数族分布， 其概率质量函数为 \[\begin{aligned} f(y; p) =& \binom{n}{y} p^y (1-p)^{n-y}, \ y=0,1,\dots,n \\ =& \binom{n}{y} (1-p)^n \left( \frac{p}{1-p} \right)^y \\ =& \exp\left\{ y \ln \frac{p}{1-p} + n \ln(1-p) + \ln \binom{n}{y} \right\} \\ =& \exp\{ a(y) b(\theta) + c(\theta) + d(y) \}, \end{aligned}\] 这里 \(\theta = p = \frac{1}{n} EY\) , 其中 \(n\) 已知。 \(b(p) = \ln \frac{p}{1-p}\) ， 称 \(\ln \frac{p}{1-p}\) 为 \(\text{logit}(p)\) ， 是二项分布的经典联系函数。 \[\begin{aligned} \text{logit}(p) =& \ln \frac{p}{1-p}, \ p \in (0, 1), \\ \quad \text{logit}^{-1}(u) =& \frac{e^u}{1 + e^u} = \frac{1}{1 + e^{-u}} \in (0, 1), \ u \in (-\infty, \infty) . \end{aligned}\] 一元逻辑斯谛回归模型可以写成 \[\begin{aligned} Y \sim & \text{B}(n, \theta), \\ g(\theta) =& \text{logit}(\theta) = a + b x . \end{aligned}\]

对于以上的三个分布， 使用经典联系函数将因变量分布参数与自变量线性组合联系起来可以给出一些粗浅的解释。

设随机变量 \(Y\) 的分布受到自变量 \(x_1, \dots, x_p\) 的影响， 在给定 \(x_1, \dots, x_p\) 的条件下 \(Y\) 服从泊松分布， 则可以用泊松回归描述 \(Y\) 与 \(x_1, \dots, x_p\) 的关系： \[\begin{aligned} Y \sim& \text{Poisson}(\lambda), \\ g(\lambda) =& \beta_0 + \beta_1 x_1 + \dots + \beta_p x_p . \end{aligned}\] 其中 \(g(\lambda)\) 经常取为经典联系函数 \(\ln\lambda\) 。

泊松回归的假定：

考虑单个自变量的泊松回归， 设 \(Y_i \sim \text{Poisson}(\lambda_i)\) , \(\ln \lambda_i = \beta_0 + \beta_1 x_i\) ， 于是单个观测 \((x_i, y_i)\) 对于似然函数的贡献为 \frac{\lambda_i^{y_i}}{y_i!} e^{-\lambda_i}, 对于对数似然函数的贡献为 y_i \ln\lambda_i - \ln y_i! - \lambda_i = y_i (\beta_0 + \beta_1 x_i) - \ln y_i! - e^{\beta_0 + \beta_1 x_i}, 其中 \(-\ln y_i!\) 在对数似然函数中是常数项， 所以对数似然函数为 \ell (\beta_0, \beta_1) = \text{常数} + \sum_{i=1}^n \left\{ y_i (\beta_0 + \beta_1 x_i) - e^{\beta_0 + \beta_1 x_i} \right\} .

可以用Newton-Raphson迭代方法进行最大似然估计的数值求解， 并且可以使用Fisher得分方法改造Newton-Raphson方法， 具体方法是在迭代中用负信息阵的当前估计值代替对数似然函数的海色阵。 这样还可以从最大值点处的信息阵的逆矩阵给出参数估计的协方差阵估计。 这种利用Fisher得分的Newton-Raphson迭代方法， 等价于迭代再加权最小二乘方法， 每步迭代是一个加权最小二乘计算， 加权按照当前的方差倒数来计算。

线性回归模型使用残差平方和以及方差分析评估模型的拟合优劣。 因为广义线性模型是对因变量分布参数的拟合而非直接对因变量值的拟合， 所以没有原来的残差以及残差平方和的概念， 而是用偏差(deviance)代替。

设要评估的模型的对数似然函数在最大似然估计处的值为 \(\ell_M\) ， 为了比较， 考虑一个参数最复杂的“饱和模型”， 使得因变量得到最大程度的解释。 在方差分析中对因素的每种组合都单独使用一个参数。 在饱和模型中， 应该每一组不同的自变量组合都单独设置一个因变量独立参数值， 但是在R统计软件 glm() 计算时实际是对每个观测单独使用参数。 计算饱和模型的对数似然函数的最大值 \(\ell_S\) 。 饱和模型的自由参数（独立参数）个数是最多的。

要评估的模型的偏差(deviance)定义为 \text{Dev} = -2[\ell_M - \ell_S], 偏差一定是非负的， 对于某些广义线性模型， 如果模型M很好地解释了数据， 则模型M与饱和模型S应该效果相近， 在一定条件下偏差近似服从卡方分布， 自由度等于两个模型独立参数的个数差， 当偏差不超过卡方的 \(1-\alpha\) 分位数时可以认为模型M是充分的。

类似于线性模型中的平方和分解， 偏差也可以进行分解。 有两个模型 \(M_0\) 和 \(M_1\) ， \(M_0\) 是 \(M_1\) 的特例， 即 \(M_0\) 是 \(M_1\) 的参数取特定值或满足特定约束后的模型， 称 \(M_1\) 为完全模型， 称 \(M_0\) 为简化模型， 称这两个模型有嵌套关系。 为了在 \(M_1\) 的假设下检验 \(M_0\) 是否成立， 计算似然比统计量 -2[\ell_0 - \ell_1] = -2[\ell_0 - \ell_S] + 2[\ell_1 - \ell_S] = \text{Dev}_0 - \text{Dev}_1, 当 \(M_0\) 与 \(M_1\) 效果相同的零假设成立时统计量服从卡方分布， 自由度等于两个模型的独立参数个数差， 这个统计量可以用两个模型的偏差之差计算。

对于泊松回归， 也可以定义残差， 但是不像线性回归那样只有观测值减去预测值这一种， 因为广义线性模型中自变量预测的是因变量的期望值而非因变量本身。

以一元泊松回归为例。 响应残差(response residual)定义为 \(e_i = y_i - \hat\mu_i\) ， 而 \(\hat\mu_i\) 估计为 \(g^{-1}(\hat\beta_0 + \hat\beta_1 x_i)\) ， 其中 \(g(\cdot)\) 是联系函数， 对泊松回归一般取为 \(\ln\) ， 所以 \(\hat\mu_i = \exp(\hat\beta_0 + \hat\beta_1 x_i)\) 。

定义标准化残差或Pearson残差为 e_i^{(P)} = \frac{y_i - \hat\mu_i}{\sqrt{\hat\mu_i}} .

定义偏差残差(deviance residual)为 e_i^{(D)} = \text{sgn}(y_i - \hat\mu_i) \left\{ 2 \left[ y_i (\ln y_i - \ln\hat\mu_i) - (y_i - \hat\mu_i) \right] \right\}^{1/2} . 其平方和恰好等于模型的残差偏差(residual deviance)。 偏差残差一般与Pearson残差近似相等。

对 glm() 结果， 可以用 residuals() 函数提取各种残差， 指定 type="response" 或 "pearson" 、 "deviance" ， 缺省为 type="deviance" 。

基本R软件的配套（默认安装）扩展包stats提供了 glm() 函数， 用来进行广义线性模型建模。 对家庭人口数问题， 以 total 为因变量， 以户主年龄 age 为自变量， 先建立仅有一个自变量的模型：

可以建立一个没有自变量的模型作为基准模型， 称这个模型为空模型(null model)：

从探索性数据分析结果发现年龄对平均人口数的影响并不是单调的， 而是在中间达到最大值。 所以有可能需要增加年龄的二次项。 平均年龄约为53， 为了避免共线问题，增加 \((\text{age}-53)^2\) 项：

其它自变量对家庭人口数也可以有影响。 先增加location变量， 这是一个有5个水平的因子， 建模时会被转化为4个二值哑变量：

流行病学研究中， 经常要研究某种病的发病率如何受到自变量的影响。 虽然这属于二项分布， 但是因为发病率都比较低， 可以用泊松分布作为近似模型。 设 \(T\) 为总人数和累计时间（如人年单位）， \(\lambda\) 为 \(T\) 对应的平均发病人数， 则针对参数 \rho = \lambda / T 这样才可以对人口数不同的地区使用统一的模型。 利用泊松回归， \log(\lambda) = \beta_0 + \beta_1 x + \log(T) , glm() 公式中可以用 offset() 表示系数固定等于1的自变量。

例36.1 ISwR包的eba1977包含类丹麦的4个城市肺癌发病率数据， 以年龄和城市为自变量建模。

在泊松回归这样的因变量为计数变量问题中， 可能会出现如下的问题：

逻辑斯谛回归模型有如下的模型假定：

数据也可以是 \((x_i, y_i, z_i), i=1,2,\dots,n\) ， 其中 \(y_i\) 代表第 \(i\) 个观测的成功数， \(z_i\) 代表失败数， \[\begin{aligned} y_i \sim& \text{B}(m_i, p_i), \\ \text{logit}(p_i) =& a + b x_i, \ i=1,2,\dots, n . \end{aligned}\] 其中 \(m_i = y_i + z_i\) 是第 \(i\) 个观测的试验次数， 是非随机的正整数值。

先考虑观测的因变量是二值因变量情形。 数据为 \((x_i, y_i)\) , \(i=1,\dots, n\) ， 一个观测对似然函数的贡献为 p_i^{y_i} (1 - p_i)^{1 - y_i}, p_i = \text{logit}^{-1}(a + b x_i) = \frac{e^{a + b x_i}}{1 + e^{a + b x_i}}, 对数似然函数为 \[\begin{aligned} \ell(a, b) =& \sum_{i=1}^n \left\{ y_i \ln p_i + (1-y_i) \ln(1-p_i) \right\} \\ =& \sum_{i=1}^n \left\{ y_i \ln \frac{p_i}{1-p_i} + \ln(1-p_i) \right\} \\ =& \sum_{i=1}^n \left\{ y_i (a + b x_i) - \ln(1 + e^{a + b x_i}) \right\} . \end{aligned}\]

若观测是二项数据 \((x_i, m_i, y_i)\) ， 其中 \(y_i \sim \text{B}(m_i, p_i)\) ， p_i = \text{logit}^{-1}(a + b x_i), 则第 \(i\) 个观测对似然函数的贡献为 \binom{m_i}{y_i} p_i^{y_i} (1 - p_i)^{m_i - y_i}, 对于对数似然函数的贡献为 \[\begin{aligned} & \ln \binom{m_i}{y_i} + y_i \ln p_i + (m_i - y_i) \ln(1 - p_i) \\ =& \ln \binom{m_i}{y_i} + y_i \ln \frac{p_i}{1 - p_i} + m_i \ln(1 - p_i) \\ =& \ln \binom{m_i}{y_i} + y_i \ln(a + b x_i) - m_i \ln(1 + e^{a + b x_i}) , \end{aligned}\] 忽略与待估参数 \(a, b\) 无关的常数项后， 对数似然函数为 因变量的对数似然函数为 \[\begin{align} \ell(a, b) =& \sum_{i=1}^n \left[ y_i \ln p_i + (m_i - y_i) \ln(1 - p_i) \right] \tag{36.1} \\ =& \sum_{i=1}^n \left[ y_i \ln(a + b x_i) - m_i \ln(1 + e^{a + b x_i}) \right] . \end{align}\]

需要用数值迭代算法求最大似然估计。 若样本为满足模型的随机样本， 最大似然估计是相合估计， 渐近有效估计， 具有渐近正态分布。

得到最大似然估计 \((\hat a, \hat b)\) 后， 一般用信息阵的逆矩阵估计其协方差阵， 信息阵为： I(a,b) = - E \left[ \frac{\partial^2 \ell(a,b)}{\partial (a, b) \partial (a, b)^T} \right] , 在用数值迭代算法计算最大似然估计时， 一般会得到最大值点处的对数似然函数的海色阵， 加上负号并求逆矩阵， 就可以作为参数协方差阵估计。

有了参数估计的协方差阵估计， 就得到了参数估计的标准误差。 设参数估计为 \(\hat\theta\) ， 估计的方差开平方根作为其抽样分布标准差估计， 称为标准误差，记为 \(\text{SE}(\hat\theta)\) 。 Z = \frac{\hat\theta}{\text{SE}(\hat\theta)}, 可以用 \(Z\) 统计量检验 \(H_0: \theta=0\) ， 在 \(H_0\) 成立、大样本且满足适当正则性条件时 \(Z\) 近似服从标准正态分布。 也可以用 \(Z^2\) 近似服从 \(\chi^2(1)\) 分布进行检验。 这样对回归系数进行的检验称为Wald类型的检验。 此检验对中小样本有一定缺陷， 当参数的实际数值绝对值较大时， 容易错判成不显著。 另一种可用的方法是下面所述的利用偏差统计量的方法。

当因变量是多次试验的二项分布结果时， 虽然可以计算比例与模型预测概率的差， 但残差平方和与线性正态情形的回归模型的残差平方和不同， 所以用偏差(deviance)代替了残差平方和。

对 \(y_i \sim \text{B}(m_i, p_i)\) ， 设 \(\hat p_i = \text{logit}^{-1}(\hat a + \hat b x_i)\) , \(\hat y_i = m_i \hat p_i\) ， 定义逻辑斯谛回归的偏差为 = 2 \sum_{i=1}^n \left[ y_i \log \left( \frac{y_i}{\hat y_i} \right) + (m_i - y_i) \log \left( \frac{m_i - y_i}{m_i - \hat y_i} \right) \right] .

当因变量观测值 \(y_i\) 与预测值 \(\hat y_i\) 都很接近时， 偏差 \(G^2\) 的值比较小； 当模型不准确时，偏差的值比较大。 记此模型为 \(M\) ， 设模型 \(M\) 有 \(p\) 个回归系数， 模型 \(M\) 拟合充分时一定条件下偏差服从自由度为 \(n-p\) 的卡方分布。 下面解释上述偏差公式。

在R软件的 glm() 函数计算时， “饱和模型”（称为模型S）就是对每个观测估计一个不同的 \(\hat y_i\) ， 即 \(\hat p_i = y_i / m_i\) , \(\hat y_i = y_i\) 。这样有 \(n\) 个参数。 但是要注意， 饱和模型使用 \(n\) 个参数， 默认假定了各个观测的自变量组合都是互不相同的， 因为某一个特定的自变量组合如果在多个观测中重复出现应该给出相同的预测值， R中并没有自动进行各个观测的自变量组合互不相同的验证。

考虑足球比赛罚球得分的概率， 研究得分概率是否与罚球球队当前是否比分领先有关系。 若干次罚球的汇总表格数据如下：

定义自变量 X = \begin{cases} 1, & \text{罚球球队比分领先}, \\ 0, & \text{比分不领先}, \end{cases} 则问题用逻辑回归表示， \[\begin{aligned} Y \sim& \text{B}(m, p), \\ \ln\frac{p}{1-p} =& \beta_0 + \beta_1 X . \end{aligned}\] 问题的数据比较简单， 可以看成只有 \(X=0\) 和 \(X=1\) 两个观测， 于是当 \(X=0\) (罚球球队比分不领先时)， 记这时罚球成功概率为 \(p_0\) ， \ln\frac{p_0}{1-p_0} = \beta_0; 当 \(X=1\) ，即罚球球队比分领先时， 记这时罚球成功概率为 \(p_1\) ，则 \ln\frac{p_1}{1-p_1} = \beta_0 + \beta_1; 由此可以解释 \(\beta_1\) 的含义为： \[\begin{aligned} \beta_1 =& (\beta_0 + \beta_1) - \beta_0 = \ln\frac{p_1}{1-p_1} - \ln\frac{p_0}{1-p_0} \\ =& \ln \frac{\frac{p_1}{1-p_1}}{\frac{p_0}{1-p_0}}, \end{aligned}\] 其中 \(\frac{\frac{p_1}{1-p_1}}{\frac{p_0}{1-p_0}}\) 是比分领先相对于比分不领先时的优势比(odds ratio)， 所以 \(e^{\beta_1}\) 等于优势比。

用二项分布的分布概率可以写出模型的似然函数为 \[\begin{aligned} L(\beta_0, \beta_1) =& \binom{24}{22} p_1^{22} (1 - p_1)^{2} \ \binom{180}{141} p_0^{141} (1 - p_1)^{39} \\ \propto & \left( \frac{e^{\beta_0 + \beta_1}}{1 + e^{\beta_0 + \beta_1}}\right)^{22} \left(1 - \frac{e^{\beta_0 + \beta_1}}{1 + e^{\beta_0 + \beta_1}}\right)^{2} \left(\frac{e^{\beta_0}}{1 + e^{\beta_0}}\right)^{141} \left(1 - \frac{e^{\beta_0}}{1 + e^{\beta_0}}\right)^{39} . \end{aligned}\] 可以用数值优化方法求得最大似然估计 \(\hat\beta_0 = 1.2852\) , \(\hat\beta_1 = 1.1127\) 。 从而优势比的估计为 \(e^{1.1127} = 3.04\) ， 与直接计算优势比完全相同。

使用R的 glm() 函数来计算这个例子。 将因变量输入为由成功数与失败数两列组成的矩阵。

数据框Michelin中包含了Zagat Survey 2006考察的164家法餐馆的顾客评分和是否被米其林推荐收录。 数据中每一行是一个评分值对应的数据， 变量 Food 是Zagat评分，取值从15到28， InMichelin 是某一评分被米其林收录家数， NotInMichelin 是未被米其林收录家数， mi 是该评分的总家数， proportion 是收录比例。

以是否收录作为因变量， 以Zagat评分作为自变量，建立逻辑斯谛回归：

考虑在某地的4个不同区域进行关于是否修建铁路的意见征询投票的数据分析， 赞同修建铁路的比例可能受到所在区域到要修建的铁路的距离(单位：英里)的影响， 距离越远越容易反对修建， 可能受到投票人中黑人比例的影响。 数据如下：

文件“data/remiss.csv”中保存了癌症病人康复的数据， 变量remiss为康复与否(1为康复，0为未康复)， 另外的6个变量为可能影响康复概率的自变量。 研究这些自变量对康复概率的影响， 并筛选最有解释能力的自变量。

读入数据:

逻辑斯谛回归是因变量为二值变量时的广义线性模型。 在有些应用问题中， 因变量是有序变量， 比如满意程度取“低”、“中”、“高”， 这时不应该分别对“中对低”、“高对低”用二值的逻辑斯谛回归建模， 因为“高”和“中”之间也是可比的， 分别建模的结果不一定能保证这个次序。

采用如下的理论模型。 设某个潜在的连续取值的因变量为 \(Z\) ， 服从如下的逻辑斯谛分布： F(x; \eta) = \text{logit}^{-1}(x-\eta), 易见 \(E Z = \eta\) 。 设观测到的取有序因子值的因变量 \(Y\) 是将 \(Z\) 的值按照如下的分点 \zeta_0 = -\infty < \zeta_1 < \dots < \zeta_{K-1} < \zeta_K = +\infty 将实数轴分为 \(K\) 段，各段分别编码为 \(1,2,\dots, K\) ， \[\begin{aligned} P(Y \leq k) =& P(Z \leq \zeta_k) = \text{logit}^{-1}(\zeta_k - \eta), \\ \text{logit}[P(Y \leq k)] =& \zeta_k - \eta . \end{aligned}\]

设潜在变量 \(Z\) 的均值 \(\eta\) 依赖于自变量的线性组合， 则有如下的有序因子值的因变量的“比例发生比逻辑斯谛回归”(POLR)模型： \[\begin{align} \text{logit}[P(Y \leq k \;|\; x)] = \zeta_k - \beta x, k=2,3,\dots, K . \tag{36.2} \end{align}\] 这里 \(k=1\) 作为一个基线值。 因为 \(\zeta_k\) 未知，所以回归的截距项也归入到 \(\{ \zeta_k \}\) 中。 模型左边是因变量类别为 \(k\) 以及低于 \(k\) 的概率的对数发生比， 如果自变量的值不变， 对于不同的 \(k\) 这些对数发生比只差一个常数， 发生比只差一个倍数， 所以称为“比例发生比逻辑斯谛回归”模型； 因为是对 \(P(Y \leq k)\) 从 \(k=1,2,\dots,K\) 依次建模， 所以也称为累积逻辑斯谛回归模型。