R/python-回归分析

最近用R和python做了一些回归分析的项目，发现这两个工具各有特点，如果目的是分析研究用R即可，毕竟R语言有很多方便的函数；如果目的是为了建模预测，用python的sklearn（集成了很多机器学习算法和模型检验函数），工具只是手段，我们对结果负责就好。

帖子目录（也是回归分析步骤）：
1、作图，观察变量分布，均值和中位数等统计学指标，宏观了解数据
2、使用全子集回归选择变量，生成模型，评估模型效果
3、回归诊断
4、对异常值，强影响点和高杠杆指进行处理，重新生成模型，再次评估模型效果
5、改进措施(变量类型转换，删除变量，因子分析）
6、非线性回归问题（python:logistic ,ridge,lasso,弹性网络 ） 测试集和训练集，模型评估 K折交叉验证等等）

注：本例数据集来源于R语言基础包中的state.x77数据集，我们想探究一个州的犯罪率和其他因素的关系，自变量为人口、文盲率、平均收入和结霜天数（温度在冰点以下的平均天数），因变量为谋杀率

多元线性回归步骤：
一、画图观察变量分布以及变量间关系

library(car)
states <- as.data.frame(state.x77[,c("Murder", "Population",
                                     "Illiteracy", "Income", "Frost")])
#绘制变量间的散点图，并添加平滑和线性拟合曲线
#对角线区域绘制每个变量的密度图和轴须图
scatterplotMatrix(states, spread=FALSE, smoother.args=list(lty=2),
                  main="Scatter Plot Matrix")
#cor()函数提供了二变量之间的相关系数， 
cor(states)

R-bianliang.png

二、拟合回归曲线
1、全子集回归

library(leaps)
leaps <-regsubsets(Murder ~ Population + Illiteracy + Income +
                     Frost, data=states, nbest=4)
plot(leaps, scale="adjr2")

R-quan.png

从图中可看出，双预测变量模型（ Population和Illiteracy）是最佳模型（R^2最高）

2、拟合回归曲线
通过全子集回归，可以看出双预测变量模型为最佳模型，但是为了对模型概要（summary(fit)）参数进行解释，同时为了和《R语言实战》回归部分保持一致，仍然用Population + Illiteracy + Income + Frost四个变量拟合曲线
注：如果有二次项I(Income^2)

fit <- lm(Murder ~ Population + Illiteracy + Income + Frost,
          data=states)
summary(fit)

运行结果如下：

summary(fit).png

结果解释：
1）文盲率的回归系数为4.14，表示控制人口、收入和温度不变时，文盲率上升1%，谋杀率将会上升4.14%，
2）文盲率的回归系数在p<0.001的水平下显著不为0（p=0.000022）。
3）Frost的系数没有显著不为0（ p=0.954），表明当控制其他变量不变时， Frost与Murder不呈线性相关。
4）总体来看，所有的预测变量解释了各州谋杀率57%(Multiple R-squared=0.567）的方差。
5）F的值是回归方程的显著性检验，表示的是模型中被解释变量与所有解释变量之间的线性关系在总体上是否显著做出推断

三、回归诊断
1、综合性检验

library(gvlma)
gvmodel <- gvlma(fit)
summary(gvmodel)

运行结果如下：

gvmodel.jpg

结果解释：
1）从输出项（ Global Stat中的文字栏）我们可以看到数据满足OLS回归模型所有的统计假设（ p=0.597）。
2）若Decision下的文字表明违反了假设条件（比如p<0.05），你可以使用下面的方法（正态性，误差独立性，线性和同方差性）来判断哪些假设没有被满足

2、正态性：当预测变量值固定时，因变量成正态分布，则残差值也应该是一个均值为0的正态分布
1)QQ图

qqPlot(fit, labels=row.names(states), id.method="identify",
       simulate=TRUE, main="Q-Q Plot")
#除输出图形之外，还会输出利群点位置

QQ.png

2）、学生化残差直方图

residplot <- function(fit, nbreaks=10) {
  z <- rstudent(fit)
  hist(z, breaks=nbreaks, freq=FALSE,
       xlab="Studentized Residual",
       main="Distribution of Errors")
  rug(jitter(z), col="brown")
  curve(dnorm(x, mean=mean(z), sd=sd(z)),
        add=TRUE, col="blue", lwd=2)
  lines(density(z)$x, density(z)$y,
        col="red", lwd=2, lty=2)
  legend("topright",
         legend = c( "Normal Curve", "Kernel Density Curve"),
         lty=1:2, col=c("blue","red"), cex=.7)
}
residplot(fit)

Studentized Residual.png

3、误差的独立性：判断因变量值（或残差）是否相互独立

#该检验适用于时间独立的数据
#结果p值不显著说明无自相关性，误差项之间独立
durbinWatsonTest(fit)

durbinWatsonTest.png

结果显示：
1）p值不显著（ p=0.28）说明无自相关性，误差项之间独立。
2）滞后项（ lag=1）表明数据集中每个数据都是与其后一个数据进行比较的。

4、线性：若因变量与自变量线性相关，那么残差值与预测值就没有任何系统关联
绘制变量的成分残差图，如果两条线（直线和拟合曲线）接近，说明可以进行线性拟合，即模型是线性的

crPlots(fit)

component plus residual plot.png

5、同方差性：指总体回归函数中的随机误差项在解释变量条件下具有不变的方差
若p值显著，则说明存在异方差性

ncvTest(fit)

ncvTest.png

结果显示：
a）计分检验不显著（ p=0.2），说明满足方差不变假设

6、多重共线性
1）与统计假设没有直接关联，但是对于解释多元回归的结果非常重要
2）多重共线性可用统计量VIF（ Variance Inflation Factor， 方差膨胀因子）进行检测，一般原则下， vif >2就表明存在多重共线性问题（自变量之间相关性很大）

vif(fit)
sqrt(vif(fit)) > 2

vif.png

结果解释：
a）vif值都小于2，不存在共线性

3. 解决方案：

四、异常值，高杠杆指，强影响点观测
0、一图绘制异常值，高杠杆指，强影响点

par(mfrow=c(1,1))
influencePlot(fit, id.method="identify", main="Influence Plot",
              sub="Circle size is proportional to Cook's distance")

yichang.png

结果显示：
1）上图反映出Nevada和Rhode Island是离群点（纵坐标超过+2或小于–2）
2）New York、 California、 Hawaii和Washington有高杠杆值（水平轴超过0.2或0.3）
3）Nevada、 Alaska和Hawaii为强影响点（圆圈大的为强影响点，Hawaii图中未显示）
1、离群点：
1）QQ图
2）标准化残差值大于2或者小于–2的点可能是离群点

2、高杠杆值点：是一个异常的预测变量值的组合
1）高杠杆值的观测点可通过帽子统计量（ hat statistic）判断。
2）对于一个给定的数据集，帽子均值为p/n，其中p是模型估计的参数数目（包含截距项）， n是样本量。一般来说，若观测点的帽子值大于帽子均值的2或3倍，就可以认定为高杠杆值点

hat.plot <- function(fit) {
  p <- length(coefficients(fit))
  n <- length(fitted(fit))
  plot(hatvalues(fit), main="Index Plot of Hat Values")
  abline(h=c(2,3)*p/n, col="red", lty=2)
  identify(1:n, hatvalues(fit), names(hatvalues(fit)))
}
hat.plot(fit)

hat.png

图形解释：
1)水平线标注的即帽子均值2倍和3倍的位置
2）此图中，可以看到Alaska和California非常异常，查看它们的预测变量值，与其他48个州进行比较发现： Alaska收入比其他州高得多，而人口和温度却很低； California人口比其他州府多得多，但收入和温度也很高。
3）高杠杆值点可能是强影响点，也可能不是，这要看它们是否是离群点。

