Vierde beweging · Verschil onthechten

Oefening 11.1 · Mediatie: direct en indirect

Mediatie: direct, indirect en totaal · mediation

Data bij deze oefening — luchtkwaliteit.csv luchtmetingen in New York · diamantjes.csv de twaalf steentjes uit het boek · wat elke kolom betekent

prikkel

Eén helling van 0.13 — maar langs welke weg?

zon, temperatuur en ozon · de omweg

Een zomer boven New York, honderdelf dagen lang gemeten.¹ Op de dagen met veel zon ligt er meer ozon in de lucht; op de bewolkte dagen minder. Reken je die band uit, dan komt er een helling van 0.13 uit: elke extra streepje zonnestraling gaat, in déze zomer, samen met zo’n 0.13 ppb meer ozon.

Zon maakt ozon — het klinkt als een rechte lijn van het één naar het ander. Maar ozon ontstaat niet zómaar uit licht; het vormt zich fotochemisch, en dat gaat sneller als de lucht warm is. En zon maakt de lucht warm. Dus misschien loopt het zon-effect niet kaarsrecht naar de ozon, maar neemt het onderweg een omweg: zon → warmere lucht → meer ozon. De vraag van dit blok is niet óf zon en ozon samenhangen — dat weten we al — maar langs welke weg. Hoeveel van die 0.13 gaat rechtstreeks, en hoeveel gaat via de temperatuur?

speel het

Knip het effect in tweeën in R

paden a en b · direct en indirect effect

Eerst zoals in het boek — de twaalf diamantjes

d <- read.csv("diamantjes.csv")

# drie modellen, en let op de derde: die levert er twee paden tegelijk
# karaat → gladheid
model.a      <- lm(gladheid ~ karaat, data = d)

# karaat → glans, kaal
model.c      <- lm(glans ~ karaat, data = d)

# hierin zitten b én c'
model.beide  <- lm(glans ~ karaat + gladheid, data = d)

# stap 1: voorspelt karaat de gladheid?
# a: coef() → 3.93
coef(model.a)[2]

# stap 2: en wat doet gladheid met glans, naast karaat?
# b: coef() → 7.93
coef(model.beide)[3]

# het pad langs gladheid is a maal b
# ab, indirect: 3.93 * 7.93 → 31.16
3.93 * 7.93

# leg dat naast het totale effect en wat er direct overblijft
# c, totaal: coef() → 32.14
coef(model.c)[2]

# c', direct: coef() → 0.98
coef(model.beide)[2]

Het pad karaat → gladheid → glans draagt a × b = 31.16 van het totale effect van 32.14: 97% loopt via de slijper, en wat rechtstreeks overblijft is die 0.98 uit het blok over meervoudige regressie.

Nu jij — met de luchtmetingen die je in oefening 2.1 al zag. Bekende data, nieuwe vraag: die helling van 0.13, langs welke weg loopt die eigenlijk?

We knippen dat ene getal van 0.13 in twee stukken. Draai daarvoor drie modellen: de hele band tussen zon en ozon, dan de band tussen zon en temperatuur, en dan het model waarin zon én temperatuur samen de ozon voorspellen.

In R · trek eraan

lucht <- read.csv("luchtkwaliteit.csv")

# n = 111
lucht <- lucht[complete.cases(lucht[c("ozon","zon","temp_F")]), ]

# het totale effect
model.c     <- lm(ozon ~ zon, data = lucht)

# zon → temperatuur
model.a     <- lm(temp_F ~ zon, data = lucht)

# hierin zitten b én c'
model.beide <- lm(ozon ~ zon + temp_F, data = lucht)

# het totale effect (c): zon → ozon, kaal
# coef() → 0.127
coef(model.c)["zon"]

# pad a: zon → temperatuur
# coef() → 0.0307
coef(model.a)["zon"]

# het model met beide: hierin zitten pad b én het directe effect c'
# zon is het directe effect c', temp_F is pad b: coef() → 0.05710959  2.27846684
coef(model.beide)

