回归分析篇

非线性回归

在许多实际问题中，变量之间的关系并不都是线性的，通常我们会碰到某些现象的被解释变量与解释变量之间呈现某种曲线关系。对于曲线形式的回归问题，我们就不能照搬之前所说的多元线性回归的建模方法了

通常我们是对自变量或因变量进行函数变换，让变换后的变量能应用多元线性回归方法来建模

---

可化为线性回归的曲线方程

实际问题中，有许多回归模型的被解释变量 $y$与解释变量 $x$之间的关系都不是线性的，其中一些回归模型可以用函数变换来化为线性关系，而有些却不能。一般我们得想办法利用函数变换和变量替换来将原方程转换为线性方程表示.

---

比如：

$$y = \beta_0 + \beta_1 e^x + \epsilon$$ 只需令 $x' = e^x$即可化为 $y = \beta_0 + \beta_1 x' + \epsilon$

---

$$y = \beta_0 + \beta_1 x + \beta_2 x^2 + \cdots + \beta_p x^p + \epsilon$$ 可以令 $x_1 = x, x_2 = x^2, x_3 = x^3,\cdots,x_p = x^p$，于是可得 $y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \cdots + \beta_p x_p + \epsilon$

---

$$y = a e^{bx} e^{\epsilon}$$ 可以两边取自然对数得 $ln(y) = ln(a) + bx + \epsilon$，于是令 $y' = ln(y), \beta_0 = ln(a), \beta_1 = b$，所以原方程变换成了 $y' = \beta_0 + \beta_1 x + \epsilon$

---

上述方程都是可以转换为线性回归的那种情况，但是像方程

$$y = a e^{bx} + \epsilon$$ 就无法转换为线性，这种就只能用非线性最小二乘法来求解了.

---

上述最后两个回归模型有相同的回归函数 $a e^{bx}$，只是误差项 $\epsilon$的形式不一样，$a e^{bx} e^{\epsilon}$的误差项被称为乘性误差项，$a e^{bx} + \epsilon$的误差项被称为加性误差项。因而一个非线性回归模型能否线性化，不仅与回归函数的形式有个，还与误差项的形式有关，当然，误差项的形式还可以有其他多种形式.

---

非线性最小二乘

非线性回归模型一般可表示为

$$y_i = f(\boldsymbol{x}_i,\boldsymbol{\theta}) + \epsilon_i,\ \ i = 1,2,\cdots,n$$ 其中，$y_i$为因变量；非随机变量 $\boldsymbol{x}_i = (x_{i1},x_{i2},\cdots,x_{ik})'$ 是自变量；$\boldsymbol{\theta} = (\theta_0,\theta_1,\cdots,\theta_p)'$ 为未知参数向量；$\epsilon_i$ 为随机误差项并且满足独立同分布假设，即 $$\begin{cases} E(\boldsymbol{\epsilon}_i) = 0 \\ cov(\epsilon_i,\epsilon_j) = \begin{cases} \sigma^2, i = j \\ 0, i \neq j \end{cases} \end{cases}$$ 其中，$i,j = 1,2,\cdots,n$

如果 $f(\boldsymbol{x}_i,\boldsymbol{\theta}) = \theta_0 + x_1 \theta_1 + x_2 \theta_2 + \cdots + x_p \theta_p$，那么就是前面我们讲的多元线性回归模型，而且必然有 $k = p$；对于一般情况的非线性模型，参数的数目与自变量的数目并没有一定的对应关系，不要求 $k = p$

---

对于非线性回归模型，我们任使用最小二乘法估计参数 $\boldsymbol{theta}$，即求使得

$$Q(\boldsymbol{\theta}) = \sum_{i=1}^n{(y_i - f(x_i,\boldsymbol{\theta}))^2}$$ 达到最小的 $\hat{\boldsymbol{\theta}}$，称为 $\boldsymbol{\theta}$的非线性最小二乘估计。在假定 $f$函数对参数 $\boldsymbol{\theta}$连续可微时，可以利用微分法，建立正规方程组，求解使 $Q(\boldsymbol{\theta})$达最小的 $\hat{\boldsymbol{\theta}}$。将 $Q$函数对参数 $\boldsymbol{\theta}_j$求偏导，并令其为 $0$，得 $p+1$个方程 $$\left.\frac{\partial Q}{\partial \theta_j}\right|_{\theta_j = \hat{\theta}_j} = -2 \sum_{i=1}^n{(y_i - f(x_i,\hat{\theta}))\left.\frac{\partial f}{\partial \theta_j}\right|_{\theta_j = \hat{\theta}_j}} = 0$$ 其中，$j = 0,1,2,\cdots,p$， 非线性最小二乘估计 $\hat{\boldsymbol{\theta}}$就是该微分方程的解，称为非线性最小二乘估计的正规方程组，它是未知参数的非线性方程组。一般用Newton迭代法求解此正规方程组.

也可以直接极小化残差平方和 $Q(\boldsymbol{\theta})$，求出未知参数 $\boldsymbol{\theta}$的非线性最小二乘估计 $\hat{\boldsymbol{\theta}}$

---

