Hoofdstuk 36 Multilevel-analyse

In dit hoofdstuk wordt besproken:
  • multilevel-modellen met meerdere voorspellers
  • verschil tussen fixed en random effecten
  • vergelijken van modellen
  • multilevel-modellen met meerdere voorspellers
  • cross-level interactie
Deze stof wordt behandeld in de volgende Open Universiteitscursus(sen):
  • Onderzoekspractium longitudinaal onderzoek (PB1712)
Dit hoofdstuk bouwt voort op deze andere hoofdstukken:
  • Multiple-regressie

36.1 Inleiding

Als we een continue variabele willen voorspellen of verklaren uit andere continue variabelen, dan is regressieanalyse de meest toegepaste analysetechniek. Bij regressieanalyse was de aanname dat ieder datapunt volkomen onafhankelijk is. Maar wat te doen wanneer dit niet het geval is? Als onderzoek gedaan wordt onder leerlingen, zijn leerlingen binnen dezelfde klas, dezelfde school en misschien wel hetzelfde district niet op zijn minst een beetje afhankelijk van elkaar? Dit zijn databestanden die beter met een multilevelanalyse (MLA) zijn te analyseren dan met een klassieke regressieanalyse, omdat MLA is ontworpen om rekening te houden met het hiërarchische karakter van de data. Multilevelanalyse wordt dan ook steeds vaker toegepast in psychologisch onderzoek. In dit hoofdstuk geven we een beknopte uitleg over wat MLA is, en hoe we de analyse resultaten ervan moeten interpreteren.

36.2 Wat zijn multileveldata?

Binnen de psychologie vormen mensen vrijwel altijd de onderzoekseenheden. Om het verband tussen variabelen te schatten worden deze variabelen daarom gemeten bij mensen, bijvoorbeeld via vragenlijsten. De onderzochte mensen vormen meestal een steekproef uit een grotere populatie, bijvoorbeeld alle mensen in Groningen met een betaalde baan of alle gepensioneerde mensen met een vaste partner. Daarbij gaan we ervan uit dat de mensen in de steekproef onafhankelijk van elkaar zijn, dat wil zeggen dat de antwoorden of observaties van de ene persoon niet afhangen van de antwoorden van een ander persoon. Deze onafhankelijkheid is een van de voorwaarden bij regressieanalyse. Als we deze voorwaarde schenden, dan levert een regressieanalyse resultaten op die niet zuiver zijn.

Er zijn echter veel voorbeelden te bedenken waarbij de antwoorden of observaties niet onafhankelijk zijn. Stel, we willen onderzoeken of bij werknemers in verschillende bedrijven leiderschapsstijl invloed heeft op werkprestaties. Als er meerdere werknemers uit een bedrijf in de steekproef zitten, dan zullen hun antwoorden van elkaar afhangen omdat ze met dezelfde kenmerken van het bedrijf te maken hebben. De onderzoekseenheden (werknemers) zijn hier dus gegroepeerd binnen bedrijven. De gegevens hebben dus verschillende niveaus (levels): werknemers vormen in dit voorbeeld het eerste niveau en bedrijven het tweede.

Het meest klassieke voorbeeld van multileveldata is afkomstig uit onderzoek naar kinderen op scholen. Prestaties van kinderen (niveau 1) hangen mogelijk af van de klas (niveau 2) waarin ze zitten en de school (niveau 3) waarop ze zitten. Sommige klassen scoren mogelijk hoger dan andere vanwege een goede docent of een positief klimaat. Ook is bekend dat sommige scholen beter scoren dan andere.

Een ander voorbeeld is crosscultureel onderzoek naar de opvattingen van mensen. De mensen zijn hier gegroepeerd binnen landen. Ook in dit voorbeeld kunnen antwoorden van mensen op bepaalde vragen binnen een land meer op elkaar lijken dan de antwoorden van mensen uit verschillende landen. De cultuur van een land zorgt ervoor dat de antwoorden van landgenoten met elkaar gaan correleren.

In deze voorbeelden zijn de multileveldata hiërarchisch geordend. Een andere term voor multileveldata is dan ook hiërarchische data. Bij hiërarchische data zijn de basisonderzoekseenheden gegroepeerd binnen een variabele op een hoger niveau, die ook weer gegroepeerd kan zijn binnen een nog hoger niveau. Het voorbeeld van drie niveaus werd al eerder gegeven, namelijk kinderen die gegroepeerd zijn binnen klassen die weer gegroepeerd zijn binnen scholen. Een ander voorbeeld van drie niveaus is dat van werknemers die gegroepeerd zijn binnen afdelingen die weer gegroepeerd zijn binnen bedrijven. In beide voorbeelden zijn de individuen (de kinderen of de werknemers) het eerste niveau. Ieder individu heeft een eigen rij in het databestand. Voor alle hogere niveaus (zoals de klas of de school) moeten we een variabele aanmaken die aangeeft in welke categorie het individu zit.

Een voorbeeld staat in tabel 1. De eerste twee kinderen komen uit de eerste klas van de eerste school, kind 3 en 4 komen uit dezelfde school maar uit de tweede klas, enzovoort. Bij elke kind meten we een bepaalde variabele zoals in dit voorbeeld taalvaardigheid en rekenvaardigheid.

Tabel 36.1: Voorbeelddata.
id klas school lezen rekenen
1 1 Beatrix 6 7
2 1 Beatrix 8 9
3 2 Beatrix 7 6
4 2 Beatrix 8 7
5 3 De Rank 6 4
6 3 De Rank 4 5
7 3 De Rank 7 9
8 4 De Rank 8 6

