Basishandelingen in R

de gereedschapskist — base R en dplyr naast elkaar

Waarvoor deze pagina is: de handelingen die je in élke analyse nodig hebt, in twee talen naast elkaar. Geen les — naslag. Hij staat achteraan omdat je hem niet vooraf hoeft te leren; je slaat hem open op het moment dat je hem nodig hebt. Op de plekken in de oefeningen waar een handeling voor het eerst langskomt, staat een pijltje hierheen.

Alles hieronder draait op de diamantjes: de twaalf steentjes uit het boek. Twaalf rijen die je in één oogopslag overziet, met genoeg kolommen om alles op te oefenen. Waar het verschil pas bij échte aantallen zichtbaar wordt, pakken we de 342 pinguïns erbij.

⬇ Download de data — diamantjes.csv

kolom wat
naam Anna, Bram, Cees … Leen
karaat gewicht van het steentje
glans 0–100
gladheid 0–10, hoe goed geslepen
merk golfje of sterretje
D 1 als het een sterretje is, anders 0
DxG D × gladheid — het product, voor moderatie
gepakt 1 als de kraai het steentje die nacht oppakte, anders 0
mooi 1 als de glans boven de 50 ligt, anders 0

Dat is het kortste lijstje, genoeg om verder te kunnen. Wil je van dit bestand of van een van de andere weten waar het vandaan komt, wat de eenheid is en wat er precies in staat, sla dan De datasets op.

Twee talen, één taak. Links base R — dat is R zoals hij uit de doos komt. Rechts dplyr — een pakket dat dezelfde dingen doet met werkwoorden en een pipe. De pipe leest als “en dan”: data |> filter(...) |> arrange(...) is neem de data, en dan filteren, en dan sorteren. R heeft sinds versie 4.1 zijn eigen pipe |>; %>% uit dplyr doet vrijwel hetzelfde en zie je in ouder werk.

OpmerkingEenmalig — de pakketten installeren 
install.packages("tidyverse")   # duurt een paar minuten, één keer

library(dplyr)      # rijen, kolommen, groepen
library(tidyr)      # wide ↔ long
library(ggplot2)    # plotten

Inlezen en wegschrijven

OpmerkingIn R 
d <- read.csv("diamantjes.csv")
dim(d)                            # → 12 9: twaalf steentjes, negen kolommen

# row.names = FALSE, anders schrijft R een kolom met rijnummers mee
write.csv(d, "mijn_versie.csv", row.names = FALSE)
library(readr)
d <- read_csv("diamantjes.csv")
dim(d)                            # → 12 9

write_csv(d, "mijn_versie.csv")

dim() is de goedkoopste controle die er is: twee getallen, en je weet meteen of je het goede bestand te pakken hebt. Klopt het aantal rijen niet, dan is er iets met je bestand of met je werkmap — kijk daar eerst, en ga niet verder rekenen. Een kolom te veel betekent meestal dat er een rijnummer-kolom is meegereisd van de vorige keer dat je iets wegschreef.

De valkuil van dag één: R zoekt het bestand in zijn werkmap. Weet je niet waar die is, dan getwd(). Staat het bestand ergens anders, dan geef je het hele pad — of makkelijker: zet je script en je data in dezelfde map en open het script via het bestand zelf, dan klopt de werkmap vanzelf.

Nog een verschil dat je een keer moet zien: read.csv geeft je een data.frame, read_csv een tibble. Een tibble print netter (tien rijen, met het type erbij) en zeurt eerder als je iets vaags vraagt. Verder doen ze hetzelfde, en alles op deze pagina werkt op allebei.

Kijken wat je hebt

OpmerkingIn R 
head(d)          # eerste zes rijen
nrow(d)          # → 12
names(d)         # de kolomnamen
str(d)           # kolommen + type
summary(d)       # per kolom een samenvatting
head(d)
nrow(d)
names(d)
glimpse(d)       # str, maar leesbaarder
summary(d)

Doe dit altijd eerst. Niet omdat het hoort, maar omdat de helft van alle fouten hier al zichtbaar is: een kolom die chr is waar je een getal verwachtte, een nrow die niet klopt met wat je dacht in te lezen.

Even checken

En dan de gewoonte die eronder ligt, want die is meer waard dan alle handelingen op deze bladzij bij elkaar.

Bij de spreiding reken je met afwijkingen van het gemiddelde, en die tellen altijd op tot nul. Dat is geen weetje maar een cadeautje: komt er bij jou geen nul uit, dan klopt er iets niet. Zoek voor jezelf van die momentjes om te controleren — het zou nul moeten zijn, klopt dat?

Zo’n momentje kost één regel. Je hebt net iets gedaan, je weet wat eruit zou moeten komen, en je kijkt of dat er staat:

d$karaat_c <- d$karaat - mean(d$karaat)
mean(d$karaat_c)         # → 0: dat moet, je hebt er net het gemiddelde afgehaald

Let op de vorm, want die komt op deze hele bladzij terug. Achter de regel staat een pijltje, en achter dat pijltje staat wat je hoort te zien — niet wat er toevallig uitkwam. Dat verschil is de hele truc. Schrijf je het pas op nadat je gekeken hebt, dan is elk getal goed en heb je niets gecontroleerd. Schrijf je het ervóór op, dan kan het botsen, en botsen is precies wat je wilt.

Verwacht van een controle niet dat hij bewijst dat het klopt. Hij vangt de domme helft, en de domme helft is de grote helft: een dummy die precies andersom staat, een hercodering waar drie steentjes uit vielen, een kolom die stiekem tekst is.

Vanaf hier krijgt zo’n momentje ook een eigen vorm in de tekst. Ziet hij er zo uit, dan is het geen uitleg maar een opdracht: draai dit, en kijk of je dat ziet.

Draai nrow(d) — hier hoort 12 te staan, want je hebt twaalf steentjes ingelezen. Staat er iets anders, dan lees je niet het bestand dat je denkt.

De rest van deze bladzij is bijna alleen maar zulke momentjes. Bij elke handeling die stil mis kan gaan staat er één regel bij met wat je hoort te zien.

Kolommen — kiezen, hernoemen, maken

OpmerkingIn R 
d$glans                                   # één kolom eruit
d[, c("naam", "glans")]                   # twee kolommen houden
d[, !(names(d) %in% "DxG")]               # één kolom weggooien

anders <- d                                        # hernoemen op een kópie,
names(anders)[names(anders) == "glans"] <- "schittering"   # zie de waarschuwing hieronder
names(anders)                             # → schittering staat er, glans niet meer

