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?
źródło
Odpowiedzi:
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).
Zauważ, że
R
s ? 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: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):
Jeśli jesteśmy z tego zadowoleni, możemy teraz wykorzystać ten proces jako dodatek do symulacji danych.
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,
źródło
Musisz modelować heteroskedastyczność. Jednym podejściem jest pakiet R (CRAN)
dglm
, uogólniony model dyspersyjny. Jest to rozszerzenie glm, które, oprócz zwykłegoglm
, 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: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
R
pakietugamlss
, który wykorzystuje inne podejście do modelowania wariancji jako funkcji zmiennych zmiennych).źródło