LASSO逻辑回归模型

主要用于对样本分类，假设样本个数是 $n$ ，变量个数是 $p$ ， $p ≫ n$ ，因变量 $Y$ 是0-1变量，表示样本的两种分类，比如取值1表示样本是diseased，取值0表示样本是正常样本。
数据
- 因变量 $Y = (Y_{1}, Y_{2}, \dots, Y_{n})^{T}$
- 特征矩阵 $X = (X_{1}, X_{2}, \dots, X_{p}) = (X_{(1)}, X_{(2)}, \dots, X_{(n)})^{T}$
假设 $Y \sim B e r n o u l l i (π)$
建立逻辑回归模型

$l n \frac{π_{i}}{1 - π_{i}} = X_{(i)}^{T} β, i = 1, 2, \dots, n$

1. LASSO逻辑回归模型的求解过程

使用极大似然估计法估计参数 $β$ ，似然函数为：

$\prod_{i = 1}^{n} π_{i}^{Y_{i}} (1 - π_{i})^{(1 - Y_{i})}$
对似然函数取负对数得到损失函数：

$l (β) = - \frac{1}{n} \sum_{i = 1}^{n} [Y_{i} l n (π_{i}) + (1 - Y_{i}) l n (1 - π_{i})]$
通过牛顿-拉普森算法将最小化上面的损失函数 $l (β)$ 转化为求解线性回归的重加权最小二乘问题
- 将 $l (β)$ 在 $β$ 的当前迭代值 $\tilde{β}$ 处进行二阶泰勒展开，得到：
  
  $l_{Q} (β) = l (\tilde{β}) + (β - \tilde{β})^{T} \nabla l (\tilde{β}) + \frac{1}{2} (β - \tilde{β})^{T} H (β - \tilde{β})$
  
  其中， $\nabla l (\tilde{β}) = \frac{\partial l (β)}{\partial β} |_{β = \tilde{β}} = - \frac{1}{n} X^{T} (Y - \tilde{Π})$ 是 $l (β)$ 在 $\tilde{β}$ 处的一阶导数， $H = \frac{\partial^{2} l (β)}{\partial β \partial β^{T}} |_{β = \tilde{β}} = \frac{1}{n} X^{T} Q X$ 是 $l (β)$ 在 $\tilde{β}$ 处的Hessian矩阵。
  
  $\tilde{Π} = [\begin{matrix} {\tilde{π}}_{1} \\ {\tilde{π}}_{2} \\ \dots \\ {\tilde{π}}_{n} \end{matrix}] = [\begin{matrix} \frac{e x p {X_{(1)}^{T} \tilde{β}}}{1 + e x p {X_{(1)}^{T} \tilde{β}}} \\ \frac{e x p {X_{(2)}^{T} \tilde{β}}}{1 + e x p {X_{(2)}^{T} \tilde{β}}} \\ \dots \\ \frac{e x p {X_{(n)}^{T} \tilde{β}}}{1 + e x p {X_{(n)}^{T} \tilde{β}}} \end{matrix}]$
  $Q = [\begin{matrix} q_{11} & 0 & \dots & 0 \\ 0 & q_{22} & \dots & 0 \\ ⋮ & ⋮ & ⋱ & ⋮ \\ 0 & 0 & \dots & q_{n n} \end{matrix}] = [\begin{matrix} {\tilde{π}}_{1} (1 - {\tilde{π}}_{1}) & 0 & \dots & 0 \\ 0 & {\tilde{π}}_{2} (1 - {\tilde{π}}_{2}) & \dots & 0 \\ ⋮ & ⋮ & ⋱ & ⋮ \\ 0 & 0 & \dots & {\tilde{π}}_{n} (1 - {\tilde{π}}_{n}) \end{matrix}] = [\begin{matrix} \frac{e x p {X_{(1)}^{T} \tilde{β}}}{{(1 + e x p {X_{(1)}^{T} \tilde{β}})}^{2}} & 0 & \dots & 0 \\ 0 & \frac{e x p {X_{(2)}^{T} \tilde{β}}}{{(1 + e x p {X_{(2)}^{T} \tilde{β}})}^{2}} & \dots & 0 \\ ⋮ & ⋮ & ⋱ & ⋮ \\ 0 & 0 & \dots & \frac{e x p {X_{(n)}^{T} \tilde{β}}}{{(1 + e x p {X_{(n)}^{T} \tilde{β}})}^{2}} \end{matrix}]$
