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
model.a      <- lm(gladheid ~ karaat, data = d)             # karaat → gladheid
model.c      <- lm(glans ~ karaat, data = d)                # karaat → glans, kaal
model.beide  <- lm(glans ~ karaat + gladheid, data = d)     # hierin zitten b én c'

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

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

# het pad langs gladheid is a maal b
3.93 * 7.93               # ab = 31.16  indirect

# leg dat naast het totale effect en wat er direct overblijft
coef(model.c)[2]          # c  = 32.14  totaal
coef(model.beide)[2]      # c' = 0.98   direct

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")
lucht <- lucht[complete.cases(lucht[c("ozon","zon","temp_F")]), ]  # n = 111

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

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

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

# het model met beide: hierin zitten pad b én het directe effect c'
coef(model.beide)
#   zon      0.057   ← direct effect (c')
#   temp_F   2.278   ← pad b

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, {
  s  <- lucht[sample(nrow(lucht), replace = TRUE), ]   # nieuwe steekproef
  a  <- coef(lm(temp_F ~ zon, data = s))["zon"]
  b  <- coef(lm(ozon ~ zon + temp_F, data = s))["temp_F"]
  a * b                                                 # herbereken a·b
})
hist(indirect)                          # een scheef bergje, rechts uitlopend
quantile(indirect, c(.025, .975))       # 0.023  tot  0.113

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.)

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. Hardop voorlezen mag.

◠ 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.

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

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.↩︎