Blandede modeller

En blandet modell er en statistisk modell som inkluderer både faste og tilfeldige effekter. Konvensjonelle lineære regresjonsmodeller antar at observasjonene er uavhengige fra hverandre og identisk fordelte, mens blandede modeller tillater korrelasjon innad i klynger. Blandede modeller håndterer manglende data på en robust måte, og foretrekkes derfor ofte fremfor tradisjonelle metoder som ANOVA for repeterte målinger. Det kan være vanskelig å trekke en tydelig skillelinje mellom hva som regnes som faste og tilfeldige effekter, og i praksis kan dette avhenge av forskningsspørsmålet. En vanlig brukt definisjon er at faste effekter er karakteristikker ved en hel populasjon som en gjør vil inferens på, mens nivåene til en tilfeldig effekt ses på som et tilfeldig utvalg fra en større populasjon av verdier.

Motivasjon

Effekter på klyngenivå kan i prinsippet også estimeres ved bruk av faste effekter i en konvensjonell lineær modell, men denne framgangsmåten har flere svakheter. En slik modell vil for det første estimere parametre for hver klynge, som sannsynligvis ikke er interessante - både fordi de ikke vil kunne generaliseres til å si noe om klynger som ikke var med i studien, og fordi det potensielt store antallet parametre "forbruker" frihetsgrader i modellen. I motsetning til dette vil en blandet modell kun estimere varianskomponentene, som, avhengig av kovariansstruktur og antall klynger, kan være betraktelig færre. Bruken av faste effekter for å estimere klyngeffekter vil ofte også møte problemer med identifiserbarhet - hvis klyngene er konfundert med en kategorisk variabel (individer i en studie har for eksempel kun fått en type behandling), er det mulig å legge til og trekke fra et vilkårlig tall fra regresjonskoeffisientene og ende opp med identiske prediksjoner. For å unngå denne situasjonen, må man introdusere noen avgrensninger (constraints).

Definisjon

En vanlig parametrisering av en blandet modell med matrisenotasjon er

$\mathbf {y} =\mathbf {X} {\boldsymbol {\beta }}+\mathbf {Zb} +{\boldsymbol {\epsilon }}$

Der

• $\mathbf {y}$ er den observerte responsvariabelen, med forventning $E(\mathbf {y} )=\mathbf {X} {\boldsymbol {\beta }}$
• $\mathbf {X}$ er designmatrisen for de faste effektene
• ${\boldsymbol {\beta }}$ er en vektor med faste regresjonskoeffisientene
• $\mathbf {Z}$ er designmatrisen for de tilfeldige effektene
• $\mathbf {b}$ er en vektor med de tilfeldige koeffisientene, med forventning $E(\mathbf {b} )=\mathbf {0}$, og varians-kovarians-matrise $\mathrm {var} (\mathbf {b} )={\boldsymbol {\psi _{\theta }}}$
• ${\boldsymbol {\epsilon }}$ er en vektor med residualene, med forventning $E({\boldsymbol {\epsilon }})=\mathbf {0}$, og varians-matrise $\mathrm {var} ({\boldsymbol {\epsilon }})=\mathbf {\Lambda } \sigma ^{2}$

Som et alternativ til å skrive modellen som én likning på denne måten, kan en også skrive modellen på klyngenivå, der indeks $i$ refererer til klyngen, og hver klynge har $n_{i}$ observasjoner:

$\mathbf {y} _{i}=\mathbf {X} _{i}{\boldsymbol {\beta }}+\mathbf {Z} _{i}\mathbf {b} _{i}+{\boldsymbol {\epsilon }}_{i}$

Dette fordi observasjoner på tvers av $i$ skal være uavhengige, og fordi $\mathbf {Z} =\mathrm {diag} (\mathbf {Z_{1},...,Z_{N}} )$ og dermed kan bli svært stor.

Korrelasjon og varians

Et mål på hvor korrelerte observasjonene innad i en klynge er, er intraklasse-korrelasjon. For en enkel blandet modell med tilfeldige konstantledd, er observasjoner $j$ fra den samme klyngen $i$ korrelert med følgende korrelasjonskoeffisient:

$\rho ={\frac {var(b_{i})}{var(b_{i}+\epsilon _{i,j})}}={\frac {\sigma _{b}^{2}}{\sigma _{b}^{2}+\sigma ^{2}}}$

Der

• $\sigma _{b}^{2}$ refererer til variansen mellom ulike klynger
• $\sigma ^{2}$ refererer til variansen innad i klyngene

