我正在尝试模拟与我拥有的经验数据相匹配的数据集，但不确定如何估计原始数据中的错误。经验数据包括异方差性，但我对将其转换掉不感兴趣，而是使用带有误差项的线性模型来重现经验数据的模拟。
例如，假设我有一些经验数据集和模型：
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，我们用线性函数对其进行近似。如果您确实先验地知道真实的数据生成过程（在这种情况下，因为您模拟了原始数据），那么您不妨使用它。您可以决定此处的近似值是否足以满足您的目的。然而，我们通常不知道真实的数据生成过程，并且基于奥卡姆剃刀，使用最简单的函数来充分拟合我们给出的可用信息量的数据。如果您愿意，也可以尝试样条曲线或更高级的方法。双变量分布看起来与我相当相似，
您需要对异方差进行建模。一种方法是通过 R 包 (CRAN) dglm，即色散广义线性模型。这是 glm 的扩展，除了通常的 之外glm，它还适合第二个 glm，用于从第一个 glm 的残差中分散。我没有使用此类模型的经验，但它们似乎很有希望......这是一些代码：
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)
library(dglm) ### double glm's
mod2 <- dglm(y ~ n, ~ n, gaussian,ykeep=TRUE,xkeep=TRUE,zkeep=TRUE)
### This uses log link for the dispersion part, should also try identity link ..
y2 <- simulate(mod2)
plot(n, y2$sim_1)
mod3 <- dglm(y ~ n, ~ n, gaussian, dlink="identity", ykeep=TRUE,xkeep=TRUE,zkeep=TRUE) ### This do not work because it leads to negative weights!
模拟图如下所示：
该图确实看起来模拟使用了估计的方差，但我不确定，因为模拟（）函数没有用于 dglm 的方法......
（另一种可能的研究方法是使用R包gamlss，它使用另一种方法将方差建模为协变量的函数。）