R語言用線性模型進(jìn)行臭氧預(yù)測： 加權(quán)泊松回歸，普通最小二乘，加權(quán)負(fù)二項(xiàng)式模型，多重
?原文鏈接：http://tecdat.cn/?p=11386
?在這篇文章中，我將從一個(gè)基本的線性模型開始，然后從那里嘗試找到一個(gè)更合適的線性模型。
數(shù)據(jù)預(yù)處理
由于空氣質(zhì)量數(shù)據(jù)集包含一些缺失值，因此我們將在開始擬合模型之前將其刪除，并選擇70％的樣本進(jìn)行訓(xùn)練并將其余樣本用于測試：
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.traintrainset <- sample(seq_len(nrow(ozone)), N.train)testset <- setdiff(seq_len(nrow(ozone)), trainset)
普通最小二乘模型
作為基準(zhǔn)模型，我們將使用普通的最小二乘（OLS）模型。在定義模型之前，我們定義一個(gè)用于繪制線性模型的函數(shù)：
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 <- NULLif (!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 datatest.residuals <- test.labels - test.predstest.res.df <- data.frame("x" = test.labels, "y" = test.preds,"x1" = test.labels, "y2" = test.preds + test.residuals,"DataSet" = "test")# append to existing dataplot.df <- rbind(plot.df, test.df)plot.res.df <- rbind(plot.res.df, test.res.df)# annotate model with R^2 valuer.squared <- rsquared(test.preds, test.labels)}#######library(ggplot2)p <- ggplot()return(p)}
現(xiàn)在，我們使用lm并研究特征估計(jì)的置信區(qū)間來建立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
我們看到模型似乎對截距的設(shè)置不太確定。讓我們看看模型是否仍然表現(xiàn)良好：
?查看模型的擬合度，有兩個(gè)主要觀察結(jié)果：
高臭氧水平被低估
預(yù)計(jì)臭氧含量為負(fù)
下面讓我們更詳細(xì)地研究這兩個(gè)問題。
高臭氧水平被低估
從圖中可以看出，當(dāng)臭氧在[0,100]范圍內(nèi)時(shí)，線性模型非常適合結(jié)果。但是，當(dāng)實(shí)際觀察到的臭氧濃度高于100時(shí)，該模型會(huì)大大低估該值。
我們應(yīng)該問一個(gè)問題，這些高臭氧含量是否不是測量誤差的結(jié)果。考慮到典型的臭氧水平，測量值似乎是合理的。最高臭氧濃度為168 ppb（十億分之一），美國城市的典型峰值濃度為150至510 ppb。這意味著我們確實(shí)應(yīng)該關(guān)注離群值。低估高臭氧含量將特別有害，因?yàn)楦吆康某粞鯐?huì)危害健康。讓我們調(diào)查數(shù)據(jù)以確定模型為何存在這些異常值的問題。
?
?直方圖表明殘差分布右尾的值確實(shí)存在問題。由于殘差不是真正的正態(tài)分布，因此線性模型不是最佳模型。實(shí)際上，殘差似乎遵循某種形式的泊松分布。為了找出最小二乘模型的擬合對離群值如此之差的原因，我們再來看一下數(shù)據(jù)。
## ? ? ?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.00
summary(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
從兩組觀測值的分布來看，我們看不到高臭氧觀測值與其他樣本之間的巨大差異。但是，我們可以使用上面的模型預(yù)測圖找到罪魁禍?zhǔn)?。在該圖中，我們看到大多數(shù)數(shù)據(jù)點(diǎn)都以[0，50]臭氧范圍為中心。為了很好地?cái)M合這些觀察值，截距的負(fù)值為-65.77，這就是為什么該模型低估了較大臭氧值的臭氧水平的原因，在訓(xùn)練數(shù)據(jù)中臭氧值不足。
該模型預(yù)測負(fù)臭氧水平
如果觀察到的臭氧濃度接近于0，則該模型通常會(huì)預(yù)測負(fù)臭氧水平。當(dāng)然，這不可能是因?yàn)闈舛炔荒艿陀?。再次，我們調(diào)查數(shù)據(jù)以找出為什么模型仍然做出這些預(yù)測。
為此，我們將選擇臭氧水平在第5個(gè)百分位數(shù)的所有觀測值，并調(diào)查其特征值：
# compare observations with low ozone with whole data setsummary(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.0
summary(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
我們發(fā)現(xiàn)，在低臭氧水平下，平均太陽輻射要低得多，而平均風(fēng)速要高得多。要了解為什么我們會(huì)有負(fù)面的預(yù)測，現(xiàn)在讓我們看一下模型系數(shù)：
coefficients(model)
## ?(Intercept) ? ? ?Solar.R ? ? ? ? Temp ? ? ? ? Wind## -65.76603538 ? 0.05308965 ? 1.56320267 ?-2.61904128
因此，對于較低的臭氧水平，的正系數(shù)Solar.R不能彌補(bǔ)截距，Wind因?yàn)閷τ谳^低的臭氧水平，的值Solar.R較低，而的值Wind較高。
處理負(fù)面的臭氧水平預(yù)測
讓我們首先處理預(yù)測負(fù)臭氧水平的問題。
截短的最小二乘模型
處理負(fù)面預(yù)測的一種簡單方法是將其替換為盡可能小的值。這樣，如果我們將模型交給客戶，他就不會(huì)開始懷疑模型有問題。我們可以使用以下功能來做到這一點(diǎn)：
現(xiàn)在讓我們驗(yàn)證這將如何改善我們對測試數(shù)據(jù)的預(yù)測。請記住，[R2[R2?最初的模型是?0.6040.604。
nonnegative.preds <- predict.nonnegative(model, ozone[testset,])# plot only the test predictions to verify the resultsplot.linear.model(model, nonnegative.preds, ozone$Ozone[testset], test.only = TRUE)
?如我們所見，此hack可以抑制問題并增加?[R2[R2?至?0.6460.646。但是，以這種方式校正負(fù)值不會(huì)改變我們的模型錯(cuò)誤的事實(shí)，因?yàn)閿M合過程并未考慮到負(fù)值應(yīng)該是不可能的。
泊松回歸
為了防止出現(xiàn)負(fù)估計(jì)，我們可以使用假定為泊松分布而非正態(tài)分布的廣義線性模型（GLM）：
plot.linear.model(pois.model, pois.preds, ozone$Ozone[testset])
?的?[R2[R2值0.616表示泊松回歸比普通最小二乘（0.604）稍好。但是，其性能并不優(yōu)于將負(fù)值為0（0.646）的模型。這可能是因?yàn)槌粞跛降姆讲畋炔此赡Ｐ图僭O(shè)的要大得多：
# mean and variance should be the same for Poisson distributionmean(ozone$Ozone)
## [1] 42.0991var(ozone$Ozone)## [1] 1107.29
對數(shù)轉(zhuǎn)換
處理負(fù)面預(yù)測的另一種方法是取結(jié)果的對數(shù)：
print(rsquared(log.preds, test.labels))
## [1] 0.616
請注意，盡管結(jié)果與通過Poisson回歸得出的結(jié)果相同，但這兩種方法通常并不相同。
應(yīng)對高估臭氧水平的低估
理想情況下，我們將在臭氧水平較高的情況下更好地進(jìn)行測量。但是，由于我們無法收集更多數(shù)據(jù)，因此我們需要利用已有的資源。應(yīng)對低估高臭氧水平的一種方法是調(diào)整損失函數(shù)。
加權(quán)回歸
使用加權(quán)回歸，我們可以影響離群值殘差的影響。為此，我們將計(jì)算臭氧水平的z得分，然后將其指數(shù)用作模型的權(quán)重，從而增加異常值的影響。
plot.linear.model(weight.model, weight.preds, ozone$Ozone[testset])
?
?該模型絕對比普通的最小二乘模型更合適，因?yàn)樗梢愿玫靥幚黼x群值。
采樣
讓我們從訓(xùn)練數(shù)據(jù)中進(jìn)行采樣，以確保不再出現(xiàn)臭氧含量過高的情況。這類似于進(jìn)行加權(quán)回歸。但是，我們沒有為低臭氧水平的觀測值設(shè)置較小的權(quán)重，而是將其權(quán)重設(shè)置為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"
現(xiàn)在，讓我們基于采樣數(shù)據(jù)構(gòu)建一個(gè)新模型：
rsquared(sampled.preds, test.labels)
## [1] 0.612
如我們所見，基于采樣數(shù)據(jù)的模型的性能并不比使用權(quán)重的模型更好。
結(jié)合
看到泊松回歸可用于防止負(fù)估計(jì)，加權(quán)是改善離群值預(yù)測的成功策略，我們應(yīng)該嘗試將兩種方法結(jié)合起來，從而得出加權(quán)泊松回歸。
加權(quán)泊松回歸
p.w.pois
?
?如我們所見，該模型結(jié)合了使用泊松回歸（非負(fù)預(yù)測）和使用權(quán)重（低估離群值）的優(yōu)勢。確實(shí)，[R2[R2該模型的最低價(jià)（截?cái)嗑€性模型為0.652 vs 0.646）。讓我們研究模型系數(shù)：
coefficients(w.pois.model)
## ?(Intercept) ? ? ?Solar.R ? ? ? ? Temp ? ? ? ? Wind## ?2.069357230 ?0.002226422 ?0.029252172 -0.104778731
該模型仍然由截距控制，但現(xiàn)在是正數(shù)。因此，如果所有其他特征的值為0，則模型的預(yù)測仍將為正。
但是，假設(shè)均值應(yīng)等于泊松回歸的方差呢？
print(paste0(c("Var: ", "Mean: "), c(round(var(ozone$Ozone), 2),round(mean(ozone$Ozone), 2))))
## [1] "Var: 1107.29" "Mean: 42.1"
該模型的假設(shè)絕對不滿足，并且由于方差大于該模型假設(shè)，因此我們遭受了過度分散的困擾。
加權(quán)負(fù)二項(xiàng)式模型
因此，我們應(yīng)該嘗試選擇一個(gè)更適合過度分散的模型，例如負(fù)二項(xiàng)式模型：
plot.linear.model(model.nb, preds.nb, test.labels)
?
?因此，就測試集的性能而言，加權(quán)負(fù)二項(xiàng)式模型并不比加權(quán)泊松模型更好。但是，在進(jìn)行推斷時(shí)，該值應(yīng)該更好，因?yàn)槠浼僭O(shè)沒有被破壞。查看這兩個(gè)模型，很明顯它們的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-59
coef(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
雖然泊松模型聲稱所有系數(shù)都非常顯著，但負(fù)二項(xiàng)式模型表明截距并不顯著。 負(fù)二項(xiàng)式的置信帶可以通過以下方式找到：
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))
使用包含測試集中的特征值以及帶有其置信帶的預(yù)測的構(gòu)造數(shù)據(jù)框，我們可以繪制估計(jì)值如何根據(jù)獨(dú)立變量而波動(dòng)：
# plot estimates vs individual features in the test setplots <- list()grid.arrange(plots[[1]], plots[[2]], plots[[3]])
?這些圖說明了兩件事：
有清晰的線性關(guān)系
Wind和Temperature。估計(jì)的臭氧水平Wind隨增加而下降，而估計(jì)的臭氧水平隨增加而Temp增加。該模型對低臭氧水平最有信心，但對高臭氧水平不太有信心
數(shù)據(jù)集擴(kuò)充
優(yōu)化模型后，我們現(xiàn)在返回初始數(shù)據(jù)集。還記得我們在分析開始時(shí)就刪除了所有缺失值的觀察結(jié)果嗎？好吧，這是不理想的，因?yàn)槲覀円呀?jīng)舍棄了有價(jià)值的信息，這些信息可以用來獲得更好的模型。
調(diào)查缺失值
讓我們首先調(diào)查缺失的值：
# ratio of missing valuesratio.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
調(diào)查顯示，由于缺少值，以前排除了相當(dāng)多的觀察值。更具體地說，唯一缺少的功能是Ozone（37次）和Solar.R（7次）。通常，只有一項(xiàng)功能缺失（40次），很少有兩項(xiàng)功能缺失（2次）。
調(diào)整訓(xùn)練和測試指標(biāo)
為了確保與以前使用相同的觀測值進(jìn)行測試，我們必須 映射到完整的空氣質(zhì)量數(shù)據(jù)集：
trainset <- c(trainset, na.idx)testset <- setdiff(seq_len(nrow(ozone)), trainset)
估算缺失值
為了獲得缺失值的估計(jì)值，我們可以使用插補(bǔ)。這種方法的想法是使用已知特征來形成預(yù)測模型，以便估計(jì)缺失的特征。?
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使用不同的引導(dǎo)程序樣本進(jìn)行多個(gè)插補(bǔ)，可以使用n.impute參數(shù)指定。由于我們要使用所有運(yùn)行的推算而不是單個(gè)運(yùn)行，因此我們將使用fit.mult.impute函數(shù)定義模型：
# compute new weightsplot.linear.model(fmi, fmi.preds, ozone$Ozone[testset])
?
?讓我們僅使用一個(gè)插補(bǔ)以便再次指定權(quán)重：
rsquared(w.pois.preds.imputed, imputed.data$Ozone[testset])
## [1] 0.431
在這種情況下，基于估算數(shù)據(jù)的加權(quán)泊松模型的性能不會(huì)比僅排除丟失數(shù)據(jù)的模型更好。這表明對缺失值的估算比將噪聲引入數(shù)據(jù)中要多得多，而不是我們可以使用的信號(hào)。可能的解釋是，具有缺失值的樣本具有不同于所有測量可用值的分布。
摘要
我們從OLS回歸模型開始（[R2=?0.604[R2=0.604），并試圖找到一個(gè)更合適的線性模型。第一個(gè)想法是將模型的預(yù)測截?cái)酁?（[R2=?0.646[R2=0.646）。為了更準(zhǔn)確地預(yù)測離群值，我們訓(xùn)練了加權(quán)線性回歸模型（[R2=?0.621[R2=0.621）。接下來，為了僅預(yù)測正值，我們訓(xùn)練了加權(quán)Poisson回歸模型（[R2=?0.652[R2=0.652）。為了解決泊松模型中的過度分散問題，我們制定了加權(quán)負(fù)二項(xiàng)式模型。盡管此模型的表現(xiàn)不如加權(quán)Poisson模型（[R2=?0.638 ），則在進(jìn)行推理時(shí)可能會(huì)更好。
此后，我們嘗試通過使用Hmisc包估算缺失值來進(jìn)一步改進(jìn)模型。盡管生成的模型比初始OLS模型要好，但是它們沒有獲得比以前更高的性能（[R2=?0.627[R2=0.627）。
那么，最好的模型到底是什么？就模型假設(shè)的正確性而言，這是加權(quán)負(fù)二項(xiàng)式模型。就決定系數(shù)而言，[R2[R2，這是加權(quán)Poisson回歸模型。因此，出于預(yù)測臭氧水平的目的，我將選擇加權(quán)Poisson回歸模型。
您可能會(huì)問：所有這些工作值得嗎？實(shí)際上，初始模型和加權(quán)泊松模型的預(yù)測在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
當(dāng)我們 比較它們時(shí)，模型之間的差異變得顯而易見：
?總之，我們從預(yù)測負(fù)值和低估高臭氧水平的模型（左側(cè)顯示OLS模型）到?jīng)]有此類明顯缺陷的模型（右側(cè)加權(quán)Poisson模型）。
?