d$per_karaat <- d$glans / d$karaat        # nieuwe kolom → Anna: 20
d$per_karaat[1]                           # → 20: Anna heeft glans 20 en karaat 1.0
d |> pull(glans)
d |> select(naam, glans)
d |> select(-DxG)

d |> rename(schittering = glans)
d |> mutate(per_karaat = glans / karaat)

Let op het verschil. Base R verándert d als je d$per_karaat <- ... schrijft. Een dplyr-regel geeft alleen een nieuw resultaat terug; wil je dat bewaren, dan d <- d |> mutate(...). Dat is de meestgemaakte fout van week één: je draait de regel, je ziet het goede resultaat op je scherm, en even later is het weg.

Na elke mutate() waarvan je iets wilde bewaren: names(d). Staat je nieuwe kolom er niet bij, dan ben je de toewijzing vergeten. Het scherm heeft je hem wél laten zien, en dat is precies wat deze fout zo taai maakt — je hebt het resultaat met eigen ogen gezien, dus je gaat er niet meer naar zoeken.

En hernoemen verandert d wél. names(d)[...] <- "schittering" schrijft in d zelf, en daarna bestáát de kolom glans niet meer. Vraag je er verderop toch om, dan zegt base R daar niets van: d$glans is dan NULL, NULL > 50 selecteert niets, en d[d$glans > 50, ] geeft je nul steentjes waar er vijf hoorden te staan. rename() uit dplyr laat d juist met rust — vandaar de kopie hierboven. Ben je glans toch kwijt: lees d gewoon opnieuw in.

Rijen — filteren en sorteren

OpmerkingIn R 
d[d$glans > 50, ]                    # → 5 steentjes
d[d$merk == "golfje", ]              # → 6
d[d$merk == "golfje" & d$glans >= 50, ]   # twee eisen tegelijk → 3
d[d$merk == "golfje" | d$glans >= 50, ]   # of de een, of de ander

d[order(d$glans), ]                  # oplopend
d[order(-d$glans), ]                 # aflopend → Geer bovenaan
head(d[order(-d$glans), ], 3)        # de top drie: Geer, Joop, Faas
d |> filter(glans > 50)
d |> filter(merk == "golfje")
d |> filter(merk == "golfje", glans >= 50)   # komma = EN
d |> filter(merk == "golfje" | glans >= 50)  # | = OF

d |> arrange(glans)
d |> arrange(desc(glans))
d |> slice_max(glans, n = 3)         # let op: → 4 rijen, zie hieronder

“De top drie” bestaat niet altijd. Vraag je hier om de drie glanzendste, dan geeft head(d[order(-d$glans), ], 3) er netjes drie — Geer, Joop, Faas — en slice_max(glans, n = 3) er vier. Geen van beide heeft ongelijk: Faas en Ida hebben allebei een glans van 70, dus er is geen derde plaats. slice_max houdt ze allebei en zegt daarmee de waarheid; order legt ze op volgorde van hun rijnummer en head hakt er dan willekeurig één af — zonder een kik. Wil je er per se drie, zeg dan slice_max(glans, n = 3, with_ties = FALSE), zodat het jouw keuze is en niet die van het toeval in je bestand.

De komma achter de haak. In d[d$glans > 50, ] staat een komma en dan niets. Dat lege plekje betekent alle kolommen: vóór de komma kies je rijen, erachter kolommen. Vergeet je die komma, dan denkt R dat je kolommen bedoelt en krijg je een foutmelding die nergens over lijkt te gaan.

En de dubbele ==. Eén = is toewijzen, twee is vergelijken. merk = "golfje" is een opdracht, merk == "golfje" is een vraag. Je wilt de vraag.

Nieuwe variabelen maken en hercoderen

Hercoderen is: van een bestaande kolom een nieuwe maken die anders is ingedeeld. Een getal in groepen hakken, twee categorieën tot een 0/1-kolom maken, de volgorde van je categorieën vastleggen.

OpmerkingIn R · rekenen en omzetten 
# tweedeling op een grens
d$hoog <- ifelse(d$glans >= 50, "hoog", "laag")     # → 7 hoog, 5 laag
table(d$hoog, d$glans >= 50)       # → hoog alleen bij TRUE, laag alleen bij FALSE

# een getal in klassen hakken
d$klasse <- cut(d$glans,
                breaks = c(0, 40, 60, 100),
                labels = c("mat", "middel", "helder"))
table(d$klasse, useNA = "ifany")   # → mat 5, middel 3, helder 4 — en géén <NA>

# een tekstkolom tot dummy: 1 voor sterretje, 0 voor golfje
d$D2 <- ifelse(d$merk == "sterretje", 1, 0)         # gelijk aan kolom D
table(d$D2, d$merk)                # → de enen staan in de sterretje-kolom
d <- d |> mutate(
  hoog   = if_else(glans >= 50, "hoog", "laag"),
  klasse = cut(glans, breaks = c(0, 40, 60, 100),
               labels = c("mat", "middel", "helder")),
  D2     = if_else(merk == "sterretje", 1, 0)
)

cut is links-open. breaks = c(0, 40, 60, 100) maakt de vakjes (0, 40], (40, 60] en (60, 100]: de grens hóórt bij het vakje eronder. Een glans van precies 40 valt dus in “mat”, niet in “middel”.

Tel je klassen na met table(d$klasse, useNA = "ifany") en kijk of het totaal nog 12 is. Die useNA is het hele punt: valt een waarde buiten je uiterste grenzen, dan maakt cut er zwijgend een NA van, en table() laat lege plekken standaard niet zien. Zonder dat ene argument telt je tabel netjes tot elf en merk je nooit dat er een steentje uit is gevallen.

En dezelfde beweging bij een tweedeling: table(d$hoog, d$glans >= 50) zet je nieuwe kolom naast de vraag waar hij uit komt. Zeven en vijf zijn goed te verwarren met vijf en zeven; twee kolommen naast elkaar zijn dat niet.

Meer dan twee gevallen? Dan is een stapel ifelse-en onleesbaar. Er is voor allebei een nette vorm:

OpmerkingIn R · meer dan twee gevallen 
d$oordeel <- with(d, ifelse(glans >= 70, "top",
                     ifelse(glans >= 40, "prima", "matig")))
d <- d |> mutate(
  oordeel = case_when(
    glans >= 70 ~ "top",
    glans >= 40 ~ "prima",
    .default    = "matig"          # alles wat overblijft
  )
)