3强影响点：它对模型参数的估计产生的影响过大，非常不成比例
有两种方法可以检测强影响点：
1）Cook距离：Cook’s D值大于4/(n–k–1)，则表明它是强影响点，其中n为样本量大小， k是预测变量数目

cutoff <- 4/(nrow(states)-length(fit$coefficients)-2)
par(mfrow=c(1,1))
plot(fit, which=4, cook.levels=cutoff)
abline(h=cutoff, lty=2, col="red")

cutoff.png

结果解释：
a)通过图形可以判断Alaska、 Hawaii和Nevada是强影响点。若删除这些点，将会导致回归模型截距项和斜率发生显著变化。
b)注意，虽然该图对搜寻强影响点很有用，但我逐渐发现以1为分割点比4/(n–k–1)更具一般性

2)变量添加图：对于每个预测变量Xk，绘制Xk在其他k–1个预测变量上回归的残差值相对于响应变量在其他k–1个预测变量上回归的残差值的关系图
注：Cook’s D图有助于鉴别强影响点，但是并不提供关于这些点如何影响模型的信息。变量添加图弥补了这个缺陷。对于一个响应变量和k个预测变量，你可以如下图创建k个变量添加图

avPlots.png

结果解释：

五、改进措施
1、当模型不符合正态性、线性或者同方差性假设时，一个或多个变量的变换通常可以改善或调整模型效果
1）当模型违反正态假设时，通常可以对响应变量尝试某种变换

