statistical-computation.qmd

# 统计计算 {#sec-statistical-computation}

```{r}
#| echo: false

source("_common.R")
```

::: hidden
$$
 \def\bm#1{{\boldsymbol #1}}
$$
:::

每一个统计模型的背后都有一个优化问题，统计计算的任务就是求解优化问题。

## 回归问题与优化问题 {#sec-regression-optimization}

1996 年出现 Lasso （Least Absolute Selection and Shrinkage Operator，简称 Lasso）[@lasso1996]，由于缺少高效的求解算法，Lasso 在高维小样本特征选择研究中没有广泛流行，最小角回归（Least Angle Regression，简称 LAR）算法 [@lar2004] 的出现有力促进了 Lasso 在高维小样本数据中的应用。为了解决 Lasso 的有偏估计问题，自适应 Lasso、松弛 Lasso， SCAD （Smoothly Clipped Absolute Deviation，简称 SCAD）[@scad2008]，MCP (Minimax Concave Penalty，简称 MCP)[@mcp2010] 陆续出现。经典的普通最小二乘、广义最小二乘、岭回归、逐步回归、Lasso 回归、最优子集回归都可转化为优化问题。具体地，一个带 L1 正则项的线性回归模型，其对应的优化问题如下：

$$
\arg \min_{\boldsymbol{\beta},\lambda} ~~ \frac{1}{2} || \bm{y} - X \boldsymbol{\beta} ||_2^2 +  \lambda ||\boldsymbol{\beta}||_1
$$

其中，$X \in \mathbb{R}^{n\times k}$， $\bm{y} \in \mathbb{R}^n$，$\boldsymbol{\beta} \in \mathbb{R}^k$， $0 < \lambda \in \mathbb{R}$ 。下面以逻辑回归模型为例，介绍 R 语言中求解此类优化问题的方法。

## 对数似然与损失函数 {#sec-log-likelihood}

### Logistic 分布 {#sec-logistic-distribution}

在介绍逻辑回归之前，先了解一下 Logistic 分布。一个均值为 $m$ ，方差为 $\frac{\pi^2}{3}s^2$ 的 Logistic 分布函数的形式为

$$
F(x) = \frac{1}{1 + \exp(-\frac{x - m}{s})}
$$

密度函数的形式为

$$
f(x) = \frac{\exp(-\frac{x - m}{s})}{s(1 + \exp(-\frac{x-m}{s}))^2} = \frac{\exp(\frac{x - m}{s})}{s(1 + \exp(\frac{x-m}{s}))^2}
$$

密度函数与分布函数的关系如下：

$$
\frac{dF(x)}{dx} = f(x) = sF(x)(1 - F(x))
$$

也就是说 Logistic 分布是上述微分方程的解。

```{r}
#| label: fig-logistic
#| echo: false
#| fig-cap: "逻辑斯谛分布"
#| fig-subcap:
#| - 概率密度函数
#| - 概率分布函数
#| fig-width: 4
#| fig-height: 3
#| layout-ncol: 2
#| fig-showtext: true

library(ggplot2)
ggplot() +
  geom_function(
    fun = dlogis, args = list(location = 0, scale = 0.5),
    colour = "#E41A1C", linewidth = 1.2, xlim = c(-6, 6)
  ) +
  geom_function(
    fun = dlogis, args = list(location = 0, scale = 1),
    colour = "#377EB8", linewidth = 1.2, xlim = c(-6, 6)
  ) +
  geom_function(
    fun = dlogis, args = list(location = 0, scale = 2),
    colour = "#4DAF4A", linewidth = 1.2, xlim = c(-6, 6)
  ) +
  theme_classic() +
  labs(x = expression(x), y = expression(f(x)))

ggplot() +
  geom_function(
    fun = plogis, args = list(location = 0, scale = 0.5),
    colour = "#E41A1C", linewidth = 1.2, xlim = c(-6, 6)
  ) +
  geom_function(
    fun = plogis, args = list(location = 0, scale = 1),
    colour = "#377EB8", linewidth = 1.2, xlim = c(-6, 6)
  ) +
  geom_function(
    fun = plogis, args = list(location = 0, scale = 2),
    colour = "#4DAF4A", linewidth = 1.2, xlim = c(-6, 6)
  ) +
  theme_classic() +
  labs(x = expression(x), y = expression(bolditalic(F)(x)))
```

