Wizualizacja wyników modeli mieszanych

Jednym z problemów, które zawsze miałem z modelami mieszanymi, jest wymyślanie wizualizacji danych - takich, które mogłyby skończyć się na papierze lub plakacie - gdy tylko uzyska się wyniki.

Obecnie pracuję nad modelem efektów mieszanych Poissona z formułą, która wygląda mniej więcej tak:

a <- glmer(counts ~ X + Y + Time + (Y + Time | Site) + offset(log(people))

Z czymś dopasowanym w glm () można łatwo użyć predykcji (), aby uzyskać prognozy dla nowego zestawu danych i zbudować coś z tego. Ale przy takim wyniku - jak skonstruowałbyś coś w rodzaju wykresu prędkości w czasie z przesunięciami z X (i prawdopodobnie z ustaloną wartością Y)? Myślę, że można dobrze przewidzieć dopasowanie tylko na podstawie oszacowań efektów stałych, ale co z 95% CI?

Czy jest coś jeszcze, co ktoś może wymyślić, aby pomóc w wizualizacji wyników? Wyniki modelu są poniżej:

Random effects:
 Groups     Name        Variance   Std.Dev.  Corr          
 Site       (Intercept) 5.3678e-01 0.7326513               
            time        2.4173e-05 0.0049167  0.250        
            Y           4.9378e-05 0.0070270 -0.911  0.172 

Fixed effects:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept) -8.1679391  0.1479849  -55.19  < 2e-16
X            0.4130639  0.1013899    4.07 4.62e-05
time         0.0009053  0.0012980    0.70    0.486    
Y            0.0187977  0.0023531    7.99 1.37e-15

Correlation of Fixed Effects:
         (Intr) Y    time
X      -0.178              
time    0.387 -0.305       
Y      -0.589  0.009  0.085

r data-visualization mixed-model

— Fomite
źródło

(+1) @EpiGrad: Dlaczego martwisz się o CI (tj. O błąd standardowy) prognoz z części Twojego modelu o stałym efekcie?

— boscovich

@ andrea Odpowiedź intelektualna i praktyczna: pod względem intelektualnym ogólnie wolę kwantyfikację i wizualizację niepewności, kiedy tylko mogę. Praktycznie, bo jestem pewien, że recenzent poprosi o to.

— Fomite,

Tak, tak, jasne, ale miałem na myśli coś innego. Przepraszam, mój komentarz nie był wystarczająco jasny. W swoim pytaniu piszesz „ale co z 95% CI?”. Mój komentarz brzmi: dlaczego nie obliczasz błędu standardowego prognozy na podstawie części modelu o stałym efekcie? Jeśli jesteś w stanie obliczyć przewidywane wartości z części o stałym efekcie, to możesz również obliczyć SE, a tym samym CI. @EpiGrad

— boscovich

@andrea Ah. Problem polega na tym, że jedna z rzeczy, które chciałbym przewidzieć, czas, ma również losowy efekt, z którym nie mam pojęcia, co z tym zrobić.

— Fomite,

Cóż, chcesz do przewidzenia counts, nie time. Naprawić wartości X, Ya timei stosując stałe-efektów część modelu przewidzieć counts. To prawda, że timejest uwzględniony w twoim modelu również jako efekt losowy (podobnie jak punkt przechwytywania i Y), ale nie ma to tutaj znaczenia, ponieważ użycie tylko części modelu z efektem stałym do prognozowania jest jak ustawienie efektów losowych na 0 @EpiGrad

— boscovich

Prognozowanie countsprzy użyciu części modelu z efektami stałymi oznacza, że ustawiłeś zero (tj. Ich średnią) efektów losowych. Oznacza to, że można „zapomnieć” o nich i użyć standardowych maszyn do obliczenia prognoz i standardowych błędów prognoz (za pomocą których można obliczyć przedziały ufności).

To jest przykład użycia Staty, ale przypuszczam, że można ją łatwo „przetłumaczyć” na język R:

webuse epilepsy, clear
xtmepoisson seizures treat visit || subject: visit
predict log_seiz, xb
gen pred_seiz = exp(log_seiz)
predict std_log_seiz, stdp
gen ub = exp(log_seiz+invnorm(.975)*std_log_seiz)
gen lb = exp(log_seiz-invnorm(.975)*std_log_seiz)

tw (line pred_seiz ub lb visit if treat == 0, sort lc(black black black) ///
 lp(l - -)), scheme(s1mono) legend(off) ytitle("Predicted Seizures") ///
 xtitle("Visit")

Wykres odnosi się treat == 0i ma być przykładem ( visitnie jest to naprawdę zmienna ciągła, ale tylko po to, aby uzyskać pomysł). Linie przerywane to 95% przedziały ufności.

wprowadź opis zdjęcia tutaj

— Boskowicz
źródło