在非线性回归中，平方和分解式 $SST = SSR + SSE$不再成立，类似于线性回归中的复决定系数，定义非线性回归的相关系数：

$$R^2 = 1 - \frac{SSE}{SST}$$

---

两个重要的非线性回归模型

非线性回归模型比较多，这里介绍两个比较重要的非线性回归模型，即分段回归和逻辑回归；而对于像多项式回归等比较简单的我们这里就不作细致的讲解，一般都是用函数变换和变量替换就可以转换为线性回归模型了

分段回归

在实际问题中，我们会碰到某些变量在不同的影响因素范围内变化趋势截然不同，例如经济问题涉及到经济政策有较大调整时，调整前与调整后的变化幅度会有很大不同。对于这种问题，我们有时用单一的一种曲线拟合效果任不能令人满意。

如果做残差分析，会发现残差不是随机的，而具有一定的系统性。对于这样一类问题，人们自然考虑到利用分段回归的方法来做处理.

例如，在观察样本散点图时如果发现，自变量在一段区域内符合一种模型，而在另一段区域时模型发生很大改变，又符合另一种模型，这个时候就可以利用分段回归了，在分段点的左右两边分别用不同的曲线去拟合；而这个问题的一个重点是如何找到分段点，因为有时候分段点的界限不是那么明显。

分段点确定的那种情况比较简单，我们这里就不进行说明，当分段点事先不知道时，我们只能通过样本数据本身来识别.

---

设样本满足分段回归方程

$$f(x) = \begin{cases} g(x), x \leq k \\ h(x), x > k \end{cases}$$ 而分段点 $k$ 我们事先不知道，于是我们假设分段点为 $k'$，然后分别对 $x \leq k' 和 x > k'$ 进行拟合，则总体误差函数为 $$Q(\boldsymbol{x}) = \sum_{i=1}^k'{(g'(x_i) - g(x_i))^2} + \sum_{i=k'+1}^n{(h'(x_i) - h(x_i))^2}$$ 寻找参数 $k'$使得总体误差函数 $Q(\boldsymbol{x})$最小，可以想象，若 $k'$与$k$越接近，则分段拟合函数的总体误差值就越小，就可以用 $k'$ 来估计 $k$

在后面的R程序演示中，笔者会用一个实例来介绍如何用改进逼近算法来寻找最佳的分段点.

---

逻辑回归 (Logistic回归模型)

设 $y$是 $0~1$ 型变量，$x_1,x_2,\cdots,x_p$ 是与 $y$相关的确定性变量，$n$组观测数据为 $(x_{i1},x_{i2},\cdots,x_{ip};y_i)(i=1,2,\cdots,n)$，其中，$y_1,y_2,\cdots,y_n$是取值$0$或$1$的随机变量，$y_i$与$x_{i1},x_{i2},\cdots,x_{ip}$的关系如下：

$$E(y_i) = \pi_i = f(\beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \cdots + \beta_p x_{ip})$$ 其中，函数 $f(x)$ 是值域在 $[0,1]$区间内的单调递增函数.

对于 Logistic回归

$$f(x) = \frac{e^x}{1 + e^x}$$ 于是 $y_i$ 是均值为 $\pi_i = f(\beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \cdots + \beta_p x_{ip})$的 0-1型分布，概率函数为 $$P(y_i = 1) = \pi_i \\ P(y_i = 0) = 1 - \pi_i$$ 可以把 $y_i$的概率函数合写为 $$P(y_i) = \pi_j^{y_i} (1 - \pi_i)^{1-y_i} , \ \ y_i = 0,1;i = 1,2,\cdots,n$$ 于是，$y_1,y_2,\cdots,y_n$的似然函数为 $$L = \prod_{i=1}^n{P(yi)} = \prod_{i=1}^n{\pi_j^{y_i} (1 - \pi_i)^{1-y_i}}$$ 对似然函数取自然对数，得 $$lnL = \sum_{i=1}^n{[y_i ln{\pi_i} + (1-y_i)ln(1-\pi_i)]} \\ = \sum_{i=1}^n{[y_i ln{\frac{\pi_i}{(1 - \pi_i)}} + ln{(1 - \pi_i)}]}$$ 对于Logistic回归，将 $$\pi_i = \frac{exp(\beta_0 + \beta_1 x_{i1} + \cdots + \beta_p x_{ip})}{1 + exp(\beta_0 + \beta_1 x_{i1} + \cdots + \beta_p x_{ip})}$$ 带入得 $$lnL = \sum_{i=1}^n{[y_i (\beta_0 + \beta_1 x_{i1} + \cdots + \beta_p x_{ip}) - ln(1 + exp(\beta_0 + \beta_1 x_{i1} + \cdots + \beta_p x_{ip}))]}$$ 极大似然估计就是选取 $\beta_0,\beta_1,\beta_2,\cdots,\beta_p$ 的估计值 $\hat{\beta}_0,\hat{\beta}_1,\hat{\beta}_2,\cdots,\hat{\beta}_p$，使上式达到极大。而求解过程需要用到数值计算，后面我们会用R程序演示.

---

R程序演示