原文链接：http://tecdat.cn/?p=11386
原文出处：拓端数据部落公众号
最近咱们被客户要求撰写对于线性模型的钻研报告，包含一些图形和统计输入。
在这篇文章中，我将从一个根本的线性模型开始，而后从那里尝试找到一个更适合的线性模型。
数据预处理
因为空气质量数据集蕴含一些缺失值，因而咱们将在开始拟合模型之前将其删除，并抉择 70%的样本进行训练并将其余样本用于测试：
data(airquality)
ozone <- subset(na.omit(airquality),
select = c("Ozone", "Solar.R", "Wind", "Temp"))
set.seed(123)
N.train <- ceiling(0.7 * nrow(ozone))
N.test <- nrow(ozone) - N.train
trainset <- sample(seq_len(nrow(ozone)), N.train)
testset <- setdiff(seq_len(nrow(ozone)), trainset)一般最小二乘模型
作为基准模型，咱们将应用一般的最小二乘（OLS）模型。在定义模型之前，咱们定义一个用于绘制线性模型的函数：
rsquared <- function(test.preds, test.labels) {return(round(cor(test.preds, test.labels)^2, 3))
}
plot.linear.model <- function(model, test.preds = NULL, test.labels = NULL,
test.only = FALSE) {
r.squared <- NULL
if (!is.null(test.preds) && !is.null(test.labels)) {
# store predicted points:
test.df <- data.frame("Prediction" = test.preds,
"Outcome" = test.labels, "DataSet" = "test")
# store residuals for predictions on the test data
test.residuals <- test.labels - test.preds
test.res.df <- data.frame("x" = test.labels, "y" = test.preds,
"x1" = test.labels, "y2" = test.preds + test.residuals,
"DataSet" = "test")
# append to existing data
plot.df <- rbind(plot.df, test.df)
plot.res.df <- rbind(plot.res.df, test.res.df)
# annotate model with R^2 value
r.squared <- rsquared(test.preds, test.labels)
}
#######
library(ggplot2)
p <- ggplot()
return(p)
}当初，咱们应用 lm 并钻研特色预计的置信区间来建设 OLS 模型：
confint(model)## 2.5 % 97.5 %
## (Intercept) -1.106457e+02 -20.88636548
## Solar.R 7.153968e-03 0.09902534
## Temp 1.054497e+00 2.07190804
## Wind -3.992315e+00 -1.24576713咱们看到模型仿佛对截距的设置不太确定。让咱们看看模型是否依然体现良好：
查看模型的拟合度，有两个次要察看后果：
- 高臭氧程度被低估
- 预计臭氧含量为负
上面让咱们更具体地钻研这两个问题。
高臭氧程度被低估
从图中能够看出，当臭氧在 [0,100] 范畴内时，线性模型非常适合后果。然而，当理论察看到的臭氧浓度高于 100 时，该模型会大大低估该值。
咱们应该问一个问题，这些高臭氧含量是否不是测量误差的后果。思考到[典型的臭氧程度]()，测量值仿佛是正当的。最高臭氧浓度为 168 ppb（十亿分之一），美国城市的典型峰值浓度为 150 至 510 ppb。这意味着咱们的确应该关注离群值。低估高臭氧含量将特地无害，因为高含量的臭氧会危害衰弱。让咱们考察数据以确定模型为何存在这些异样值的问题。
直方图表明残差散布右尾的值的确存在问题。因为残差不是真正的正态分布，因而线性模型不是最佳模型。实际上，残差仿佛遵循某种模式的泊松散布。为了找出最小二乘模型的拟合对离群值如此之差的起因，咱们再来看一下数据。
## Ozone Solar.R Wind Temp
## Min. :110.0 Min. :207.0 Min. :2.300 Min. :79.00
## 1st Qu.:115.8 1st Qu.:223.5 1st Qu.:3.550 1st Qu.:81.75
## Median :120.0 Median :231.5 Median :4.050 Median :86.50
## Mean :128.0 Mean :236.2 Mean :4.583 Mean :86.17
## 3rd Qu.:131.8 3rd Qu.:250.8 3rd Qu.:5.300 3rd Qu.:89.75
## Max. :168.0 Max. :269.0 Max. :8.000 Max. :94.00summary(ozone)## Ozone Solar.R Wind Temp
## Min. : 1.0 Min. : 7.0 Min. : 2.30 Min. :57.00
## 1st Qu.: 18.0 1st Qu.:113.5 1st Qu.: 7.40 1st Qu.:71.00
## Median : 31.0 Median :207.0 Median : 9.70 Median :79.00
## Mean : 42.1 Mean :184.8 Mean : 9.94 Mean :77.79
## 3rd Qu.: 62.0 3rd Qu.:255.5 3rd Qu.:11.50 3rd Qu.:84.50
## Max. :168.0 Max. :334.0 Max. :20.70 Max. :97.00从两组观测值的散布来看，咱们看不到高臭氧观测值与其余样本之间的微小差别。然而，咱们能够应用下面的模型预测图找到罪魁祸首。在该图中，咱们看到大多数数据点都以 [0，50] 臭氧范畴为核心。为了很好地拟合这些察看值，截距的负值为 -65.77，这就是为什么该模型低估了较大臭氧值的臭氧程度的起因，在训练数据中臭氧值有余。
该模型预测负臭氧程度
如果察看到的臭氧浓度靠近于 0，则该模型通常会预测负臭氧程度。当然，这不可能是因为浓度不能低于 0。再次，咱们考察数据以找出为什么模型依然做出这些预测。
为此，咱们将抉择臭氧程度在第 5 个百分位数的所有观测值，并考察其特征值：
# compare observations with low ozone with whole data set
summary(ozone[idx,])## Ozone Solar.R Wind Temp
## Min. :1.0 Min. : 8.00 Min. : 9.70 Min. :57.0
## 1st Qu.:4.5 1st Qu.:20.50 1st Qu.: 9.85 1st Qu.:59.5
## Median :6.5 Median :36.50 Median :12.30 Median :61.0
## Mean :5.5 Mean :37.83 Mean :13.75 Mean :64.5
## 3rd Qu.:7.0 3rd Qu.:48.75 3rd Qu.:17.38 3rd Qu.:67.0
## Max. :8.0 Max. :78.00 Max. :20.10 Max. :80.0summary(ozone)## Ozone Solar.R Wind Temp
## Min. : 1.0 Min. : 7.0 Min. : 2.30 Min. :57.00
## 1st Qu.: 18.0 1st Qu.:113.5 1st Qu.: 7.40 1st Qu.:71.00
## Median : 31.0 Median :207.0 Median : 9.70 Median :79.00
## Mean : 42.1 Mean :184.8 Mean : 9.94 Mean :77.79
## 3rd Qu.: 62.0 3rd Qu.:255.5 3rd Qu.:11.50 3rd Qu.:84.50
## Max. :168.0 Max. :334.0 Max. :20.70 Max. :97.00咱们发现，在低臭氧程度下，均匀太阳辐射要低得多，而均匀风速要高得多。要理解为什么咱们会有负面的预测，当初让咱们看一下模型系数：
coefficients(model)## (Intercept) Solar.R Temp Wind
## -65.76603538 0.05308965 1.56320267 -2.61904128因而，对于较低的臭氧程度，的正系数 Solar.R 不能补救截距，Wind因为对于较低的臭氧程度，的值 Solar.R 较低，而的值 Wind 较高。
解决负面的臭氧程度预测
让咱们首先解决预测负臭氧程度的问题。
截短的最小二乘模型
解决负面预测的一种简略办法是将其替换为尽可能小的值。这样，如果咱们将模型交给客户，他就不会开始狐疑模型有问题。咱们能够应用以下性能来做到这一点：
当初让咱们验证这将如何改善咱们对测试数据的预测。请记住，[R2[R2 最后的模型是 0.6040.604。
nonnegative.preds <- predict.nonnegative(model, ozone[testset,])
# plot only the test predictions to verify the results
plot.linear.model(model, nonnegative.preds, ozone$Ozone[testset], test.only = TRUE)如咱们所见，此 hack 能够克制问题并减少 [R2[R2 至 0.6460.646。然而，以这种形式校对负值不会扭转咱们的模型谬误的事实，因为拟合过程并未思考到负值应该是不可能的。
泊松回归
为了防止出现负预计，咱们能够应用假设为泊松散布而非正态分布的狭义线性模型（GLM）：
plot.linear.model(pois.model, pois.preds, ozone$Ozone[testset])的 [R2[R2 值 0.616 示意泊松回归比一般最小二乘（0.604）稍好。然而，其性能并不优于将负值为 0（0.646）的模型。这可能是因为臭氧程度的方差比泊松模型假如的要大得多：
# mean and variance should be the same for Poisson distribution
mean(ozone$Ozone)## [1] 42.0991var(ozone$Ozone)## [1] 1107.29对数转换
解决负面预测的另一种办法是取后果的对数：
print(rsquared(log.preds, test.labels))## [1] 0.616请留神，只管后果与通过 Poisson 回归得出的后果雷同，但这两种办法通常并不相同。
应答高估臭氧程度的低估
现实状况下，咱们将在臭氧程度较高的状况下更好地进行测量。然而，因为咱们无奈收集更多数据，因而咱们须要利用已有的资源。应答低估高臭氧程度的一种办法是调整损失函数。
加权回归
应用加权回归，咱们能够影响离群值残差的影响。为此，咱们将计算臭氧程度的 z 得分，而后将其指数用作模型的权重，从而减少异样值的影响。
plot.linear.model(weight.model, weight.preds, ozone$Ozone[testset])
该模型相对比一般的最小二乘模型更适合，因为它能够更好地解决离群值。
采样
让咱们从训练数据中进行采样，以确保不再呈现臭氧含量过高的状况。这相似于进行加权回归。然而，咱们没有为低臭氧程度的观测值设置较小的权重，而是将其权重设置为 0。
print(paste0("N (trainset before):", length(trainset)))## [1] "N (trainset before): 78"
print(paste0("N (trainset after):", length(trainset.sampled)))## [1] "N (trainset after): 48"当初，让咱们基于采样数据构建一个新模型：
rsquared(sampled.preds, test.labels)## [1] 0.612如咱们所见，基于采样数据的模型的性能并不比使用权重的模型更好。
联合
看到泊松回归可用于避免负预计，加权是改善离群值预测的胜利策略，咱们应该尝试将两种办法联合起来，从而得出加权泊松回归。
加权泊松回归
p.w.pois
如咱们所见，该模型联合了应用泊松回归（非负预测）和使用权重（低估离群值）的劣势。的确，[R2[R2 该模型的最低价（截断线性模型为 0.652 vs 0.646）。让咱们钻研模型系数：
coefficients(w.pois.model)## (Intercept) Solar.R Temp Wind
## 2.069357230 0.002226422 0.029252172 -0.104778731该模型依然由截距管制，但当初是负数。因而，如果所有其余特色的值为 0，则模型的预测仍将为正。
然而，假如均值应等于泊松回归的方差呢？
print(paste0(c("Var:", "Mean:"), c(round(var(ozone$Ozone), 2),
round(mean(ozone$Ozone), 2))))## [1] "Var: 1107.29" "Mean: 42.1"该模型的假如相对不满足，并且因为方差大于该模型假如，因而咱们蒙受了适度扩散的困扰。
加权负二项式模型
因而，咱们应该尝试抉择一个更适宜适度扩散的模型，例如[负二项式模型]()：
plot.linear.model(model.nb, preds.nb, test.labels)
因而，就测试集的性能而言，加权负二项式模型并不比加权泊松模型更好。然而，在进行推断时，该值应该更好，因为其假如没有被毁坏。查看这两个模型，很显著它们的 p 值相差很大：
coef(summary(w.pois.model))## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.069357230 0.2536660583 8.157801 3.411790e-16
## Solar.R 0.002226422 0.0003373846 6.599061 4.137701e-11
## Temp 0.029252172 0.0027619436 10.591155 3.275269e-26
## Wind -0.104778731 0.0064637151 -16.210295 4.265135e-59coef(summary(model.nb))## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.241627650 0.640878750 1.937383 5.269853e-02
## Solar.R 0.002202194 0.000778691 2.828071 4.682941e-03
## Temp 0.037756464 0.007139521 5.288375 1.234078e-07
## Wind -0.088389583 0.016333237 -5.411639 6.245051e-08尽管泊松模型宣称所有系数都十分显着，但负二项式模型表明截距并不显着。负二项式的相信带能够通过以下形式找到：
CI.int <- 0.95
# calculate CIs:
df <- data.frame(ozone[testset,], "PredictedOzone" = ilink(preds.nb.ci$fit), "Lower" = ilink(preds.nb.ci$fit - ci.factor * preds.nb.ci$se.fit),
"Upper" = ilink(preds.nb.ci$fit + ci.factor * preds.nb.ci$se.fit))应用蕴含测试集中的特征值以及带有其相信带的预测的结构数据框，咱们能够绘制估计值如何依据独立变量而稳定：
# plot estimates vs individual features in the test set
plots <- list()
grid.arrange(plots[[1]], plots[[2]], plots[[3]])这些图阐明了两件事：
- 有清晰的线性关系
Wind和Temperature。预计的臭氧程度Wind随减少而降落，而预计的臭氧程度随减少而Temp减少。 - 该模型对低臭氧程度最有信念，但对高臭氧程度不太有信念
数据集裁减
优化模型后，咱们当初返回初始数据集。还记得咱们在剖析开始时就删除了所有缺失值的察看后果吗？好吧，这是不现实的，因为咱们曾经舍弃了有价值的信息，这些信息能够用来取得更好的模型。
考察缺失值
让咱们首先考察缺失的值：
# ratio of missing values
ratio.missing <- length(na.idx) / nrow(ozone)
print(paste0(round(ratio.missing * 100, 3), "%"))## [1] "27.451%"# which features are missing most often?
nbr.missing <- apply(ozone, 2, function(x) length(which(is.na(x))))
print(nbr.missing)## Ozone Solar.R Wind Temp
## 37 7 0 0# multiple features missing in one observation?
nbr.missing <- apply(ozone, 1, function(x) length(which(is.na(x))))
table(nbr.missing)## nbr.missing
## 0 1 2
## 111 40 2考察显示，因为短少值，以前排除了相当多的察看值。更具体地说，惟一短少的性能是Ozone（37 次）和Solar.R（7 次）。通常，只有一项性能缺失（40 次），很少有两项性能缺失（2 次）。
调整训练和测试指标
为了确保与以前应用雷同的观测值进行测试，咱们必须 映射到残缺的空气质量数据集：
trainset <- c(trainset, na.idx)
testset <- setdiff(seq_len(nrow(ozone)), trainset)估算缺失值
为了取得缺失值的估计值，咱们能够应用插补。这种办法的想法是应用已知特色来造成预测模型，以便预计缺失的特色。
summary(as.numeric(imputed.data$Ozone))## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.00 16.00 30.00 41.66 59.00 168.00请留神，aregImpute应用不同的疏导程序样本进行多个插补，能够应用 n.impute 参数指定。因为咱们要应用所有运行的推算而不是单个运行，因而咱们将应用 fit.mult.impute 函数定义模型：
# compute new weights
plot.linear.model(fmi, fmi.preds, ozone$Ozone[testset])
让咱们仅应用一个插补以便再次指定权重：
rsquared(w.pois.preds.imputed, imputed.data$Ozone[testset])## [1] 0.431在这种状况下，基于估算数据的加权泊松模型的性能不会比仅排除失落数据的模型更好。这表明对缺失值的估算比将噪声引入数据中要多得多，而不是咱们能够应用的信号。可能的解释是，具备缺失值的样本具备不同于所有测量可用值的散布。
摘要
咱们从 OLS 回归模型开始（[R2= 0.604[R2=0.604），并试图找到一个更适合的线性模型。第一个想法是将模型的预测截断为 0（[R2= 0.646[R2=0.646）。为了更精确地预测离群值，咱们训练了加权线性回归模型（[R2= 0.621[R2=0.621）。接下来，为了仅预测正值，咱们训练了加权 Poisson 回归模型（[R2= 0.652[R2=0.652）。为了解决泊松模型中的适度扩散问题，咱们制订了加权负二项式模型。只管此模型的体现不如加权 Poisson 模型（[R2= 0.638），则在进行推理时可能会更好。
尔后，咱们尝试通过应用 Hmisc 包估算缺失值来进一步改良模型。只管生成的模型比初始 OLS 模型要好，然而它们没有取得比以前更高的性能（[R2= 0.627[R2=0.627）。
那么，最好的模型到底是什么？就模型假如的正确性而言，这是加权负二项式模型。就决定系数而言，[R2[R2，这是加权 Poisson 回归模型。因而，出于预测臭氧程度的目标，我将抉择加权 Poisson 回归模型。
您可能会问：所有这些工作值得吗？实际上，初始模型和加权泊松模型的预测在 5%的程度上存在显着差别：
##
## Wilcoxon signed rank test
##
## data: test.preds and w.pois.preds
## V = 57, p-value = 1.605e-05
## alternative hypothesis: true location shift is not equal to 0当咱们 比拟它们时，模型之间的差别变得不言而喻：
总之，咱们从预测负值和低估高臭氧程度的模型（左侧显示 OLS 模型）到没有此类显著缺点的模型（右侧加权 Poisson 模型）。