Als we zijn geïnteresseerd in een algemene vraag waarbij klassen of afdelingen uit het onderzoek een steekproef vormen (en we dus niet zijn geïnteresseerd in het specifieke resultaat per klas) dan is het uitvoeren van gewone regressieanalyse op het hele bestand niet correct omdat dan de eventuele afhankelijkheid van de residuen of error-termen binnen een groep (klas, afdeling) wordt genegeerd. Een belangrijke aanname bij veel statistische procedures zoals regressieanalyse is namelijk dat de gegevens onafhankelijk zijn. Als de steekproef afkomstig is uit bepaalde groepen van deelnemers (klassen, scholen, afdelingen en bedrijven), dan hebben we een andere analysetechniek nodig, namelijk multilevelanalyse. Bij multilevelanalyse zijn we vooral geïnteresseerd in een verband tussen variabelen over alle groepen heen, waarbij we rekening kunnen houden met de afhankelijkheden tussen de personen binnen de groepen. Een ander voordeel van multilevel analyse is dat deze techniek heel geschikt is om met missing data om te gaan. Met name in experimentele designs met herhaalde metingen maakt een multilevel aanpak veel beter gebruik van de beschikbare data dan de variantieanalyse voor herhaalde metingen.

36.3 Een illustratie van multileveldata met gewone regressieanalyse

Aan de hand van een verzonnen voorbeeld illustreren we het probleem van multileveldata. We gaan een eenvoudig regressiemodel testen op een databestand van \(N = 80\) respondenten, namelijk werknemers bij vier verschillende bedrijven. Elk van de vier bedrijven heeft \(20\) respondenten. Het bestand bevat een afhankelijke variabele (Stress) een gestandaardiseerde onafhankelijke variabele (Werkdruk) en een groepsindeling (Bedrijf). Het bestand is zodanig geconstrueerd dat voor de vier bedrijven geldt:

\[\begin{equation} stress = b_0 + b_1werkdruk + \epsilon \end{equation}\]

Echter, voor ieder bedrijf geldt een andere waarde voor het intercept \(b_0\) en een andere voor de richtingscoefficient \(b_1\).

Binnen elk bedrijf geldt dat er een positief verband is tussen werkdruk en stress; dus globaal is het zo dat hoe meer werkdruk de werknemers ervaren hoe meer stress ze zullen hebben. Daarnaast spelen natuurlijk ook individuele factoren een rol in de beleving van stress, die we met de error term \(\epsilon\) aanduiden.

Wanneer we een regressieanalyse uitvoeren op deze data, dan geeft dat de volgende resultaten: \(b_0 =\) 4.19 en \(b_1 =\) 0.38. De verklaarde variantie van stress door werkdruk is 10.9%, want de \(R^2 =\) 0.109. Als we vier aparte regressies doen op elk van de vier bedrijven, dan zijn de resultaten als volgt:

Tabel 36.2: Regressiecoefficienten per bedrijf in voorbeelddata.
Bedrijf b0 b1 p R2
1 3.07 0.61 0.00 0.53
2 4.93 0.39 0.00 0.42
3 4.07 0.71 0.00 0.38
4 4.81 0.06 0.55 0.02

Het berekende intercept van de totale groep ligt ongeveer in het midden van de intercepten van de vier groepen. De regressiecoëfficiënt ligt ook bij het gemiddelde van de coëfficiënten van elke afzonderlijk bedrijf. Met name de bedrijven met de hoge b-waarden laten een vrij hoge verklaarde variantie zien. In Figuur 36.1 staat deze analyse grafisch weergegeven. De vier afzonderlijke regressielijnen zoals gepresenteerd in bovenstaande tabel en de figuur passen beter bij de data dan de gemiddelde regressielijn zoals hierboven beschreven.

Voorbeeld data met vier bedrijven

Figuur 36.1: Voorbeeld data met vier bedrijven

De lijnen behorend bij bedrijf 1 en 3 liggen lager dan de andere lijnen (bedrijf 2 en 4), vanwege de lagere waarde van het intercept. De punten behorend bij het eerste bedrijf liggen allemaal dicht bij hun lijn (dus relatief weinig error); er is hier een hoge mate van verklaarde variantie. Bij de lijn van bedrijf 4 is de variatie van de punten rondom de lijn het grootst, deze lijn past daarom niet zo goed bij de punten.

Het doel bij multilevelanalyse is om een intercept en een regressiecoëfficiënt (of meerdere coëfficiënten) te schatten voor de hele groep, waarbij men rekening houdt met de fluctuaties over de bedrijven. Het gaat daarbij om fluctuaties van de hoogte van de lijnen (intercepts) en van de richting van de lijnen (regressiecoëfficiënten).

36.4 Fixed en random effecten

Om de werking van een multilevelanalyse helder te maken, kijken we eerst nog een keer terug naar het meest bekende (basis)regressiemodel. Dit ziet er als volgt uit:

\[\begin{equation} Y = b_{0} + b_{1}X + \epsilon \end{equation}\]

Hierbij is Y de afhankelijke variabele, \(b_0\) het intercept (oftewel, het snijpunt met de y-as, en dus de waarde van Y als X nul is), \(b_1\) is de richtingscoëfficiënt die de relatie weergeeft tussen de onafhankelijke variabele X en de afhankelijke variabele Y, en \(\epsilon\) is de individuele variantie (de afwijking van de regressielijn) die niet door het model verklaard wordt.

De \(b_0\) en de \(b_1\) in dit model gelden voor de gehele onderzoekspopulatie en noemen we daarom ook wel fixed. De \(\epsilon\) is een afwijking per individu, en varieert dus binnen de onderzoekspopulatie en noemen we daarom ook wel random. Het basis regressiemodel bestaat dus uit een fixed intercept en een fixed richtingscoëfficiënt die gelden voor de gehele onderzoekspopulatie, en een random error die varieert voor de individuen binnen de onderzoekspopulatie.