R 语言中分别表示逻辑斯谛分布的密度函数、分布函数、分位函数和随机数生成函数如下：

``` r
dlogis(x, location = 0, scale = 1, log = FALSE)
plogis(q, location = 0, scale = 1, lower.tail = TRUE, log.p = FALSE)
qlogis(p, location = 0, scale = 1, lower.tail = TRUE, log.p = FALSE)
rlogis(n, location = 0, scale = 1)
```

如果函数参数 `location` 或 `scale` 没有指定，则分别取默认值 0 和 1，就是标准的逻辑斯谛分布。位置参数（类似正态分布中的均值 $\mu$）为 `location = m` ，尺度参数（类似正态分布中的标准差 $\sigma$）为 `scale = s`，逻辑斯谛分布是一个长尾分布。

### 逻辑回归 {#sec-logistic-regression}

响应变量 $Y$ 服从伯努利分布 $\mathrm{Bernoulli}(p)$，取值是 0 或 1，对线性预测 $X\boldsymbol{\beta}$ 做 Logistic 变换

$$
\bm{p} = \mathsf{E}Y = \mathrm{Logistic}(X\boldsymbol{\beta}) = \frac{1}{1 + e^{-(\alpha + X\boldsymbol{\beta})}} = \frac{e^{\alpha + X\boldsymbol{\beta}}}{1 + e^{\alpha + X\boldsymbol{\beta}}}
$$

Logistic 的逆变换

$$
\mathrm{Logistic}^{-1}(\bm{p})= \ln\big(\frac{\bm{p}}{1 - \bm{p}}\big) = \alpha + X\boldsymbol{\beta}
$$

记数据矩阵 $X$ 为

$$
X = \begin{bmatrix}
    x_{11} & x_{12} & x_{13} & \dots  & x_{1k} \\
    x_{21} & x_{22} & x_{23} & \dots  & x_{2k} \\
    \vdots & \vdots & \vdots & \ddots & \vdots \\
    x_{n1} & x_{n2} & x_{n3} & \dots  & x_{nk}
\end{bmatrix}
= 
\begin{bmatrix}
\bm{x}_1^{\top} \\
\bm{x}_2^{\top} \\
\vdots \\
\bm{x}_n^{\top}
\end{bmatrix}
$$

每一行表示一次观测，每一列表示一个变量的 $n$ 次观测，记 $X = (X_1, X_2, \cdots, X_k)$ 是一个 $n \times k$ 数据矩阵，其中 $\bm{x}_i^{\top}$ 表示矩阵 $X$ 的第 $i$ 行，一共有 $n$ 行，可以看作是 $1 \times k$ 的矩阵，$X_j, j = 1,2, \cdots, k$ 表示矩阵 $X$ 的第 $j$ 列，一共有 $k$ 列。类似地， $\boldsymbol{\beta} = (\beta_1, \beta_2, \cdots, \beta_k)^{\top}$ 是一个列向量，可以看作是 $k \times 1$ 的矩阵，$\beta_j$ 表示第 $j$ 个变量 $X_j$ 的系数。对第 $i$ 次观测

$$
\mathrm{Logistic}^{-1}(p_i)= \ln\big(\frac{p_i}{1-p_i}\big) = \alpha + \bm{x}_i^{\top}\boldsymbol{\beta}
$$

关于参数 $\alpha,\boldsymbol{\beta}$ 的似然函数如下：