case_when loopt van boven naar beneden en pakt de eerste regel die klopt. Daarom hoeft de tweede regel niet “tussen 40 en 70” te zeggen: wie boven de 70 zat, is al weg.

Van tekst naar factor

Een kolom met woorden erin is voor R gewoon tekst — chr, zoals str() het noemt. Een factor is diezelfde kolom, maar dan met de mededeling erbij: dit zijn categorieën, en dit is hun volgorde. Die volgorde is niet cosmetisch. Hij bepaalt in welke volgorde je balken en boxplots staan, en welke categorie de referentie wordt zodra je gaat modelleren — dat is de groep waar je model straks vandaan meet, en waar alle andere groepen mee vergeleken worden. Verderop komt hij nog terug; voor nu is het genoeg om te weten dat de eerste categorie in je rijtje die rol krijgt.

Twaalf steentjes met twee merken zijn hier te klein om het te voelen, dus we pakken de 342 pinguïns erbij. Hun kolom soort heeft drie categorieën.

OpmerkingIn R · van chr naar factor 
peng <- read.csv("penguins.csv")

class(peng$soort)                # → "character": nog gewoon tekst

peng$soort <- factor(peng$soort)
class(peng$soort)                # → "factor", niet "character"
levels(peng$soort)               # → "Adelie" "Chinstrap" "Gentoo": Adelie is de referentie
table(peng$soort)                # → 151 / 68 / 123, samen 342

Die levels-regel is de belangrijkste van de drie, en hij is één woord lang. Zonder eigen opdracht zet R de categorieën alfabetisch neer, en de eerste in dat rijtje wordt je referentie. Hier is dat Adélie, en dat is toeval: de A ging voor. Bij merk in de diamantjes zou het golfje zijn, ook alfabetisch.

Wil je de volgorde zelf bepalen, geef hem dan mee:

OpmerkingIn R · de volgorde zelf zetten 
d$merk <- factor(d$merk, levels = c("golfje", "sterretje"))
levels(d$merk)          # → "golfje" "sterretje": golfje is nu de referentie

Draai levels() na élke factor(). Het is de enige plek waar je te zien krijgt wat R van je volgorde gemaakt heeft, en na een tikfout in een categorienaam is table() de tweede: een niveau met nul waarnemingen erin is een verschrijving, geen lege groep.

En terug naar tekst of naar getal:

OpmerkingIn R 
as.character(d$merk)
as.numeric(as.character(d$D))   # NOOIT as.numeric(factor) rechtstreeks!

De valkuil die stil misgaat: as.numeric() op een factor geeft je de nummers van de niveaus, niet je eigen getallen. Een factor met waarden 10, 20, 30 wordt dan 1, 2, 3 — geen foutmelding, gewoon andere data. Vandaar de omweg via as.character().

Het ijkpunt verzetten — relevel

De referentiecategorie is de groep waar je model vandaan meet. Alles wat er straks over de andere groepen in je tabel staat, staat er ten opzichte van díé groep. Dat maakt het alfabet een rare baas: die weet niet wat jouw onderzoeksvraag is.

Verzetten kost één regel.

OpmerkingIn R · een andere referentie kiezen 
peng$soort <- relevel(peng$soort, ref = "Gentoo")
levels(peng$soort)[1]            # → "Gentoo": nu is Gentoo het ijkpunt
levels(peng$soort)               # → "Gentoo" "Adelie" "Chinstrap"

relevel() werkt alleen op een factor. Staat de kolom nog als tekst in je bestand, dan krijg je een foutmelding — eerst factor(), dan relevel().

levels(peng$soort)[1] is de kortste controle die er is: de eerste in het rijtje ís je referentie. En let op de rest van het rijtje — relevel haalt Gentoo naar voren en laat de andere twee in hun oude volgorde staan. Hij sorteert niet, hij verhuist er één.

Zelf een dummy bouwen

Een dummy is een kolom met alleen nullen en enen: hoort deze pinguïn bij die groep, ja of nee. Drie categorieën worden zo twee dummy’s, en niet drie. De groep die overblijft is je ijkpunt, en die heeft geen eigen kolom nodig — want wie in beide dummy’s een nul heeft, is vanzelf die derde. Aan het eind van dit hoofdstukje, bij contrasts(), zie je dat zwart op wit staan.

Zelf bouwen is één vergelijking:

OpmerkingIn R · een dummy met de hand 
peng$adelie <- as.numeric(peng$soort == "Adelie")
table(peng$adelie, peng$soort)   # → 151 enen, en alleen in de Adelie-kolom
peng <- peng |> mutate(adelie = as.numeric(soort == "Adelie"))
peng |> count(adelie, soort)     # → 151 enen, en alleen bij Adelie

Die table()-regel is de beste controle op deze hele bladzij, en het is het enige stukje van dit hoofdstukje dat je echt moet onthouden. Hij zet je nieuwe kolom náást de kolom waar hij uit komt, en dan zie je in één blik of de enen op de goede plek staan.

Dat is het verschil tussen tellen en kijken, en het verschil wordt pas duidelijk als je je vergist. Stel dat je een dummy voor de Chinstraps wilde maken en per ongeluk "Gentoo" typte. Dan geeft sum() je 123 — een keurig getal, waar niets geks aan te zien is. Om die fout te vángen zou je uit je hoofd moeten weten dat er 68 Chinstraps zijn. De kruistabel vraagt dat niet: daar staan de enen gewoon in de verkeerde kolom, en dat zie je zonder iets te onthouden.

sum() vertelt je hoeveel. table() vertelt je welke. Bij het bouwen van een dummy is het tweede de vraag.

Of gewoon de factor in het model gooien

Want dat kan ook, en meestal is het de betere weg. Zet je een factor in lm(), dan maakt R de dummy’s zelf — je hoeft er geen kolom voor bij te bouwen.

Eerst even opruimen. Die relevel() van hierboven zit nog als een vlaggetje “Gentoo eerst” op je kolom geplakt, en factor() leest dat vlaggetje gewoon opnieuw. as.character() trekt het eraf.

OpmerkingIn R · twee wegen, hetzelfde model 
# vlaggetje eraf, en dan pakt factor() weer het alfabet
peng$soort <- factor(as.character(peng$soort))
levels(peng$soort)[1]                     # → "Adelie": even natrekken

model.soort <- lm(gewicht ~ soort, data = peng)
coef(model.soort)
# → intercept 3700.66 | soortChinstrap 32.43 | soortGentoo 1375.35

summary(model.soort)$r.squared             # → 0.67: soort verklaart tweederde
                                           #   van de verschillen in gewicht