Altså øker korrelasjonen innad i klyngene med økende variasjon mellom klyngene. Hvis $\sigma _{b}^{2}=0$, er det ingen variasjon mellom ulike klynger og modellen tilsvarer en enkel lineær regresjonsmodell. Typisk er en mer interessert i variansen til en tilfeldige effekten på tvers av klynger, enn størrelsen/retningen på de enkelte klyngeeffektene. Er variansen mellom klynger stor for en fast effekt, kan det tyde på at det er riktig å inkludere en tilfeldig effekt. Hvis korrelasjonen $\rho$ er stor, vil det også bety at den effektive utvalgsstørrelsen er mindre enn $\sum _{i}n_{i}$.

Mer generelt kan korrelasjonen mellom to observasjoner $M1$ og $M2$ i samme klynge skrives som

$\rho ={\frac {{\hat {\mathrm {Cov} }}(\epsilon _{M1},\epsilon _{M2})}{{\sqrt {{\hat {\mathrm {Var} }}(\epsilon _{M1})}}{\sqrt {{\hat {\mathrm {Var} }}(\epsilon _{M2})}}}}$

Både ${\boldsymbol {\psi _{\theta }}}$ og ${\boldsymbol {\Lambda }}$ kan bestemmes for å håndtere ulike korrelasjonsstrukturer i data. Dette vil si at det er en struktur i korrelasjonen mellom observasjonene i en klynge, for eksempel at observasjoner som ligger nært hverandre i tid eller rom er mer korrelerte enn observasjoner som er fjernt fra hverandre. En ustrukturert korrelasjonsmatrise legger ingen begrensninger på parametrene, men øker dermed antallet parametre som må estimeres. Spesielt for ${\boldsymbol {\Lambda }}$ brukes derfor ofte en av følgende mer avgrensede korrelasjonsmatriser:

• Diagonal: Impliserer at det ikke er noen korrelasjon mellom residualene i samme klynge.

$\sigma ^{2}{\begin{pmatrix}1&0&0&\cdots &0\\0&1&0&\cdots &0\\0&0&1&\cdots &0\\\vdots &\vdots &\vdots &\ddots &\vdots \\0&0&\cdots &\cdots &1\end{pmatrix}}$

• Sammensatt symmetrisk ("compound symmetric"): Impliserer at det er en konstant korrelasjon mellom residualene.

$\sigma ^{2}{\begin{pmatrix}1&\rho &\rho &\cdots &\rho \\\rho &1&\rho &\cdots &\rho \\\rho &\rho &1&\cdots &\rho \\\vdots &\vdots &\vdots &\ddots &\vdots \\\rho &\rho &\cdots &\cdots &1\end{pmatrix}}$

• Ulike former autoregressive (AR) strukturer, der korrelasjon modelleres som en funksjon av avstanden mellom observasjonene. Denne strukturen antar at avstanden mellom naboobservasjoner er konstant. En første-ordens (AR(1)) matrise vil se slik ut:

$\sigma ^{2}{\begin{pmatrix}1&\rho &\rho ^{2}&\cdots &\rho ^{p}\\\rho &1&\rho &\cdots &\rho ^{p-1}\\\rho ^{2}&\rho &1&\cdots &\rho ^{p-2}\\\vdots &\vdots &\vdots &\ddots &\vdots \\\rho ^{p}&\rho ^{p-1}&\cdots &\cdots &1\end{pmatrix}}$

Estimering

Maximum likelihood

En blandet modell kan forenkles ved å omformuleres til en marginal modell

$y=\mathbf {X} {\boldsymbol {\beta }}+\mathbf {e} ,\mathbf {e} \sim N(0,\mathbf {\Sigma } _{\theta }\sigma ^{2})$

Der

• $\mathbf {e} =\mathbf {Zb} +{\boldsymbol {\epsilon }}$
• $\mathbf {\Sigma } _{\theta }=\mathbf {Z} {\boldsymbol {\psi _{\theta }}}\mathbf {Z} ^{T}/\sigma ^{2}+\mathbf {I}$.

Siden parameter-vektoren ${\boldsymbol {\theta }}$ ikke er kjent, må en bruke Maximum likelihood for å estimere denne. Likelihood-funksjonen er

$L({\boldsymbol {\beta }},{\boldsymbol {\theta }},\sigma ^{2})={\frac {1}{\sqrt {(2\pi \sigma ^{2})^{n}|\mathbf {\Sigma _{\theta }} |}}}\mathrm {exp} \left[-(\mathbf {y} -\mathbf {X} {\boldsymbol {\beta }})^{T}\mathbf {\Sigma _{\theta }} ^{-1}(\mathbf {y} -\mathbf {X} {\boldsymbol {\beta }})/2\sigma ^{2}\right]$

