广义线性模型

Y1,Y2,...,YNY_1,Y_2,...,Y_NY1,Y2,...,YN是服从指数分布族某一分布的被解释变量，并且EYi=μiEY_i=\mu_iEYi=μi，存在某个函数ggg使得解释变量与g(μi)g(\mu_i)g(μi)之间具有线性关系
g(μi)=Xiβg(\mu_i) = X_i \beta g(μi)=Xiβ
这样的回归模型叫广义线性回归模型。显然当g(μi)=μig(\mu_i)=\mu_ig(μi)=μi时，回归模型是多元线性回归，当ggg是Logistics函数的反函数时，是Logistics回归。

二值被解释变量

回归模型
Yi=Xiβ+ϵiY_i = X_i \beta + \epsilon_i Yi=Xiβ+ϵi
中，有时Yi=0,1Y_i = 0,1Yi=0,1，这类解释变量叫二值被解释变量，这种回归可以用来做两分类问题。如果把YiY_iYi视为Bernoulli随机变量，则
pi=P(Yi=1)=E[Yi]=Xiβ^p_i =P(Y_i = 1)= E[Y_i] = X_i \hat{\beta} pi=P(Yi=1)=E[Yi]=Xiβ^
表示的是成功概率。这个模型比较直白，问题也比较多。

残差项不满足正态假设
给定样本时，残差只有两个可能的取值，当Yi=0Y_i=0Yi=0时，ϵi=−Xiβ^\epsilon_i = -X_i \hat{\beta}ϵi=−Xiβ^，当Yi=1Y_i=1Yi=1时，ϵi=1−Xiβ^\epsilon_i = 1-X_i \hat{\beta}ϵi=1−Xiβ^，显然这不服从正态分布。
同方差假设不成立
残差与YiY_iYi同分布，σ2(ϵi)=−Xiβ^(1−Xiβ^)\sigma^2(\epsilon_i)=-X_i \hat{\beta}(1-X_i \hat{\beta})σ2(ϵi)=−Xiβ^(1−Xiβ^)，显然与XiX_iXi有关，同方差假设不成立。
回归方程取值受限
由于拟合值的含义是概率，因此E[Yi]=Xiβ^∈[0,1]E[Y_i] = X_i \hat{\beta} \in [0,1]E[Yi]=Xiβ^∈[0,1]，否则拟合值无意义。

所以二值被解释变量一般做如下处理。假设YicY_i^cYic是一个连续型的变量，但在观测时只取0和1，假设取值的规律为
Yi=1,ifYic≤Y0Yi=0,ifYic>Y0Y_i = 1,\ if\ Y_i^c\le Y_0 \\Y_i = 0,\ if\ Y_i^c > Y_0 Yi=1, if Yic≤Y0Yi=0, if Yic>Y0
则
P(Yi=1)=P(Yic≤Y0)P(Y_i=1) = P(Y_i^c \le Y_0) P(Yi=1)=P(Yic≤Y0)
假设YicY_i^cYic满足线性回归模型
Yic=Xiβ+ϵi,ϵi∼N(0,σc2)Y_i^c = X_i \beta + \epsilon_i,\ \epsilon_i \sim N(0,\sigma^2_c) Yic=Xiβ+ϵi, ϵi∼N(0,σc2)
则
P(Yi=1)=P(Yic≤Y0)=P(Xiβ+ϵi≤Y0)=P(ϵi≤Y0−Xiβ)=P(ϵiσc≤Y0σc−Xiβσc)P(Y_i=1) = P(Y_i^c \le Y_0) = P(X_i \beta + \epsilon_i \le Y_0 )\\ =P(\epsilon_i \le Y_0 - X_i \beta) = P(\frac{\epsilon_i}{\sigma_c} \le \frac{Y_0}{\sigma_c} - X_i \frac{\beta}{\sigma_c}) P(Yi=1)=P(Yic≤Y0)=P(Xiβ+ϵi≤Y0)=P(ϵi≤Y0−Xiβ)=P(σcϵi≤σcY0−Xiσcβ)

Probit模型