contrasts(peng$soort)                      # → welke dummy's R voor je bouwde

Lees die uitvoer eens langzaam. Het intercept, 3700.66 gram, is het gemiddelde gewicht van de Adélie’s — de groep die géén eigen rij heeft. De twee rijen die er wél staan zijn afstanden vanaf dat punt: een Chinstrap weegt gemiddeld 32.43 gram méér, een Gentoo 1375.35 gram méér. Wil je het gemiddelde gewicht van de Gentoo’s weten, dan tel je op: 3700.66 + 1375.35.

Tel het na met tapply(peng$gewicht, peng$soort, mean). Er hoort 3700.66 / 3733.09 / 5076.02 uit te komen, en die derde is je intercept plus de Gentoo-rij. Tel je de twee afgeronde getallen op je rekenmachine op, dan krijg je 5076.01 — één cent verschil, want R telt door achter de komma en rondt pas daarna af. Dat soort verschil in de laatste decimaal hoort erbij. Een ánder getal hoort er niet bij: dan kijk je niet naar de referentie die je denkt.

En contrasts() laat zien wat R achter je rug om gebouwd heeft. Twee dummy’s voor drie categorieën, precies zoals je ze zelf gemaakt zou hebben:

Chinstrap Gentoo
Adelie 0 0
Chinstrap 1 0
Gentoo 0 1

Dáár is het antwoord op de vraag waarom de derde groep geen eigen kolom heeft. Adélie krijgt in beide dummy’s een nul, en dat is genoeg om hem aan te wijzen — “niet Chinstrap en niet Gentoo” laat er maar één over. De rij met alleen nullen is dus je referentie. Dezelfde informatie als levels()[1], maar dan zoals het model hem ziet.

Welke van de twee wegen je neemt, is een keuze en geen smaak. Ben zegt het in de les zo: je wilt het liever altijd zelf coderen, zodat je het kunt controleren. Let op waar de klem ligt — op dat controleren, niet op het met de hand bouwen van kolommen. Wat je niet wilt, is dat het alfabet stilletjes voor jou beslist. Er zijn dus geen twee wegen maar drie, en alleen de eerste is de verkeerde:

wie kiest de referentie wanneer je hem neemt
lm(gewicht ~ soort) zónder meer het alfabet, buiten jou om eigenlijk nooit bewust — dit is wat er gebeurt als je er niet over nadenkt
relevel() en dan lm() jij, in één regel bijna altijd: R bouwt de dummy’s, jij zegt vanaf waar er gemeten wordt
zelf dummy’s bouwen jij, door te kiezen welke kolom je niet maakt als je precies één groep tegenover de rest wilt zetten

Het verschil tussen die eerste en tweede rij is één regel code en het hele punt van dit hoofdstukje. lm(gewicht ~ soort) is niet fout — hij is alleen stil. Er stáát nergens dat Adélie je ijkpunt is; dat volgt uit de A, en dat kan morgen anders zijn als er een soort bij komt. Zet je er relevel() boven, dan staat je keuze zwart op wit in je eigen script, en dát is wat Ben bedoelt met controleren.

Verzet je referentie één keer met relevel() en draai hetzelfde model opnieuw. Alle b-waarden veranderen: met Gentoo als ijkpunt lees je 5076.02 (het gemiddelde van de Gentoo’s), −1375.35 (Adélie’s wegen zoveel mínder) en −1342.93 (Chinstraps ook). Maar R² blijft .67 — hetzelfde model, vanaf een ander nulpunt bekeken. Verandert R² wél, dan heb je per ongeluk je model veranderd en niet alleen je ijkpunt.

Centreren

Centreren is één aftrekking: je haalt het gemiddelde van elke waarde af. Het verandert niets aan je gegevens en alles aan wat je aflezen kunt.

OpmerkingIn R 
d$karaat_c <- d$karaat - mean(d$karaat)     # mean(karaat) = 1.6
mean(d$karaat_c)                            # → 0, per constructie

# hetzelfde, met de ingebouwde functie
d$karaat_c <- as.numeric(scale(d$karaat, scale = FALSE))
d <- d |> mutate(karaat_c = karaat - mean(karaat))

# meerdere kolommen tegelijk
d <- d |> mutate(across(c(karaat, gladheid), \(x) x - mean(x),
                        .names = "{.col}_c"))

Let op de scale = FALSE. scale(x) doet standaard twéé dingen: centreren én delen door de standaardafwijking. Dat tweede maakt er een z-score van, en dan lees je je helling in standaardafwijkingen in plaats van in karaat. Wil je alleen centreren, dan moet die scale = FALSE erbij. En scale() geeft een matrix terug, geen kolom — vandaar as.numeric() eromheen.

Waarom je het doet

Het intercept van een regressie is de voorspelling bij x = 0. Zolang x rauw is, is dat vaak een steentje van nul karaat — een ding dat niet bestaat.

OpmerkingIn R · kijk wat er verschuift 
lm(glans ~ karaat,   data = d)    # → intercept −1.43, helling 32.14
lm(glans ~ karaat_c, data = d)    # → intercept  50.00, helling 32.14

De helling is exact hetzelfde — 32.14 glanspunten per karaat, of je nu centreert of niet. Alleen het intercept verhuist: van −1.43 (de glans van een onbestaanbaar steentje van nul karaat, en nog negatief ook) naar 50.00 — en dat is precies mean(d$glans), de glans van een gemiddeld steentje. Van een onzin-getal naar een getal dat je kunt voorlezen.

Bij een interactie wordt het meer dan cosmetiek:

OpmerkingIn R · en bij een interactie 
d$gladheid_c <- d$gladheid - mean(d$gladheid)      # mean = 5

lm(glans ~ gladheid   * D, data = d)   # D-effect →  1.84
lm(glans ~ gladheid_c * D, data = d)   # D-effect → 15.03

Dezelfde data, hetzelfde model, hetzelfde R² (.97) en dezelfde interactieterm (2.64). Maar het merkverschil leest de ene keer als 1.84 en de andere keer als 15.03. Allebei waar — ze beantwoorden alleen een andere vraag. Zonder centreren is het “het merkverschil bij gladheid = 0”, en gladheid nul komt in deze twaalf steentjes niet voor: dat is een uitspraak over een plek buiten je data. Gecentreerd is het “het merkverschil bij gemiddelde gladheid”, en dat is een steentje dat je zó kunt aanwijzen.

