全文鏈接：http://tecdat.cn/?p=29277

原文出處：拓端數(shù)據(jù)部落公眾號

一些標準的圖形工具可以極大地幫助理解數(shù)據(jù)集并評估所建議模型的質量。

學生考試成績

例如，數(shù)據(jù)集包含600個觀察結果，用于國家統(tǒng)計教育中心對學生進行的一項非常大的研究。數(shù)據(jù)集中的一些變量包括：

?性別：性別男性或女性。

?種族：種族或民族，具有西班牙裔，亞洲人，非洲裔美國人，白人的水平。

?學校類型，公立和私立。

?軌跡：控制位點，一個連續(xù)的協(xié)變量，指示受試者對影響他們的事件的自我感知控制程度（更高=更感知的控制）

?概念：自我概念，反映受試者自我感知的學術能力的連續(xù)協(xié)變量（更高=更多的自我感知能力）。

?動機，一個連續(xù)的協(xié)變量（更高=更多的自我感知的動機）。

?數(shù)學：標準化的數(shù)學考試成績。

? 閱讀：標準化閱讀測試成績。

要基本了解數(shù)據(jù)集中的變量，我們可以使用摘要命令：

1. 


2. summary(hsb)

?

然而，由于變量太少，我們實際上可以用成對命令把所有的連續(xù)變量繪制在一起。

?

圖中顯示，雖然動機是一個連續(xù)變量，但它只具有四個類別。這很可能是因為學生被要求在4個點的量表上評估他們的自我激勵（0，0.33，0.66，和1）。另外，數(shù)學和任何一個連續(xù)預測因素之間最強的關聯(lián)似乎是與閱讀。

1. anova(M1, M2)

2. 


?

有兩個標準圖來檢查M1是否似乎違反了線性模型的假設。第一種是繪制殘差ei = yi - x ′ iβ?作為預測值的函數(shù)?yi = x ′ iβ?。除了通常的殘差之外，還有兩種類型的殘差經常被考慮。

1. 標準化殘差：ri = ei σ?。這些殘差是單位的

2. 研究性殘差：di = ei σ? √ 1 - hi?

這樣，如果模型是真實的，殘差大約是N（0，1）。

圖1顯示了殘差與M1的擬合值之間的關系，使用的是通常的殘差和 studentized殘差的一個版本，ei/ √ 1 - hi 。其他統(tǒng)計學家更喜歡在標準化的尺度上繪制殘差（除以σ?），這樣Y軸上的單位對每張殘差圖都有相同的解釋。

1. 


2. plot(predict(M1), res, pch = 21, bg = "black", cex = cex, cex.axis = cex,

3. ex)

?由于σ2（1-hi）是一個方差，我們必須有（1-hi）>0=?hi<1=?（1-hi）<1。因此，studentized殘差ei/ √ 1 - hi在規(guī)模上大于原始的ei。

通過繪制殘差與理論上的正態(tài)分布的直方圖，可以更容易地看到這一點，如圖所示。

1. 