用Φ\PhiΦ表示标准正态分布的分布函数，假设
P(Yi=1)=P(ϵiσc≤Y0σc−Xiβσc)=Φ(Xiβ∗)P(Y_i=1) = P(\frac{\epsilon_i}{\sigma_c} \le \frac{Y_0}{\sigma_c} - X_i \frac{\beta}{\sigma_c}) = \Phi(X_i \beta^*) P(Yi=1)=P(σcϵi≤σcY0−Xiσcβ)=Φ(Xiβ∗)
其中
β0∗=Y0σc−β0σcβi∗=−βiσc\beta_0^* = \frac{Y_0}{\sigma_c} - \frac{\beta_0}{\sigma_c} \\ \beta_i^* = - \frac{\beta_i}{\sigma_c} β0∗=σcY0−σcβ0βi∗=−σcβi
这个模型叫Probit模型。

Logit模型

Logit模型又叫Logistics回归。对于
P(Yi=1)=P(ϵiσc≤Y0σc−Xiβσc)=P(ϵiσc≤Xiβ∗)P(Y_i=1) = P(\frac{\epsilon_i}{\sigma_c} \le \frac{Y_0}{\sigma_c} - X_i \frac{\beta}{\sigma_c}) = P(\frac{\epsilon_i}{\sigma_c} \le X_i \beta^*) P(Yi=1)=P(σcϵi≤σcY0−Xiσcβ)=P(σcϵi≤Xiβ∗)
定义
ϵL=π3ϵiσc,β=π3β∗\epsilon_L = \frac{\pi}{\sqrt{3}} \frac{\epsilon_i}{\sigma_c}, \beta = \frac{\pi}{\sqrt{3}} \beta^* ϵL=3πσcϵi,β=3πβ∗
这里的π3\frac{\pi}{\sqrt{3}}3π没有别的意思，就是随便给的。假设ϵL\epsilon_LϵL的分布函数为Logistics函数
F(ϵL)=exp⁡(ϵL)1+exp⁡(ϵL)F(\epsilon_L) = \frac{\exp{(\epsilon_L)}}{1+\exp{(\epsilon_L)}} F(ϵL)=1+exp(ϵL)exp(ϵL)
从而
P(Yi=1)=P(ϵiσc≤Xiβ∗)=exp⁡(Xiβ)1+exp⁡(Xiβ)=11+exp⁡(−Xiβ)P(Y_i=1) = P(\frac{\epsilon_i}{\sigma_c} \le X_i \beta^*) = \frac{\exp{(X_i \beta)}}{1+\exp{(X_i \beta)}} = \frac{1}{1+\exp{(-X_i\beta)}} P(Yi=1)=P(σcϵi≤Xiβ∗)=1+exp(Xiβ)exp(Xiβ)=1+exp(−Xiβ)1
这个模型就是Logit模型。相比Probit模型，Logit模型的系数更好解释。
exp⁡(Xiβ)=P(Yi=1)P(Yi=0)⟺Xiβ=ln(P(Yi=1)P(Yi=0))=ln(oddi)\exp{(X_i\beta)} = \frac{P(Y_i=1)}{P(Y_i=0)} \Longleftrightarrow X_i\beta = ln(\frac{P(Y_i=1)}{P(Y_i=0)}) = ln(odd_i) exp(Xiβ)=P(Yi=0)P(Yi=1)⟺Xiβ=ln(P(Yi=0)P(Yi=1))=ln(oddi)
其中P(Yi=1)P(Yi=0)\frac{P(Y_i=1)}{P(Y_i=0)}P(Yi=0)P(Yi=1)叫做输赢比（odd ratio），因此系数的含义是解释变量对对数输赢比的单位影响。
β=∂ln(oddi)∂Xi\beta = \frac{\partial ln(odd_i)}{\partial X_i} β=∂Xi∂ln(oddi)

系数的最大似然估计