De regel om te onthouden: in een model mét interactie zijn de hoofdeffecten altijd effecten bij nul van de ander. Centreer je voorspellers en die nul wordt het gemiddelde — dan pas is de tabel te lezen zonder erbij te rekenen. De interactie zelf verandert nooit.

Ontbrekende waarden — NA

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

sum(is.na(peng$sekse))            # → 9
colSums(is.na(peng))              # per kolom, in één blik

mean(c(1, NA, 3))                 # → NA
mean(c(1, NA, 3), na.rm = TRUE)   # → 2

peng[!is.na(peng$sekse), ]        # → 333 rijen
peng |> summarise(across(everything(), \(x) sum(is.na(x))))

peng |> filter(!is.na(sekse))     # → 333 rijen
peng |> drop_na(sekse)            # zelfde, uit tidyr

NA is besmettelijk, en dat is expres. Eén ontbrekende waarde maakt het hele gemiddelde NA — R weigert te doen alsof hij iets weet wat hij niet weet. Met na.rm = TRUE zeg je: laat weg wat er niet is. Doe dat bewust, en noteer hoeveel je weglaat, want het is geen nul: negen pinguïns zonder sekse zijn geen negen pinguïns zonder gewicht. is.na(x) is de enige manier om erop te toetsen — x == NA geeft altijd NA terug, ook als het antwoord “ja” had moeten zijn.

Samenvatten, en per groep

OpmerkingIn R 
mean(d$glans)                                  # → 50
sd(d$glans)                                    # → 21.74

aggregate(glans ~ merk, data = d, FUN = mean)  # → golfje 40, sterretje 60
tapply(d$glans, d$merk, mean)                  # zelfde, andere vorm
table(d$merk)                                  # → golfje 6, sterretje 6

# twee groeperende kolommen: 3 soorten × 2 sekses → 6 rijen
aggregate(gewicht ~ soort + sekse, data = peng, FUN = mean)
d |> summarise(gem = mean(glans), sd = sd(glans))

d |> group_by(merk) |>
  summarise(gem = mean(glans), n = n())

d |> count(merk)

peng |> group_by(soort, sekse) |>
  summarise(gem = mean(gewicht), n = n(), .groups = "drop")

Dit is waar de pipe zijn geld verdient. “Neem de diamantjes, en dan per merk, en dan het gemiddelde” staat er in die volgorde — precies zoals je het zegt.

Zet er altijd n = n() bij. Een gemiddelde zonder aantal is een half getal: 40 uit zes steentjes en 40 uit één steentje zien er in de tabel identiek uit. En vergeet .groups = "drop" niet als je met twee kolommen groepeert, anders blijft het resultaat gegroepeerd en gedraagt de volgende regel zich raar.

Alles tegen alles — de correlatiematrix

cor(x, y) geeft je één samenhang tussen twee variabelen. Wil je ze allemaal tegelijk zien, dan geef je cor() een hele tabel in plaats van twee kolommen. Dat klinkt als gemak, en het is ook gemak — maar het is de handeling op deze bladzij waar de meeste onzin uit komt, en dat is precies waarom hij hier staat.

Eerst wat er gebeurt als je het naïef probeert:

OpmerkingIn R · de matrix 
peng <- read.csv("penguins.csv")   # vers, zonder de kolommen van hierboven

cor(peng)                          # → Error: 'x' must be numeric

# dus eerst de numerieke kolommen eruit
num <- peng[sapply(peng, is.numeric)]
names(num)                         # → id, snavellengte, snaveldiepte, vinlengte, gewicht
round(cor(num, use = "complete.obs", method = "pearson"), 2)
peng <- read_csv("penguins.csv")   # vers, zonder de kolommen van hierboven

num <- peng |> select(where(is.numeric))
names(num)                         # → id, snavellengte, snaveldiepte, vinlengte, gewicht
round(cor(num, use = "complete.obs", method = "pearson"), 2)

Die foutmelding is een vriend. soort en sekse zijn woorden, en een correlatie tussen woorden bestaat niet — R weigert netjes. Wat je terugkrijgt is dit:

id snavellengte snaveldiepte vinlengte gewicht
id 1.00 0.77 −0.26 0.46 0.28
snavellengte 0.77 1.00 −0.24 0.66 0.60
snaveldiepte −0.26 −0.24 1.00 −0.58 −0.47
vinlengte 0.46 0.66 −0.58 1.00 0.87
gewicht 0.28 0.60 −0.47 0.87 1.00

De diagonaal is overal 1.00 — elke variabele met zichzelf — en de matrix is symmetrisch: linksonder staat hetzelfde als rechtsboven. Je hoeft dus maar één helft te lezen.

Kijk eerst naar die diagonaal en naar de spiegeling. Staat er ergens op de diagonaal iets anders dan 1.00, of is de matrix niet symmetrisch, dan kijk je niet naar een correlatiematrix maar naar iets anders. En tel de rijen: hier horen er 5 te zijn, evenveel als er numerieke kolommen zijn.

Heb je die regel gedraaid zonder peng opnieuw in te lezen, dan zijn het er zes: de dummy-kolom adelie van een paar hoofdstukjes terug is een getal, dus is.numeric pakt hem mee. Dat is geen fout van R — het is de eerste keer dat je ziet wat het probleem van deze handeling is. Daar gaat het volgende stukje over.

Waarom je je kolommen met de hand kiest

Kijk nog eens naar de eerste rij. id correleert .77 met snavellengte — sterker dan snavellengte met gewicht (.60), en bijna zo sterk als vinlengte met gewicht (.87). Als dat waar was, zou het volgnummer van een pinguïn iets zeggen over zijn snavel.

Het is niet waar. id is het rijnummer in het bestand, en het bestand staat gesorteerd op soort: de Adélie’s staan op 1 tot en met 151, de Gentoo’s op 152 tot en met 274, de Chinstraps op 275 tot en met 342. En die drie soorten hebben toevallig een oplopende gemiddelde snavellengte — 38.79, 47.50 en 48.83 millimeter. Loop je van boven naar beneden door het bestand, dan lopen de snavels mee omhoog. Meer is er niet aan de hand.

Drie platte wolkjes, samen een trap: Adélie de rondjes, Gentoo de driehoekjes, Chinstrap de vierkantjes, en het gestippelde lijntje per soort is zijn gemiddelde — de vlakke trede. Binnen elke soort zegt het volgnummer vrijwel niets over de snavellengte, maar het bestand staat op soort gesorteerd en de soorten verschillen — dus klimt er één lijn door alles heen naar r = .77. (Voorlopige werkfiguur — vervangt Bens eigen tekening.)