2. hist(res/sigma.hat, breaks = 50, freq = FALSE, cex.axis = cex,

??對于這個殘差圖，我個人更喜歡使用標準化的殘差，即除以?σ。這是因為我更容易在N（0，1）上發(fā)現(xiàn)大的殘差值（標準正態(tài)在±2之間有95%的機會，在±3之間有99.7%的機會）。圖2證實了殘差略微偏向負值，表明需要進一步建模。

線性預測模型

回顧一下滿意度數(shù)據(jù)集，即病人自我報告對在醫(yī)院接受的護理質量的滿意度，由年齡、壓力和病情的嚴重程度來預測。下面是潛在模型之間的一些比較。

1. 


2. anova(lm(Satisfaction ~ Age

?

anova(lm(Satisfaction ~ Age, data = satis),

?

anova(lm(Satisfaction ~ Age + Severity, data = satis),

?

?

根據(jù)上面的測試結果，我們可以決定考慮使用年齡和壓力的模型，因為它的F統(tǒng)計量比年齡和嚴重程度稍顯重要。然而，這種方法會增加選擇錯誤模型的機會，正如下面的模擬實驗所顯示的那樣。

?

1. 


2. pval <- 2*pt(abs(beta.hat)/sqrt(diag(vcov(M))[-1]),

3. df = n-3, lower.tail = FALSE) # pvalues

4. sim.out[ii,] <- round(c(beta.hat[1], pval[1], beta.hat[2], pval[2]),2)

5. }

6. sim.out

?

在10次實驗運行中，只有兩次對兩個協(xié)變量都有顯著的P值，但對這些協(xié)變量的估計卻有很大偏差。事實上，其中兩個復制甚至把協(xié)變量的符號弄錯了。雖然這個例子并不現(xiàn)實（協(xié)變量之間的相關性幾乎不可能達到0.999），但它說明了這樣一個事實：當協(xié)變量高度相關時，回歸很難弄清y的變化是由一個協(xié)變量還是另一個協(xié)變量引起的。換句話說，β?的估計值從一個隨機樣本到另一個隨機樣本會有很大的變化。這種現(xiàn)象被稱為方差膨脹

衡量預測因子之間的共線性的一個簡單方法是所謂的方差膨脹因子（VIF）。

其中C是協(xié)變量之間的樣本相關矩陣。

下面的R代碼說明了VIF的兩種形式之間的平等。

1. VIF <- diag(solve(cor(X)))

2. VIF

?

1. 


2. names(satis.star) <- c("Satisfaction", paste0(names(satis)[-1], ".star"))

3. head(satis.star)

?

1. 


2. diag(vcov(M.orig))["Severity"]/diag(vcov(M.trans))["Severity.star"]

?

VIF["Severity"]

?

嚴重程度和壓力的方差膨脹系數(shù)為2.00，這兩個變量之間的相關度為0.67。相比之下，模擬實驗中的方差膨脹系數(shù)是1/（1-.9992）≈500。作為一條經驗法則，VIF大于10是值得關注的。

離群值

簡單地說，離群值是指與其他觀測值相比具有異常大的殘差。然而，為了提高模型的擬合度而刪除這些觀測值的情況比較少見，因為出錯的幾乎總是模型，而不是數(shù)據(jù)。也就是說，異常值通常是模型中遺漏了重要協(xié)變量的結果。例如，數(shù)據(jù)集包含了105個世界國家的嬰兒死亡率和人均收入的數(shù)據(jù)。數(shù)據(jù)集中的變量是

- 收入：以元計的人均收入。

- 嬰兒：每1000個活產的嬰兒死亡率。

- 地區(qū)：大陸，有非洲、美洲、亞洲、歐洲等級別。

- 石油：石油出口國，級別為無、有。?

對數(shù)尺度的線性模型適合于數(shù)據(jù)，殘差與擬合值繪制在圖3中。

1. 


2. coef(M) # 系數(shù)

3. 


1. 


2. as.character(formula(M))[c(2,1,3)]), collapse = " "), # 模型

3. xlab = "Predicted Log-Infant Mortality",

這張圖上的異常點是最上面的三角形，它對應的是沙特阿拉伯。目前還不清楚為什么沙特阿拉伯的嬰兒死亡率比模型預測的要高很多。一種可能性是，它的石油產量太高，以至于歪曲了平均收入。在這種情況下，我們有以下選擇。

1. 將沙特阿拉伯排除在分析之外，承認我們并不完全知道原因（不推薦）。

2. 將所有9個產油國排除在分析之外，并說明我們的預測只針對非產油國。

3.將沙特阿拉伯排除在外，并評估這對我們的分析有何影響。

樹木的針葉產生的陰涼面積

數(shù)據(jù)集包含以下關于35種針葉樹的變量。- 針葉面積。樹木的針葉產生的陰涼面積 - 高度：樹木的高度 - 樹干尺寸。樹的直徑 以下模型適合于該數(shù)據(jù)。

1. 


2. 


3. X <- model.matrix(M)

4. H <- X %*% solve(crossprod(X), t(X)) # HAT matrix

5. head(diag(H))

?

1. range(h - diag(H))

2. 


?

在這四種類型的殘差中，每一種都把大殘差放大得越來越多。為了衡量每個觀測點的總體影響，圖中的庫克距離與杠桿率作了對比

1. 


2. 


3. n <- nrow(lforest)

4. hbar <- p/n # 平均杠桿率

5. abline(v = 2*hbar, col = "gray60", lty = 2) # 2x平均杠桿率

其中一個觀測值的庫克距離幾乎是其他觀測值的3倍以上（紅色），而其中的e個觀測值的平均杠桿率是兩倍（藍色）。為了了解為什么這些觀測值有很高的杠桿率，我們可以看一下所有協(xié)變量的配對圖（在這種情況下，這只是高度和樹干尺寸）。

我們可以看到，其中三個標記點（有影響力的點和兩個高杠桿點）的TrunkSize值非常大（記得這些值是按對數(shù)比例繪制的）。剩下的那個高杠桿點的樹干尺寸相對于其高度來說非常小。至于為什么只有三棵樹中樹干尺寸最大的一棵被標記為離群點，我們可以將NeedleArea與每個協(xié)變量（包括實際值和預測值）作對比。為了進行比較，預測是在所有觀測值和省略一個觀測值的情況下進行的：要么是有影響力的觀測值，要么是有最高杠桿的觀測值。

1. # 用省略的值進行預測

2. omit.ind <- c(which(inf.ind), # 最有影響力的

3. which.max(h)) # 最高杠桿率

4. which(infl.ind)中的錯誤："which "的參數(shù)不符合邏輯

5. names(omit.ind) <- c("inf", "lev")

6. omit.ind

?

1. yobs <- lforest[, "NeedleArea"] # 觀察值

2. Mf <- lm(NeedleArea ~ Height + TrunkSize, data = lforest) # 所有觀測值

3. 


4. for(jj in 1:length(opt.ind)) {

5. # 用省略的觀察值建立模型

6. Mo <- lm(NeedleArea ~ Height + TrunkSize,

7. data = lforest, subset = -omit.ind[jj])

8. yfito <- predict(Mo, newdata = lforest) # 所有觀測值的擬合值

9. for(ii in c("Height", "TrunkSize")) {

我們可以看到，當移除有影響的觀測值時，有無遺漏的預測值之間的差異更加明顯。這在針孔面積與樹干尺寸的圖中最為明顯，這是對具有高杠桿作用的三個點和具有高影響的點的預測。在這個特殊的案例中，我們確定具有最大樹干尺寸的三棵樹的測量是不正確的，它們可以從分析中移除。

最受歡迎的見解

1.R語言多元Logistic邏輯回歸 應用案例

2.面板平滑轉移回歸(PSTR)分析案例實現(xiàn)

3.matlab中的偏最小二乘回歸（PLSR）和主成分回歸（PCR）

4.R語言泊松Poisson回歸模型分析案例

5.R語言回歸中的Hosmer-Lemeshow擬合優(yōu)度檢驗

6.r語言中對LASSO回歸，Ridge嶺回歸和Elastic Net模型實現(xiàn)

7.在R語言中實現(xiàn)Logistic邏輯回歸

8.python用線性回歸預測股票價格

9.R語言如何在生存分析與Cox回歸中計算IDI，NRI指標