Bij multilevelanalyse schat men, evenals bij regressieanalyse, regressiecoëfficiënten die het verband tussen twee variabelen weergeven. Net zoals de onderzoekseenheden van het eerste niveau (bijvoorbeeld mensen) een random steekproef vormen uit een grotere populatie, vormen de eenheden van het tweede niveau (bijvoorbeeld klassen, bedrijven) ook een random steekproef uit een grotere populatie van die eenheden. Bij zowel regressieanalyse als multilevelanalyse zijn we vooral geïnteresseerd in de gemiddelde verbanden of effecten, dat wil zeggen gemiddeld over alle kinderen of over alle klassen. Bij multilevelanalyse houdt men daarbij rekening met de mogelijkheid dat deze effecten kunnen variëren over de verschillende groepen (bijvoorbeeld klassen) binnen de onderzoekspopulatie. Dit noemen we ook wel random effecten. Het belangrijke verschil tussen multilevelanalyse en regressieanalyse is de aanwezigheid van deze random effecten. Bij een standaard regressieanalyse gaan we ervan uit dat de effecten niet variëren over groepen maar voor de hele steekproef vastliggen. Zoals eerder beschreven, worden dit ook wel fixed effecten genoemd. Het multilevelmodel schat naast de fixed effecten (die dus gelden voor de gehele onderzoekspopulatie) ook nog random effecten die variëren tussen de groepen binnen de onderzoekspopulatie. De random effecten kunnen we dus zien als een extra variantie binnen het model (naast de individuele variantie e van het basisregressiemodel), waarbij een afwijking van de algemene (fixed) regressielijn veroorzaakt wordt door de groep waar iemand zich in bevindt (de niveau 2-variabele). Als we bijvoorbeeld willen kijken hoe de wiskundescore (Y) afhangt van het IQ van de leerling (X), kan het zo zijn dat dit mede afhankelijk is van de klas waarin de leerling zich bevindt, wat kan leiden tot een afwijking van het algemene regressiemodel. Deze afwijking kan zich bevinden in:

  • het intercept (\(b_0\)): een bepaalde klas heeft bijvoorbeeld een hoger algemeen wiskunde-niveau dan de andere klas. In dat geval varieert dus het intercept van het model per klas. Deze variatie tussen de intercepten van de verschillende klassen noemen we ook wel het random intercept. Als alleen het intercept van het model varieert over de klassen, bestaat het regressiemodel dus uit een aantal lijnen (1 per groep) die parallel lopen aan elkaar (zie figuur 36.2, links).

  • de richtingscoëfficiënt (\(b_1\)): de relatie tussen de wiskundescore en IQ kan bijvoorbeeld in een klas met een betere leraar minder sterk zijn doordat een goede leraar een leerling met een lager IQ ook wiskunde kan aanleren. Deze variantie in de richtingscoëfficiënt noemen we de random slope. Als ook de richtingscoëfficiënt beïnvloed wordt door de groep waar iemand zich in bevindt, betekent dit dus dat de regressielijnen per groep niet langer parallel lopen aan elkaar (zie figuur 36.2, midden). Merk op dat wanneer de klas de relatie tussen X en Y versterkt of verzwakt (zoals het geval is bij een random slope) dat dit te vergelijken is met het moderatiemodel dat we eerder tegenkwamen in de cursus. Het standaard moderatiemodel is echter op dit moment niet geschikt, omdat deze uitgaat van onafhankelijke metingen, terwijl het multilevelmodel uit gaat van metingen die gecorreleerd zijn binnen de klas.

  • zowel het intercept als in de richtingscoëfficiënt (zie figuur 36.2, rechts).

Illustratie random richtingscoëfficiënten en random intercepten

Figuur 36.2: Illustratie random richtingscoëfficiënten en random intercepten

Voor zowel het intercept als voor de richtingscoëfficiënt geldt dat die in de beschrijving van het multilevelmodel een extra index \(j\) krijgt als deze afhankelijk is van de groep waarin iemand zich bevindt. Een multilevelmodel met zowel een random intercept als een random slope ziet er dus als volgt uit:

\[\begin{equation} Y = b_{0j} + b_{1j}X + \epsilon \end{equation}\]

De \(b_{0j}\) kunnen we hierbij opsplitsen in de \(b_0\) die we eerder terugzagen in het basisregressiemodel (het fixed intercept), en de variantie eromheen (het random intercept) die veroorzaakt wordt door de variabele op het tweede niveau (aangeduid met index \(j\)). Hetzelfde geldt voor de richtingscoëfficiënt die in dit geval ook bestaat uit een fixed deel (geldende voor de hele onderzoekspopulatie) en uit een random deel (de variatie in de richtingscoëfficiënt afhankelijk van de groep waarin iemand zich bevindt). Oftewel: voor de onderzoekspopulatie als geheel is er een algemeen verband tussen wiskundescore en IQ (fixed richtingscoëfficiënt) waarbij de minimumwaarde van de wiskundescore afhankelijk is van het algemene wiskundeniveau van de populatie (fixed intercept). Per klas kan worden afgeweken van dit gemiddelde in zowel het algemene wiskundeniveau (random intercept) als in de sterkte van de relatie tussen IQ en wiskundescore (random slope). De variantie in het model die afhankelijk is van de groep (zowel in het intercept als in de richtingscoëfficiënt) wordt in de literatuur vaak uitgedrukt met de term G.

Een vaak berekende maat bij het uitvoeren van multilevelanalyses is de intraclasscorrelatie (ICC). Met deze maat berekenen we hoe groot de bijdrage van de groepen (variantie van het tweede niveau, uitgedrukt in G) is aan de totale variantie binnen het volledige model (oftewel, de \(\sigma^2(\epsilon)\) de variantie van de error component (\(\epsilon\)), die we zagen in het basis model, plus de G de variantie afkomstig van het tweede niveau). Uitgedrukt in een formule:

\[\begin{equation} ICC = \frac{G}{G + \sigma^2(\epsilon)} \end{equation}\]

Een hoge ICC betekent dat een groot deel van de variantie toe te schrijven is aan de groepsindeling. Bij een heel lage ICC (bijvoorbeeld 2%) is maar een klein deel van de variantie toe te schrijven aan de groepsindeling, en zijn de gegevens blijkbaar toch onafhankelijk van elkaar. In dat geval kan zonder probleem gewone regressieanalyse worden toegepast. De ICC kan het meest zuiver worden uitgerekend bij het zogenaamde intercept-only (IO) model. Dit is een model waarin alleen een intercept zit als fixed en random effect, dus waarin verder geen predictoren zijn opgenomen.

36.5 Eenvoudige multilevel modellen

36.5.1 Voorbeeld van multileveldata

