Symuluj regresję liniową z heteroscedastycznością

Próbuję symulować zestaw danych, który pasuje do posiadanych danych empirycznych, ale nie jestem pewien, jak oszacować błędy w oryginalnych danych. Dane empiryczne obejmują heteroscedastyczność, ale nie jestem zainteresowany jej przekształceniem, ale raczej stosuję model liniowy ze składnikiem błędu do odtworzenia symulacji danych empirycznych.

Załóżmy na przykład, że mam jakiś empiryczny zestaw danych i model:

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)

za pomocą plot(n,y)otrzymujemy następujące.

Jeśli jednak spróbuję zasymulować dane, simulate(mod)heteroscedastyczność zostanie usunięta i nie zostanie przechwycona przez model.

Mogę użyć uogólnionego modelu najmniejszych kwadratów

VMat <- varFixed(~n)
mod2 = gls(y ~ n, weights = VMat)

zapewnia to lepsze dopasowanie modelu na podstawie AIC, ale nie wiem, jak symulować dane przy użyciu danych wyjściowych.

Moje pytanie brzmi: jak stworzyć model, który pozwoli mi symulować dane w celu dopasowania do oryginalnych danych empirycznych (n i y powyżej). W szczególności potrzebuję sposobu oszacowania sigma2, czyli błędu, przy użyciu albo przy użyciu modelu?

Aby symulować dane ze zmienną wariancją błędu, należy określić proces generowania danych dla wariancji błędu. Jak zauważono w komentarzach, zrobiłeś to podczas generowania oryginalnych danych. Jeśli masz rzeczywiste dane i chcesz tego spróbować, wystarczy zidentyfikować funkcję, która określa, w jaki sposób rezydualna wariancja zależy od zmiennych towarzyszących. Standardowym sposobem na to jest dopasowanie modelu, sprawdzenie, czy jest to uzasadnione (inne niż heteroscedastyczność) i zapisanie resztek. Te reszty stają się zmienną Y nowego modelu. Poniżej zrobiłem to dla twojego procesu generowania danych. (Nie widzę, gdzie ustawiłeś losowe ziarno, więc nie będą to dosłownie te same dane, ale powinny być podobne, i możesz odtworzyć moje za pomocą mojego ziarna).

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)

Zauważ, że Rs ? Plot.lm da ci wykres (por. Tutaj ) pierwiastka kwadratowego z bezwzględnych wartości reszt, pomocnie nałożony z dopasowaniem lowess, co jest właśnie tym, czego potrzebujesz. (Jeśli masz wiele zmiennych towarzyszących, możesz chcieć to ocenić osobno dla każdej zmiennej towarzyszącej). Jest najmniejszy ślad krzywej, ale wygląda na to, że linia prosta dobrze dopasowuje dane. Dopasujmy więc wyraźnie ten model:

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)

Nie musimy się obawiać, że wariancja rezydualna wydaje się również zwiększać na wykresie lokalizacji skali dla tego modelu - to w zasadzie musi się zdarzyć. Znowu jest najdelikatniejszy ślad krzywej, więc możemy spróbować dopasować kwadrat do kwadratu i sprawdzić, czy to pomaga (ale nie pomaga):

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

Jeśli jesteśmy z tego zadowoleni, możemy teraz wykorzystać ten proces jako dodatek do symulacji danych.

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

Należy pamiętać, że proces ten nie gwarantuje dokładniejszego znalezienia prawdziwego procesu generowania danych niż jakakolwiek inna metoda statystyczna. Użyłeś funkcji nieliniowej do wygenerowania błędów SD, a my przybliżyliśmy ją funkcją liniową. Jeśli faktycznie znasz prawdziwy proces generowania danych a-priori (jak w tym przypadku, ponieważ symulowałeś oryginalne dane), równie dobrze możesz go użyć. Możesz zdecydować, czy przybliżenie tutaj jest wystarczające dla twoich celów. Zazwyczaj jednak nie znamy prawdziwego procesu generowania danych i na podstawie brzytwy Ockhama zastosowaliśmy najprostszą funkcję, która odpowiednio pasuje do danych, które podaliśmy, o ilości dostępnych informacji. Możesz również wypróbować splajny lub bardziej wyszukane podejścia, jeśli wolisz. Dwuwymiarowe rozkłady wyglądają dość podobnie do mnie,

Musisz modelować heteroskedastyczność. Jednym podejściem jest pakiet R (CRAN) dglm, uogólniony model dyspersyjny. Jest to rozszerzenie glm, które, oprócz zwykłego glm, pasuje do drugiego glm w celu zdyspergowania resztek z pierwszego glm. Nie mam doświadczenia z takimi modelami, ale wydają się obiecujące ... Oto kod:

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!

Symulowany wykres pokazano poniżej:

Wykres wygląda na to, że symulacja wykorzystała oszacowaną wariancję, ale nie jestem pewien, ponieważ funkcja symulacji () nie ma metod dla dglm ...

(Inną możliwością zbadania jest użycie Rpakietu gamlss, który wykorzystuje inne podejście do modelowania wariancji jako funkcji zmiennych zmiennych).