Vier getallen, en ze passen op elkaar als puzzelstukjes. Het totale effect van zon op ozon is 0.127. Zet je de temperatuur ernaast, dan zakt de zon-helling naar 0.057 — dát stuk loopt rechtstreeks, buiten de warmte om. Wat ertussenuit viel, ging via de omweg: zon warmt de lucht op met 0.0307 graad per streepje (pad a), en elke graad warmer levert 2.278 ppb meer ozon op bij gelijke zon (pad b). Vermenigvuldig die twee — 0.0307 × 2.278 — en je krijgt 0.070: het stuk zon-effect dat door de temperatuur is doorgegeven. En kijk: 0.057 + 0.070 = 0.127. Direct plus indirect is precies het totaal.

trek het opnieuw

Bootstrap het indirecte pad

bootstrappen · trekken met teruglegging · 95%-interval

Dat indirecte pad van 0.070 — is dat echt, of kan zo’n omweg-getal ook uit toeval opborrelen? We willen het toetsen, maar er is een addertje: a·b is een product van twee geschatte getallen, en zulke producten hebben een scheve verdeling. Een nette symmetrische foutmarge klopt hier niet. Dus we doen het anders: we laten de data ons zélf vertellen hoe wiebelig die 0.070 is, door hem duizend keer opnieuw te trekken.

Het heet bootstrappen. Je trekt uit je 111 dagen een nieuwe steekproef van 111 dagen — mét teruglegging, dus sommige dagen komen dubbel voorbij en andere niet. Op die nieuwe steekproef herbereken je a·b. Duizend keer, duizend nét iets andere werelden, duizend waarden voor het indirecte pad. Samen vormen ze een verdeling, en in het middelste stuk daarvan — tussen de 2.5% en de 97.5% — ligt je 95%-interval.

In R · trek duizend keer opnieuw

set.seed(42)
indirect <- replicate(1000, {
  # nieuwe steekproef
  s  <- lucht[sample(nrow(lucht), replace = TRUE), ]
  a  <- coef(lm(temp_F ~ zon, data = s))["zon"]
  b  <- coef(lm(ozon ~ zon + temp_F, data = s))["temp_F"]
  # herbereken a·b
  a * b
})

# een scheef bergje, rechts uitlopend
hist(indirect)

# quantile() → 0.023  0.113
quantile(indirect, c(.025, .975))

Kijk waar dat bergje ligt: het 95%-interval loopt van 0.023 tot 0.113 — en het komt met geen enkele voet aan de nul. Duizend hertrekkingen, en het indirecte pad blijft staan. De omweg via de temperatuur is niet toeval; hij is echt.

snap het

Wat mediatie wél en niet bewijst

mediator · richting van de pijlen · gemeenschappelijke oorzaak

Wat je hier hebt gedaan heet mediatie, en de kern past in één zin: het totale effect valt uiteen in een direct deel en een indirect deel.