定义pi=P(Yi=1)p_i=P(Y_i=1)pi=P(Yi=1)，显然YiY_iYi服从Bernoulli分布Ber(pi)Ber(p_i)Ber(pi)，因此
f(Yi)=piYi(1−pi)1−Yif(Y_i) = p_i^{Y_i}(1-p_i)^{1-Y_i} f(Yi)=piYi(1−pi)1−Yi
样本的对数似然函数为
l(β)=∑i=1N[Yilnpi+(1−Yi)ln(1−pi)]=∑i=1NYiln⁡(pi1−pi)+∑i=1Nln(1−pi)l(\beta) = \sum_{i=1}^N [Y_ilnp_i+(1-Y_i)ln(1-p_i)] \\= \sum_{i=1}^N Y_i \ln{(\frac{p_i}{1-p_i})} + \sum_{i=1}^N ln(1-p_i) l(β)=i=1∑N[Yilnpi+(1−Yi)ln(1−pi)]=i=1∑NYiln(1−pipi)+i=1∑Nln(1−pi)
其中pip_ipi是β\betaβ的函数
pi=11+exp⁡(−Xiβ),1−pi=11+exp⁡(Xiβ)p_i = \frac{1}{1+\exp{(-X_i\beta)}},1-p_i=\frac{1}{1+\exp{(X_i\beta)}} pi=1+exp(−Xiβ)1,1−pi=1+exp(Xiβ)1
ln⁡(pi1−pi)=Xiβ\ln{(\frac{p_i}{1-p_i})} =X_i\beta ln(1−pipi)=Xiβ
因此
l(β)=∑i=1NYiXiβ−∑i=1Nln(1+exp⁡(Xiβ))l(\beta) = \sum_{i=1}^N Y_i X_i \beta- \sum_{i=1}^N ln(1+\exp{(X_i\beta)}) l(β)=i=1∑NYiXiβ−i=1∑Nln(1+exp(Xiβ))
这个对数似然函数最大化的解析解没法计算出来，通常考虑数值方法求解，比较常用的是上一篇博文中提到的Iterative Reweighted Least Square.

系数的推断

定义对数似然函数的Hessian Matrix为
G=∂2l(β)∂β2G = \frac{\partial^2 l(\beta) }{\partial \beta^2} G=∂β2∂2l(β)
则系数最大似然估计的方差估计为
s2(β^)=(−G(β^))−1s^2(\hat{\beta}) = (-G(\hat{\beta}))^{-1} s2(β^)=(−G(β^))−1

Wald检验

在Logistics回归中，对单个系数的检验需要Wald检验。
H0:β1=0Ha:β1≠0H_0: \beta_1 = 0 \\ H_a: \beta_1 \ne 0 H0:β1=0Ha:β1=0
当样本数足够大时，
Z=β^1−β1s(β^1)∼N(0,1)Z=\frac{\hat{\beta}_1 - \beta_1}{s(\hat{\beta}_1)} \sim N(0,1) Z=s(β^1)β^1−β1∼N(0,1)
因此Wald检验其实是一种Z检验。

似然比检验

似然比检验可以用来对部分或者所有系数进行检验。
H0:βq=βq+1=...=βp−1=0Ha:notallcoefficientsinH0equaltozeroH_0: \beta_q = \beta_{q+1} = ... =\beta_{p-1} = 0 \\ H_a: not\ all\ coefficients\ in\ H_0\ equal\ to\ zero H0:βq=βq+1=...=βp−1=0Ha:not all coefficients in H0 equal to zero
Suppose L(R)L(R)L(R) to be likelihood of reduced model and L(F)L(F)L(F) to be likelihood of full model. When sample size is large enough
G2=−ln(L(R)L(F))∼χ2(p−q)G^2 = -ln(\frac{L(R)}{L(F)}) \sim \chi^2(p-q) G2=−ln(L(F)L(R))∼χ2(p−q)

二项回归

假设对解释变量的取值Xj,j=1,2,...,cX_j, j= 1,2,...,cXj,j=1,2,...,c做了njn_jnj次重复观察，得到的被解释变量的取值为Yij,j=1,2,...,njY_{ij},j=1,2,...,n_jYij,j=1,2,...,nj，则解释变量取值为XjX_jXj时，观察到被解释变量等于1的次数为
Y.j=∑i=1njYij∼Binom(nj,pj)Y_{.j} = \sum_{i=1}^{n_j} Y_{ij} \sim Binom(n_j,p_j) Y.j=i=1∑njYij∼Binom(nj,pj)
pj=11+exp⁡(−Xjβ)p_j = \frac{1}{1+\exp{(-X_j\beta)}} pj=1+exp(−Xjβ)1
这个模型叫二项回归（Binomial Regression）。