- 令 $l_{Q} (β)$ 的梯度为0：
  
  $\frac{\partial l_{Q} (β)}{\partial β} = \nabla l (\tilde{β}) + H (β - \tilde{β}) = 0$
- 得到 $β$ 的迭代公式：
  
  $β = (X^{T} Q X)^{- 1} X^{T} (Y - \tilde{Π}) + (X^{T} Q X)^{- 1} X^{T} Q X \tilde{β} = (X^{T} Q X)^{- 1} X^{T} Q [Q^{- 1} (Y - \tilde{Π}) + X \tilde{β}]$
- 令 $Z = Q^{- 1} (Y - \tilde{Π}) + X \tilde{β}$ ，则 $β = (X^{T} Q X)^{- 1} X^{T} Q Z$ ， $β$ 可以看作新的因变量 $Z$ 与数据矩阵 $X$ 的线性回归加权最小二乘的解，因此，可以通过坐标下降法求解。
通过坐标下降法求解具有LASSO惩罚项的加权最小二乘问题
- 此时的目标函数为：
  
  $L (β) = \frac{1}{2 n} \sum_{i = 1}^{n} q_{i i} (Z_{i} - X_{(i)}^{T} β)^{2} + λ \sum_{j = 1}^{p} | β_{j} |$
- 令梯度为0：
  
  $\frac{\partial L (β)}{\partial β_{j}} = - \frac{1}{n} \sum_{i = 1}^{n} q_{i i} x_{i j} (Z_{i} - \sum_{l \neq j}^{p} x_{i l} β_{l}) + \frac{1}{n} \sum_{i = 1}^{n} q_{i i} x_{i j}^{2} β_{j} + λ s i g n (β_{j}) = 0$
- 分情况讨论 $β_{j}$ 的迭代公式
  
  如果 $β_{j} > 0$ ，有 $β_{j} = \frac{\sum_{i = 1}^{n} q_{i i} x_{i j} (Z_{i} - \sum_{l \neq j}^{p} x_{i l} β_{l}) - n λ}{\sum_{i = 1}^{n} q_{i i} x_{i j}^{2}}$
  
  如果 $β_{j} < 0$ ，有 $β_{j} = \frac{\sum_{i = 1}^{n} q_{i i} x_{i j} (Z_{i} - \sum_{l \neq j}^{p} x_{i l} β_{l}) + n λ}{\sum_{i = 1}^{n} q_{i i} x_{i j}^{2}}$
  
  $β_{j} = 0$ 当且仅当 $n λ \geq | \sum_{i = 1}^{n} q_{i i} x_{i j} (Z_{i} - \sum_{l \neq j}^{p} x_{i l} β_{l}) |$

2. MATLAB代码

通过glmnet包求解

clear;
clc;
load data;   #加载数据，包括数据矩阵x，因变量y
x=zscore(x);
fold = 10;    #10折交叉验证
m = 10;       #参数lambda的个数
options = glmnetSet;
options.alpha = 1;
options.nlambda = m;
fit = glmnet(x,y+1,'binomial',options);
CVerr = cvglmnet(x,y+1,fold,[],'class','binomial',glmnetSet,0);
beta_path = fit.beta;    #得到对应于所有lambda的系数矩阵
intercept_path = fit.a0;   #截距项
lambda_path = fit.lambda;   #lambda的取值
opt_index = find(CVerr.glmnetOptions.lambda == CVerr.lambda_min); #最优lambda所在位置
beta = beta_path(:,opt_index);
selectedgene = find(beta~=0);   #模型筛选出来的变量所在位置