\[ \underbrace{0.127}_{\text{totaal } (c)} \;=\; \underbrace{0.057}_{\text{direct } (c')} \;+\; \underbrace{0.070}_{\text{indirect } (a\cdot b)} \]

Ruim de helft — 55% — van het zon-effect op ozon loopt niet rechtstreeks, maar via de temperatuur. Zon warmt de lucht, warme lucht maakt ozon. De temperatuur is de mediator: de tussenstap die een stuk van het effect doorgeeft. Dat is de vraag die mediatie stelt — hoeveel van het effect van A op C loopt via B? — en het antwoord is hier: meer dan de helft.

Maar nu het addertje, en het is een grote. Die pijlen heb jíj erin gelegd. De data tonen je drie variabelen die samenhangen — meer niet. Dat zon vóór temperatuur komt, en temperatuur vóór ozon, dat weet je niet uit de getallen; dat weet je uit de natuurkunde. De fotochemie zegt: zon warmt op, warmte maakt ozon. Die kennis leg je bovenop de cijfers, en pas dán mag je het een mediatiemodel noemen. Draai de pijlen om en de wiskunde loopt net zo hard door — de getallen protesteren niet. Een mediatiemodel is een hypothese over hoe het werkt, geen bewijs dát het zo werkt.

Waarom dat scherp houden? Neem een tegenvoorbeeld. Stel je zou “pinguïnsoort → vinlengte → gewicht” als mediatie opschrijven: soort werkt op gewicht via de vin. De rekensom levert keurige getallen op — maar de pijl klopt niet, want vinlengte veroorzaakt geen gewicht; het zijn twee dingen die de soort samen groot maakt. Daar is de “mediator” in werkelijkheid een gemeenschappelijke oorzaak, en het model liegt met een uitgestreken gezicht. Bij zon → temperatuur → ozon durven we de pijlen wél te tekenen, niet omdat de cijfers mooier zijn, maar omdat de fysica de richting kent. De statistiek levert de grootte; de richting moet ergens anders vandaan komen.

jouw beurt

Nu de omweg van de wind

negatief indirect pad · tekens van de paden

Jouw beurt

Jouw beurt — dezelfde ontleding, maar nu met de wind als startpunt. In luchtkwaliteit.csv zit ook de kolom wind (windsnelheid in mijl/uur). Meer wind gaat samen met mínder ozon; de vraag is weer: loopt dat effect rechtstreeks, of deels via de temperatuur? (Wind koelt de lucht af, koelere lucht maakt minder ozon — dat is de omweg die je gaat opmeten.)

Onder Tentamenvragen, aan het eind van deze oefening, staan dezelfde onderwerpen nog een keer in de vorm waarin een tentamen ze stelt.

Bereken eerst het totale effect van wind op ozon: bouw model.c <- lm(ozon ~ wind, data = lucht) en lees coef(model.c)["wind"] af. Ter controle kom je uit op ongeveer −5.73: elke mijl/uur meer wind gaat samen met bijna zes ppb minder ozon.
Draai nu de twee andere modellen — model.a <- lm(temp_F ~ wind, data = lucht) voor pad a, en model.beide <- lm(ozon ~ wind + temp_F, data = lucht) voor het directe effect (c′) én pad b. Reken het indirecte pad a·b uit. Ter controle: het directe effect is ongeveer −3.29 en het indirecte pad ongeveer −2.43 — samen precies de −5.73 van stap 1 (op afronding na). Ruim veertig procent van het wind-effect loopt dus via de temperatuur.
Denk-vraag: alle getallen zijn hier negatíef, en tóch klopt de optelsom direct + indirect = totaal net zo netjes als bij de zon. Waaróm mag een indirect pad negatief zijn — en wat betekent een negatief a·b eigenlijk voor de omweg die de wind neemt? (Hint: kijk naar het teken van pad a — wat doet wind mét de temperatuur? — en het teken van pad b.)

rapportage

Schrijf het op als een onderzoeker

rapportagezin · c, c′ en ab · bootstrap-CI

Rapportage

Ook een mediatie mag je in gewone-mensentaal opschrijven — het blijft dat boodschappenlijstje in zinsvorm, alleen staan er nu drie getallen op in plaats van twee: wat het totale effect was, hoeveel daarvan rechtstreeks liep, en hoeveel de omweg nam. Kijk maar:

Om te onderzoeken of de zomerzon boven New York zijn effect op de ozonconcentratie deels via de luchttemperatuur doorgaf, ontleedden we het totale zon-effect (n = 111) in een direct en een indirect deel. Het totale effect bedroeg c = 0.13; daarvan liep c′ = 0.06 rechtstreeks en ab = 0.07 via de temperatuur (95% bootstrap-CI [0.02, 0.11], 1000 hertrekkingen) — ruim de helft van het effect nam dus de omweg.

Zie hoe de zin op eigen benen staat: hij noemt de steekproef, de variabelen én de drie paden, en je kunt hem zó uit de opdracht tillen. Twee vormafspraken lopen weer mee: alléén de symbolen c, c′, ab en n zijn cursief — de gewone woorden blijven rechtop — en de getallen krijgen een punt (0.13, geen komma), netjes op twee decimalen.

Jouw beurt: schrijf zelf zo’n zelfstandige zin voor de omweg van de wind (c = −5.73, c′ = −3.29, ab = −2.43). Let op de min-tekens — een omweg mag best de andere kant op wijzen. En dat −3.29 + −2.43 op −5.72 uitkomt en niet op −5.73, is afronding, geen fout: het model rekent met alle decimalen. Hardop voorlezen mag.

tentamen

Tentamenvragen

direct, indirect en totaal · richting van de pijlen · teken van a·b

Deze zes hebben de vorm waarin een tentamen dit onderwerp stelt: vier opties, en precies één is de beste — de andere drie kunnen best waar zijn en toch niet het beste antwoord. Bij vraag 1 reken je; bij de rest kies je.

Tentamenvragen

1. Bij de zomerzon boven New York is het totale effect van zon op ozon c = 0.127 en het directe effect c′ = 0.057. Hoe groot is het indirecte effect ab?

0.031
0.070
0.184
2.278

2. In het model met zon én temperatuur staat de temperatuur op 2.278. Welke omschrijving van dat getal is de beste?

Elke graad warmere lucht gaat samen met 2.278 ppb meer ozon, bij gelijke zonnestraling.
Elk streepje zonnestraling erbij gaat samen met 2.278 graad warmere lucht.
De temperatuur geeft 2.278 keer zoveel effect op de ozon door als de zon rechtstreeks doet.
Van het totale zon-effect op de ozon loopt 2.278 ppb via de temperatuur.

3. Bij Bootstrap het indirecte pad geeft dit blok voor ab een bootstrap-interval en geen gewone foutmarge met een p. Wat is de beste reden?

Met 111 dagen is de steekproef te klein voor een toets op ab.
Een bootstrap maakt geen enkele aanname over de data, en zo’n toets wel.
ab is een product van twee schattingen, en dat ligt scheef.
Het interval laat zien dat de nul er niet in valt, en dat kan een p niet.

4. Iemand schrijft “pinguïnsoort → vinlengte → gewicht” op als mediatiemodel: de soort werkt op het gewicht via de vin. Waarom is dat geen mediatie?

Omdat vinlengte en gewicht te sterk samenhangen (r = .87) om hun paden nog uit elkaar te trekken.
Omdat pinguïnsoort een categorie is, en mediatie drie getallen nodig heeft.
Omdat de rekensom op deze data geen bruikbare getallen oplevert.
Omdat de vin het gewicht niet veroorzaakt: de soort maakt vin én gewicht groot.

5. Bij de wind als startpunt komt het indirecte pad op ab = −2.43 uit: een negatief getal. Pad a is daar −1.33 en pad b is +1.83 — die twee reken je zelf uit bij Jouw beurt. Waar komt dat minteken vandaan?

Pad a is negatief en pad b positief; min keer plus is min.
Pad a en pad b zijn allebei negatief, en min keer min is min.
Zodra het totale effect negatief is, is het indirecte pad dat ook.
Er kwam minder ozon uit dan het model voorspelde, en dat verschil draagt een minteken.

6. Van het effect van zon op ozon loopt 55% via de temperatuur. Wat mag je daaruit besluiten?

Dat de temperatuur voor de ozon belangrijker is dan de zon.
Dat de temperatuur meer dan de helft van het effect van zon op ozon veroorzaakt.
Dat 55% van het zon-effect op ozon met de temperatuur meebeweegt — meer niet.
Dat de pijlen deze kant op móéten lopen, want andersom klopt de optelsom c′ + ab = c niet meer.

Antwoord — open na je eigen poging

1 b · 2 a · 3 c · 4 d · 5 a · 6 c

Hieronder staat per vraag eerst waaróm het goede antwoord goed is, en daaronder een korte toelichting per foute optie. Heb je er een goed, dan kun je die overslaan.

1. b. Direct plus indirect is het totaal, dus ab = 0.127 − 0.057 = 0.070.

a — dat is pad a: hoeveel warmer de lucht wordt per streepje zon.
c — de twee opgeteld in plaats van afgetrokken; dan wordt de omweg groter dan het hele effect.
d — dat is pad b, en het staat in ppb per graad in plaats van per streepje zon.

2. a. Pad b is de helling van ozon op temperatuur in het model waar de zon óók in zit — dus bij gelijke zonnestraling.

b — dat is pad a (0.0307), met graden als uitkomst in plaats van als ingang.
c — 2.278 is geen keer-getal maar een aantal ppb per graad.
d — het stuk dat via de temperatuur loopt is ab = 0.070, niet pad b.

3. c. Een product van twee schattingen trekt scheef, en een symmetrische foutmarge past daar niet omheen.

a — 111 dagen is ruim genoeg; het zit in de vórm van de verdeling, niet in de omvang.
b — een bootstrap neemt nog altijd aan dat je steekproef de populatie vertegenwoordigt.
d — waar over dit interval, en toch de verkeerde reden: je kiest een bootstrap vóórdat je de uitkomst ziet, en wel omdat de verdeling scheef is. Een p-waarde kan ook prima “niet nul” zeggen.

4. d. De pijl klopt niet: de soort maakt vin en gewicht samen groot. Wat hier “mediator” heet, is een gemeenschappelijke oorzaak.

a — die r van .87 komt uit Oefening 3.2 en maakt de schatting onzeker, niet de richting fout.
b — een categorie mag prima in een mediatiemodel staan; deze pijl mag niet.
c — de rekensom levert juist keurige getallen op, en dát is het gevaar.

5. a. Min keer plus is min. Wind koelt de lucht en koelere lucht maakt minder ozon: twee stappen, en een omweg die daardoor omlaag loopt.

b — twee mintekens vóelen hier kloppend, maar min keer min is plus; dan was ab positief geweest. Pad b is +1.83: warmere lucht geeft méér ozon, of je nu met wind of met zon begint.
c — zo’n regel bestaat niet; een totaal en een omweg mogen tegengesteld van teken zijn.
d — dat beschrijft een residu, niet een pad.

6. c. De ontleding meet hoe het effect verdeeld is. Dát de pijlen die kant op lopen weet je uit de fotochemie, niet uit deze getallen.

a — pad b en het totale zon-effect op de ozon staan in verschillende eenheden en zijn niet te vergelijken.
b — het antwoord dat bijna iedereen geeft, en precies één woord te sterk: veroorzaakt.
d — draai de pijlen om en de wiskunde loopt net zo hard door; de som komt nog steeds uit.

Wil je de getallen uit deze zes vragen zelf zien ontstaan — geen ervan is overgetikt:

lucht <- read.csv("luchtkwaliteit.csv")
L <- lucht[complete.cases(lucht[c("ozon","zon","temp_F","wind")]), ]

# vraag 1 — het totale effect c → 0.127
coef(lm(ozon ~ zon, data = L))[["zon"]]

# vraag 1 — het directe effect c' → 0.057
coef(lm(ozon ~ zon + temp_F, data = L))[["zon"]]

# vraag 1 — en dus het indirecte effect ab → 0.070
coef(lm(ozon ~ zon, data = L))[["zon"]] -
  coef(lm(ozon ~ zon + temp_F, data = L))[["zon"]]

# vraag 1a — pad a, hoeveel warmer per streepje zon → 0.031
coef(lm(temp_F ~ zon, data = L))[["zon"]]

# vraag 1c — de afleider: de twee opgeteld in plaats van afgetrokken → 0.184
coef(lm(ozon ~ zon, data = L))[["zon"]] +
  coef(lm(ozon ~ zon + temp_F, data = L))[["zon"]]

# vraag 2 — pad b, ppb ozon per graad bij gelijke zon → 2.278
coef(lm(ozon ~ zon + temp_F, data = L))[["temp_F"]]

peng <- read.csv("penguins.csv")

# vraag 4d — hoe sterk vinlengte en gewicht samenhangen → 0.87
cor(peng$vinlengte, peng$gewicht)

# vraag 5 — pad a bij de wind: koeler, dus negatief → -1.33
coef(lm(temp_F ~ wind, data = L))[["wind"]]

# vraag 5 — pad b bij de wind: warmer geeft meer ozon, dus positief → 1.83
coef(lm(ozon ~ wind + temp_F, data = L))[["temp_F"]]

# vraag 5 — min keer plus: het indirecte pad → -2.43
coef(lm(temp_F ~ wind, data = L))[["wind"]] *
  coef(lm(ozon ~ wind + temp_F, data = L))[["temp_F"]]

# vraag 6 — welk deel van het zon-effect de omweg neemt, in procenten → 55.1
100 * (coef(lm(temp_F ~ zon, data = L))[["zon"]] *
       coef(lm(ozon ~ zon + temp_F, data = L))[["temp_F"]]) /
      coef(lm(ozon ~ zon, data = L))[["zon"]]

◠ Zelfde vorm in het boek

Hier: één totaal effect dat uiteenvalt in een rechtstreeks deel en een deel dat via een tussenstap wordt doorgegeven.

In het boek: W11 — Het doorgegeven verschil

Daar draait het om dezelfde ontleding — direct, indirect, totaal — en om de waarschuwing eronder: de pijlen tussen de variabelen leg jij erin, de data tonen alleen dat ze samen bewegen.

Hayes in dit blok

Eén naam gaat hier een paar keer vallen: Andrew F. Hayes. Hij is de methodoloog achter PROCESS, het gereedschap waarmee zo ongeveer iedereen mediatie en moderatie draait, en hij is de man van de modelnummers: zeg waar dan ook “model 4” en een onderzoeker weet dat je mediatie bedoelt, zeg “model 1” en het is moderatie. Maar dat is niet waarom hij hier staat. Hayes is interessant omdat hij zijn eigen gereedschap durft af te danken. De Sobel-toets, jarenlang dé manier om een indirect effect te toetsen — en die jij straks gewoon met de hand uitrekent — heeft hij zelf ten grave gedragen toen er iets beters kwam. Zo hoort een methodoloog te werken: niet trouw aan zijn gereedschap, trouw aan de vraag.

De Sobel-toets op een kladje

Je hebt het indirecte pad: a × b = 31.16. Maar twaalf steentjes is een kleine wereld — kan zo’n omweg-getal ook uit toeval opborrelen? Dat is precies wat de Sobel-toets vraagt. Hij zet het product ab af tegen zijn eigen standaardfout: dat geeft een z, en die z geeft een p. Meer is het niet — één breuk. Je hebt er geen macro voor nodig en geen menu; een kladje en de getallen die je net zelf hebt uitgerekend zijn genoeg. Pak een pen.

\[ z \;=\; \frac{a \times b}{\sqrt{\,b^2 \cdot SE_a^2 \;+\; a^2 \cdot SE_b^2\,}} \]

Invullen met wat je hebt — a = 3.93 (SE = 1.54) en b = 7.93 (SE = 1.62):

onder het wortelteken: 7.93² × 1.54² = 150.0, plus 3.93² × 1.62² = 40.4, samen 190.4
de wortel daarvan: SE van ab = 13.80
en dus z = 31.16 / 13.80 = 2.26, met p = .024

Van kladje naar PROCESS

Daar staat het: z = 2.26, p = .024. Significant, zegt Sobel. En toch is dit het moment waarop Hayes zelf zou ingrijpen. De toets leunt namelijk op een aanname: dat ab netjes normaal verdeeld is. Maar een product van twee geschatte getallen is dat niet — de verdeling trekt scheef, en bij twaalf steentjes al helemaal. Daarom geeft PROCESS voor het indirecte effect geen z en geen p, maar een bootstrap-betrouwbaarheidsinterval: duizenden keren opnieuw trekken en kijken waar ab landt. Waarom dan toch eerst dat kladje? Omdat je nu weet wát PROCESS voor je toetst. Een macro die je aanzet voordat je snapt wat hij doet, is een knop. Nu is het een verrekijker.

Extra — hetzelfde met PROCESS (model 4)

PROCESS bestaat óók voor R: één bestand van Hayes’ site, geen pakketten nodig. Staat het nog niet naast je script, dan haal je het op zoals bij PROCESS installeren beschreven staat. Inlezen met source() en dan is het één regel.

# eenmalig per sessie inladen
source("process.R")

process(data = d, y = "glans", x = "karaat", m = "gladheid",
        model = 4, boot = 5000, seed = 31216)

Je krijgt terug wat je net met de hand deed — pad a, pad b, het totale en het directe effect — plus het indirecte effect met een bootstrap-interval in plaats van een Sobel-z.

Voetnoten

airquality — dagelijkse metingen, New York, mei–september 1973 (New York State Department of Conservation & National Weather Service); via het R-basispakket datasets. Complete gevallen n = 111. Alle getallen in dit blok zijn op die échte data nagerekend.↩︎