如果我的線性回歸數據包含幾個混合的線性關係怎麼辦?
假設我正在研究水仙花對各種土壤條件的反應。我收集了關於土壤 pH 值與水仙花成熟高度的數據。我期待一個線性關係,所以我開始運行線性回歸。
然而,當我開始研究時,我並沒有意識到人口中實際上包含兩種水仙花,每種水仙花對土壤 pH 值的反應都非常不同。所以該圖包含兩個不同的線性關係:
當然,我可以觀察它並手動將其分開。但我想知道是否有更嚴格的方法。
問題:
- 是否有統計測試來確定數據集是否更適合單行或 N 行?
- 我將如何運行線性回歸來擬合 N 條線?換句話說,我如何解開混合混合的數據?
我可以想到一些組合方法,但它們的計算成本似乎很高。
說明:
- 在收集數據時,尚不清楚是否存在兩個品種。每個水仙花的品種沒有被觀察,沒有被記錄,也沒有被記錄。
- 無法恢復此信息。自收集數據以來,水仙花已經死亡。
我的印像是這個問題類似於應用聚類算法,因為你幾乎需要在開始之前知道聚類的數量。我相信對於任何數據集,增加行數都會減少總 rms 誤差。在極端情況下,您可以將數據集劃分為任意對,並在每對之間簡單地畫一條線。(例如,如果您有 1000 個數據點,您可以將它們分成 500 個任意對並在每對之間畫一條線。)擬合將是精確的,rms 誤差將完全為零。但這不是我們想要的。我們想要“正確”的行數。
如果我們假設您擁有不同品種的標籤,我認為 Demetri 的答案是一個很好的答案。當我閱讀您的問題時,對我來說似乎並非如此。我們可以使用基於 EM 算法的方法來基本擬合 Demetri 建議的模型,但不知道品種的標籤。幸運的是,R 中的 mixtools 包為我們提供了這個功能。由於您的數據是完全分離的,而且您似乎有相當多的數據,因此應該相當成功。
library(mixtools) # Generate some fake data that looks kind of like yours n1 <- 150 ph1 = runif(n1, 5.1, 7.8) y1 <- 41.55 + 5.185*ph1 + rnorm(n1, 0, .25) n2 <- 150 ph2 <- runif(n2, 5.3, 8) y2 <- 65.14 + 1.48148*ph2 + rnorm(n2, 0, 0.25) # There are definitely better ways to do all of this but oh well dat <- data.frame(ph = c(ph1, ph2), y = c(y1, y2), group = rep(c(1,2), times = c(n1, n2))) # Looks about right plot(dat$ph, dat$y) # Fit the regression. One line for each component. This defaults # to assuming there are two underlying groups/components in the data out <- regmixEM(y = dat$y, x = dat$ph, addintercept = T)
我們可以檢查結果
> summary(out) summary of regmixEM object: comp 1 comp 2 lambda 0.497393 0.502607 sigma 0.248649 0.231388 beta1 64.655578 41.514342 beta2 1.557906 5.190076 loglik at estimate: -182.4186
所以它適合兩個回歸,它估計 49.7% 的觀察結果落入組件 1 的回歸中,50.2% 落入組件 2 的回歸中。我模擬數據的方式是 50-50 拆分,所以這很好。
我用於模擬的“真實”值應該給出以下幾行:
y = 41.55 + 5.185ph 和 y = 65.14 + 1.48148ph
(我從您的圖中“手動”估計,以便我創建的數據看起來與您的相似)以及 EM 算法在這種情況下給出的行是:
y = 41.514 + 5.19ph 和 y = 64.655 + 1.55ph
非常接近實際值。
我們可以將擬合線與數據一起繪製
plot(dat$ph, dat$y, xlab = "Soil Ph", ylab = "Flower Height (cm)") abline(out$beta[,1], col = "blue") # plot the first fitted line abline(out$beta[,2], col = "red") # plot the second fitted line