拟合优度检验

拟合优度检验（Goodness of fit test）的作用是检验模型整体的拟合质量，对小部分拟合不好的地方不会很敏感。如果样本数据存在replication，可以使用卡方拟合优度检验或者Deviance拟合优度检验，如果样本数据不存在replication，可以使用Hosmer-Lemeshow拟合优度检验。Logistics回归的拟合优度检验假设为
H0:P(Yi=1)=11+exp⁡(−Xiβ)Ha:P(Yi=1)≠11+exp⁡(−Xiβ)H_0: P(Y_i=1) = \frac{1}{1+\exp{(-X_i\beta)}} \\ H_a: P(Y_i=1) \ne \frac{1}{1+\exp{(-X_i\beta)}} H0:P(Yi=1)=1+exp(−Xiβ)1Ha:P(Yi=1)=1+exp(−Xiβ)1

Pearson卡方拟合优度检验

使用二项回归关于样本replication的设定，
O1j=Y.j,O0j=nj−Y.jO_{1j} = Y_{.j}, O_{0j} = n_j - Y_{.j} O1j=Y.j,O0j=nj−Y.j
在原假设下，假设拟合值为
p^j=11+exp⁡(−Xiβ^)\hat{p}_j = \frac{1}{1+\exp{(-X_i\hat{\beta})}} p^j=1+exp(−Xiβ^)1
被解释变量取1和0的理论样本数为
E1j=njp^j,E0j=nj(1−p^j)E_{1j} = n_j \hat{p}_j, E_{0j} = n_j (1 - \hat{p}_j) E1j=njp^j,E0j=nj(1−p^j)
当p<cp<cp<c且replication数量njn_jnj足够大时
χ2=∑j=1c[(O0j−E0j)2E0j+(O1j−E1j)2E1j]∼χ2(c−p)\chi^2 = \sum_{j=1}^c [ \frac{(O_{0j}-E_{0j})^2}{E_{0j}} + \frac{(O_{1j}-E_{1j})^2}{E_{1j}} ] \sim \chi^2(c-p) χ2=j=1∑c[E0j(O0j−E0j)2+E1j(O1j−E1j)2]∼χ2(c−p)

Deviance拟合优度检验

Deviance拟合优度检验用的是似然比检验的思想。其reduced model是Logistics回归，full model是
E(Yij)=pjE(Y_{ij}) = p_j E(Yij)=pj
这时似然比检验的统计量又被叫做Deviance，记作Dev(X)Dev(X)Dev(X)，它代表这两个模型之间的deviation。

Hosmer-Lemeshow拟合优度检验

Hosmer-Lemeshow拟合优度检验是卡方拟合优度检验的直接推广，nj=1n_j=1nj=1代表样本数据没有replication，nj>1n_j>1nj>1代表样本数据有replication。

多值被解释变量

