模擬具有異方差性的線性回歸

January 27, 2017

我正在嘗試模擬與我擁有的經驗數據相匹配的數據集，但不確定如何估計原始數據中的錯誤。經驗數據包括異方差性，但我對將其轉換掉不感興趣，而是使用帶有誤差項的線性模型來重現經驗數據的模擬。

例如，假設我有一些經驗數據集和模型：
n=rep(1:100,2)
a=0
b = 1
sigma2 = n^1.3
eps = rnorm(n,mean=0,sd=sqrt(sigma2))
y=a+b*n + eps
mod <- lm(y ~ n)
使用plot(n,y)我們得到以下信息。

但是，如果我嘗試模擬數據，simulate(mod)則異方差會被移除，並且不會被模型捕獲。

我可以使用廣義最小二乘模型
VMat <- varFixed(~n)
mod2 = gls(y ~ n, weights = VMat)
它提供了基於 AIC 的更好的模型擬合，但我不知道如何使用輸出來模擬數據。

我的問題是，如何創建一個模型，讓我能夠模擬數據以匹配原始的經驗數據（上面的 n 和 y）。具體來說，我需要一種使用模型來估計 sigma2 誤差的方法？

要模擬具有不同誤差方差的數據，您需要指定誤差方差的數據生成過程。正如評論中指出的那樣，您在生成原始數據時就這樣做了。如果您有真實數據並想嘗試這個，您只需要確定指定殘差如何取決於協變量的函數。執行此操作的標準方法是擬合您的模型，檢查它是否合理（異方差除外），並保存殘差。這些殘差成為新模型的 Y 變量。下面我已經為您的數據生成過程做到了這一點。（我看不到你在哪裡設置了隨機種子，所以這些數據實際上不是相同的數據，但應該是相似的，你可以使用我的種子精確地複制我的數據。）
set.seed(568)  # this makes the example exactly reproducible

n      = rep(1:100,2)
a      = 0
b      = 1
sigma2 = n^1.3
eps    = rnorm(n,mean=0,sd=sqrt(sigma2))
y      = a+b*n + eps
mod    = lm(y ~ n)
res    = residuals(mod)

windows()
 layout(matrix(1:2, nrow=2))
 plot(n,y)
 abline(coef(mod), col="red")
 plot(mod, which=3)
請注意，R? plot.lm將為您提供殘差絕對值的平方根的圖（參見此處），有助於覆蓋低擬合，這正是您所需要的。（如果您有多個協變量，您可能希望分別針對每個協變量進行評估。）有最輕微的曲線暗示，但這看起來像一條直線可以很好地擬合數據。因此，讓我們明確地擬合該模型：
res.mod = lm(sqrt(abs(res))~fitted(mod))
summary(res.mod)
# Call:
# lm(formula = sqrt(abs(res)) ~ fitted(mod))
# 
# Residuals:
# Min 1Q Median 3Q Max 
# -3.3912 -0.7640 0.0794 0.8764 3.2726 
# 
# Coefficients:
# Estimate Std. Error t value Pr(>|t|) 
# (Intercept) 1.669571 0.181361 9.206 < 2e-16 ***
# fitted(mod) 0.023558 0.003157 7.461 2.64e-12 ***
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 1.285 on 198 degrees of freedom
# Multiple R-squared: 0.2195, Adjusted R-squared: 0.2155 
# F-statistic: 55.67 on 1 and 198 DF, p-value: 2.641e-12
windows()
 layout(matrix(1:4, nrow=2, ncol=2, byrow=TRUE))
 plot(res.mod, which=1)
 plot(res.mod, which=2)
 plot(res.mod, which=3)
 plot(res.mod, which=5)
我們不必擔心這個模型的尺度位置圖中的剩餘方差似乎也在增加——這基本上是必鬚髮生的。再次有最輕微的曲線暗示，所以我們可以嘗試擬合一個平方項，看看是否有幫助（但沒有幫助）：
res.mod2 = lm(sqrt(abs(res))~poly(fitted(mod), 2))
summary(res.mod2)
# output omitted
anova(res.mod, res.mod2)
# Analysis of Variance Table
# 
# Model 1: sqrt(abs(res)) ~ fitted(mod)
# Model 2: sqrt(abs(res)) ~ poly(fitted(mod), 2)
# Res.Df RSS Df Sum of Sq F Pr(>F)
# 1 198 326.87 
# 2 197 326.85 1 0.011564 0.007 0.9336
如果我們對此感到滿意，我們現在可以將此過程用作附加組件來模擬數據。
set.seed(4396)  # this makes the example exactly reproducible
x = n
expected.y = coef(mod)[1] + coef(mod)[2]*x
sim.errors = rnorm(length(x), mean=0,
                  sd=(coef(res.mod)[1] + coef(res.mod)[2]*expected.y)^2)
observed.y = expected.y + sim.errors
請注意，與任何其他統計方法相比，此過程無法保證找到真實的數據生成過程。您使用非線性函數來生成誤差 SD，我們用線性函數對其進行近似。如果您確實先驗地知道真實的數據生成過程（在這種情況下，因為您模擬了原始數據），那麼您不妨使用它。您可以決定此處的近似值是否足以滿足您的目的。然而，我們通常不知道真實的數據生成過程，並且基於奧卡姆剃刀，使用最簡單的函數來充分擬合我們給出的可用信息量的數據。如果您願意，也可以嘗試樣條曲線或更高級的方法。雙變量分佈看起來與我相當相似，