StdDev eines Random Effect aus summary(lmerMod) aus gamm4

Modelle zur Korrelations- und Regressionsanalyse

Moderator: EDi

Antworten
wepos

StdDev eines Random Effect aus summary(lmerMod) aus gamm4

Beitrag von wepos »

Liebes Forum,

ich arbeite seit kurzem mit gam, gamm und gamm4 und bitte euch um Unterstützung beim Verständnis dieser Funktionen:

Ich habe ein multilineares Modell mit einer GAM komponente an einen Datensatz anpassen lassen:

Code: Alles auswählen

rawmod = gamm4(data = dfscaled, formula = bm ~ mnh_d_ar + sdi + bon_100 + maxh + api + s(DG, bs='cr', k = -1), random = ~(1|Tnr), REML=REML, weights = weig_lm_r)
rawmod ist eine Liste mit zwei modellen (1) einem lmerMod (modl im folgenden) und (2) einem gam.

Ich habe mir summary(modl) ausgeben lassen

Code: Alles auswählen

Random effects:
 Groups   Name        Variance      Std.Dev. 
 Tnr      (Intercept) 10376213.2738 3221.2130
 Xr       s(DG)       10092728.7380 3176.9055
 Residual                    0.7849    0.8859
Number of obs: 7784, groups:  Tnr, 4858; Xr, 8
Da die Modellresiduen heteroskedastisch sind, wollte ich gerne wissen, ob der Zufallseffekt der Gruppe Tnr auch mit dem geschätzten Erwartungswert variiert. Dafür wollte ich mir die Varianz dieses Effektes gegen die Vorhersage des Modelles betrachten.

Dafür habe ich mir die Werte des Effektes erzeugt wie folgt:

Code: Alles auswählen

p_l_ti0 = predict(modl, re.form = ~0) # no random effects 
p_l_tnr = predict(modl, re.form = ~(1|Tnr)) # prediction on level of group Tnr
zu meiner Überraschung ergibt

Code: Alles auswählen

sd(p_l_tnr - p_l_ti0)
eine Standardabweichung des Effektes von 1330.03. Das ist viel weniger als unter Std.Dev von Summary ausgegeben. Offenbar gehe ich an diesem Punkt bereits von einer falschen Vorstellung aus: Weshalb sind nicht beide identisch?

Vielen Dank
Werner
Benutzeravatar
EDi
Beiträge: 1599
Registriert: Sa Okt 08, 2016 3:39 pm

Re: StdDev eines Random Effect aus summary(lmerMod) aus gamm4

Beitrag von EDi »

Wieso sollten diese beiden Sachen gleich sein?

Die RE-SD ist die SD im Intercept.
Was dein anderer Ausdruck sein soll weiß ich nicht (außer das es die SD der differenz zwischen den predictions einmal mit RE und einmal mit RE=0).
Wie die beiden verbunden sind sehe ich derzeit (noch) nicht - ist aber auch gerade spät....
Bitte immer ein reproduzierbares Minimalbeispiel angeben. Meinungen gehören mir und geben nicht die meines Brötchengebers wieder.

Dieser Beitrag ist lizensiert unter einer CC BY 4.0 Lizenz
Bild.
wepos

Re: StdDev eines Random Effect aus summary(lmerMod) aus gamm4

Beitrag von wepos »

Hallo EDi,

ich würde gerne angeben, welche Abweichung der Anwender des Modelles maximal zu erwarten hat, wenn er das Modell mit einem beliebigen Wert von Tnr anwendet. Das ist zwar für die eigentliche Zielstellung meiner Arbeit unerheblich, scheint aber in meinem Bereich ein Standard zu sein, der zu erfüllen ist.

Wenn man dann also ein gemischtes Modell anwendet

yij = Summe(aij * xij) + bj + eij

wobei Indizes i und j stehen für Index der Beobachtung bzw. der Gruppe b und

bj der Effekt der Gruppe b ist und eij die Residuenverteilung mit bij~V1 und eij~V2
dann werden in der Tabelle mit den Modellparametern auch die Parameter der Verteilungen V1 und V2 erwartet, bei Normalveilung (~N(0, sigma)) dann sigma.

Das macht (oder?) wohl nur Sinn, wenn die Verteilung des Gruppeneffektes homoskedastisch ist, sonst hängt diese Abweichung (s.o.) ja von dem geschätzten Erwartungswert ab. Außerdem würde ich gerne sehen, ob die Verteilung des Effektes eine Gaussverteilung ist, sonst reicht sigma als Verteilungsparameter nicht.

Also wollte ich mir ein Sample von Gruppeneffekten erzeugen als {y^i - y^ij}
also der Differenz zwischem bedingten Erwartungswert ohne Zufallseffekte und bedingtem Erwartungswert wenn Effekt der Gruppe bj berücksichtigt

und habe (vielleicht irrtümlich) erwartet dass die Standardabweichung dieser Menge der des Intercept unter Voraussetzung dass Gruppe Tnr berücksichtigt ist entsprechen müsste.

Ciao
Werner
Benutzeravatar
EDi
Beiträge: 1599
Registriert: Sa Okt 08, 2016 3:39 pm

Re: StdDev eines Random Effect aus summary(lmerMod) aus gamm4

Beitrag von EDi »

Außerdem würde ich gerne sehen, ob die Verteilung des Effektes eine Gaussverteilung ist, sonst reicht sigma als Verteilungsparameter nicht.
Reicht da nicht einfach ein QQ-Plot / histogramm deiner Gruppeneffekte?
Analog, kann man auch die homoskedastizizäz graphisch überprüfen.

Wenn es nicht N(0, sigma) ist, müsste du wohl weg von nlme/lme4 (hin zu STAN oder BUGS, wo das vermutlich einfacher ist) .
Bitte immer ein reproduzierbares Minimalbeispiel angeben. Meinungen gehören mir und geben nicht die meines Brötchengebers wieder.

Dieser Beitrag ist lizensiert unter einer CC BY 4.0 Lizenz
Bild.
Antworten