Snijd het per soort door en er blijft bijna niets van over:

OpmerkingIn R · dezelfde samenhang, binnen de soorten 
ad <- peng[peng$soort == "Adelie", ]
ch <- peng[peng$soort == "Chinstrap", ]
ge <- peng[peng$soort == "Gentoo", ]

cor(ad$id, ad$snavellengte, method = "pearson")   # → 0.01
cor(ch$id, ch$snavellengte, method = "pearson")   # → 0.02
cor(ge$id, ge$snavellengte, method = "pearson")   # → 0.21

Nul komma nul één, nul komma nul twee, nul komma twee één. Die laatste is niet precies nul — maar .21 tegenover .77 is het verschil tussen “misschien iets” en “dit is een van de sterkste verbanden in mijn tabel”. De .77 zat niet in de pinguïns maar in de volgorde van het bestand, en een correlatiematrix zegt daar niets over. Hij rekent gehoorzaam elke kolom tegen elke andere, ook de kolommen die geen meting zijn.

Daarom is sapply(peng, is.numeric) een begin en geen antwoord. Die test vindt alles wat een getal is, en een volgnummer is een getal. Zo ook een postcode, een respondentnummer, een jaartal en een gecodeerd geslacht. Kies je kolommen dus met de hand:

OpmerkingIn R · de kolommen die je bedoelt 
metingen <- peng[, c("snavellengte", "snaveldiepte", "vinlengte", "gewicht")]
round(cor(metingen, use = "complete.obs", method = "pearson"), 2)

Loop je matrix één keer af en vraag bij elke rij: is dit een meting? Bij id is het antwoord nee, en dan hoort die rij er niet te staan. Dit is de enige controle op deze bladzij die je niet kunt automatiseren, en daarom de enige die je echt zelf moet doen.

use = — wat je met de ontbrekende waarden doet

Zodra er ergens een NA in je tabel zit, moet cor() weten wat hij ermee aan moet. Er zijn drie antwoorden en ze kosten alle drie iets anders.

Op de pinguïns zie je het verschil niet: de negen ontbrekende waarden zitten in sekse, en dat is tekst, dus die kolom staat niet eens in de matrix. De vier metingen zijn compleet, en dan geven de drie opties precies hetzelfde. Voor het verschil pakken we luchtkwaliteit.csv erbij — 116 dagen, waarvan er 5 geen zonmeting hebben.

OpmerkingIn R · de drie opties 
lucht <- read.csv("luchtkwaliteit.csv")
nl    <- lucht[sapply(lucht, is.numeric)]
sum(!complete.cases(nl))              # → 5 dagen missen iets

round(cor(nl, use = "everything"), 2)             # de standaard
round(cor(nl, use = "complete.obs"), 2)
round(cor(nl, use = "pairwise.complete.obs"), 2)

Vijf numerieke kolommen heeft dat bestand — ozon, zon, temp_F, wind en maand — dus zon heeft vier partners om mee te correleren. Onthoud die vier, want ze komen in de eerste regel van de tabel terug.

use = wat hij doet wat het kost
"everything" rekent gewoon door, met NA erin de hele rij en kolom zon wordt NA — vijf ontbrekende metingen maken alle vier de correlaties van zon onleesbaar
"complete.obs" gooit elke dag weg die érgens iets mist alles rust op dezelfde 111 dagen, maar je gooit ook vijf complete ozon-metingen weg
"pairwise.complete.obs" per paar wat er is ozon × temperatuur op 116 dagen, alles met zon op 111 — de cellen rusten dus op verschillende aantallen

Die laatste is verleidelijk, want hij lijkt het meeste te bewaren. Kijk dan naar wat er gebeurt met de samenhang tussen ozon en maand: onder complete.obs leest hij .14, onder pairwise.complete.obs .16. Dezelfde data, dezelfde vraag, en toch twee getallen — want het zijn twee verschillende groepen dagen.

Dat is nog te overzien bij vijf ontbrekende metingen. Bij een vragenlijst waar iedereen ergens anders een vraag oversloeg, staan de cellen van je matrix op steekproeven die soms honderden mensen schelen, en dan kun je twee getallen uit dezelfde tabel niet meer met elkaar vergelijken.

Het datablad als rooster met gaten: 116 dagen, en vijf ervan — waarvan drie aaneengesloten — missen hun zonmeting. complete.obs gooit die vijf hele dagen weg — de lichte banden — zodat elke cel van de matrix op dezelfde 111 dagen rust. (Voorlopige werkfiguur — vervangt Bens eigen tekening.)

Kies dus complete.obs, en meld hoeveel je weglaat. Dan staat je hele matrix op dezelfde mensen en zijn twee cellen weer met elkaar te vergelijken. Dat is één regel in je verslag — “vijf van de 116 dagen misten een meting en zijn buiten de correlaties gelaten” — en die regel is meer waard dan de vijf dagen die je ermee koopt.

sum(!complete.cases(nl)) zegt hoeveel rijen dat je kost. Hier 5, van de 116.

Wil je zien waar de gaten precies zitten, dan is er nog één regel:

OpmerkingIn R · waar rust elke cel op? 
crossprod(!is.na(nl))              # → per paar het aantal dagen dat meetelde

Die geeft een tabelletje in de vorm van je correlatiematrix, maar met aantallen in plaats van correlaties: ozon × temperatuur op 116 dagen, alles met zon op 111. Je hoeft crossprod niet te snappen om hem te gebruiken — draai hem, en kijk of alle getallen ongeveer gelijk zijn. Lopen ze ver uiteen, dan rust je matrix op verschillende steekproeven en is complete.obs geen luxe meer.

Twee tabellen aan elkaar — merge en join

Stel dat je er een tweede tabelletje bij hebt: waar elk merk geslepen wordt.

OpmerkingIn R 
slijperijen <- data.frame(
  merk      = c("golfje", "sterretje"),
  slijperij = c("De Vijl", "Sterrenstof"),
  plaats    = c("Antwerpen", "Amsterdam")
)

# by = de kolom die in beide tabellen staat
samen <- merge(d, slijperijen, by = "merk")
nrow(samen)    # → 12
slijperijen <- tibble(
  merk      = c("golfje", "sterretje"),
  slijperij = c("De Vijl", "Sterrenstof"),
  plaats    = c("Antwerpen", "Amsterdam")
)

samen <- d |> left_join(slijperijen, by = "merk")
nrow(samen)    # → 12