summary(powerTransform(states$Murder))

zhengtai.png

结果表明:
a)可以用Murder^0.6来正态化变量Murder。由于0.61很接近0.5，你可以尝试用平方根变换来提高模型正态性的符合程度。但在本例中， λ=1的假设也无法拒绝（ p=0.15），因此没有强有力的证据表明本例需要变量变换

2)当违反了线性假设时，对预测变量进行变换常常会比较有用

boxTidwell(Murder~Population+Illiteracy,data=states)

xianxing.png
结果显示：
a)使用变换Population^{0.87和Illiteracy}1.36能够大大改善线性关系。但是对
Population（ p=0.7）和Illiteracy（ p=0.5）的计分检验又表明变量并不需要变换

3)响应变量变换还能改善异方差性

ncvTest(fit)
spreadLevelPlot(fit)

fangcha.png

结果显示：
a)计分检验不显著（ p=0.19），说明满足方差不变假设
b)经过p次幂（ Y^p）（ suggested power transformation）变换，非恒定的误差方差将会平稳,本例建议的p=1.21，约等于1

2、增删变量
1）多重共线性时，删去vif^0.5<2的变量

3、因子分析
因为
1）因子分析通过寻找一组更小的、潜在的或隐藏的结构来解释已观测到的、显式的变量间的关系。
2）分析步骤：
a)数据预处理
可以输入原始数据矩阵或者相关系数矩阵到principal()和fa()函数中。
若输入初始数据，相关系数矩阵将会被自动计算，
在计算前请确保数据中没有缺失值

library(psych)
covariances <- ability.cov$cov
correlations <- cov2cor(covariances)

b)判断要选择的因子数目:
判断规则：
i) Kaiser-Harris准则建议保留特征值大于1的主成分（PCA：y=1直线上部分因子，EFA:y=0直线上部分因子)
ii) 若基于真实数据的某个特征值大于一组随机数据矩阵相应的平均特征值，
那么该主成分可以保留。该方法称作平行分析（虚线上部分因子）

#fa="fa",fa="pc",fa="both"
#码展示了基于观测特征值的碎石检验、
#根据100个随机数据矩阵推导出来的特征值均值（虚线），
#以及大于1的特征值准则（ y=1的水平线）
fa.parallel(correlations, n.obs=112, fa="fa", n.iter=100,
              main="Scree plots with parallel analysis") 