In het vervolg van deze tekst wordt uitgelegd hoe we zo’n multilevelmodel (en variaties hierop) kunnen analyseren en hoe we de output hiervan kunnen interpreteren. Hierbij gebruiken we een voorbeeldbestand (\(N = 196\)) dat bestaat uit 10 klassen met per klas circa 20 leerlingen. De leerlingen zijn hier genest binnen de klassen. De data bestaan uit de volgende variabelen: klasnummer, wiskunde, intelligentie, ervaringLeerkracht, aantalLeerlingen, sekseLeerling, sekseLeerkracht. De wiskundescore, de intelligentie (IQ), en het geslacht van de leerling zijn per leerling gemeten (niveau 1). Daarnaast bevatten de data gegevens over het aantal jaren ervaring van de leerkracht, de klassengrootte en het geslacht van de leerkracht, die gemeten zijn per klas (niveau 2). Via het klasnummer kunnen we zien in welke klas een leerling zit. De afhankelijke variabele in dit bestand is de wiskundescore (\(M = 5.8, SD = 1.7\)).

Het aantal klassen (10) in dit voorbeeld is eigenlijk te laag voor een multilevelanalyse. Dat heeft te maken met de power (een effect vinden in de steekproef dat aanwezig is de populatie). Hoewel bij een normaal regressiemodel de power afhankelijk is van het aantal personen in de steekproef, is de power voor het schatten van de random effecten vooral afhankelijk van het aantal groepen (niveau 2) in de steekproef, in dit voorbeeld de klassen. Hierbij geldt als vuistregel dat er minimaal 20 groepen moeten zijn voor het uitvoeren van multilevelanalyses (Kreft & de Leeuw, 1998). Maar de precieze berekening van de power hangt van meerdere factoren af. In dit voorbeeld is de power (met slechts 10 klassen) dus aan de lage kant.

36.5.2 Centreren van de variabelen

Voordat we het uitvoeren van de multilevelanalyses in de praktijk brengen, is het van belang om te weten hoe u een variabele moet centreren. Het centreren van de predictorvariabelen en eventuele covariaten speelt een belangrijke rol bij multilevelanalyse, vooral om de uitkomsten beter te kunnen interpreteren. Als er twee niveaus zijn, kunnen we op verschillende manieren centreren. In de eerste plaats kunnen we het algemene gemiddelde (van alle observaties) gebruiken om te centreren. De gecentreerde variabele maken we door dit gemiddelde van alle scores af te trekken. In de tweede plaats kunnen we het gemiddelde van iedere groep (bijvoorbeeld klas) gebruiken om te centreren. De gecentreerde variabele wordt dan gemaakt door het desbetreffende groepsgemiddelde van de scores van de desbetreffende groep af te trekken. Anders gezegd: we kunnen centreren voor de dataset als geheel (“grand mean centering”), of per groep (“group mean centering”). In het eerste geval geven de individuele scores de afwijking van het algemene gemiddelde weer (net als bij data met maar een enkel niveau) en in het tweede geval de afwijking van hun groepsgemiddelde. Bij centreren rondom de groepsgemiddelden, halen we de groepsgemiddelden uit de desbetreffende variabele. Na deze vorm van centreren verschillen de groepsgemiddelden niet meer van elkaar, want deze zijn allemaal 0 geworden. De groepsgemiddelden kunnen dan wel als extra variabele (op niveau 2) aan het model worden toegevoegd. De keuze van de manier waarop we centreren heeft gevolgen voor de interpretatie van de uitkomsten. Het hangt van de onderzoeksvraag af welke van de twee het meest geschikt is, maar in het algemeen heeft centreren rondom de groepsgemiddelden de voorkeur. Voor meer details hierover zie Enders en Tofighi (2007).

36.5.3 Model A: een predictor en een groepsvariabele

In het eerste model wordt verondersteld dat de voorspellende variabele IQ van invloed is op de afhankelijke variabele wiskundescore. Hierbij kunnen we veronderstellen dat het effect op de wiskundescores afhangt van de klas waarin iemand zit, waarbij de klas een groepsvariabele (klasnummer) is die het tweede niveau definieert. De klas kan hierbij zowel het intercept (\(b_0\)) als de richtingscoëfficiënt (\(b_1\)) beïnvloeden, wat dus leidt tot een random intercept of een random slope (of beide). Of er een multilevelmodel met een random intercept of met een random slope (of beiden) geanalyseerd moet worden, is afhankelijk van iemands eigen verwachtingen. Verwacht de onderzoeker dat de klas van invloed is op het intercept, of op de richtingscoëfficiënt?

Hieronder bespreken we eerst de analyse met alleen een random intercept, en vervolgens de analyse met zowel een random intercept als een random slope.

36.5.4 Model A1: predictor IQ en random intercept

Het model met een fixed predictor en alleen een random intercept kunnen we uitdrukken in de volgende formule:

\[\begin{equation} wiskunde = b_{0j} + b_{1}intelligentie + \epsilon_{ij} \end{equation}\]

Merk op dat alleen de \(b_{0j}\) een index-\(j\) bezit, en dus een variantie heeft op het tweede niveau en dat we voor de volledigheid ook de \(\epsilon\) hebben voorzien van twee indices, want de error is ook volledig random. De \(b_{0j}\) kunnen we op de volgende manier opsplitsen in twee componenten:

\[\begin{equation} b_{0j} = \gamma_{00} + u_{0j} \end{equation}\]

Hierbij staat \(\gamma_{00}\) voor het fixed deel van het intercept, een geschat gemiddeld intercept voor de hele populatie, en is \(u_{0j}\) de afwijkingen rond het algemene intercept die veroorzaakt worden door de klas waar iemand in zit. In het algemeen zijn we niet in de \(u_{0j}\) geinterseerd, maar hooguit in de variantie van deze waarden, geschreven als \(\sigma^2_{u_{0j}}\). Wanneer we deze uitgesplitste \(b_{0j}\) weer invullen in de voorgaande formule, leidt dit tot de volgende formule :

\[\begin{equation} wiskunde = \gamma_{00} + u_{0j} + b{_1}intelligentie + \epsilon_{ij} \end{equation}\]