De sleutel is de kolom die in beide tabellen staat — hier merk. Wat je moet weten: merge houdt standaard alleen rijen die in allebei voorkomen, left_join houdt álle rijen van links en zet NA waar rechts niets had.

Laat het jezelf één keer zien met een half tabelletje, waarin alleen golfje staat:

OpmerkingIn R · het verschil, zichtbaar 
half <- slijperijen[1, ]                    # alleen golfje

nrow(merge(d, half, by = "merk"))           # →  6   de sterretjes zijn wég
nrow(d |> left_join(half, by = "merk"))     # → 12   met NA bij de sterretjes

Zes rijen verdwenen zonder één waarschuwing. Bij twaalf steentjes zie je dat meteen; bij tienduizend rijen niet, en dan is het het verschil tussen kloppen en stilletjes je halve steekproef kwijtraken. Tel dus altijd je rijen vóór en na een koppeling.

Wide en long — dezelfde data, andere vorm

Nu staat elk steentje op één rij met drie eigenschappen naast elkaar: dat is wide. Eén rij per steentje, één kolom per eigenschap. Prima om naar te kijken.

Maar zodra je wilt plotten, wil ggplot het andersom. En daar zit de hele reden voor dit hoofdstukje: ggplot zet één kolom op de x-as en één op de y-as. Wil je drie eigenschappen náást elkaar in één plaatje, dan moeten die drie kolommen eerst één kolom worden — met ernaast een kolom die zegt wélke eigenschap het is. Dat is long: één rij per meting in plaats van per steentje.

OpmerkingIn R 
# reshape wil precies weten wat er kantelt en hoe de nieuwe kolommen heten.
# as.data.frame() omdat reshape base R is en een data.frame wil: heb je
# hierboven het dplyr-tabblad gedraaid, dan is `d` een tibble en geeft dit
# blok er stil drie in plaats van zesendertig terug.
lang <- reshape(
  as.data.frame(d)[, c("naam", "karaat", "glans", "gladheid")],
  direction = "long",
  varying   = c("karaat", "glans", "gladheid"),
  v.names   = "waarde",
  timevar   = "eigenschap",
  times     = c("karaat", "glans", "gladheid"),
  idvar     = "naam"
)
nrow(lang)    # → 36  (12 steentjes × 3 eigenschappen)

breed <- reshape(lang, direction = "wide",
                 idvar = "naam", timevar = "eigenschap")
nrow(breed)   # → 12  weer terug
library(tidyr)

# drie kolommen worden twee: een naam-kolom en een waarde-kolom
lang <- d |>
  select(naam, karaat, glans, gladheid) |>
  pivot_longer(cols      = c(karaat, glans, gladheid),
               names_to  = "eigenschap",
               values_to = "waarde")
nrow(lang)    # → 36

breed <- lang |>
  pivot_wider(names_from = eigenschap, values_from = waarde)
nrow(breed)   # → 12

Zo ziet de kanteling eruit — links wide, rechts long:

naam karaat glans gladheid
Anna 1.0 20 1
Bram 1.1 30 2

wordt

naam eigenschap waarde
Anna karaat 1.0
Anna glans 20
Anna gladheid 1
Bram karaat 1.1
… … …

Waar iedereen op vastloopt: je moet drie dingen benoemen — wélke kolommen gaan kantelen, hoe de kolom gaat heten die hun namen opvangt, en hoe die met hun waarden. Teken het één keer op papier: twaalf rijen worden zesendertig, vier kolommen worden er drie. Dan klopt het beeld voorgoed.

En nu waarvoor je het deed

Dit kán niet in de wide-vorm, en in de long-vorm is het één regel:

OpmerkingIn R · drie eigenschappen, drie vakjes 
# elk vakje zijn eigen verdeling — let op scales = "free_x",
# want karaat loopt tot 2 en glans tot 80
ggplot(lang, aes(waarde)) +
  geom_histogram(bins = 6) +
  facet_wrap(~ eigenschap, scales = "free_x")

# of elk steentje als een lijntje langs zijn drie eigenschappen
ggplot(lang, aes(eigenschap, waarde, group = naam)) +
  geom_line(alpha = .5) +
  geom_point()

Dát is het motief: long maak je niet omdat het netter is, maar omdat de plot erom vraagt. Elke keer dat je in ggplot een kleur, een facet of een lijn per groep wilt en die “groep” nu nog een rij kolomnamen is — dan is het tijd om te kantelen. Klaar met plotten? Dan mag je gerust terug naar wide om naar te kijken.

Plotten

OpmerkingIn R · de drie basisplaatjes 
model.karaat <- lm(glans ~ karaat, data = d)   # eerst het model, dan pas tekenen

plot(d$karaat, d$glans)                    # puntenwolk
abline(model.karaat)                       # de lijn erdoor

hist(d$glans)
boxplot(glans ~ merk, data = d)
library(ggplot2)

ggplot(d, aes(karaat, glans)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE)

ggplot(d, aes(glans)) + geom_histogram(bins = 6)
ggplot(d, aes(merk, glans)) + geom_boxplot()

Het verschil in één zin: base tekent, ggplot bouwt op. In ggplot zeg je eerst wélke data en wélke kolom op welke as (aes), en daarna stapel je lagen met +. Dat voelt omslachtig bij één plaatje en wint zodra je er iets bij wilt — en dat “iets bij” is de rest van dit hoofdstukje.

Onthoud de vorm, dan onthoud je ggplot:

ggplot(DATA, aes(WELKE KOLOM WAAR)) +   # wat waar
  geom_...()                        +   # welke vorm punten/balken/lijnen
  labs(...)                         +   # hoe het heet
  scale_...(...)                    +   # hoe de assen en kleuren lopen
  facet_...(...)                    +   # in hoeveel vakjes
  theme_...()                           # hoe het eruitziet

Punten, lijnen, balken

OpmerkingIn R 
ggplot(d, aes(karaat, glans)) + geom_point()                 # stippen
ggplot(d, aes(karaat, glans)) + geom_point() + geom_line()   # stippen én lijn
ggplot(d, aes(karaat, glans)) + geom_point() +
  geom_smooth(method = "lm", se = FALSE)                     # de rechte erdoor
ggplot(d, aes(karaat, glans)) + geom_point() + geom_smooth() # een buigende trend

ggplot(d, aes(merk)) + geom_bar()                            # tellen
ggplot(d, aes(merk, glans)) + geom_boxplot()
ggplot(d, aes(merk, glans)) + geom_boxplot() +
  geom_jitter(width = .1)                                    # de punten erbovenop

