逻辑回归

1、问题提???

　　实际中，往往会出现响应变量为二分类变量（或多分类）的情况，此时若继续使用OLS方法对数据建立回归模型，会出现预测值与实际值取值范围不一致的情况，为此我们引出逻辑回归模型???

2、Logistic模型

???1）非线性变???
　　考虑到线性模???\(\sum \beta_k x_k\)的取值范围为\((-\infty,\infty)\),我们希望通过某种变换将其取值范围映射到区间(0,1)上，而logit变换和probit变换正好符合此条件。下面分别介绍logit模型和probit模型???
???1-1）logit模型???
　　logit函数的形式如下： \[logit(p)=ln\frac{p}{1-p},p\in(0,1)，logit(p)\in(-\infty,\infty)\] 其反函数称作logistic函数??? \[logistic(\eta)=\frac{exp(\eta)}{1+exp(\eta)},\eta \in (-\infty,\infty),logistic(\eta)\in(0,1)\] ???\(\eta=ln\frac{p}{1-p}\)时，???\(logit(p)=\eta\),\(logistic(\eta)=p\)，下面做出二者的函数图像???

par(mfrow=c(1,2))
p=seq(0,1,by=0.0001)
plot(p,log(p/(1-p)),type="l",ylab="logit(p)","main"="Logit")
n=seq(-5,5,0.1)
plot(n,exp(n)/(1+exp(n)),type="l",ylab="logistic(n)",main="Logistic")

???\(\eta=\sum \beta_k x_k\),???\(logistic(\eta)\in (0,1)\)，解决了区间不一致的问题???
???1-2）probit模型
　　和logit模型相同，probit模型选取标准正态分布的分布函数作为变换的“机制”，标准正态分布的反函数即probit函数,???\(\eta=\sum \beta_k x_k\),有： \[probit(p)=\Phi^{-1}(p)=\eta=\sum\beta_k x_k\] \[p=\Phi(\sum\beta_k x_k)\] 其中\(\Phi\)表示标准正态分布的分布函数???

无论logit模型还是probit模型，在解决这一问题时都在利用分布函数的值域是（0,1)这一性质。logit模型利用logistic分布，probit模型利用标准正态分布???

???2）潜变量模型

　　对逻辑回归模型也可以从潜变量模型的角度来解释：认为存在一个未被观测到的或潜在的连续变???\(y^*\),它与样本i是否出现\(y_i=1\)的的关系表示??? \[y_i= \begin{cases}1,\quad\ y_i^* \geq 0\\ 0,\quad\ y_i^* \leq 0 \end{cases} \] 潜在变量同时可以表示???\(x_{ik}\)与误???\(e_i\)的线性函数： \[y^*=\sum\beta_k x_{ik}+e_i\] \(y^*\)相当???\(y_i\)???\(\sum\beta_k x_{ik}+e_i\)之间建立的一座“桥梁。于???: \[P(y_i=1)=P(y_i^*\geq 0)=P(e_i\geq -\eta_i)=P(e_i\leq\eta_i)\]
???\(e_i\)~\(N(0,1)\), \(p_i=\Phi(\eta_i)\);???\(e_i\)~\(logistic(0,1)\),\(p_i=\frac{exp(\eta)}{1+exp(\eta)}\)

3、参数估???

???1）极大似估计

　　在样本独立同分布的情况下，可以通过构造似然函数来求模型参数。似然函数的形式如下??? \[ L(\theta)= \begin{cases}\prod f(x_i,\theta),(连续)\\ \prod p(x_i,\theta),(离散) \end{cases} \] eg：已知独立同分布的样本：\(x_1,x_2,\cdots,x_i\)假设服从正态分???/二项分布/泊松分布，分别求总体参数\(\mu、p、\lambda\)???
　　对于极大似然估计，应该注意两个事实：???1）在正态分布的情况下，极大似然估计和最小二乘估计等价；???2）极大似然估计是M估计的一个特例。下面对这两个事实作简要说明： 设：\(x_1,x_2,\cdots,x_n \sim f(x_i,\theta)\),???:\(\theta_{MLE}=argmax\sum lnf(x_i,\theta)=argmin\sum -lnf(x_i,\theta)=argmin\sum P(x)\)???(2)成立。当\(x_1,x_2,\cdots,x_n \sim N(\mu,\sigma^2)\),\(\theta_{MLE}=argmin\sum (x_i-\mu)^2K\),等价于OLS，（1）成立???

???2）参数估???