EFA.png

结果解释：
i)虚线上部分，y=0直线的上部分都表明应该保留两个因子

c)提取公因子
提取公共因子的方法很多，包括最大似然法（ ml）、主轴迭代法（ pa）、
加权最小二乘法（ wls）、广义加权最小二乘法（ gls）和最小残差法（ minres）。

#nfactors=2,提取两个因子
#scores设定是否计算因子得分（默认不计算
#fm设定因子化方法（默认极小残差法）
#主轴迭代法提取未旋转的因子fm="pa"
fa <- fa(correlations, nfactors=2, rotate="none", fm="pa") 

fa.png

结果解释：
i)PA1，PA2栏（提取的公因子）包含了成分载荷(观测变量与主成分的相关系数)
ii)h2栏指成分公因子方差，即主成分对每个变量的方差解释度。
iii)u2栏指成分唯一性，即方差无法被主成分解释的比例（ 1–h2）。
iv)SS loadings行包含了与主成分相关联的特征值，指的是与特定主成分相关联的标准化后的方差值（2.75,0.83）
v)Proportion Var行表示的是每个主成分对整个数据集的解释程度，可以看到，两个因子解释了六个心理学测验60%的方差。

d)旋转因子：得到因子与变量相关系数
1）用正交旋转提取因子:人为地强制两个因子不相关
fa.varimax <- fa(correlations, nfactors=2, rotate="varimax", fm="pa")

fa.varimax.jpg

i)结果显示阅读(reading)和词汇(vocab)在第一因子上载荷较大，画图(picture)、积木图案(blocks)和迷宫(maze)在第二因子上载荷较大，
ii)非语言的普通智力测量(general)在两个因子上载荷较为平均，这表明存在一个语言智力因子和一个非语言智力因子。

2）斜交转轴法提取因子:允许两个因子相关（但模型将更符合真实数据）

#而对于斜交旋转，因子分析会考虑三个矩阵：
#1）因子结构矩阵：变量与因子的相关系数
#2）因子模式矩阵：标准化的回归系数矩阵。它列出了因子预测变量的权重
#3）因子关联矩阵：即因子相关系数矩阵。

library(GPArotation)
fa.promax <- fa(correlations, nfactors=2, rotate="promax", fm="pa")

fa.promax.png

结果解释：
i)PA1和PA2栏中的值组成了因子模式矩阵。它们是标准化的回归系数，而不是相关系数
ii)因子关联矩阵显示两个因子的相关系数为0.55，相关性很大。如果因子间的关联性很低，你可能需要重新使用正交旋转来简化问题
iii)但你可以使用公式F = P*Phi得到因子结构矩阵（或称因子载荷阵），其中F是因子载荷阵， P为因子模式矩阵， Phi为因子关联矩阵

fsm <- function(oblique) {
  if (class(oblique)[2]=="fa" & is.null(oblique$Phi)) {
    warning("Object doesn't look like oblique EFA")
  } else {
    P <- unclass(oblique$loading)
    F <- P %*% oblique$Phi
    colnames(F) <- c("PA1", "PA2")
    return(F)
  }
}
fsm(fa.promax)

fsm(fa.promax).png

结果解释：
i)此图显示了变量与因子间的相关系数。与正交旋转所得因子载荷阵相比，该载荷阵列的噪音比较大（因为允许潜在因子相关）。
ii)虽然斜交方法更为复杂，但模型将更符合真实数据

e）作图显示正交或斜交结果的图形

factor.plot(fa.promax, labels=rownames(fa.promax$loadings))

fa.promax-plot.png

结果解释：
i)词汇（vocab）和阅读(read)在第一个因子（ PA1）上载荷较大，
ii)而积木图案、画图和迷宫在第二个因子（ PA2）上载荷较大。
iii)普通智力测验在两个因子上较为平均

可视化显示相关系数：

fa.diagram(fa.promax, simple=FALSE)

xishu.png