Het uitsplitsen van de formule in deze vier componenten (parameters) is van belang om later de output goed te kunnen interpreteren. In totaal worden er in dit model twee fixed parameters (\(\gamma_{00}\) en \(b_1\)) en twee random parameters (\(\sigma^2_{u_{0j}}\) en \(\sigma^2_{\epsilon_{ij}}\)) geschat, dit zijn varianties van de random termen.

Om dit model te analyseren, kunnen diverse softwarepakketten worden gebruikt. In elk programma dienen in ieder geval de volgende onderdelen te worden gespecificeerd. Ten eerste de afhankelijke variabele, ten tweede de predictor variabelen die de fixed effecten weergegeven, en ten derde de predictor variabele die de random effecten weergegeven. Een schematische korte weergave van het te analyseren model kan als volgt worden genoteerd:

ModelA1: Afhankelijk = wiskunde
      / Fixed = intercept intelligentie
      / Random = intercept (klasnummer)

We geven het model een naam, in dit geval “modelA1” en daarna worden de drie onderdelen gespecificeerd. De afhankelijke variabele is de variable wiskunde. In het fixed deel van het model wordt de predictor intelligentie opgegeven. Het intercept staat hier ook bij om dat deze standaard in het model zit. In veel programma’s hoef je deze term daarom niet expliciet op te geven. Het intercept wordt hier als ook random term opgenomen, met tussen haakjes de variabele waarover deze term random varieert. Iedere klas krijgt hier dus zijn eigen intercept. Na analyse van dit model zijn de resultaten als volgt. Allereerst zien we de fixed effecten.

Tabel 36.3: Schattingen fixed effecten model A1
Value Std.Error DF t-value p-value lower upper
intercept 5.84 0.19 185 31.52 0 5.48 6.21
intelligentie 1.13 0.07 185 15.47 0 0.99 1.28

Het fixed intercept is 5.84, hetgeen een schatting is van de gemiddelde wiskunde score over alle klassen heen. Deze schatting ligt heel dicht bij het eerder getoonde gemiddelde van wiskunde in de data. Intelligentie heeft een regressieparameter van 1.13. Eén standaardafwijking stijging van de intelligentie voorspelt ongeveer een punt hoger voor wiskunde. Hierna worden de random effecten getoond.

Tabel 36.4: Schattingen random effecten model A1
random effect
(Intercept) 0.291
Residual 1.035

Het tweede deel van de output laat zien dat de geschatte variantie van het random deel van het intercept (\(\sigma^2_{u_{0j}}\)) 0.291 is. De residuele variantie ( \(\sigma^2_{\epsilon_{ij}}\)) is 1.035.

Met bovenstaande gegevens kunnen we de intraclasscorrelatie uitrekenen, die de verhouding weergeeft tussen die variantie tussen de groepen ten opzichte van de totale variantie:

\[\begin{equation} ICC = \frac{0.291}{0.291 + 1.035} = 0.219 \end{equation}\]

Dit betekent dat circa 21.9 procent van de totale variantie van de wiskundescores is veroorzaakt door de variantie tussen de klassen. Overigens is dit een ICC, die gecorrigeerd is voor IQ in het model. Voor een ongecorrigeerde ICC zou het intercept-only (IO) model moeten worden gebruikt, waarin geen predictoren zijn opgenomen.

Ook staan bij de analyseresultaten een aantal schattingen van hoe goed dit model bij de data past. Een van de maten die men hiervoor gebruikt is de zogenaamde -2 Restricted Log Likelihood (-2LL), deze maat noemen we \(deviance\). De deviance van model A1 is \(584.64\). Hoe lager deze waarde is, hoe beter het model past bij de data. Er is geen criterium dat aangeeft wat een goede deviance is, omdat de waarde afhangt van het aantal variabelen in het model en het aantal waarnemingen. De deviancewaarden zeggen niet zo veel, maar worden vooral gebruikt om twee modellen met elkaar te vergelijken. Het model met de laagste deviancewaarde past beter bij de data dan een model met een hogere waarde. Op basis van deze maat kunnen we dus niet zeggen of een bepaald model goed is, maar wel of een model beter past dan een ander model. We kunnen toetsen of een verschil in deviance significant is met een \(\chi^2\)-toets als twee modellen genest zijn. Uitleg over wat geneste modellen zijn staat in de verdieping: toetsen van geneste modellen.

Met behulp van een \(\chi^2\)-toets kunnen we nagaan of het verschil tussen twee deviancewaarden significant is, en het ene model dus significant beter bij de data past dan het andere. U kunt hiervoor gebruik maken van een tabel met \(\chi^2\)-waarden, die op diverse plekken op het internet te vinden is en ook voor de vrijheidsgraden achterin dit hoofdstuk. Het aantal vrijheidsgraden (\(df\) = degrees of freedom) is in dit geval het verschil in geschatte parameters van de twee modellen. Als het verschil tussen de twee deviancewaarden van de modellen groter is dan de waarde die achter de df in de tabel staat, dan is het verschil tussen de twee modellen significant. Het model met de laagste deviance past dan beter bij de data dan het andere model.

36.5.5 Model A2: predictor IQ en random intercept en random slope

Het model met een random intercept en een random slope kunnen we uitdrukken in de volgende formule:

\[\begin{equation} wiskunde = b_{0j} + b_{1j}intelligentie + \epsilon_{ij} \end{equation}\]

Merk op dat alleen de \(b_{0j}\) een index-\(j\) bezit, en dus een variantie heeft op het tweede niveau en dat we voor de volledigheid ook de \(\epsilon\) hebben voorzien van twee indices, want de error is ook volledig random. De \(b_{0j}\) kunnen we op de volgende manier opsplitsen in twee componenten:

\[\begin{align*} b_{0j} = \gamma_{00} + u_{0j} \\ b_{1j} = \gamma_{10} + u_{1j} \end{align*}\]

De \(\gamma\)’s hebben dus in beide gevallen betrekking op het fixed deel, en de \(u\)’s op het random deel (oftewel, de variantie veroorzaakt door de klassen). Als we beide uitgesplitste \(b\)’s invullen in de eerdere formule, leidt dit tot de volgende formule:

