Regression

如果我的線性回歸數據包含幾個混合的線性關係怎麼辦?

  • January 28, 2019

假設我正在研究水仙花對各種土壤條件的反應。我收集了關於土壤 pH 值與水仙花成熟高度的數據。我期待一個線性關係,所以我開始運行線性回歸。

然而,當我開始研究時,我並沒有意識到人口中實際上包含兩種水仙花,每種水仙花對土壤 pH 值的反應都非常不同。所以該圖包含兩個不同的線性關係:

土壤 pH 值與花高 (cm)

當然,我可以觀察它並手動將其分開。但我想知道是否有更嚴格的方法。

問題:

  1. 是否有統計測試來確定數據集是否更適合單行或 N 行?
  2. 我將如何運行線性回歸來擬合 N 條線?換句話說,我如何解開混合混合的數據?

我可以想到一些組合方法,但它們的計算成本似乎很高。


說明:

  1. 在收集數據時,尚不清楚是否存在兩個品種。每個水仙花的品種沒有被觀察,沒有被記錄,也沒有被記錄。
  2. 無法恢復此信息。自收集數據以來,水仙花已經死亡。

我的印像是這個問題類似於應用聚類算法,因為你幾乎需要在開始之前知道聚類的數量。我相信對於任何數據集,增加行數都會減少總 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

通過 EM 擬合線

引用自:https://stats.stackexchange.com/questions/389545

comments powered by Disqus