假设YiY_iYi有mmm个可能的取值，不失一般性，假设为1,2,...,m1,2,...,m1,2,...,m。假设YiY_iYi服从Boltzmann分布
P(Yi=j)=11+∑k=1m−1exp⁡(−Xiβk),j=1P(Yi=j)=exp⁡(Xiβj)1+∑k=1m−1exp⁡(−Xiβk),j>1P(Y_i=j) = \frac{1}{1+\sum_{k=1}^{m-1}\exp{(-X_i \beta_k)}}, j = 1 \\ P(Y_i=j) = \frac{\exp{(X_i \beta_j)}}{1+\sum_{k=1}^{m-1}\exp{(-X_i \beta_k)}}, j > 1 P(Yi=j)=1+∑k=1m−1exp(−Xiβk)1,j=1P(Yi=j)=1+∑k=1m−1exp(−Xiβk)exp(Xiβj),j>1
这个模型叫多值Logit模型（MLogit）。Boltzmann分布用来描述第jjj个能级的粒子数，假设共有mmm个能级，则第jjj个能级的粒子数占系统总粒子数的比值为
pj=exp(−ϵj/kT)∑k=1mexp(−ϵk/kT)=exp((ϵ1−ϵj)/kT)∑k=1mexp((ϵ1−ϵk)/kT)p_j = \frac{exp(-\epsilon_j/kT)}{\sum_{k=1}^m exp(-\epsilon_k/kT)} = \frac{exp((\epsilon_1-\epsilon_j)/kT)}{\sum_{k=1}^m exp((\epsilon_1-\epsilon_k)/kT)} pj=∑k=1mexp(−ϵk/kT)exp(−ϵj/kT)=∑k=1mexp((ϵ1−ϵk)/kT)exp((ϵ1−ϵj)/kT)
Boltzmann分布用来描述多分类任务的想法还是比较直接。系数βk\beta_kβk的大小就相当于是ϵk−ϵ1\epsilon_k-\epsilon_1ϵk−ϵ1，即第kkk能级与第1能级的能量差。能级的能量越大，能级的粒子数就越少，因此当解释变量的取值范围为正实数时，系数βk\beta_kβk越大，Yi=kY_i=kYi=k的可能性就越小。第1能级包含的粒子数占比被标准化为
11+∑k=2mexp⁡(ϵ1−ϵk)/kT),\frac{1}{1+\sum_{k=2}^{m}\exp{(\epsilon_1-\epsilon_k)/kT)}}, 1+∑k=2mexp(ϵ1−ϵk)/kT)1,
对应在mLogit中就是Yi=1Y_i=1Yi=1的概率。

UA MATH571A 多元线性回归IV 广义线性模型相关推荐

UA MATH571A 一元线性回归IV 模型诊断
UA MATH571A 一元线性回归IV 模型诊断解释变量解释变量的可视化残差残差的性质 Semistudentized Residual 残差的可视化残差关于解释变量的图残差关于拟合值的 ...
UA MATH571A 多元线性回归V 自相关与非线性模型简介
UA MATH571A 多元线性回归V 自相关与非线性模型简介一阶误差自相关模型 Durbin-Watson检验一阶自相关的消去 Cochrane-Orcutt方法 Hildreth-Lu方法非 ...
UA MATH571A 多元线性回归I 模型设定与推断
UA MATH571A 多元线性回归I 模型设定与推断模型设定最小二乘法(Method of Least Square) 系数 Mean Response and Residual 多元回归的AN ...
UA MATH571A 多元线性回归II 变量选择
UA MATH571A 多元线性回归II 变量选择多项式回归与交互项回归阶数的确定含质量型变量的回归含质量型变量的交互项二值变量与二值变量的交互项二值变量与数量型变量的交互项变量选择的准 ...
UA MATH571A 一元线性回归III 方差分析与相关性分析
UA MATH571A 一元线性回归III 方差分析与相关性分析 ANOVA Table F检验回归系数的F检验 F检验与t检验等价广义线性检验方法 R2R^2R2 数值例子:女性肌肉量与年龄的关 ...
UA MATH571A 一元线性回归II 统计推断2
UA MATH571A 一元线性回归II 统计推断2 β0\beta_0β0的分布拟合与预测拟合值的区间估计预测值的区间估计数值例子:女性肌肉量与年龄的关系 β0\beta_0β0的分布 ...
UA MATH571A 一元线性回归II 统计推断1
UA MATH571A 一元线性回归II 统计推断1 β1\beta_1β1的假设检验与置信区间 Gauss-Markov定理检验的势双边检验,单边检验与置信区间置信区间双边检验单边检验 ...
UA MATH571A 一元线性回归I 模型设定与估计
UA MATH571A 一元线性回归I 模型设定与估计模型设定最小二乘法(Method of Least Square) Coefficients Mean Response and Residu ...
PySpark线性回归与广义线性模型
PySpark线性回归与广义线性模型 1.线性回归 2.岭回归(Ridge Regression)与LASSO回归(LASSO Regression) 3.广义线性模型 (GLM) 本文为销量预测第7 ...

UA MATH571A 多元线性回归IV 广义线性模型