20180404（从有道迁移）
回归
-
回归的多面性
回归分析的各种变体
回归类型 用 途 简单线性 用一个量化的解释变量预测一个量化的响应变量 多项式 用一个量化的解释变量预测一个量化的响应变量，模型的关系是n阶多项式 多元线性 用两个或多个量化的解释变量预测一个量化的响应变量 多变量 用一个或多个解释变量预测多个响应变量 Logistic 用一个或多个解释变量预测一个类别型响应变量 泊松 用一个或多个解释变量预测一个代表频数的响应变量 Cox比例风险 用一个或多个解释变量预测一个事件（死亡、失败或旧病复发）发生的时间 时间序列 对误差项相关的时间序列数据建模 非线性 用一个或多个量化的解释变量预测一个量化的响应变量，不过模型是非线性的 非参数 用一个或多个量化的解释变量预测一个量化的响应变量，模型的形式源自数据形式，不事先设定 稳健 用一个或多个量化的解释变量预测一个量化的响应变量，能抵御强影响点的干扰 -
OLS 回归
-
普通最小二乘（OLS）回归法的适用情境
OLS回归是通过预测变量的加权和来预测量化的因变量，其中权重是通过数据估计而得的参数
-
OLS回归拟合模型的形式,其中，n 为观测的数目，k 为预测变量的数目。
\widehat{y}_{i}=\widehat{\beta}_{0}+\widehat{\beta}_{1}X_{1i}+...+\widehat{\beta}_{k}X_{ki} i=1...nimage 为了能够恰当地解释OLS模型的系数，数据必须满足以下统计假设。
正态性 对于固定的自变量值，因变量值成正态分布。
独立性 Yi值之间相互独立。
线性 因变量与自变量之间为线性相关。
-
同方差性 因变量的方差不随自变量的水平不同而变化。也可称作不变方差，但是说同方差性感觉上更犀利。
如果违背了以上假设，你的统计显著性检验结果和所得的置信区间很可能就不精确。注意，OLS回归还假定自变量是固定的且测量无误差，但在实践中通常都放松了这个假设
-
用lm()拟合回归模型
-
在R中，拟合线性模型最基本的函数就是lm()，格式为：
myfit <- lm(formula,data),其中:- formula指要拟合的模型形式
- data是一个数据框，包含了用于拟合模型的数据。
- 结果对象（本例中是myfit）存储在一个列表中，包含了所拟合模型的大量信息。
-
表达式（formula）形式如下：
y ~ x1+x2+x3+...+xk- ~左边为响应变量
- 右边为各个预测变量，预测变量之间用+符号分隔
-
R表达式中常用的符号
符号 用 途 ~ 分隔符号，左边为响应变量，右边为解释变量。</br>例如，要通过x、z和w预测y，代码为y ~ x + z + w + 分隔预测变量 : 表示预测变量的交互项。</br>例如，要通过x、z及x与z的交互项预测y，代码为y ~ x + z + x:z * 表示所有可能交互项的简洁方式。</br>代码y~ x * z * w可展开为y ~ x + z + w + x:z + x:w + z:w +x:z:w ^ 表示交互项达到某个次数。</br>代码y ~ (x + z + w)^2可展开为y ~ x + z + w + x:z + x:w + z:w . 表示包含除因变量外的所有变量。</br>例如，若一个数据框包含变量x、y、z和w，代码y ~ .可展开为y ~ x +z + w - 减号，表示从等式中移除某个变量。</br>例如，y ~ (x + z + w)^2 – x:w可展开为y ~ x + z + w + x:z + z:w -1 删除截距项。</br>例如，表达式y ~ x - 1拟合y在x上的回归，并强制直线通过原点 I() 从算术的角度来解释括号中的元素。</br>例如，y ~ x + (z + w)^2将展开为y ~ x + z + w + z:w。</br>相反, 代码y~ x + I((z + w)^2)将展开为y ~ x + h, h是一个由z和w的平方和创建的新变量 function 可以在表达式中用的数学函数。</br>例如，log(y) ~ x + z + w表示通过x、z和w来预测log(y) -
对拟合线性模型非常有用的其他函数
函 数 用 途 summary() 展示拟合模型的详细结果 coefficients() 列出拟合模型的模型参数（截距项和斜率） confint() 提供模型参数的置信区间（默认95%） fitted() 列出拟合模型的预测值 residuals() 列出拟合模型的残差值 anova() 生成一个拟合模型的方差分析表，或者比较两个或更多拟合模型的方差分析表 vcov() 列出模型参数的协方差矩阵 AIC() 输出赤池信息统计量 plot() 生成评价拟合模型的诊断图 predict() 用拟合模型对新的数据集预测响应变量值
-
-
简单线性回归：当回归模型包含一个因变量和一个自变量时，我们称为简单线性回归
- 示例：基础安装中的数据集women提供了15个年龄在30~39岁间女性的身高和体重信息，我们想通过身高来预测体重，获得一个等式可以帮助我们分辨出那些过重或过瘦的个体
- 代码：
## 拟合回归模型 fit <- lm(weight ~ height, data=women) ## 查看拟合后的详细结果 summary(fit) Call: lm(formula = weight ~ height, data = women) Residuals: Min 1Q Median 3Q Max -1.7333 -1.1333 -0.3833 0.7417 3.1167 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -87.51667 5.93694 -14.74 1.71e-09 *** height 3.45000 0.09114 37.85 1.09e-14 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.525 on 13 degrees of freedom Multiple R-squared: 0.991, Adjusted R-squared: 0.9903 F-statistic: 1433 on 1 and 13 DF, p-value: 1.091e-14 ## 在Pr(>|t|)栏，可以看到回归系数（3.45）显著不为0（p<0.001），表明身高每增高1英寸，体重将预期增加3.45磅 ## R平方项（0.991）表明模型可以解释体重99.1%的方差，它也是实际和预测值之间的相关系数（R2 = r2ŶY）。 ## 残差标准误（1.53 lbs）则可认为是模型用身高预测体重的平均误差。 ## F统计量检验所有的预测变量预测响应变量是否都在某个几率水平之上。由于简单回归只有一个预测变量，此处F检验等同于身高回归系数的t检验 women$weight ## 列出拟合模型的预测值 fitted(fit) ## 列出拟合模型的残差值 residuals(fit) ## 提供模型参数的置信区间（默认95%） confint(fit) plot(women$height,women$weight, main="Women Age 30-39", xlab="Height (in inches)", ylab="Weight (in pounds)") # add the line of best fit abline(fit)
-
多项式回归,。当只有一个预测变量，但同时包含变量的幂（比如，X、X2、X3）时，我们称之为多项式回归
针对上诉身高问题，你可以通过添加一个二次项（即
X2）来提高回归的预测精度，因此可以拟合含二次项的等式：fit2 <- lm(weight ~ height + I(height^2),data=women)-
代码
fit2 <- lm(weight ~ height + I(height^2), data=women) summary(fit2) Call: lm(formula = weight ~ height + I(height^2), data = women) Residuals: Min 1Q Median 3Q Max -0.50941 -0.29611 -0.00941 0.28615 0.59706 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 261.87818 25.19677 10.393 2.36e-07 *** height -7.34832 0.77769 -9.449 6.58e-07 *** I(height^2) 0.08306 0.00598 13.891 9.32e-09 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.3841 on 12 degrees of freedom Multiple R-squared: 0.9995, Adjusted R-squared: 0.9994 F-statistic: 1.139e+04 on 2 and 12 DF, p-value: < 2.2e-16 ## 在p<0.001水平下，回归系数都非常显著。模型的方差解释率已经增加到了99.9%。 ## 二次项的显著性（t = 13.89，p<0.001）表明包含二次项提高了模型的拟合度。 plot(women$height,women$weight, main="Women Age 30-39", xlab="Height (in inches)", ylab="Weight (in lbs)") lines(women$height,fitted(fit2)) note
image- 在car包中的scatterplot()函数，它可以很容易、方便地绘制二元关系图,==这个功能加强的图形，既提供了身高与体重的散点图、线性拟合曲线和平滑拟合（loess）曲线，还在相应边界展示了每个变量的箱线图==
- spread=FALSE选项删除了残差正负均方根在平滑曲线上的展开和非对称信息。
- lty.smooth=2选项设置loess拟合曲线为虚线。
- pch=19选项设置点为实心圆（默认为空心圆）
library(car) scatterplot(weight ~ height, data=women, spread=FALSE, smoother.args=list(lty=2), pch=19, main="Women Age 30-39", xlab="Height (inches)", ylab="Weight (lbs.)")
-
多元线性回归,当有不止一个预测变量时，则称为多元线性回归
多元回归分析中，第一步最好检查一下变量间的相关性。cor()函数提供了二变量之间的相关系数，car包中scatterplotMatrix()函数则会生成散点图矩阵
-
使用lm()函数拟合多元线性回归模型
- 当预测变量不止一个时，回归系数的含义为，一个预测变量增加一个单位，其他预测变量保持不变时，因变量将要增加的数量。
cor(states) library(car) scatterplotMatrix(states,spread = FLASE,smooth = 2,main="Scatter Plot Matrix") fit <- lm(Murder ~ Population+Illiteracy+Income+Frost,data=states) summary(fit)
-
有交互项的多元线性回归
在初步判断多个变量之间存在互相作用时，可以将变量之间的交互项加入拟合模型中，通过对交互项在模型中的pvalue验证是否相关性显著
-
可以通过effects包中的effect()函数，你可以用图形展示交互项的结果
- 格式：plot(effect(term,mod,xlevels),multiple=TRUE)
- term即模型要画的项，mod为通过lm()拟合的模型，xlevels是一个列表，指定变量要设定的常量值，multiline=TRUE选项表示添加相应直线
library(effects) plot(effect("hp:wt", fitcar, xlevels=list(wt=c(2.2, 3.2, 4.2))), multiline=TRUE)
-
-
回归诊断
-
标准方法：最常见的方法就是对lm()函数返回的对象使用plot()函数，可以生成评价模型拟合情况的四幅图形
image正态性：当预测变量值固定时，因变量成正态分布，则残差值也应该是一个均值为0的正态分布。正态Q-Q图（Normal Q-Q，右上）是在正态分布对应的值下，标准化残差的概率图。若满足正态假设，那么图上的点应该落在呈45度角的直线上；若不是如此，那么就违反了正态性的假设。
独立性：从这些图中分辨出因变量值是否相互独立，只能从收集的数据中来验证。
线性：若因变量与自变量线性相关，那么残差值与预测（拟合）值就没有任何系统关联。换句话说，除了白噪声，模型应该包含数据中所有的系统方差。在“残差图与拟合图”（Residuals vs Fitted，左上）中可以清楚的看到一个曲线关系，这暗示着你可能需要对回归模型加上一个二次项。
同方差性：若满足不变方差假设，那么在位置尺度图（Scale-Location Graph，左下）中，水平线周围的点应该随机分布。该图似乎满足此假设。
-
最后一幅“残差与杠杆图”（Residuals vs Leverage，右下）提供了你可能关注的单个观测点的信息。从图形可以鉴别出离群点、高杠杆值点和强影响点。
- 一个观测点是离群点，表明拟合回归模型对其预测效果不佳（产生了巨大的或正或负的残差）。
- 一个观测点有很高的杠杆值，表明它是一个异常的预测变量值的组合。也就是说，在预测变量空间中，它是一个离群点。因变量值不参与计算一个观测点的杠杆值。
- 一个观测点是强影响点（influential observation），表明它对模型参数的估计产生的影响过大，非常不成比例。强影响点可以通过Cook距离即Cook’s D统计量来鉴别。
-
改进的方法
-
car包提供了大量函数，大大增强了拟合和评价回归模型的能力
函 数 目 的 qqPlot() 分位数比较图 durbinWatsonTest() 对误差自相关性做Durbin-Watson检验 crPlots() 成分与残差图 ncvTest() 对非恒定的误差方差做得分检验 spreadLevelPlot() 分散水平检验 outlierTest() Bonferroni离群点检验 avPlots() 添加的变量图形 inluencePlot() 回归影响图 scatterplot() 增强的散点图 scatterplotMatrix() 增强的散点图矩阵 vif() 方差膨胀因子 - 正态性：与基础包中的plot()函数相比，qqPlot()函数提供了更为精确的正态假设检验方法，它画出了在n-p-1个自由度的t分布下的学生化残差（studentized residual，也称学生化删除残差或折叠化残差）图形，其中n是样本大小，p是回归参数的数目（包括截距项）。
- 独立性：Durbin-Watson检验函数，能够检测误差的序列相关性，该检验适用于时间独立的数据，对于非聚集型的数据并不适用
- 线性：通过成分残差图（component plus residual plot）也称偏残差图（partial residual plot），可以看看因变量与自变量之间是否呈非线性关系，也可以看看是否有不同于已设定线性模型的系统偏差，图形可用car包中的crPlots()函数绘制。若图形存在非线性，则说明你可能对预测变量的函数形式建模不够充分，那么就需要添加一些曲线成分，比如多项式项，或对一个或多个变量进行变换（如用log(X)代替X），或用其他回归变体形式而不是线性回归
- 同方差性：car包提供了两个有用的函数，可以判断误差方差是否恒定。ncvTest()函数生成一个计分检验，零假设为误差方差不变，备择假设为误差方差随着拟合值水平的变化而变化。若检验显著，则说明存在异方差性（误差方差不恒定）。spreadLevelPlot()函数创建一个添加了最佳拟合曲线的散点图，展示标准化残差绝对值与拟合值的关系；通过分布水平图看到这一点，其中的点在水平的最佳拟合曲线周围呈水平随机分布。若违反了该假设你将会看到一个非水平的曲线。代码结果建议幂次变换（suggested power transformation）的含义是，经过p次幂（Y p）变换，非恒定的误差方差将会平稳。例如，若图形显示出了非水平趋势，建议幂次转换为0.5，在回归等式中用 Y 代替Y，可能会使模型满足同方差性。若建议幂次为0，则使用对数变换。
-
-
线性模型假设的综合验证
library(gvlma) gvmodel <- gvlma(fit2) summary(gvmodel) 多重共线性
多重共线性可用统计量VIF（Variance Inflation Factor，方差膨胀因子）进行检测。VIF的平方根表示变量回归参数的置信区间能膨胀为与模型无关的预测变量的程度（因此而得名）。car包中的vif()函数提供VIF值。一般原则下， vif的平方根>2就表明存在多重共线性问题
-
-
异常观测值
一个全面的回归分析要覆盖对异常值的分析，包括离群点、高杠杆值点和强影响点
-
离群点
- 离群点是指那些模型预测效果不佳的观测点。它们通常有很大的、或正或负的残差。正的残差说明模型低估了响应值，负的残差则说明高估了响应值。
- 鉴别离群点的方法：的Q-Q图，落在置信区间带外的点即可被认为是离群点。另外一个粗糙的判断准则：标准化残差值大于2或者小于2的点可能是离群点，需要特别关注。
- car包提供了一种离群点的统计检验方法。outlierTest()函数可以求得最大标准化残差绝对值Bonferroni调整后的p值.该函数只是根据单个最大（或正或负）残差值的显著性来判断是否有离群点。若不显著，则说明数据集中没有离群点；若显著，则你必须删除该离群点，然后再检验是否还有其他离群点存在
-
高杠杆值点
- 高杠杆值观测点，即是与其他预测变量有关的离群点。换句话说，它们是由许多异常的预测变量值组合起来的，与响应变量值没有关系。
- 高杠杆值的观测点可通过帽子统计量（hat statistic）判断。对于一个给定的数据集，帽子均值为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) -
强影响点
强影响点，即对模型参数估计值影响有些比例失衡的点。
有两种方法可以检测强影响点：Cook距离，或称D统计量，以及变量添加图（added variableplot）。
-
Cook’s D值大于4/(n-k-1)，则表明它是强影响点，其中n 为样本量大小，k 是预测变量数目
# Cooks Distance D # identify D values > 4/(n-k-1) cutoff <- 4/(nrow(states)-length(fit$coefficients)-2) plot(fit, which=4, cook.levels=cutoff) abline(h=cutoff, lty=2, col="red") -
Cook’s D图有助于鉴别强影响点，但是并不提供关于这些点如何影响模型的信息。变量添加图弥补了这个缺陷。对于一个响应变量和k个预测变量，你可以创建k个变量添加图。所谓变量添加图，即对于每个预测变量Xk，绘制Xk 在其他k-1个预测变量上回归的残差值相对于响应变量在其他k-1个预测变量上回归的残差值的关系图。car包中的avPlots()函数可提供变量添加图
avPlots(fit, ask=FALSE, id.method="identify")
-
利用car包中的influencePlot()函数，你还可以将离群点、杠杆值和强影响点的信息整合到一幅图形中
influencePlot(fit, id.method="identify", main="Influence Plot", sub="Circle size is proportial to Cook's Distance" )
-
-
改进措施
有四种方法可以处理违背回归假设的问题：
-
删除观测点:
删除离群点通常可以提高数据集对于正态假设的拟合度，而强影响点会干扰结果，通常也会被删除。删除最大的离群点或者强影响点后，模型需要重新拟合。若离群点或强影响点仍然存在，重复以上过程直至获得比较满意的拟合。
-
变量变换；
当模型不符合正态性、线性或者同方差性假设时，一个或多个变量的变换通常可以改善或调整模型效果
-
常见的变换
imagesummary(powerTransform(states$Murder)) bcPower Transformation to Normality Est Power Rounded Pwr Wald Lwr bnd Wald Upr Bnd states$Murder 0.6055 1 0.0884 1.1227 Likelihood ratio tests about transformation parameters LRT df pval LR test, lambda = (0) 5.665991 1 0.01729694 LR test, lambda = (1) 2.122763 1 0.14512456 ## 结果表明，你可以用Murder0.6来正态化变量Murder。 ## 由于0.6很接近0.5，你可以尝试用平方根变换来提高模型正态性的符合程度。 ## 但在本例中，λ= 1的假设也无法拒绝（p=0.145），因此没有强有力的证据表明本例需要变量变换 -
当违反了线性假设时，对预测变量进行变换常常会比较有用。car包中的boxTidwell()函数通过获得预测变量幂数的最大似然估计来改善线性关系
boxTidwell(Murder~Population+Illiteracy,data=states) Score Statistic p-value MLE of lambda Population -0.3228003 0.7468465 0.8693882 Illiteracy 0.6193814 0.5356651 1.3581188 iterations = 19 ## 结果显示，使用变换Population0.87和Illiteracy1.36能够大大改善线性关系 ## 但是对Population（p=0.75）和Illiteracy（p=0.54）的计分检验又表明变量并不需要变换。
-
添加或删除变量
- 最常见的方法就是删除某个存在多重共线性的变量（某个变量 sqrt(vif)> 2 ）。
- 另外一个可用的方法便是岭回归——多元回归的变体，专门用来处理多重共线性问题
使用其他回归方法。
-
-
选择“最佳”的回归模型
- 模型比较
- 用基础安装中的anova()函数可以比较两个嵌套模型的拟合优度。所谓嵌套模型，即它的一些项完全包含在另一个模型中
states <- as.data.frame(state.x77[,c("Murder", "Population", "Illiteracy", "Income", "Frost")]) fit1 <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states) fit2 <- lm(Murder ~ Population + Illiteracy, data=states) anova(fit2, fit1) - AIC（Akaike Information Criterion，赤池信息准则）也可以用来比较模型，它考虑了模型的统计拟合度以及用来拟合的参数数目。AIC值越小的模型要优先选择，它说明模型用较少的参数获得了足够的拟合度
fit1 <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states) fit2 <- lm(Murder ~ Population + Illiteracy, data=states) AIC(fit1,fit2)
- 用基础安装中的anova()函数可以比较两个嵌套模型的拟合优度。所谓嵌套模型，即它的一些项完全包含在另一个模型中
- 变量选择
-
逐步回归法（stepwise method）：逐步回归中，模型会一次添加或者删除一个变量，直到达到某个判停准则为止。
- 向前逐步回归（forward stepwise）每次添加一个预测变量到模型中，直到添加变量不会使模型有所改进为止。
- 向后逐步回归（backward stepwise）从模型包含所有预测变量开始，一次删除一个变量直到会降低模型质量为止。
- 向前向后逐步回归（stepwise stepwise，通常称作逐步回归），结合了向前逐步回归和向后逐步回归的方法，变量每次进入一个，但是每一步中，变量都会被重新评价，对模型没有贡献的变量将会被删除，预测变量可能会被添加、删除好几次，直到获得最优模型为止。
- MASS包中的stepAIC()函数可以实现逐步回归模型（向前、向后和向前向后），依据的是精确AIC准则
library(MASS) states <- as.data.frame(state.x77[,c("Murder", "Population", "Illiteracy", "Income", "Frost")]) fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states) stepAIC(fit, direction="backward") Start: AIC=97.75 Murder ~ Population + Illiteracy + Income + Frost Df Sum of Sq RSS AIC - Frost 1 0.021 289.19 95.753 - Income 1 0.057 289.22 95.759 <none> 289.17 97.749 - Population 1 39.238 328.41 102.111 - Illiteracy 1 144.264 433.43 115.986 Step: AIC=95.75 Murder ~ Population + Illiteracy + Income Df Sum of Sq RSS AIC - Income 1 0.057 289.25 93.763 <none> 289.19 95.753 - Population 1 43.658 332.85 100.783 - Illiteracy 1 236.196 525.38 123.605 Step: AIC=93.76 Murder ~ Population + Illiteracy Df Sum of Sq RSS AIC <none> 289.25 93.763 - Population 1 48.517 337.76 99.516 - Illiteracy 1 299.646 588.89 127.311 Call: lm(formula = Murder ~ Population + Illiteracy, data = states) Coefficients: (Intercept) Population Illiteracy 1.6515497 0.0002242 4.0807366 ## 逐步回归法其实存在争议，虽然它可能会找到一个好的模型，但是不能保证模型就是最佳模型，因为不是每一个可能的模型都被评价了 -
全子集回归（all-subsets regression），即所有可能的模型都会被检验，可以选择展示所有可能的结果，也可以展示n 个不同子集大小（一个、两个或多个预测变量）的最佳模型
- 全子集回归可用leaps包中的regsubsets()函数实现。你能通过R平方、调整R平方或Mallows Cp统计量等准则来选择“最佳”模型
- R平方含义是预测变量解释响应变量的程度
- Mallows Cp统计量也用来作为逐步回归的判停规则。广泛研究表明，对于一个好的模型，它的Cp统计量非常接近于模型的参数数目（包括截距项）
- 可用leaps包中的plot()函数绘制，或者用car包中的subsets()函数绘制全子集回归图
library(leaps) states <- as.data.frame(state.x77[,c("Murder", "Population", "Illiteracy", "Income", "Frost")]) leaps <-regsubsets(Murder ~ Population + Illiteracy + Income + Frost, data=states, nbest=4) plot(leaps, scale="adjr2") library(car) subsets(leaps, statistic="cp", main="Cp Plot for All Subsets Regression") abline(1,1,lty=2,col="red")
-
- 模型比较
-
深层次分析
-
交叉验证，评价回归方程的泛化能力
- 所谓交叉验证，即将一定比例的数据挑选出来作为训练样本，另外的样本作保留样本，先在训练样本上获取回归方程，然后在保留样本上做预测。由于保留样本不涉及模型参数的选择，该样本可获得比新数据更为精确的估计。
- k 重交叉验证，在k 重交叉验证中，样本被分为k个子样本，轮流将k1个子样本组合作为训练集，另外1个子样本作为保留集。这样会获得k 个预测方程，记录k个保留样本的预测表现结果，然后求其平均值。
- bootstrap 包中的 crossval() 函数可以实现 k 重交叉验证。
- 通过选择有更好泛化能力的模型，还可以用交叉验证来挑选变量。
- 下方示例，通过定义shrinkage()函数对模型的R平方统计量做k 重交叉验证
shrinkage <- function(fit,k=10){ require(bootstrap) # define functions theta.fit <- function(x,y){lsfit(x,y)} theta.predict <- function(fit,x){cbind(1,x)%*%fit$coef} # matrix of predictors x <- fit$model[,2:ncol(fit$model)] # vector of predicted values y <- fit$model[,1] results <- crossval(x,y,theta.fit,theta.predict,ngroup=k) r2 <- cor(y, fit$fitted.values)**2 # raw R2 r2cv <- cor(y,results$cv.fit)**2 # cross-validated R2 cat("Original R-square =", r2, "\n") cat(k, "Fold Cross-Validated R-square =", r2cv, "\n") cat("Change =", r2-r2cv, "\n") } # using it states <- as.data.frame(state.x77[,c("Murder", "Population", "Illiteracy", "Income", "Frost")]) fit <- lm(Murder ~ Population + Income + Illiteracy + Frost, data=states) shrinkage(fit) fit2 <- lm(Murder~Population+Illiteracy,data=states) shrinkage(fit2) -
相对重要性
- 最简单的莫过于比较标准化的回归系数，它表示当其他预测变量不变时，该预测变量一个标准差的变化可引起的响应变量的预期变化（以标准差单位度量）。在进行回归分析前，可用scale()函数将数据标准化为均值为0、标准差为1的数据集，这样用R回归即可获得标准化的回归系数。（注意，scale()函数返回的是一个矩阵，而lm()函数要求一个数据框，你需要用一个中间步骤来转换一下。）
states <- as.data.frame(state.x77[,c("Murder", "Population", "Illiteracy", "Income", "Frost")]) zstates <- as.data.frame(scale(states)) zfit <- lm(Murder~Population + Income + Illiteracy + Frost, data=zstates) coef(zfit) (Intercept) Population Income Illiteracy Frost -2.054026e-16 2.705095e-01 1.072372e-02 6.840496e-01 8.185407e-03 ## 当其他因素不变时，文盲率一个标准差的变化将增加0.68个标准差的谋杀率。根据标准化的回归系数，我们可认为Illiteracy是最重要的预测变量，而Frost是最不重要的- 相对权重（relative weight）是一种比较有前景的新方法，它是对所有可能子模型添加一个预测变量引起的R平方平均增加量的一个近似值
- 下方示例提供了一个生成相对权重的函数
relweights <- function(fit,...){ R <- cor(fit$model) nvar <- ncol(R) rxx <- R[2:nvar, 2:nvar] rxy <- R[2:nvar, 1] svd <- eigen(rxx) evec <- svd$vectors ev <- svd$values delta <- diag(sqrt(ev)) lambda <- evec %*% delta %*% t(evec) lambdasq <- lambda ^ 2 beta <- solve(lambda) %*% rxy rsquare <- colSums(beta ^ 2) rawwgt <- lambdasq %*% beta ^ 2 import <- (rawwgt / rsquare) * 100 import <- as.data.frame(import) row.names(import) <- names(fit$model[2:nvar]) names(import) <- "Weights" import <- import[order(import),1, drop=FALSE] dotchart(import$Weights, labels=row.names(import), xlab="% of R-Square", pch=19, main="Relative Importance of Predictor Variables", sub=paste("Total R-Square=", round(rsquare, digits=3)), ...) return(import) } # Applying the relweights function states <- as.data.frame(state.x77[,c("Murder", "Population", "Illiteracy", "Income", "Frost")]) fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states) relweights(fit, col="blue") ## 对结果进行查看，可以看到各个预测变量对模型方差的解释程度（R平方=0.567），Illiteracy解释了59%的R平方，Frost解释了20.79%。根据相对权重法，Illiteracy有最大的相对重要性，余下依次是Frost、Population和Income
-