\[\begin{equation} wiskunde = (\gamma_{00} + u_{0j}) + (\gamma_{10} + u_{1j})intelligentie + \epsilon_{ij} \end{equation}\]

Dit model bevat dus in totaal vijf parameters, die we later in de analyseresultaten weer zullen tegenkomen.

Een schematische weergave van dit model kan als volgt worden genoteerd:

ModelA2: Afhankelijk = wiskunde
      / Fixed = intercept intelligentie
      / Random = intercept intelligentie (klasnummer)

Het random deel van het model heeft dus een extra term (intelligentie) vergeleken met model A2.

Wanneer we dit model gebruiken om de data te analyseren, dan zijn de resultaten weer gesplitst in een fixed en een random deel. Er zijn twee fixed effecten: ten eerste het algemene intercept (\(\gamma_{00}\)) en ten tweede de regressiecoëfficiënt behorend bij het effect van IQ (\(\gamma_{10}\)). Er zijn in totaal drie random effecten: de bekende errorvariantie op niveau 1: \(\sigma^2_{\epsilon_{ij}}\) en twee errorvarianties op niveau 2: \(\sigma^2_{u_{0j}}\) behorende bij het intercept en \(\sigma^2_{u_{1j}}\) behorende bij de slope van IQ.

Tabel 36.5: Schattingen fixed effecten model A2
Value Std.Error DF t-value p-value lower upper
intercept 5.87 0.18 185 32.81 0 5.52 6.22
intelligentie 1.15 0.20 185 5.82 0 0.76 1.54

Het fixed intercept is 5.87, hetgeen ook in dit model een schatting is van de gemiddelde wiskunde score over alle klassen heen. Intelligentie heeft een regressieparameter van 1.15. Hierna worden de random effecten getoond.

Tabel 36.6: Schattingen random effecten model A2
random effect
(Intercept) 0.282
intelligentie 0.351
Residual 0.693

Het tweede deel van de output laat zien dat de geschatte variantie van het random deel van het intercept \(\sigma^2_{u_{0j}} =\) 0.282, en de geschatte variantie van het random deel van de slope \(\sigma^2_{u_{1j}} =\) 0.351. Voor de residuele variantie geldt \(\sigma^2_{\epsilon_{ij}} =\) 0.693. De residuele variantie van dit model (A2) is nu veel kleiner dan bij het voorgaande model (A1), wat aangeeft dat dit model beter past bij de data.

Zoals onderstaande tabel laat zien, waarin beide modellen worden vergeleken met betrekking tot hun fit-maten, past model A2 inderdaad beter dan model A1.

Tabel 36.7: Model vergelijking
df AIC BIC logLik Test L.Ratio p-value
Model A1 4 592.64 605.72 -292.32 NA NA
Model A2 5 540.70 557.04 -265.35 1 vs 2 53.94 0

De deviance (-2*LogLikelihood), 530.7 (met vijf geschatte parameters) is bij model A2, dus lager is dan de deviance bij model A1 van 584.64 (met vier geschatte parameters). De zogenoemde “likelihood-ratio (LR) test” geeft aan dat de verbetering ten opzichte van model A1 significant is, want de \(\chi^2(1)=\) 53.94 (dat is 584.64 - 530.7) is statistisch significant. Het toevoegen van de slope als random effect verbetert het model en het lijkt er dus op dat de effecten van IQ nogal variëren over de klassen, zowel met betrekking tot het intercept als tot de slope.

36.6 Covariaat op het eerste niveau

In de meeste modellen worden naast de voorspellende en de afhankelijke variabelen ook andere variabelen meegenomen die van invloed kunnen zijn op de afhankelijke variabelen. Deze variabelen noemen we covariaten. Denk aan de leeftijd en het geslacht van de leerling, de ervaring van de leerkracht of het wiskundeboek dat op school gebruik wordt. Door het opnemen van een covariaat in het model verwacht men de relatie tussen de predictor en de afhankelijke variabele zuiverder te kunnen schatten.

Meestal voegt men een covariaat als een fixed effect toe aan het multilevelmodel. Maar als een onderzoeker wil onderzoeken of het effect van de covariaat op de afhankelijke variabele afhankelijk is van de groepen in de onderzoekspopulatie, dan is het nodig om de covariaat ook als random effect mee te nemen.

Voordat we covariaten als random effect meenemen, moeten we bedenken dat de power voor het schatten van de random effecten afhankelijk is van het aantal groepen in de steekproef, zoals eerder besproken. Het is daarom verstandig om het aantal random effecten zo klein mogelijk te houden. Als we bepaalde variabelen dus alleen meenemen om hier in het algemeen voor te corrigeren, dan voegen we deze als fixed effect toe aan het model. Alleen als er echt gegronde reden zijn om te verwachten dat de invloed van deze variabele varieert over de verschillende groepen, neemt men deze als random effect mee. In onderstaand voorbeeld zullen we u eerst laten zien hoe u het geslacht van de student (sekseLeerling) als fixed covariaat (covariaat op het eerste niveau) mee kunt nemen in het model. Dit is een covariaat van het eerste niveau omdat hij varieert per leerling, en leerlingen vormen het eerste niveau. Deze variabele wordt toegevoegd aan het model met het random intercept en random slope over de klassen (model A2). Als we de variabele sekseLeerling toevoegen aan het model A2 dan ziet het nieuwe model (B) er als volgt uit:

\[\begin{equation} wiskunde = b_{0j} + b_{1j}intelligentie + b_2sekseLeerling + \epsilon_{ij} \end{equation}\]

Omdat het intercept en de slope van intelligentie random coëfficiënten zijn, kunnen we net als bij model A2 ook schrijven:

\[\begin{equation} wiskunde = \gamma_{00} + u_{0j} + (\gamma_{10} + u_{1j})intelligentie + b_2sekseLeerling + \epsilon_{ij} \end{equation}\]

Er is dus een extra fixed effect \(b_2\) bijgekomen in vergelijking met model A2. Deze covariaat wordt toegevoegd aan het model in het volgende modelschema:

ModelB: Afhankelijk = wiskunde
      / Fixed = intercept intelligentie sekseLeerling
      / Random = intercept intelligentie (klasnummer)