En forenkling som gjør utregningen raskere, er å i stedet bruke profil-likelihooden. Dette er en generell metode som fungerer ved å partisjonere vektoren med ukjente parametre ${\boldsymbol {\alpha }}$ i to deler $({\boldsymbol {\psi }},{\boldsymbol {\lambda }})$, og så velge en verdi for $\psi$ og maksimere likelihooden med hensyn på $\lambda$. Dette er en spesielt passende metode å bruke i denne sammenhengen, fordi maximum likelihood-estimatene for ${\boldsymbol {\beta }}$ og $\sigma ^{2}$ er enkle å regne ut hvis en kjenner ${\boldsymbol {\theta }}$. Disse kan dermed defineres som implisitte funksjoner av ${\boldsymbol {\theta }}$, slik at profil-likelihooden blir

$L_{p}({\boldsymbol {\theta }})={\frac {1}{\sqrt {(2\pi {\hat {\sigma }}_{\theta }^{2})^{n}|\mathbf {\Sigma _{\theta }} |}}}\mathrm {exp} \left[-(\mathbf {y} -\mathbf {X} {\hat {\boldsymbol {\beta }}}_{\theta })^{T}\mathbf {\Sigma _{\theta }} ^{-1}(\mathbf {y} -\mathbf {X} {\hat {\boldsymbol {\beta }}}_{\theta })/2{\hat {\sigma }}_{\theta }^{2}\right]$

der ${\hat {\boldsymbol {\beta }}}_{\theta }$ maximum likelihood-estimatet for ${\boldsymbol {\beta }}$ for en gitt ${\boldsymbol {\theta }}$, og ${\hat {\sigma }}^{2}$ er det tilhørende estimatet på $\sigma ^{2}$. En kan dermed behandle $L_{p}$ som en funksjon som kun avhenger av ${\boldsymbol {\theta }}$ og bruke numerisk optimering for å finne verdien av ${\boldsymbol {\theta }}$ som maksimerer $L_{p}$. Vanligvis brukes EM-algoritmen, Newton-Raphson-metoden eller en kombinasjon av disse. EM-algoritmen fungerer ved å ta utgangspunkt i sannsynlighetsfordelingen til $({\boldsymbol {\beta }}|\mathbf {y} )$ for det nåværende estimatet for ${\boldsymbol {\theta }}$, for å deretter finne den verdien av ${\boldsymbol {\theta }}$ som maksimerer forventningen til log-likelihooden. Newton-Raphson-metoden fungerer ved å bruke gradienten til log-likelihooden for å finne neste estimat for theta.

Begrenset sannsynlighetsmaksimering

Begrenset sannsynlighetsmaksimering (restricted maximum likelihood): Et generelt problem for maximum likelihood-estimatorer av varianskomponentene, er at de er forventningsskjeve, som kan ses ved uttrykket

$E({\hat {\sigma }}_{\scriptscriptstyle ML}^{2})={\frac {n-p}{n}}\sigma ^{2}$

Dette betyr at de vil undervurdere variansen ved økende antall parameter $p$. For å komme seg rundt dette, brukes det ofte en alternativ metode kalt REML (REstriced Maximum Likelihood). Denne finner et skalert gjennomsnitt av likelihood-funksjonen over alle verdier av ${\boldsymbol {\beta }}$. For balanserte data gir REML forventningsrette estimat for $\sigma ^{2}$ og ${\boldsymbol {\psi }}_{\theta }$.

$L_{R}({\boldsymbol {\theta }},\sigma ^{2})=\int L({\boldsymbol {\beta }},{\boldsymbol {\theta }},\sigma ^{2})d{\boldsymbol {\beta }}$

Etter en har funnet varians-komponentene i ${\boldsymbol {\theta }},\sigma ^{2}$, kan en estimere de faste parametrene i ${\boldsymbol {\beta }}$.

For store utvalg, gir ML og REML tilnærmet like resultater. For små utvalg med balanserte data, gir REML bedre estimat fordi den er forventningsrett. For små utvalg med ubalanserte data, vil både ML og REML gi upresise estimat.

Hypotesetesting/Inferens

Hypotesetesting og inferens foregår vanligvis for de faste effektene, og er presise for balanserte design. En lineær blandet modell er balansert hvis alle $n_{i}$ er like, og alle designmatrisene $\mathbf {Z} _{i}$ er like. En modell med tilfeldige koeffisienter er balansert hvis i tillegg alle designmatrisene $\mathbf {X} _{i}$ er like. Dette fordi estimatene ${\boldsymbol {\beta }}_{ML}$ og $\sigma _{ML}^{2}$ da er lik estimatene ${\boldsymbol {\beta }}_{OLS}$ og $\sigma _{OLS}^{2}$.

Frihetsgrader

