原文連結：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模型）。