$$
\begin{aligned}
\mathcal{L}(\alpha,\boldsymbol{\beta}) &= \prod_{i=1}^{n} p_i^{y_i}(1 - p_i)^{1 - y_i} \\ 
     &= \prod_{i=1}^{n} \Big(\frac{e^{\alpha + \bm{x}_i^{\top}\boldsymbol{\beta}}}{1 + e^{\alpha + \bm{x}_i^{\top}\boldsymbol{\beta}}}\Big)^{y_i}\Big(\frac{1}{e^{\alpha + \bm{x}_i^{\top}\boldsymbol{\beta}}}\Big)^{1-y_i} \\
\end{aligned}
$$ {#eq-logit-lik}

关于参数 $\alpha,\boldsymbol{\beta}$ 的对数似然函数如下：

$$
\begin{aligned}
\ell(\alpha,\boldsymbol{\beta}) &= \log \mathcal{L}(\alpha,\boldsymbol{\beta}) \\
& = \sum_{i=1}^{n} \Big[y_i \log (p_i) + (1 - y_i) \log(1-p_i)\Big] \\
&= \sum_{i=1}^{n} \Big[y_i \log \Big(\frac{e^{\alpha + \bm{x}_i\boldsymbol{\beta}}}{1 + e^{\alpha + \bm{x}_i\boldsymbol{\beta}}}\Big) + (1 - y_i) \log\Big(\frac{1}{e^{\alpha + \bm{x}_i\boldsymbol{\beta}}}\Big)\Big]
\end{aligned}
$$ {#eq-logit-lik-log}

对数似然函数 $\ell(\alpha,\boldsymbol{\beta})$ 关于参数 $\alpha,\boldsymbol{\beta}$ 的偏导数如下：

$$
\begin{aligned}
\frac{\partial \ell(\alpha,\boldsymbol{\beta})}{\partial \alpha}  &= \sum_{i=1}^{n}\Big[ \big(\frac{y_i}{p_i} -  \frac{1- y_i}{1 - p_i}\big) \frac{\partial p_i}{\partial \alpha} \Big] \\
\frac{\partial \ell(\alpha,\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} &= \sum_{i=1}^{n}\Big[\big(\frac{y_i}{p_i} -  \frac{1- y_i}{1 - p_i}\big) \frac{\partial p_i}{\partial \beta} \Big] \\
& = \sum_{i=1}^{n}\Big[\big(\frac{y_i}{p_i} -  \frac{1- y_i}{1 - p_i}\big) p_i(1- p_i) \bm{x}_i^{\top} \Big] 
\end{aligned}
$$ {#eq-logit-lik-log-partial}

其中， $p_i = \frac{e^{\alpha + \bm{x}_i\boldsymbol{\beta}}}{1 + e^{\alpha + \bm{x}_i\boldsymbol{\beta}}}$ ，要使 $\ell(\alpha,\boldsymbol{\beta})$ 取极大值，一般通过迭代加权最小二乘算法（Iteratively (Re-)Weighted Least Squares，简称 IWLS）求解此优化问题，它可以看作拟牛顿法的一种特殊情况，在 R 语言中，函数 `glm()` 是求解此类问题的办法。

## 数值优化问题求解器 {#sec-solvers}

### `optim()` {#sec-logit-optim}

从一个逻辑回归模型模拟一组样本，共 2500 条记录，即 $n = 2500$，10 个观测变量，即 $k=10$，其中，只有变量 $X_1$ 和 $X_2$ 的系数非零，参数设定为 $\alpha = 1, \beta_1 = 3,\beta_2 = -2$，而 $\beta_i = 0, i=3, \cdots, 10$ 模拟数据的代码如下：

```{r}
set.seed(2023)
n <- 2500
k <- 10
X <- matrix(rnorm(n * k), ncol = k)
y <- rbinom(n, size = 1, prob = plogis(1 + 3 * X[, 1] - 2 * X[, 2]))
```

模拟数据矩阵 X 与上述记号 $X$ 是对应的，记号 $\bm{x_i}^{\top}$ 表示数据矩阵的第 $i$ 行。$\alpha$ 是逻辑回归方程的截距，$\bm{\beta}$ 是 $k$ 维列向量，$X$ 是 $n \times k$ 维的矩阵且 $n > k$，$y$ 是 $n$ 维向量。极大化对数似然函数 @eq-logit-lik-log ，就是求解一个多维非线性无约束优化问题。方便起见，将 $\alpha$ 合并进 $\bm{\beta}$ 向量，另，函数 `optim()` 默认求极小，因此在对数似然函数前添加负号。

```{r}
# 目标函数
log_logit_lik <- function(beta) {
  p <- plogis(cbind(1, X) %*% beta)
  -sum(y * log(p) + (1 - y) * log(1 - p))
}
```

高维情形下，没法绘制似然函数图形，退化到二维，如 @fig-log-logit-lik 所示，二维情形下的逻辑回归模型的负对数似然函数曲面。

```{r}
#| label: fig-log-logit-lik
#| fig-cap: "二维情形下的逻辑回归模型的负对数似然函数曲面"
#| fig-width: 6
#| fig-height: 5
#| fig-showtext: true
#| echo: false

# 目标函数
log_logit_lik0 <- function(beta, X0 = X0, y0 = y0) {
  p <- plogis(cbind(1, X0) %*% beta)
  -sum(y0 * log(p) + (1 - y0) * log(1 - p))
}

set.seed(2023)
n <- 2500
k0 <- 1
X0 <- matrix(rnorm(n * k0), ncol = k0)
y0 <- rbinom(n, size = 1, prob = plogis(1 + 3 * X0[, 1]))

alpha <- seq(0, 2, length.out = 20)
beta <- seq(2, 4, length.out = 20)

df <- expand.grid(x = alpha, y = beta)
df$fnxy <- apply(X = df, MARGIN = 1, FUN = log_logit_lik0, X0 = X0, y0 = y0)
library(lattice)
custom_palette <- function(irr, ref, height, saturation = 0.9) {
  hsv(
    h = height, s = 1 - saturation * (1 - (1 - ref)^0.5),
    v = irr
  )
}
wireframe(
  data = df, fnxy ~ x * y,
  shade = TRUE, drape = FALSE,
  xlab = expression(alpha),
  ylab = expression(beta),
  zlab = list(expression(-loglik(alpha, beta)), rot = 90),
  shade.colors.palette = custom_palette,
  # 减少三维图形的边空
  lattice.options = list(
    layout.widths = list(
      left.padding = list(x = -.6, units = "inches"),
      right.padding = list(x = -1.0, units = "inches")
    ),
    layout.heights = list(
      bottom.padding = list(x = -.8, units = "inches"),
      top.padding = list(x = -1.0, units = "inches")
    )
  ),
  scales = list(arrows = FALSE, col = "black"),
  # 设置坐标轴字体大小
  par.settings = list(
    axis.line = list(col = "transparent")
  ),
  screen = list(z = -30, x = -70, y = 0)
)
```

当用 Base R 函数 `optim()` 来求解时，发现 Nelder-Mead 算法收敛慢，易陷入局部最优解，即使迭代 10000 次，与真值仍然相去甚远。当用 SANN （模拟退火算法）求解此 11 维非线性无约束优化问题时，迭代 10000 次后，比较接近真值。

```{r}
optim(
  par = rep(1, 11), # 初始值
  fn = log_logit_lik, # 目标函数
  method = "SANN",
  control = list(maxit = 10000)
)
```

根据目标函数计算其梯度，有了梯度信息，可以使用迭代效率更高的 L-BFGS-B 算法。

```{r}
# 梯度函数
log_logit_lik_grad <- function(beta) {
  p <- plogis(cbind(1, X) %*% beta)
  -t((y / p - (1 - y) / (1 - p)) * p * (1 - p)) %*% cbind(1, X)
}

optim(
  par = rep(1, 11), # 初始值
  fn = log_logit_lik, # 目标函数
  gr = log_logit_lik_grad, # 目标函数的梯度
  method = "L-BFGS-B"
)
```

相比于函数 `optim()`，R 包 **nloptr** 不但可以提供类似的数值优化功能，而且可以处理各类非线性约束，能力更强。仍然基于上面的优化问题， 调用 **nloptr** 包求解的代码如下：

```{r}
library(nloptr)
nlp <- nloptr(
  x0 = rep(1, 11),
  eval_f = log_logit_lik,
  eval_grad_f = log_logit_lik_grad,
  opts = list(
    "algorithm" = "NLOPT_LD_LBFGS",
    "xtol_rel" = 1.0e-8
  )
)
nlp
```

如果对数似然函数是多模态的，一般的求解器容易陷入局部最优解，推荐用 **nloptr** 包的[全局优化求解器](https://nlopt.readthedocs.io/en/latest/NLopt_Algorithms/#global-optimization)。

```{r}
#| eval: false
#| echo: false
# 加载 ROI 时不要自动加载插件
Sys.setenv(ROI_LOAD_PLUGINS = FALSE)
library(ROI)
library(ROI.plugin.nloptr)
op <- OP(
  objective = F_objective(F = log_logit_lik, n = 11L),
  types = rep("C", 11), maximum = FALSE,
  bounds = V_bound(ld = -3, ud = 3, nobj = 11L)
)
nlp <- ROI_solve(op, solver = "nloptr.directL")
nlp$solution
# L-BFGS
op <- OP(
  objective = F_objective(F = log_logit_lik, n = 11L, G = log_logit_lik_grad),
  bounds = V_bound(ld = -3, ud = 3, nobj = 11L)
)
nlp <- ROI_solve(op, solver = "nloptr.lbfgs", start = rep(1, 11))
nlp
```

### `glm()` {#sec-logit-glm}

Base R 提供的函数 `glm()` 拟合模型，指定联系函数为 logit 变换。

```{r}
fit_r <- glm(y ~ X, family = binomial(link = "logit"))
summary(fit_r)
```

或者也可以用函数 `glm.fit()`，效果类似，使用方式不同罢了。

```{r}
fit_r2 <- glm.fit(x = cbind(1, X), y = y, family = binomial(link = "logit"))
coef(fit_r2)
```

函数 `glm()` 的参数是一个公式，函数 `glm.fit()` 的参数是矩阵、向量，用函数 `glm()` 拟合模型，其内部调用的就是函数 `glm.fit()`。

### glmnet 包 {#sec-logit-glmnet}

调用 **glmnet** 包的函数 `glmnet()` 拟合模型，指定指数族的具体形式为二项分布，伯努利分布是二项分布的特殊形式，也叫两点分布或0-1分布。

```{r}
#| message: false

library(Matrix)
library(glmnet)
fit_glm <- glmnet(x = X, y = y, family = "binomial")
```

逻辑回归模型系数在 L1 正则下的迭代路径图

```{r}
#| label: fig-logit-glmnet
#| fig-cap: "回归系数的迭代路径"
#| fig-width: 5
#| fig-height: 5
#| fig-showtext: true

plot(fit_glm, ylab = "回归系数")
```

从图可见，剩余两个系数是非零的，一个是 3， 一个是 -2，其余都被压缩，而接近为 0 了。

```{r}
#| label: fig-logit-glmnet-lambda
#| fig-cap: "惩罚系数的迭代路径"
#| fig-width: 5
#| fig-height: 5
#| fig-showtext: true

plot(fit_glm$lambda,
  ylab = expression(lambda), xlab = "迭代次数",
  main = "惩罚系数的迭代路径"
)
```

随着迭代的进行，惩罚系数 $\lambda$ 越来越小，接近于 0，这也是符合预期的，因为模型本来就是简单的逻辑回归，不带惩罚项。选择一个迭代趋于稳定时的 $\lambda$ 比如 0.0005247159，此时各个参数的取值如下：

```{r}
coef(fit_glm, s = 0.0005247159)
```

截距 (Intercept) 对应 $\alpha = 0.997741857$，而 $\beta_1 = 3.076358149$ 对应 V1，$\beta_2 = -1.984018387$ 对应 V2，以此类推。

## 评估模型的分类效果 {#sec-evaluation-model-performance}

逻辑回归模型是二分类模型，评估模型的分类效果，两个办法。

1.  可以用 AUC 指标或者 ROC 曲线，**pROC** 包和 **ROCR** 包都可以绘制 ROC 曲线。
2.  可以用 Wilcoxon 检验，越显著表示分类效果越好。

### ROC 曲线和 AUC 值

ROC 是 Receiver Operating Characteristic 简写。随机抽取 2000 个样本作为训练集，余下的数据作为测试集。

```{r}
dat <- cbind.data.frame(X, y)
set.seed(20232023)
idx <- sample(x = 1:nrow(dat), size = 2000, replace = F)
# 训练集
dat_train <- dat[idx, ]
# 测试集
dat_test <- dat[-idx, ]
```

函数 `glm()` 拟合训练集数据

```{r}
fit_binom <- glm(y ~ ., data = dat_train, family = binomial(link = "logit"))
```

将训练好的模型用于测试集，调用函数 `predict()` 进行预测，`type = "response"` 获得预测概率值，它是对数几率，比值比的对数。

```{r}
dat_test$pred <- predict(fit_binom, newdata = dat_test, type = "response")
```

返回值介于 0 - 1 之间，表示预测概率。在测试集上绘制 ROC 曲线。

```{r}
#| label: fig-logit-roc
#| fig-cap: ROC 曲线
#| fig-showtext: true
#| fig-width: 4
#| fig-height: 4

pROC::plot.roc(
  y ~ pred, data = dat_test,
  col = "dodgerblue", print.auc = TRUE,
  auc.polygon = TRUE, auc.polygon.col = "#f6f6f6",
  xlab = "FPR", ylab = "TPR", main = "预测 ROC 曲线"
)
```

ROC 曲线越往左上角拱，表示预测效果越好。FPR 是 False Positive Rate 的缩写，TPR 是 True Positive Rate 的缩写。

```{r}
# 计算 AUC 值
pROC::auc(y ~ pred, data = dat_test)
```

AUC 是 area under curve 的缩写，表示 ROC 曲线下的面积，所以 AUC 指标越接近 1 越好。

### Wilcoxon 检验

对每个标签的预测概率指定服从均匀分布，相当于随机猜测，所以最后 ROC 会接近对角线，而且样本量越大越接近，AUC 会越来越接近 0.5。如果预测结果比随机猜测要好，Wilcoxon 检验会显著，预测效果越好检验会越显著，表示预测 pred 和观测 y 越接近。

```{r}
wilcox.test(pred ~ y, data = dat_test)
```

```{r}
#| eval: false
#| echo: false
# 计算 auc 的函数
# dat is a data.frame as input return AUC value
comp_auc <- function(dat, show_roc = TRUE) {
  # order label by predicted probability
  dat <- dat[order(dat$pred, dat$label, decreasing = TRUE), ]

  # total samples
  n_total <- length(dat$label)

  # number of positive label 1
  n_pos <- sum(dat$label)

  # number of negative label 0
  n_neg <- n_total - n_pos

  # calculate TPR and FPR
  tpr <- cumsum(dat$label) / n_pos
  fpr <- (1:n_total - cumsum(dat$label)) / n_neg

  # calculate auc
  auc <- 0
  for (i in 1:(n_total - 1)) {
    auc <- auc + (fpr[i + 1] - fpr[i]) * tpr[i]
  }
  # show ROC curve or not?
  if (show_roc) {
    plot(fpr, tpr, type = "l")
  }
  auc
} 
sim_dat = cbind.data.frame(pred = dat_test$pred, label =  dat_test$y) 
comp_auc(dat = sim_dat, show_roc = FALSE)
```