For T- og F-test må en regne ut frihetsgrader for å tolke observatoren. Dette blir fort komplisert for ubalanserte design, generaliserte lineære blandede modeller (GLMM) eller kryssede tilfeldige effekter, men for balanserte, hierarkiske modeller er det mulig å regne ut nøyaktig:

$denDF_{i}=m_{i}-(m_{i-1}+p_{i}),\;i=1,...,Q+1$

• Der $denDF_{i}$ er "denominator degrees of freedom", eller nevner-frihetsgrader
• $m_{i}$ er det totale antall grupper på nivå $i$. Inneholder modellen et konstantledd er $m_{0}=1$, hvis ikke er $m_{0}=0$, mens $m_{Q+1}=N$.
• $m_{i-1}$ er antallet grupper på nivået over
• $p_{i}$ er summen av "numerator degrees of freedom"/teller-frihetsgrader for koeffisienter estimert på nivå $i$.

Formelen definerer nivået en koeffisient estimeres på avhengig av om den er "innenfor" eller "utenfor" det aktuelle nivået. En variabel er innenfor hvis verdien kan forandre seg innenfor et gitt nivå av klyngefaktoren, og utenfor hvis verdien ikke kan forandre seg. Et eksempel kan illustrere formelen: Hvis man forholder seg til et datasett med $N=100$ observasjoner, $Q=2$ nivå, $p_{0}=3$ faste koeffisienter estimert på nivå $i=2$ og $m_{1}=7$ grupper, blir da frihetsgradene $denDF_{2}=100-(7+3)=90$.

T-test

T-testen tester hypotesen $H_{0}:\beta =0$ mot den alternative hypotesen $H_{A}:\beta \neq 0$. T-observatoren er definert som $T={\frac {\hat {\beta }}{SE({\hat {\beta }})}}$ Observatoren følger da Students T-fordeling med frihetsgrader bestemt av nivået av klyngevariabelen der regresjonskoeffisienten er estimert.

F-test

F-testen tar utgangspunkt i lineære hypoteser $H_{0}:\mathbf {C} {\boldsymbol {\beta }}=\mathbf {0}$ og $H_{A}:\mathbf {C} {\boldsymbol {\beta }}\neq \mathbf {0}$, der $\mathbf {C}$ er en valgt matrise for den lineære hypotesen. F-observatoren er definert som

$F={\frac {{\hat {\boldsymbol {\beta }}}\mathbf {C} ^{T}\left(\mathbf {C} (\sum _{i}^{m}\mathbf {X} _{i}^{T}\Sigma _{i}^{-1}\mathbf {X} _{i})^{-1}\mathbf {C} ^{T}\right)^{-1}\mathbf {C} {\hat {\boldsymbol {\beta }}}^{T}}{rank(\mathbf {C} )}}$

F følger da en F-fordeling med frihetsgrader bestemt av rangen til matrisen $\mathbf {C}$ og klyngenivået.

Sannsynlighetskvotetest

Hvis $L_{0}$ refererer til likelihood-funksjonen til en referansemodell med $p_{0}$ parametre og $L_{1}$ refererer til likelihood-funksjonen til en redusert modell med $p_{1}<p_{0}$ parametre, er sannsynlighetskvotetesten, eller likelihood ratio-testen, definert som

$-2\log \left({\frac {L_{1}}{L_{0}}}\right)=2\left(\log(L_{0})-\log(L_{1})\right)\sim \chi _{p_{0}-p_{1}}^{2}$

Et viktig poeng er at sannsynlighetskvotetester for å teste faste effekter forutsetter at modellene som sammenlignes er estimert med maximum likelihood. Omvendt, så forutsetter sannsynlighetskvotetester for tilfeldige effekter at modellene som sammenlignes er estimert med restricted maximum likelihood.

Sannsynlighetskvotetester der en vil teste om en tilfeldig effekt kan utelukkes fra modellen krever en ytterligere justering, da den tilfeldige effekten i denne situasjonen er på grensen til parameterrommet. Dette vil si at, i motsetning til situasjonen for faste effekter der en tester $H_{0}:\beta =0$ mot $H_{A}:\beta \neq 0$, tester en $H_{0}:\sigma ^{2}=0$ mot $H_{A}:\sigma ^{2}>0$, da $\sigma ^{2}$ per definisjon er et positivt tall. Hvis en tester om en kan utelukke én tilfeldig effekt, følger testobservatoren da ${\frac {1}{2}}(\chi _{0}^{2}+\chi _{1}^{2})={\frac {1}{2}}\chi _{1}^{2}$. I praksis anbefales det å se på p-verdier fra slike tester som unøyaktige, slik at kun svært store eller små verdier brukes til å beholde eller forkaste nullhypotesen.