　　构建似然函数:\[L(\beta)=\prod p_i^{y_i}(1-p_i)^{1-y_i}\] 左右同时取对数得到：\[ln(L(\beta))=\sum　y_i(x'\beta)-\sum ln(1+exp(x'\beta))\]求偏导可???\(\hat \beta\)???

???3）参数检验与模型检???
???4）回归系数解???

logit模型的着重点从事件的发生概率转移到时间的发生比率。注意以下几个常用指标：
* 发生比率（odds???:\(\omega=p/(1-p)\) （就???\(e^\eta\)???
* 发生比率比（odds ratio???:\(\theta =\omega_1 /\omega_2\)
* 相对风险（relative risk???:\(p=p_1/p_2\)
logit模型对分类变量和连续变量的解释完全不同???
　　

4、应用举???

???1）电击试???

　　实验：选择7头牛，六种电击强度（0,1,2,3,4,5）。每头牛被电???60下，每种强度10下，按随机次序进行。对每次电击响应变量（嘴巴运动）出现或不出现。给出每种电击强???70次实验中牛发生相应的总次数。试分析电击对牛的影响??? 根据题意，建立逻辑回归模型??? \[ln(\frac{p}{1-p})=\beta_0+\beta_1x\]

data<-read.csv("C:\\Users\\JHON\\Desktop\\应用回归分析\\逻辑回归\\lightexp.csv")
names(data)<-c("x","n","k")
mod<-glm(cbind(data[,3],data[,2]-data[,3])~data[,1],family=binomial)

根据以上结果，可以得到逻辑回归模型???
\[p=\frac{exp(-3.301+1.246x)}{1+exp(-3.301+1.246x)}\] 根据模型，我们可以预测或者控制：

pre<-predict(mod,newdata = data.frame(x=3.5),type="response")

## Warning: 'newdata' had 1 row but variables found have 6 rows

说明电流强度???3.5毫安时牛响应的概率为0.74???

-mod$coe[[1]]/mod$coe[[2]]

## [1] 2.649439

p=0.5时，\(ln\frac{p}{1-p}=0\)所???\(x=-\beta_0/\beta_1\),若要控制牛的响应概率???50%，则电流应该控制???2.65???
???2）交通事???
　　已知：y：是否出过事故；\(x_1\)：视力状况；\(x_2\)：年龄；\(x_3\)：驾车教育；分析\(x_1、x_2、x_3\)对y的影响??? 根据题意，建立逻辑回归模型??? \[ln\frac{p}{1-p}=\beta_0+\beta_1x_1+\beta_2x_2+\beta_3x_3\]

data<-read.csv("C:\\Users\\JHON\\Desktop\\应用回归分析\\逻辑回归\\driver.csv")
mod2<-glm(y~x1+x2+x3,family=binomial,data=data)
summary(mod2)$coef

##                 Estimate Std. Error     z value   Pr(>|z|)
## (Intercept)  0.597609952 0.89483074  0.66784692 0.50423131
## x1          -1.496083919 0.70486126 -2.12252254 0.03379388
## x2          -0.001595183 0.01675828 -0.09518775 0.92416570
## x3           0.315864768 0.70109346  0.45053162 0.65232716

根据检验结果，发现\(\beta_2,\beta_3\)没有通过检验，用step()筛选变量???

mod3<-step(mod2)

## Start:  AIC=65.03
## y ~ x1 + x2 + x3
## 
##        Df Deviance    AIC
## - x2    1   57.035 63.035
## - x3    1   57.232 63.232
## <none>      57.026 65.026
## - x1    1   61.936 67.936
## 
## Step:  AIC=63.03
## y ~ x1 + x3
## 
##        Df Deviance    AIC
## - x3    1   57.241 61.241
## <none>      57.035 63.035
## - x1    1   61.991 65.991
## 
## Step:  AIC=61.24
## y ~ x1
## 
##        Df Deviance    AIC
## <none>      57.241 61.241
## - x1    1   62.183 64.183

summary(mod3)$coef

##               Estimate Std. Error   z value   Pr(>|z|)
## (Intercept)  0.6190392  0.4688072  1.320456 0.18668286
## x1          -1.3728110  0.6352980 -2.160893 0.03070362

最终模型可以写成： \[p=\frac{exp(0.619-1.317x)}{1+exp(0.619-1.317x)}\]

p1<-predict(mod3,newdata=data.frame(x1=1),type="response")
p2<-predict(mod3,newdata=data.frame(x1=0),type="response")
p1/p2

##         1 
## 0.4923077