Punten of lijnen? Een lijn zegt: deze punten horen achter elkaar. Bij tijd of een oplopende dosis klopt dat; bij twaalf losse steentjes suggereert een lijn een volgorde die er niet is. geom_smooth() is iets anders dan geom_line() — die tekent geen verbinding tussen jouw punten maar een geschatte trend eróverheen. En se = FALSE haalt het grijze bandje weg; laat het staan als je de onzekerheid wilt tonen, wat meestal eerlijker is.

Assen en labels

OpmerkingIn R 
ggplot(d, aes(karaat, glans)) + geom_point() +
  labs(title    = "Zwaarder is glanzender",
       subtitle = "twaalf steentjes",
       x        = "gewicht (karaat)",
       y        = "glans (0–100)",
       colour   = "merk",              # ook de legenda-titel zet je hier
       caption  = "bron: diamantjes.csv")

# waar de as begint en eindigt, en waar de streepjes staan
ggplot(d, aes(karaat, glans)) + geom_point() +
  scale_x_continuous(limits = c(0, 2.5)) +
  scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, by = 25))

# assen omdraaien (handig bij lange categorie-namen)
ggplot(d, aes(merk, glans)) + geom_boxplot() + coord_flip()

De valkuil met limits. scale_y_continuous(limits = ...) gooit punten weg die erbuiten vallen — en dan klopt je trendlijn niet meer, want die wordt berekend over wat er overblijft. Wil je alleen in- of uitzoomen zonder data te verliezen, gebruik dan coord_cartesian(ylim = c(0, 100)). Die snijdt het plaatje bij en laat de berekening met rust.

Zet je limits, tel dan hoeveel punten er buiten vallen: sum(d$glans < 0 | d$glans > 100) hoort 0 te zijn. Is het meer dan nul, dan staan er nu minder steentjes in je plaatje dan in je data, en die zijn ook uit je trendlijn verdwenen. ggplot zegt het wel — “Removed n rows containing missing values” — maar dat is één grijze regel in je console tussen de rest, en die leest niemand.

Kleuren en vormen

OpmerkingIn R 
# één vaste kleur voor alles: BUITEN aes()
ggplot(d, aes(karaat, glans)) +
  geom_point(size = 3, colour = "steelblue", alpha = .8)

# kleur die iets betekent: BINNEN aes()
ggplot(d, aes(karaat, glans, colour = merk)) + geom_point(size = 3)
ggplot(d, aes(karaat, glans, shape  = merk)) + geom_point(size = 3)

# je eigen kleuren
ggplot(d, aes(karaat, glans, colour = merk)) + geom_point(size = 3) +
  scale_colour_manual(values = c(golfje = "#4A5D7E", sterretje = "#C2A868"))

Dit is dé ggplot-verwarring, en hij is de moeite waard om één keer goed te snappen. Staat colour = binnen aes(), dan is het een afbeelding: de kleur staat ergens voor, R maakt er een legenda bij. Staat het buiten aes(), in de geom_, dan is het gewoon verf. aes(colour = "rood") maakt daarom oranje-rode punten met een legenda waarin het woord “rood” staat — R denkt dat je een kolom met de waarde “rood” bedoelde. Zit een kleur ergens vast waar je hem niet wilt: kijk of hij aan de goede kant van dat haakje staat.

colour is de rand of de stip, fill de vulling. Bij balken, boxplots en histogrammen wil je meestal fill.

Facetten — hetzelfde plaatje, meerdere keren

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

ggplot(peng, aes(snavellengte, snaveldiepte)) +
  geom_point() +
  facet_wrap(~ soort)                     # één vakje per soort

ggplot(peng, aes(snavellengte, snaveldiepte)) +
  geom_point() +
  facet_grid(sekse ~ soort)               # rijen × kolommen

# elke as z'n eigen bereik (alleen als de schalen echt verschillen)
facet_wrap(~ eigenschap, scales = "free_x")

Facetten zijn de eerlijke versie van kleuren zodra het druk wordt. Drie soorten in drie vakjes laten zien wat drie kleuren in één wolk verstoppen — je zag dat in de oefening over correlatie, waar de wolk over alle soorten heen omláág hing en binnen elke soort omhoog.

Let op scales = "free". Standaard krijgt elk vakje dezelfde assen, en dat is precies wat je meestal wilt: alleen dan zijn de vakjes onderling te vergelijken. Zet je de assen vrij, dan lijken drie totaal verschillende wolken ineens even groot. Doe het alleen als de eenheden echt uiteenlopen — zoals bij de drie eigenschappen hierboven, waar karaat tot 2 loopt en glans tot 80.

Uiterlijk en opslaan

OpmerkingIn R 
ggplot(d, aes(karaat, glans, colour = merk)) + geom_point(size = 3) +
  theme_minimal() +                        # ook: theme_bw(), theme_classic()
  theme(legend.position = "bottom")        # of "none" om hem weg te halen

# opslaan: ggsave pakt standaard de laatste plot
ggsave("wolk.png", width = 12, height = 8, units = "cm", dpi = 300)

Maak je plaatje op het formaat waarop het gebruikt wordt. Wat op je scherm leesbaar is, is het in een verslag van 12 bij 8 centimeter meestal niet — de letters schalen niet mee met het beeld. Geef width en height mee aan ggsave en kijk naar het bestand, niet naar het RStudio-venster. Voor een verslag of een poster is dpi = 300 het minimum.

De hele boog op één dataset

Waarom juist deze twaalf steentjes: er zit precies genoeg in om alles van het boek op te draaien.

OpmerkingIn R · van gemiddelde tot moderatie 
mean(d$glans)                             # H1  het gemiddelde → 50
sd(d$glans)                               # H1  de spreiding → 21.74
cor(d$karaat, d$glans)                    # H3  samenhang → .578
lm(glans ~ karaat, data = d)              # H3  de lijn → −1.43 + 32.14 · karaat
lm(glans ~ karaat + gladheid, data = d)   # H8  en dan valt karaat weg
lm(glans ~ D, data = d)                   # H9  groepen als getallen → 40 + 20 · D
lm(glans ~ gladheid * D, data = d)        # H10 verschil in verschil

Die vijfde regel is het hele hoofdstuk over confounding in één opdracht: los gezien voorspelt karaat de glans (b = 32.14), maar zodra gladheid meedoet, zakt karaat naar bijna nul (b = 0.98) en verklaart het model ineens 82% in plaats van 33%. Niet het gewicht maakte ze mooi — de slijper deed dat.