Symuluj regresję liniową z heteroscedastycznością

9

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. wprowadź opis zdjęcia tutaj

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?

użytkownik44796
źródło
1
Zatem model liniowy nie uchwyci warunkowej heteroskedastyczności, chyba że wyraźnie spróbuje to zrobić, używając jednego z kilku podejść. Standardowe techniki ekonometryczne dostosowują standardowe błędy parametrów, aby uwzględnić heteroskedastyczność, ale nie modelują tego jawnie.
generic_user
Masz rację. Próbuję użyć modelu liniowego do uchwycenia heterogeniczności. Myślę, że powinienem używać uogólnionego modelu najmniejszych kwadratów. Jeśli są jakieś inne rekomendacje, spróbuję je.
user44796
BŁĄD W TWOIM KODIE, MUSISZ UŻYWAĆ `lm (y ~ n)`
kjetil b halvorsen 27.01.17
1
Nie rozumiem twojego pytania, ponieważ kod spełnia dokładnie to, o co w tytule prosisz: symuluje regresję liniową z błędami heteroscedastycznymi. Czy pytasz o metody szacowania jakiegoś modelu dla heteroscedastyczności? Jeśli tak, to musisz podać model!
whuber
Mam nadzieję, że wyjaśniłem moje pytanie za pomocą edycji. W powyższym pytaniu n i y reprezentują dane empiryczne. Chcę dopasować model do danych, a następnie użyć tego modelu do wygenerowania danych symulowanych, które pasują do średniej i resztek oryginalnych danych.
user44796

Odpowiedzi:

9

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)

wprowadź opis zdjęcia tutaj

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)

wprowadź opis zdjęcia tutaj

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,

wprowadź opis zdjęcia tutaj

gung - Przywróć Monikę
źródło
Był to właściwie wniosek, do którego zacząłem dochodzić, ale nigdy nie uzyskałbym tak eleganckiej odpowiedzi.
user44796
5

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:

wprowadź opis zdjęcia tutaj

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).

kjetil b halvorsen
źródło
1
podwójnie uogólniony model liniowy wydaje się odpowiednio modelować oryginalne dane. Nie jestem pewien, w jaki sposób modelowany jest błąd resztkowy za pomocą przewidywania (). Będę musiał to sprawdzić.
user44796