Wanneer we een analyse doen op dit model dan komt de interpretatie van de resultaten overeen met het voorgaande model. In de resultaten staan de schattingen van zes parameters. Drie daarvan, zijn fixed effecten, namelijk een algemeen gemiddelde (het fixed-deel van het intercept:\(\gamma_{00}\)) en twee regressiecoëfficiënten (het fixed deel van \(b_1\) (\(\gamma_{10}\)) ,en \(b_2\)) behorend bij de variabelen intelligentie en sekseLeerling. Drie zijn random parameters (\(u_{0j}\), \(u_{1j}\), en \(\epsilon_{ij}\)).

De resultaten van de fixed effecten staan in de volgende tabel.

Tabel 36.8: Schattingen fixed effecten model C
Value Std.Error DF t-value p-value lower upper
intercept 5.72 0.18 184 31.25 0.00 5.36 6.09
intelligentie 1.14 0.19 184 5.86 0.00 0.76 1.52
sekse_jongen 0.29 0.12 184 2.35 0.02 0.05 0.53

De deviance van model B is \(527.6\).

Als we kijken naar de variabele sekseLeerling, dan zien we dat deze een significant effect heeft op de wiskundescore (0.29, \(p =\) 0.02). Kijken we naar de deviance, dan blijkt dat deze lager is dan die van model A2. De LR-toets geeft echter aan dat de verbetering ten opzichte van model A2 niet significant is, want de \(\chi^2(1)=\) 3.1 (dat is 530.7 - 527.6) is statistisch niet significant.

De schattingen van de varianties van de random effecten staan hieronder.
Tabel 36.9: Schattingen random effecten model B
random effect
(Intercept) 0.259
intelligentie 0.339
Residual 0.680

36.7 Covariaat op het tweede niveau

Stel dat we naast het intelligentie van de leerling ook beschikken over informatie op het tweede niveau, dus informatie over de klassen. We weten bijvoorbeeld het aantal jaren ervaring in het lesgeven van de leerkracht. Deze informatie staat in de variabele ervaringLeerkracht, die binnen een klas altijd dezelfde waarde heeft, maar tussen klassen kan variëren. De coëfficiënt van ervaringLeerkracht varieert niet op een hoger niveau dan klas (want klas is in deze data het hoogste niveau) en kan dus in principe geen random coëfficiënt zijn in dit model.

De syntax is volledig vergelijkbaar met de syntax van Model B, waarbij sekseLeerling vervangen is door ervaringLeerkracht. In deze analyse gebruiken we de gecentreerde versie van ervaringLeerkracht (ervaringLeerkracht_c) in plaats van ervaringLeerkracht zelf.

ModelC: Afhankelijk = wiskunde 
      / Fixed = intercept intelligentie ervaringLeerkracht_c
      / Random = intercept intelligentie (klasnummer)
Tabel 36.10: Schattingen fixed effecten model C
Value Std.Error DF t-value p-value lower upper
intercept 5.87 0.16 185 36.98 0.0 5.56 6.18
intelligentie 1.14 0.20 185 5.84 0.0 0.76 1.53
ervaring 0.07 0.04 8 1.85 0.1 -0.02 0.15

De deviance van model C is \(532.32\).

De resultaten tonen eveneens de schattingen van zes parameters. Drie daarvan zijn fixed effecten, namelijk een algemeen gemiddelde (het fixed-deel van het intercept: \(\gamma_{00}\)) en twee regressiecoëfficiënten (het fixed deel van \(b_1\) (\(\gamma_{10}\)) ,en \(b_2\)) behorend bij de variabelen intelligentie en ervaringLeerkracht. Er zijn wederom drie random parameters (\(u_{0j}\), \(u_{1j}\), en \(\epsilon_{ij}\)). De resultaten staan samengevat in de volgende tabel. In de tabel is te zien dat ervaringLeerkracht geen significant effect heeft op de wiskundescore en de deviance (532.32) is hoger dan de deviance van model A2, wat betekent dat het toevoegen van ervaringLeerkracht het model niet beter maakt.

De schattingen van de varianties van de random effects staan hieronder.
Tabel 36.11: Schattingen random effecten model C
random effect
(Intercept) 0.213
intelligentie 0.344
Residual 0.694

De residuele variantie is niet omlaag gegaan ondanks het toevoegen van een extra predictor.

36.8 Interactie-effect tussen eerste en tweede niveau

Model A veronderstelde dat de wiskundescores afhingen van de intelligentie van de leerling (niveau 1). Model C voegde daar de ervaring van de leerkracht (ervaringLeerkracht_c, niveau 2) als covariaat aan toe. Hieruit bleek dat de ervaring van de leraar weinig of geen invloed had op de wiskundescore. Het is echter ook mogelijk dat niet de wiskundescore zelf, maar wel de relatie tussen intelligentie en wiskundescore beïnvloed wordt door de ervaring van de leraar. Een ervaren leerkracht slaagt er mogelijk beter in ook leerlingen met een lagere intelligentie goede scores te laten halen, terwijl bij een onervaren leerkracht de intelligentie veel bepalender is voor de wiskundescores. Anders geformuleerd: het effect van inteligentie op de wiskundescores wordt mogelijk gemodereerd door de ervaring van de leerkracht. In het algemeen geldt dat wanneer we een variabele op niveau 1 hebben waarbij we een random effect veronderstellen en we veronderstellen tevens dat dit random effect afhangt van een andere variabele dan is er sprake van moderatie. Dit wordt onderzocht met model D, waarbij dus sprake is van een zogenaamde “cross-level” interactie. Dat is dus een interactie tussen variabelen die op een verschillend niveau zijn gemeten.

Als er een interactie (moderatie) is tussen twee variabelen, dan moeten we het product van deze variabelen opnemen in het model. In dit geval moeten we dus het product ervaringLeerkracht_c x intelligentie opnemen. Zoals we in bovenstaand model ervaringLeerkracht_c toevoegden als fixed effect, voegen we bij model D de interactie tussen ervaringLeerkracht_c en intelligentie toe als fixed effect. Dit kunnen we uitdrukken in de volgende formule:

\[\begin{align*} wiskunde = b_{0j} + b_{1j}intelligentie + b_2sekseLeerling +\\ b_3ervaringLeerkracht\_c*intelligentie + \epsilon_{ij} \end{align*}\]

Het symbool “\(*\)” geeft het product van de twee variabelen aan, met andere woorden dat is de interactieterm. Dit model heeft in totaal zeven termen (bedenk dat zowel \(b_{0j}\) als \(b_{1j}\) een fixed en een random term vertegenwoordigen) waarvan we er zes al eerder zijn tegengekomen. Alleen de interactieterm (ervaringLeerkracht_c*intelligentie) is nieuw. Er zijn vier fixed effecten: het algemene intercept, de regressiecoëfficiënt van ervaring, de regressiecoëfficiënt van intelligentie en de regressiecoëfficiënt van de interactie tussen ervaringLeerkracht en intelligentie. Verder zijn er dezelfde drie random effecten als in model C.

In schematische notatie ziet het model er als volgt uit:

ModelD: Afhankelijk = wiskunde 
      / Fixed = intercept intelligentie ervaringLeerkracht_c
                intelligentie*ervaringLeerkracht_c
      / Random = intercept intelligentie (klasnummer)

De fixed effecten laten naast een significant intercept een significant effect van intelligentie zien, net als bij model C. Het effect van de interactie tussen intelligentie en ervaring is niet statistisch significant. De ervaring van de leerkracht lijkt geen modererend effect te hebben op het verband tussen het intelligentie en de wiskundescores. De resultaten van het random deel zijn vergelijkbaar met de voorgaande analyse.

Tabel 36.12: Schattingen fixed effecten model D
Value Std.Error DF t-value p-value lower upper
intercept 5.86 0.16 184 36.74 0.00 5.55 6.18
intelligentie 1.13 0.18 184 6.40 0.00 0.78 1.48
ervaring 0.07 0.04 8 1.75 0.12 -0.02 0.15
interactie 0.07 0.04 184 1.73 0.08 -0.01 0.15

De deviance van dit model is 534.05 (aantal parameters = 7). Vergeleken met model C levert dit model dus geen verbetering, hetgeen ook bevestigt dat het opnemen van de interactie geen verbetering oplevert.

De schattingen van de varianties van de random effects staan hieronder.
Tabel 36.13: Schattingen random effecten model D
random effect
(Intercept) 0.216
intelligentie 0.273
Residual 0.694

36.9 Verdieping: model D afgeleid

Een andere manier om model D af te leiden is door ervan uit te gaan dat de random parameters (\(b_{0j}\) en \(b_{1j}\)) worden voorspeld door de variabele ervaringLeerkracht_c. Voor een wat handzamere notatie vervangen we de variabelenaam intelligentie door IQ en ervaringLeerkracht_c door EL. In formules wordt dit, allereerst op niveau 1:

\[\begin{equation} wiskunde = b_{0j} + b_{1j}IQ + \epsilon_{ij} \end{equation}\]

En op niveau 2 de voorspelling van de random parameters (\(b_{0j}\) en \(b_{1j}\)):

\[\begin{align*} b_{0j} = \gamma_{00} + \gamma_{01}EL + u_{0j} \\ b_{1j} = \gamma_{10} + \gamma_{11}EL + u_{1j} \end{align*}\]

Als we de uitdrukkingen van niveau 2 invullen in de eerste formule, dan krijgen we het volgende resultaat:

\[\begin{equation} wiskunde = (\gamma_{00} + \gamma_{01}EL + u_{0j}) + (\gamma_{10} + \gamma_{11}EL+ u_{1j})IQ + \epsilon_{ij} \end{equation}\]

En dit kan verder worden uitgewerkt tot:

\[\begin{equation} wiskunde = \gamma_{00} + \gamma_{01}EL + u_{0j} + \gamma_{10}IQ + \gamma_{11}EL*IQ + u_{1j}IQ + \epsilon_{ij} \end{equation}\]

Op de naamgeving van de parameters na en de uitsplitsing van het random effect in een random en een fixed deel, is bovenstaande formule gelijk aan de in de tekst gegeven formule van model D.

36.10 Verdieping: Toetsen van geneste modellen

Modellen kunnen we vergelijken via de Likelihood-ratio (LR) test, mits de modellen genest zijn. Twee modellen noemen we genest als ze dezelfde termen bevatten en één van de modellen ten minste één extra term bevat.

Voorbeeld:

Model A: \(\hat{y} = b_1x_1 + b_2x_2\)

Model B: \(\hat{y} = b_1x_1 + b_2x_2 + b_3x_3\)

Model C: \(\hat{y} = b_1x_1 + b_3x_3\)

Model A is genest binnen model B. In model B wordt namelijk een extra parameter geschat (\(b_3\)) in vergelijking met model A. Model C is ook genest binnen model B omdat het verschil een parameter is (\(b_2\)). Model A en C zijn echter niet binnen elkaar genest.

De vraag is of het model met de extra parameter wel significant beter bij de data past. Dus: is het toevoegen van \(b_3\) in model B wel een goede keuze of levert model B nauwelijks een verbetering op ten opzichte van model A? Als twee modellen genest zijn, dan kan het verschil in de -2LL (deviance) waarden van beide modellen worden vergeleken met behulp van de \(\chi^2\)-toets. De deviances worden eerst van elkaar afgetrokken. Vervolgens trekt u ook het aantal geschatte parameters van de twee modellen van elkaar af. Het verschil in -2LL en het verschil in aantal parameters zijn dan de \(\chi^2\) en de \(df\) die gebruikt worden om de significantie van het verschil tussen de modellen vast te stellen. Is het verschil significant dan heeft het toevoegen van de extra parameter zin gehad.

Tabel 36.14: Kritieke waarden van de Chi-kwadraat verdeling
df p=0.05 p=0.01
1 3.84 6.63
2 5.99 9.21
3 7.81 11.34
4 9.49 13.28
5 11.07 15.09
6 12.59 16.81
7 14.07 18.48
8 15.51 20.09
9 16.92 21.67
10 18.31 23.21