24 Het schatten van parameters, en betrouwbaarheidsintervallen

Nu je de basis van kansrekenen onder de knie hebt, is het tijd je kennis toe te passen op de statistiek. Laten we even terugpakken waar we gebleven waren.

In Hoofdstuk 21 hebben we het gehad over schattingen. Het doel van veel onderzoeken is om een bepaalde grootheid of populatieparameter te schatten. Bijna alle getallen die je in de krant of in de wetenschappelijk literatuur tegenkomt, zijn zulke schattingen.

We hebben benadrukt dat schattingen onzeker zijn. Bij schattingen op basis van steekproeven komt dat door steekproefvariabiliteit en mogelijk door bias in de samenstelling van de steekproef. Bij metingen spelen systematische en toevallige meetfouten een rol. Daarom is het heel belangrijk om je bij elk getal dat je tegenkomt af te vragen hoe onzeker dat getal is.

Een getal zonder indicatie van de nauwkeurigheid en precisie is nutteloos.

In dit hoofdstuk leer je voor een specifieke situatie hoe je die onzekerheid kwantificeert. We doen dit door betrouwbaarheidsintervallen uit te rekenen. In de volgende paragraaf leggen we uit wat dat zijn.

24.1 Betrouwbaarheidsintervallen

Het idee van een betrouwbaarheidsinterval (BHI, confidence interval) is dat we bij het maken van een schatting niet alleen één getal berekenen (een zogenaamde puntschatting) maar ook een interval om dat punt dat aangeeft wat de onzekerheid in die schatting is. In deze cursus en daarbuiten wordt vaak een 95%-BHI gebruikt. Een 95%-BHI is een interval dat zo berekend wordt dat de ware waarde van de parameter in 95% van de gevallen binnen dat interval ligt.

Wat is een BHI?

Een 95%-BHI is een interval rondom een puntschatting, berekend op zo’n manier dat in 95% van de gevallen de ware parameter binnen dat interval ligt.

We kunnen dit concreter maken met een voorbeeld. Stel je voor dat verschillende onderzoekers geïnteresseerd zijn in dezelfde parameter, namelijk de gemiddelde chlorofylconcentratie in de bladeren van tarweplanten. Iedere onderzoeker neemt een steekproef van 10 planten en meet daarin de chlorofylconcentratie.

Laten we even aannemen dat in werkelijkheid het gemiddelde voor deze planten (de populatieparameter) gelijk is aan 2,0 mg per gram vers bladgewicht. De onderzoekers weten dat niet; die hebben alleen toegang tot hun eigen steekproef. Iedere onderzoeker berekent de gemiddelde chlorofylconcentratie van de eigen steekproef. Dat is de puntschatting. Iedere onderzoeker krijgt door steekproefvariabiliteit dus een andere schatting uit zijn experiment. Die schatting zal over het algemeen afwijken van het werkelijke gemiddelde voor deze planten.

Maar iedere onderzoeker berekent ook een 95%-BHI. Figuur 24.1 laat zien wat het resultaat zou kunnen zijn. Van alle 20 steekproeven zijn de waarnemingen in het grijs geplot. De gemiddelden van alle 20 steekproeven (de puntschattingen) zijn als paarse punten weergegeven, de ware populatieparameter als een oranje horizontale lijn. Zwarte lijnen geven voor iedere puntschatter het bijbehorende BHI aan.

Ieder BHI bevat met 95% zekerheid het ware gemiddelde. Dat betekent helaas ook dat 5% van de intervallen het gemiddelde niet bevatten. Bij 20 herhalingen komt dat gemiddeld één keer voor. In dit geval heeft onderzoeker 9 pech: hoewel deze zijn onderzoek net zo zorgvuldig heeft uitgevoerd als de andere onderzoekers, geeft zijn puntschatting en betrouwbaarheidsinterval een verkeerde indruk.

Code

# Laad ggplot2 indien nodig
if (!("ggplot2" %in% loadedNamespaces())) {
  library(ggplot2)
}

# Parameters voor de populatieverdeling
waar_gemiddelde <- 2.0
waar_sd <- 0.4
aantal_steekproeven <- 20
steekproefgrootte <- 10  # Aantal datapunten per steekproef

# Genereer 20 aselecte steekproeven
set.seed(125)  # Voor reproduceerbaarheid
steekproeven <- replicate(
  aantal_steekproeven,
  rnorm(steekproefgrootte, mean = waar_gemiddelde, sd = waar_sd),
  simplify = FALSE
)

# Maak een data frame voor de individuele datapunten
datapunten <- data.frame(
  steekproef_id = rep(1:aantal_steekproeven, each = steekproefgrootte),
  waarde = unlist(steekproeven)
)

# Bereken gemiddelden en 95%-betrouwbaarheidsintervallen met t.test()
bhi_data <- data.frame(
  steekproef_id = 1:aantal_steekproeven,
  gemiddelde = sapply(steekproeven, mean),
  ondergrens = sapply(steekproeven, function(x) t.test(x)$conf.int[1]),
  bovengrens = sapply(steekproeven, function(x) t.test(x)$conf.int[2])
)

# Plot de data
ggplot() +
  # Jitterplot voor de individuele datapunten
  geom_jitter(
    data = datapunten,
    aes(x = factor(steekproef_id), y = waarde),
    color = "gray",
    size = 1,
    width = 0.2,
    alpha = 0.5
  ) +
  # Foutbalken voor de betrouwbaarheidsintervallen
  geom_errorbar(
    data = bhi_data,
    aes(x = factor(steekproef_id), ymin = ondergrens, ymax = bovengrens),
    width = 0.2
  ) +
  # Punten voor de steekproefgemiddelden
  geom_point(
    data = bhi_data,
    aes(x = factor(steekproef_id), y = gemiddelde),
    color = "darkorchid",
    size = 3
  ) +
  # Horizontale lijn voor het ware populatiegemiddelde
  geom_hline(
    yintercept = waar_gemiddelde,
    color = "#e8601d",
    linetype = "solid",
    size = 1
  ) +
  labs(
    x = "Steekproeven",
    y = "Chlorofylconcentratie (mg/g)",
    title = NULL
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 0, hjust = 0.5))

Figuur 24.1: Illustratie van het concept van een betrouwbaarheidsinterval. Getoond worden de 95%-betrouwbaarheidsintervallen van twintig onderzoekers die hetzelfde onderzoek hebben uitgevoerd om de gemiddelde chlorofylconcentratie te bepalen in de bladeren van tarweplanten. De ruwe datapunten zijn in lichtgrijs weergegeven, en de oranje lijn geeft het ware gemiddelde aan – de populatieparameter.

In dit hoofdstuk zullen we zien hoe je een betrouwbaarheidsinterval kunt berekenen voor het populatiegemiddelde van een continue variabele. Je zult uiteindelijk zien dat de formule daarvoor helemaal niet moeilijk is. Maar, we willen graag dat je ook begrijpt waarom de formule werkt, en onder welke aannames. Om dat voor elkaar te krijgen zullen we wat dieper de theorie in moeten duiken.

Oefening 24.1 (Bloeddrukverlagende medicijnen)

Een farmaceutisch bedrijf doet onderzoek naar het effect van twee nieuwe medicijnen op de systolische bloeddruk.

Hieronder zijn voor beide medicijnen het gemiddelde effect geplot (oranje datapunten). De oranje lijnen geven de 95%-betrouwbaarheidsintervallen aan.

Geven deze gegevens de indruk dat medicijn 1 gemiddeld genomen een effect heeft op de bloeddruk?
Kun je uitsluiten dat medicijn 1 gemiddeld genomen een effect heeft op de bloeddruk?
Geven deze gegevens een indicatie dat medicijn 2 een effect heeft op de bloeddruk?
Kan het zo zijn dat medicijn 2 gemiddeld genomen géén effect heeft op de bloeddruk?

Oefening 24.2 (Waarom 95%?)

Waarom berekenen we meestal een 95%-betrouwbaarheidsinterval, en niet, zeg, een 50%-betrouwbaarheidsinterval?

24.2 Het interpreteren van onzekere gegevens vereist een statistisch model

Wanneer we conclusies proberen te trekken uit onzekere gegevens maken we, bewust of onbewust, vaak gebruik van bepaalde kennis of veronderstellingen. Dat kan ook haast niet anders, want de conclusies die je kunt trekken hangen af van de manier waarop de gegevens tot stand zijn gekomen. Bijvoorbeeld, het is onmogelijk om conclusies te trekken uit de chlorofylmetingen hierboven zonder iets te veronderstellen over de nauwkeurigheid van de meetmethode, of over de eigenschappen van de steekproef.

Daarom is iedere statistische analyse gebaseerd op een model. In het onderdeel Biologische Modellen heb je geleerd dat een model een vereenvoudigde beschrijving is van een onderdeel van de werkelijkheid, en dat een model gebaseerd is op aannames. Die aannames vloeien voort uit kennis die we hebben, veronderstellingen die we bereid zijn te maken, of hypotheses waarvan we de consequenties willen onderzoeken. In de statistiek is dat niet anders, al vragen statistische problemen typisch om een ander type modellen dan je bij het onderdee Biologische Modellen gezien hebt. Bij Biologische Modellen gingen de aannames vaak over de manier waarop variabelen elkaars dynamiek (verandering in de tijd) beïnvloeden. In de statistiek gaan de aannames met name over de kansenprocessen die ten grondslag liggen aan onze onzekere gegevens, zoals de manier waarop de steekproef is uitgevoerd of de aard van de meetfouten. Men heeft het daarom in deze context vaak over statistische modellen. Een statistisch model is dus een kansmodel dat gebruikt wordt om onzekere gegevens te interpreteren.

Dit klinkt allemaal erg abstract, maar we gaan het nu gelijk concreter maken.

24.3 Het Normale Model voor eenvoudige aselecte steekproeven uit grote populaties

In Hoofdstuk 21 hebben we het uitgebreid gehad over steekproeven en steekproevenvariabiliteit. We hebben daar gezien dat onderzoeken binnen en buiten de biologie meestal gaan over een bepaalde populatie die bestaat uit zekere eenheden. Onze interesse gaat dan uit naar een bepaalde eigenschap $X$ van die eenheden. Specifiek willen we vaak het gemiddelde van variabele $X$ weten in de volledige populatie. We noemen die populatieparameter $\mu_X$.

Als we in staat zouden zijn om de waarde van $X$ te meten voor alle eenheden in de populatie, dan zouden we $\mu_X$ precies kunnen berekenen en zouden steekproefvariabiliteit en bias geen rol spelen. In de praktijk is dat zelden mogelijk en kunnen we enkel een steekproef van $n$ eenheden bestuderen. Schatters op basis van zo’n steekproef zijn onzeker als gevolg van steekproevenvariabiliteit. Hoe kunnen we die onzekerheid kwantificeren?

Dit is waar onze kennis van kansrekenen uit Hoofdstuk 22 en Hoofdstuk 23 van pas komt.

Het idee is dat we het nemen van een steekproef gaan beschouwen als een kansexperiment, en iedere waarneming binnen de steekproef als een kansvariabele.

Om daadwerkelijk kansen te kunnen berekenen op uitkomsten van steekproeven hebben we een statistisch model nodig dat de eigenschappen van de kansvariabelen vastlegt. In de statistiek wordt vaak gebruik gemaakt van het volgende model, dat we in dit boek het Normale Model zullen noemen.

Definitie: Het Normale Model van steekproeven.

Het Normale Model is gebaseerd op de volgende aannames:

De variabele is normaal verdeeld in de populatie. We zagen eerder dat veel variabelen een klokvormige verdeling hebben. Dit model gaat ervan uit dat de verdeling benaderd kan worden door een normale verdeling.
Iedere waarde $X_i$ van de steekproef is onafhankelijk uit deze normale verdeling getrokken. Dit is een redelijke benadering als we in werkelijkheid een eenvoudige aselecte steekproef uitvoeren en bovendien de populatie veel groter is dan de steekproef, zodat de populatie niet wezenlijk verandert als we al een deel van onze steekproef hebben gekozen.

Hoewel de twee aannames van het Normale Model in de meeste toepassingen niet exact waar zijn, zijn ze vaak een goede benadering. Maar er zijn ook situaties waarvoor dit model niet geschikt is. Methoden die op het Normale Model gebaseerd zijn geven in die situaties geen betrouwbare resultaten. We komen hier in Paragraaf 24.9.5 op terug.

24.4 De kansverdeling van het steekproefgemiddelde

Goed — laten we de eigenschappen van steekproeven eens verder bestuderen, aangenomen dat het Normale Model van toepassing is. We willen dus een schatting maken van het gemiddelde (de verwachtingswaarde) van $X$ in de populatie, $\mu_X$. Daarom hebben we een steekproef van $n$ eenheden genomen. We berekenen het gemiddelde van de steekproef: \[ \overline{X} = \frac{\sum_{i=1}^n X_i}{n}. \tag{24.1}\] Het steekproefgemiddelde gebruiken we als schatter voor het populatiegemiddelde.

Het is belangrijk te beseffen dat het gemiddelde van de steekproef $\overline{X}$ zelf óók een kansvariabele is. Want als je meerdere keren op dezelfde manier een steekproef zou uitvoeren en steeds het steekproefgemiddelde zou uitrekenen, dan zou je steeds een ander gemiddelde vinden. Kunnen we iets zeggen over hoe groot de onzekerheid in de schatter is?

Om daar meer inzicht in te krijgen gaan we het nemen van steekproeven eerst eens simuleren.

Simulaties van het Normale Model

Stel dat het populatiegemiddelde gelijk is aan $\mu_X = 100$, met standaarddeviatie $\sigma_{X} = 10$. Dan ziet de verdeling van de populatie er dus zo uit:

Code

# Laad ggplot2 indien nodig
if (!("ggplot2" %in% loadedNamespaces())) { 
  library(ggplot2) 
}

# Parameters voor de normale verdeling
mu <- 100                           # Gemiddelde
sigma <- 10                         # Standaarddeviatie
aantal_sd <- 3                      # Aantal standaarddeviaties rond het gemiddelde
xmin <- mu - aantal_sd * sigma
xmax <- mu + aantal_sd * sigma
x_waarden <- seq(xmin, xmax, by = 1)  # Domein voor X
kansdichtheden <- dnorm(
  x_waarden, 
  mean = mu, 
  sd = sigma
)  # Kansdichtheden berekenen

# Maak een dataframe
data <- data.frame(
  x = x_waarden,
  kansdichtheid = kansdichtheden
)

# Plot de normale verdeling
populatie <- ggplot(data, aes(x = x, y = kansdichtheid)) +
  geom_line(color = "#e8601d", size = 1) +  # Lijn voor de curve
  labs(
    x = expression(italic(X)),                # Label voor de x-as
    y = "Kansdichtheid",                      # Label voor de y-as
    title = NULL                              # Geen titel
  ) +
  theme_minimal()                             # Minimalistisch thema

populatie

Figuur 24.2: Een normaalverdeelde populatie met gemiddelde 100 en standaarddeviatie 10.

In de praktijk weten we de waardes van $\mu_X$ en $\sigma_{X}$ niet en hebben we alleen toegang tot een steekproef uit de populatie. Daarom trekken we als voorbeeld $n = 10$ eenheden uit de populatie. In Figuur 24.3 (a) zijn de resultaten in een histogram weergegeven. In de figuur is ook het gemiddelde $\overline{X}$ van de steekproef aangegeven met een verticale lijn.

Code

# Controleer of ggplot2 geladen is, en laad indien nodig
if (!("ggplot2" %in% loadedNamespaces())) {
  library(ggplot2)
}

set.seed(130)  # Voor reproduceerbaarheid

# Stel parameters in
steekproef_grootte <- 10
staafbreedte <- 2

# Genereer drie steekproeven van grootte 10
steekproeven <- lapply(
  1:3, 
  function(i) rnorm(steekproef_grootte, mean = mu, sd = sigma)
)

# Bereken steekproefgemiddelden
gemiddelden <- lapply(steekproeven, mean)

# Functie om een histogram te maken
maak_histogram <- function(data, gemiddelde, y_label = NULL) {
  intervallen <- cut(
    data, 
    breaks = seq(xmin, xmax, by = staafbreedte), 
    include.lowest = TRUE
  )
  hoogste <- max(table(intervallen))
  
  ggplot(data.frame(x = data), aes(x = x)) +
    geom_histogram(
      binwidth = staafbreedte,
      boundary = xmin,
      fill = "#e8601d",
      alpha = 0.7,
      color = "black",
      size = 0.5
    ) +
    geom_vline(
      xintercept = gemiddelde,
      color = "DarkOrchid",
      linetype = "dashed",
      size = 1
    ) +
    xlim(xmin, xmax) +
    labs(
      x = expression(italic(X)),
      y = y_label
    ) +
    geom_curve(
      aes(
        x = gemiddelde + 1.5 * sigma,
        y = 0.85 * hoogste,
        xend = gemiddelde + sigma / 20,
        yend = 0.7 * hoogste
      ),
      curvature = -0.4,
      arrow = arrow(length = unit(0.2, "cm")),
      color = "DarkOrchid",
      size = 0.7
    ) +
    annotate(
      "text",
      x = gemiddelde + 1.5 * sigma,
      y = 0.85 * hoogste + 0.06 * hoogste,
      label = expression(bar(italic(X))),
      color = "DarkOrchid",
      hjust = 0.5,
      size = 4
    ) +
    theme_minimal() +
    theme(plot.title = element_blank())
}

# Maak histogrammen
histogrammen <- lapply(
  1:3, 
  function(i) maak_histogram(
    steekproeven[[i]], 
    gemiddelden[[i]], 
    if (i == 1) "Frequentie" else NULL
  )
)

# Toon het eerste histogram
histogrammen[[1]]
histogrammen[[2]]
histogrammen[[3]]

Het steekproefgemiddelde dat we vinden hangt natuurlijk af van de toevallige steekproef die we getrokken hebben. Als we opnieuw een steekproef uitvoeren uit dezelfde populatie vinden we net andere waarden en dus ook een ander gemiddelde. Figuur 24.3 (b) en Figuur 24.3 (c) laten twee andere steekproeven zien, ieder met hun eigen gemiddelde $\overline{X}$.

Omdat iedere steekproef een andere waarde van $\overline{X}$ geeft en dus ook een andere schatting voor $\mu_X$, rijst de vraag: Hoe ver ligt $\overline{X}$ typisch van $\mu_X$ af? Als $\overline{X}$ bij een steekproef van $n = 10$ typisch heel ver van $\mu_X$ afligt, dan is een schatting op basis van zo’n steekproef blijkbaar erg onzeker. Als $\overline{X}$ typisch juist heel dicht bij $\mu_X$ ligt, dan is die schatting blijkbaar juist heel precies en nauwkeurig.

Om daar een beeld van te krijgen trekken we nu heel vaak opnieuw een steekproef, steeds met steekproefgrootte $n = 10$, en rekenen elke keer opnieuw het gemiddelde $\overline{X}$ uit. (In werkelijkheid zouden we dat natuurlijk nooit zo doen, maar dit is een gedachtenexperiment.) Vervolgens tekenen we een histogram van alle steekproefgemiddelden die we gevonden hebben. We noemen dit de steekproevenverdeling van het gemiddelde (sampling distribution of the mean). We simuleren dit proces weer met de computer; het resultaat staat in Figuur 24.4 (b). Ter vergelijking is in Figuur 24.4 (a) weer de populatie weergegeven.

Code

# Laad ggplot2 indien nodig
if (!("ggplot2" %in% loadedNamespaces())) { 
  library(ggplot2) 
}

# Parameters instellen
aantal_herhalingen <- 10^4       # Aantal herhalingen
steekproefgroottes <- c(10, 40, 100)  # Steekproefgroottes
staafbreedte <- 1                # Binbreedte voor histogrammen

# Simuleer steekproefgemiddelden
set.seed(42)                     # Voor reproduceerbaarheid
steekproef_gemiddelden <- matrix(nrow = aantal_herhalingen, ncol = length(steekproefgroottes))
for (i in seq_along(steekproefgroottes)) {
  for (j in seq_len(aantal_herhalingen)) {
    steekproef <- rnorm(steekproefgroottes[i], mean = mu, sd = sigma)
    steekproef_gemiddelden[j, i] <- mean(steekproef)
  }
}

# Zet de steekproefgemiddelden om naar een long-formaat data frame
steekproef_gemiddelden <- data.frame(
  waarde = c(
    steekproef_gemiddelden[, 1], 
    steekproef_gemiddelden[, 2], 
    steekproef_gemiddelden[, 3]
  ),
  groep = factor(
    rep(c("n = 10", "n = 40", "n = 100"), each = aantal_herhalingen),
    levels = c("n = 10", "n = 40", "n = 100")
  )
)

# Functie om een histogram te maken
maak_verdelingsplot <- function(data, label, y_as_label) {
  plot <- ggplot(data, aes(x = waarde)) +
    geom_histogram(
      binwidth = staafbreedte, fill = "#e8601d", 
      color = "black", alpha = 0.7, size = 0.5
    ) +
    xlim(xmin, xmax) +
    labs(
      x = expression(bar(italic(X))),  # LaTeX-stijl voor X
      title = NULL                    # Titel is de groepslabel
    ) +
    theme_minimal() +
    theme(
      plot.title = element_text(hjust = 0.5, face = "bold"),
      axis.text.y = if (!y_as_label) element_blank() else element_text(),
      axis.ticks.y = if (!y_as_label) element_blank() else element_line()
    )
  
  # Voeg y-label toe indien nodig
  if (y_as_label) {
    plot <- plot + labs(y = "Frequentie")
  } else {
    plot <- plot + labs(y = NULL)
  }
  
  return(plot)
}

# Subsets van de data voor elke groep
data_n10 <- subset(steekproef_gemiddelden, groep == "n = 10")
data_n40 <- subset(steekproef_gemiddelden, groep == "n = 40")
data_n100 <- subset(steekproef_gemiddelden, groep == "n = 100")

# Maak afzonderlijke histogrammen
plot_n10 <- maak_verdelingsplot(data_n10, label = "n = 10", y_as_label = TRUE)
plot_n40 <- maak_verdelingsplot(data_n40, label = "n = 40", y_as_label = TRUE)
plot_n100 <- maak_verdelingsplot(data_n100, label = "n = 100", y_as_label = TRUE)

populatie
plot_n10
plot_n40
plot_n100

(a) Verdeling van de variabele in de populatie.

Het eerste wat opvalt aan de steekproevenverdeling in Figuur 24.4 (b) is dat het gemiddelde van de verdeling op $\overline{X} = 100$ ligt en dus gelijk is aan de populatieparameter $\mu_X$. Dat wil zeggen dat $\overline{X}$ een zuivere schatter is: gemiddeld is die gelijk aan het populatiegemiddelde. (Anders gezegd: de verwachtingswaarde van $\overline{X}$ is $\mu_X$.) Dat komt doordat het Normale Model uitgaat van een steekproef zonder bias. De schatter $\overline{X}$ heeft dus geen systematische fout, maar wel een toevallige fout.

Misschien viel je ook op dat de steekproevenverdeling van het gemiddelde (Figuur 24.4 (b)) een mooie klokvorm heeft die lijkt op een normale verdeling. Dat klopt: als de populatie normaal verdeeld is, is de steekproevenverdeling van het gemiddelde dat ook.

Vergelijk tot slot de steekproevenverdeling van het gemiddelde (Figuur 24.4 (b)) eens met de verdeling van de variabele in de populatie (Figuur 24.4 (a)). Het is duidelijk te zien dat de steekproevenverdeling een veel kleinere spreiding heeft. (De klokvorm is bij de steekproevenverdeling smaller.) Als de verschillende gemiddelden in de steekproevenverdeling sterk van elkaar verschillen, dan wil dat zeggen dat $\overline{X}$ een onzekere (weinig precieze) schatter is. Daarom is de spreiding van de steekproevenverdeling een maat van de onzekerheid.

Intuïtief zou je verwachten dat een grote steekproef een preciezere schatting zou moeten opleveren dan een kleine steekproef. Laten we daarom ook de steekproevenverdeling simuleren voor een grotere steekproef van $n = 40$. We nemen dus weer héél veel steekproeven, maar nu van $n = 40$ eenheden. Vervolgens berekenen we weer het gemiddelde van iedere steekproef, en geven al die steekproefgemiddelden weer in een histogram. Het resultaat staat in Figuur 24.4 (c). Zoals je ziet is de piek van de steekproevenverdeling bij $n = 40$ smaller dan bij $n = 10$. En de verdeling wordt nog smaller als $n = 100$; zie Figuur 24.4 (d). De conclusie is: hoe groter de steekproeven zijn, hoe dichter de gemiddelden van de steekproeven bij elkaar liggen, en hoe kleiner dus de toevallige fout is.

Conclusies op basis van de simulaties:

De spreiding van de steekproevenverdeling van het gemiddelde is een maat voor de toevallige fout als gevolg van steekproevenvariabiliteit.
Bij grotere steekproeven heeft de steekproevenverdeling een kleinere spreiding.
Grotere steekproeven geven dus preciezere schattingen.

24.5 De standaarddeviatie van het steekproefgemiddelde.

We zagen net dat de spreiding van de steekproevenverdeling van het gemiddelde aangeeft hoe ver een typische schatting van het populatiegemiddelde ligt. Het is daarom nuttig om die spreiding te kwantificeren. Een goede maat van spreiding is de standaarddeviatie. De standaarddeviatie van de steekproevenverdeling noteren we als $\sigma_{\overline{X}}$.

Met wat wiskunde (die we je hier zullen besparen) kan de volgende formule voor $\sigma_{\overline{X}}$ worden afgeleid: \[ \sigma_{\overline{X}} = \frac{\sigma_{X}}{\sqrt{n}}. \tag{24.2}\] Dit is een belangrijk resultaat! Het bevestigt wat we in Figuur 24.4 al zagen: de standaarddeviatie van de steekproevenverdeling neemt af naarmate $n$ groter wordt. Maar nu weten we precies hoe sterk die afname is: de afname komt door de $\sqrt{n}$ in de noemer van Vergelijking 24.2.

Dat wortelteken heeft grote gevolgen. Bijvoorbeeld, als we de toevallige fout twee keer zo klein willen maken, hebben we een vier keer zo grote steekproef nodig. Willen we de toevallige fout 10 keer zo klein hebben, dan moeten we zelfs een 100 keer zo grote steekproef nemen. Dat loopt snel uit hand.

Wat hebben we met al dit rekenwerk bereikt? We vatten de voorlopige conclusies samen:

Conclusies over de steekproevenverdeling van het gemiddelde

We gaan uit van het Normale Model: we trekken onafhankelijk $n$ eenheden uit een normaalverdeelde populatie met parameters $\mu_X$ en $\sigma_{X}$.

Het gemiddelde van $X$ in de steekproef, dus $\overline{X}$, is een kansvariabele.
Deze kansvariabele $\overline{X}$ is normaal verdeeld.
De verwachtingswaarde van $\overline{X}$ is het populatiegemiddelde $\mu_X$. Dat betekent dat $\overline{X}$ een zuivere schatter is voor het populatiegemiddelde.
De standaarddeviatie van $\overline{X}$ is $\sigma_{\overline{X}} = \sigma/\sqrt{n}$.

In Paragraaf 23.7.1 (vuistregels) hebben we gezien dat bij een normale verdeling 95% van de waarnemingen binnen 1,96 standaarddeviaties van het gemiddelde valt. Daaruit kunnen we concluderen dat een gemiddelde van een steekproef met 95% zekerheid binnen $1,96 \sigma_{\overline{X}}$ van de ware populatieparameter $\mu_X$ afligt.

We kunnen die redenering omdraaien. Als we in een steekproef een steekproefgemiddelde $\overline{X}$ vinden, dan betekent dit dat met 95% zekerheid $\mu_X$ binnen $1{,}96 \sigma_{\overline{X}}$ van die waarde afligt.

Dit is het idee achter het BHI van het gemiddelde. In de volgende paragraaf gaan we dat idee precies maken.

Oefening 24.3 (Steekproevenverdeling van de gemiddelde spanwijdte van buizerds)

Stel dat de gemiddelde spanwijdte van de populatie buizerds in Nederland $\mu_X = 110\textrm{\,cm}$ is, met standaarddeviatie $\sigma_{X} = 9\textrm{\,cm}$. Ga uit van het Normale Model.

Schets de verdeling van de populatie. Label het gemiddelde en de standaarddeviatie in de figuur.
Je vangt aselect 1 buizerd en meet de spanwijdte. Wat is de kans dat deze buizerd een spanwijdte heeft die meer dan 18cm afwijkt van $\mu_X$?

Arceer de kans in je schets uit onderdeel a.

Je bent van plan om eenvoudig aselect 9 vogels te vangen, hun spanwijdte te meten, en het gemiddelde van die steekproef te meten.

Wat is de verwachtingswaarde van dat steekproefgemiddelde?
Wat is de standaarddeviatie van dat steekproefgemiddelde?
Schets de kansverdeling van het steekproefgemiddelde (de steekproevenverdeling van het gemiddelde). Label het gemiddelde en de standaarddeviatie van deze verdeling.
Wat is de kans dat het steekproefgemiddelde meer dan $6\textrm{\,cm}$ afwijkt van het populatiegemiddelde $\mu_X$?

24.6 95%-BHI van het gemiddelde als de standaarddeviatie van de populatie bekend is

We hebben nu alle theorie op een rijtje om af te gaan leiden hoe je het BHI van een gemiddelde kunt berekenen.

We gaan nog steeds uit van het Normale Model. De waarde van het populatiegemiddelde $\mu_X$ is onbekend. Om die waarde te schatten hebben we een steekproef van $n$ eenheden genomen.

In deze paragraaf behandelen we een situatie die niet vaak voorkomt, maar waarbij de berekening van het BHI gemakkelijk te begrijpen is. We gaan er namelijk vanuit dat we de standaarddeviatie van de populatie, $\sigma_{X}$, precies weten. Het is moeilijk een situatie te bedenken waarbij we $\mu_X$ niet weten, maar $\sigma_{X}$ wel, dus erg realistisch is deze situatie niet. Maar dit is een opstapje naar de iets ingewikkeldere situatie waarbij we $\sigma_{X}$ niet kennen; die situatie behandelen we in de volgende paragraaf.

We hebben net gezien dat het steekproefgemiddelde $\overline{X}$ zelf een kansvariabele is. Deze kansvariabele is normaal verdeeld met verwachtingswaarde $\mu_X$ en standaarddeviatie $\sigma_{\overline{X}}$. We hebben geleerd dat we elke normaalverdeelde variabele kunnen transformeren naar een standaardnormaal verdeelde variabele $Z$ door het gemiddelde eraf te trekken en te delen door de standaarddeviatie. Oftewel, \[ Z = \frac{\overline{X} - \mu_X}{\sigma_{\overline{X}}} \tag{24.3}\] is standaardnormaal verdeeld.

Van de standaardnormale verdeling weten we dat tussen -1,96 en 1,96 een kans van 0,95 ligt: \[ \textrm{Pr}\!\left[-1{,}96 < Z < 1{,}96\right] = 0{,}95. \] We vullen Vergelijking 24.3 in en krijgen: \[ \textrm{Pr}\!\left[-1{,}96 < \frac{\overline{X} - \mu_X}{\sigma_{\overline{X}}} < 1{,}96\right] = 0{,}95. \] We gaan de ongelijkheid verder uitwerken. Eerst vermenigvuldigen we alle kanten van de ongelijkheid met $-\sigma_{\overline{X}}$: \[ \textrm{Pr}\!\left[1{,}96 \sigma_{\overline{X}} > \mu_X -\overline{X} > -1{,}96 \sigma_{\overline{X}}\right] = 0{,}95. \] Omdat we met een negatief getal vermenigvuldigd hebben, zijn de tekens van de ongelijkheden omgeklapt. Het is overzichtelijker om de volgorde van de ongelijkheden om te draaien: \[ \textrm{Pr}\!\left[-1{,}96 \sigma_{\overline{X}} < \mu_X -\overline{X} < 1{,}96 \sigma_{\overline{X}}\right] = 0{,}95. \] Ten slotte tellen we $\overline{X}$ bij alle kanten van de ongelijkheid op: \[ \textrm{Pr}\!\left[\overline{X} -1{,}96\, \sigma_{\overline{X}} < \mu_X < \overline{X} + 1{,}96\, \sigma_{\overline{X}}\right] = 0{,}95. \]

Dit is precies de vergelijking waarnaar we op zoek waren! De vergelijking zegt: elke keer dat we een steekproef nemen en $\overline{X}$ uitrekenen, hebben we een kans van 0,95 dat het populatiegemiddelde $\mu_X$ tussen de grenzen $\overline{X} -1{,}96\, \sigma_{\overline{X}}$ en $\overline{X} + 1{,}96\, \sigma_{\overline{X}}$ in valt. Dat betekent dat het volgende interval een 95%-BHI is: \[\left[\overline{X} -1{,}96\, \sigma_{\overline{X}}, \overline{X} + 1{,}96\, \sigma_{\overline{X}} \right].\] Het interval wordt ook vaak genoteerd als \[ \mu_X = \overline{X} \pm 1{,}96\, \sigma_{\overline{X}}. \] Laten we dit resultaat weer even samenvatten.

95%-BHI voor de situatie dat $\sigma_{X}$ bekend is

We gaan uit van het Normale Model, waarin variabele $X$ normaal verdeeld is in de populatie en de steekproef bestaat uit $n$ eenheden die onafhankelijk uit die normale verdeling zijn getrokken.

We weten het populatiegemiddelde $\mu_X$ niet, maar de standaarddeviatie $\sigma_{X}$ wel.

Dan wordt het 95%-BHI van het populatiegemiddelde $\mu_X$ gegeven door:

\[ \mu_X = \overline{X} \pm 1{,}96\, \sigma_{\overline{X}}, \tag{24.4}\] waarin \[ \sigma_{\overline{X}} = \frac{\sigma_{X}}{\sqrt{n}}. \]

Oefening 24.4 (BHI van de spanwijdte van buizerds)

Neem aan, net als in Oefening 24.3, dat de spanwijdte van buizerds in de populatie normaal verdeeld is met standaarddeviatie $\sigma_{X} = 9\textrm{\,cm}$. Maar in deze opgave weten we het populatiegemiddelde $\mu_X$ niet.

Je vangt 9 buizerds (eenvoudig aselect) en meet hun spanwijdte. De meetwaarden zijn, in cm:

\[117,\quad 106,\quad 123,\quad 111,\quad 124,\quad 120,\quad 108,\quad 115,\quad 119.\]

Bereken $\overline{X}$.
Bereken het betrouwbaarheidsinterval voor $\mu_X$.
Stel dat de standaarddeviatie van de populatie niet $9\textrm{\,cm}$ maar $18\textrm{\,cm}$ is. (De spreiding is dus 2 keer zo groot.)

Hoe verandert dan het betrouwbaarheidsinterval?
Ga weer uit van een standaarddeviatie van $9\textrm{\,cm}$. Stel dat je niet 9 maar 18 vogels gevangen had, maar wel hetzelfde steekproefgemiddelde had gevonden als in onderdeel a.
Hoe verandert dan het betrouwbaarheidsinterval?
Reflecteer op je antwoorden bij c en d. Welke factoren bepalen de onzekerheid in onze schatting van $\mu_X$?

24.7 BHI als de variantie van de populatie niet bekend is

In de vorige paragraaf hebben we een 95%-BHI afgeleid voor de situatie waarin de populatievariantie $\sigma_{X}$ bekend is. De afleiding voor de gebruikelijke situatie, waarin $\sigma_{X}$ niet bij voorbaat gegeven is, volgt dezelfde route — maar met een twist.

In feite gaat de afleiding van de vorige paragraaf ook op voor deze situatie, tot en met het eindresultaat van Vergelijking 24.4. Het probleem is alleen dat we nu niets aan deze formule hebben. Om de grenzen van dat interval uit te rekenen, zullen we namelijk een waarde moeten invullen voor $\sigma_{\overline{X}} = \sigma_{X}/\sqrt{n}$, en dat gaat niet als we $\sigma_{X}$ niet weten.

Om verder te kunnen komen zullen we $\sigma_{\overline{X}}$ op een of andere manier moeten schatten; daarover gaat de volgende paragraaf.

De standaardfout van het steekproefgemiddelde

Hoewel we niet bij voorbaat weten wat de waarde is van $\sigma_{X}$, de standaarddeviatie van de populatie, geeft de steekproef hier wel informatie over. De standaarddeviatie van de steekproef, $s_X$, kan namelijk gebruikt worden als een schatter voor de standaarddeviatie van de populatie $\sigma_{X}$. Dat zal je niet verrassen: als de spreiding binnen de steekproef groot is, valt te verwachten dat die van de populatie óók groot is. Uiteraard komt deze schatter met zijn eigen toevallige (en ook systematische) fouten, vooral als de steekproef klein is.

Omdat we $s_X$ kunnen gebruiken als schatter voor $\sigma_{X}$ kunnen we ook $\sigma_{\overline{X}} = \sigma_{X}/\sqrt{n}$ schatten door voor $\sigma_{X}$ de waarde van $s_X$ in te vullen. Het resultaat wordt de standaardfout van het gemiddelde genoemd: \[ \mathrm{SE}_{\overline{X}} = \frac{s_X}{\sqrt{n}}. \tag{24.5}\] Hier gebruiken we de notatie $\mathrm{SE}_{\overline{X}}$; in de literatuur kom je ook de afkorting SEM tegen, wat staat voor Standard Error of the Mean.

Definitie: Standaardfout van het gemiddelde.

\[ \mathrm{SE}_{\overline{X}} = \frac{s_X}{\sqrt{n}}. \tag{24.6}\]

Transformeren naar $t$

In de situatie dat $\sigma_{X}$ bekend was, konden we de kansvariabele $\overline{X}$ transformeren naar een standaard-normaal verdeelde variabele door middel van Vergelijking 24.3: \[ Z = \frac{\overline{X} - \mu_X}{\sigma_{\overline{X}}}. \] Om uiteindelijk een formule voor een BHI te krijgen waarin in plaats van $\sigma_{\overline{X}}$ de standaardfout $\mathrm{SE}_{\overline{X}}$ voorkomt, moeten we in de noemer $\sigma_{\overline{X}}$ vervangen door zijn schatter $\mathrm{SE}_{\overline{X}}$. Helaas heeft dat ook tot gevolg dat de getransformeerde variabele nu niet meer standaardnormaal verdeeld is. We noemen het resultaat van de transformatie daarom niet meer $Z$, maar $t$: \[ t = \frac{\overline{X} - \mu_X}{\mathrm{SE}_{\overline{X}}}. \tag{24.7}\] De keuze voor $t$ als de naam van deze getransformeerde variabele is niet willekeurig. Er kan namelijk worden aangetoond dat $t$ onder het Normale Model een $t$-verdeling volgt met $n - 1$ vrijheidsgraden. Dat wil zeggen: als je

een steekproef neemt van $n$ eenheden,
het gemiddelde $\overline{X}$ en de standaarddeviatie $s$ van die steekproef berekent
en daarmee Vergelijking 24.7 invult om $t$ te berekenen

dan is $t$ een kansvariabele met als kansverdeling de $t$-verdeling met $n-1$ vrijheidsgraden.

Helaas voert het bewijs voor deze stelling te ver voor deze cursus; je zult daarop moeten vertrouwen. Maar in de spirit van trust, but verify kun je de bewering wel controleren met behulp van simulaties.

De $t$-verdeling controleren met simulaties

Het is niet moeilijk om te verifiëren dat de $t$ uit Vergelijking 24.7 onder het Normale Model een $t$-verdeling volgt met $n - 1$ vrijheidsgraden.

De code hieronder (klik op “Code” om die zichtbaar te maken) neemt $10^4$ keer een (kleine) steekproef ($n = 4$) uit een normale verdeling, en berekent op basis van die steekproeven steeds $\overline{X}$, $s$, $\mathrm{SE}_{\overline{X}}$, en uiteindelijk $t$. Het histogram van de verkregen $t$ waarden zou dan de $t$-verdeling met $\mathrm{df}= n -1 = 3$ vrijheidsgraden moeten volgen. Om dat te verifiëren plot de code bovenop het histogram ook die $t$-verdeling (in paars).

Code

# Laad ggplot2 indien nodig
if (!("ggplot2" %in% loadedNamespaces())) { 
  library(ggplot2) 
}

# Parameters instellen
aantal_iteraties <- 10^4   # Aantal herhalingen
steekproefgrootte <- 4      # Grootte van de steekproef
mu <- 100                   # Populatiegemiddelde
sigma <- 10                 # Populatiestandaarddeviatie
vrijheidsgraden <- steekproefgrootte - 1  # Vrijheidsgraden voor t-verdeling
staafbreedte <- 0.2         # Binbreedte voor het histogram
domein <- seq(-4, 4, length.out = 400)  # Verfijnd domein voor gladde lijnen

# Genereer t-waarden
set.seed(42)  # Voor reproduceerbaarheid
t_waarden <- replicate(aantal_iteraties, {
  # Trek een steekproef
  steekproef <- rnorm(steekproefgrootte, mean = mu, sd = sigma)
  
  # Bereken steekproefgemiddelde en standaarddeviatie van de steekproef
  gemiddelde <- mean(steekproef)
  s <- sd(steekproef)
  
  # Bereken standaardfout
  SE <- s / sqrt(steekproefgrootte)
  
  # Bereken t-waarde
  t <- (gemiddelde - mu) / SE
  return(t)
})

# Zet t-waarden in een data frame
data <- data.frame(t_waarden)

# Bereid theoretische verdelingen voor
theoretische_data <- data.frame(
  x = domein,
  t_verdeling = dt(domein, df = vrijheidsgraden),
  standaard_normaal = dnorm(domein)
)

# Maak het histogram en voeg de theoretische verdelingen toe
ggplot(data, aes(x = t_waarden)) +
  geom_histogram(
    aes(y = ..density..), # schaal het histogram, totale oppervlakte 1 
    binwidth = staafbreedte, 
    fill = "#e8601d", 
    color = "black", 
    alpha = 0.5,
    show.legend = FALSE
  ) +
  geom_line(
    data = theoretische_data, 
    aes(x = x, y = t_verdeling, color = "t-verdeling"), 
    size = 1
  ) +
  geom_line(
    data = theoretische_data, 
    aes(x = x, y = standaard_normaal, color = "standaardnormaal"), 
    size = 1, linetype = "dashed"
  ) +
  scale_color_manual(
    values = c("t-verdeling" = "DarkOrchid", "standaardnormaal" = "black"),
    name = "Verdeling"
  ) +
  coord_cartesian(xlim = c(-4, 4)) +  # Beperk het domein van de plot
  labs(
    x = expression(italic(t)-waarden),
    y = "Kansdichtheid",
    title = ""
  ) +
  theme_minimal()

Figuur 24.5: Histogram van $t$-waarden vergeleken met de $t$-verdeling en standaardnormale verdeling.

Ter vergelijking is ook de standaardnormale verdeling toegevoegd (zwarte stippellijn). Zoals je geleerd hebt in Paragraaf 23.9 heeft de $t$-verdeling een lagere piek en dikkere staarten, en dus een grotere spreiding. In dit geval komt dat doordat in de noemer van Vergelijking 24.7 niet $\sigma_{\overline{X}}$ staat maar de schatter daarvan $\mathrm{SE}_{\overline{X}}$. Daardoor is nu niet alleen de teller maar ook de noemer onzeker. Het gevolg is dat $t$ (Vergelijking 24.7) onzekerder is dan dan $Z$ (Vergelijking 24.3).

We vatten de situatie weer even samen:

Belangrijk:

Onder het Normale Model is de kansverdeling van de kansvariabele \[ t = \frac{\overline{X} - \mu_X}{\mathrm{SE}_{\overline{X}}} \tag{24.8}\] is de $t$-verdeling met $n - 1$ vrijheidsgraden.

Afleiding van het 95%-BHI

Nu we dit alles weten, kunnen we de afleiding van het 95%-BHI eindelijk afmaken.

Omdat de $t$ uit Vergelijking 24.8 verdeeld is volgens de $t$-verdeling met $n-1$ vrijheidsgraden, geldt \[ \textrm{Pr}\!\left[ -t_{0{,}05 (2) \mathrm{(n-1)}} < t < t_{0{,}05 (2) \mathrm{(n-1)}} \right] = 0{,}95. \] We hebben hier de kritieke waarden gebruikt voor de $t$-verdeling (zie Vergelijking 23.15 en de paragraaf daar omheen).

Als we Vergelijking 24.8 hierin invullen, krijgen we \[ \textrm{Pr}\!\left[ -t_{0{,}05 (2) \mathrm{(n-1)}} < \frac{\overline{X} - \mu_X}{\mathrm{SE}_{\overline{X}}} < t_{0{,}05 (2) \mathrm{(n-1)}} \right] = 0{,}95. \] De afleiding is verder helemaal analoog aan die van Paragraaf 24.6. We vermenigvuldigen de ongelijkheden eerst met $-\mathrm{SE}_{\overline{X}}$: \[ \textrm{Pr}\!\left[t_{0{,}05 (2) \mathrm{(n-1)}}\, \mathrm{SE}_{\overline{X}} > \mu_X - \overline{X} > -t_{0{,}05 (2) \mathrm{(n-1)}}\,\mathrm{SE}_{\overline{X}}\right] = 0{,}95. \] Het is overzichtelijker om de ongelijkheden om te draaien: \[ \textrm{Pr}\!\left[-t_{0{,}05 (2) \mathrm{(n-1)}}\, \mathrm{SE}_{\overline{X}} < \mu_X - \overline{X} < t_{0{,}05 (2) \mathrm{(n-1)}}\,\mathrm{SE}_{\overline{X}}\right] = 0{,}95. \] Tot slot tellen we $\overline{X}$ bij de ongelijkheden op: \[ \textrm{Pr}\!\left[\overline{X} -t_{0{,}05 (2) \mathrm{(n-1)}}\, \mathrm{SE}_{\overline{X}} < \mu_X < \overline{X} + t_{0{,}05 (2) \mathrm{(n-1)}}\, \mathrm{SE}_{\overline{X}}\right] = 0{,}95. \] Dit is precies de vergelijking waarnaar we op zoek waren! De vergelijking zegt dat er iedere keer dat we een steekproef nemen een kans van 0,95 is dat de populatieparameter $\mu_X$ tussen de grenzen $\overline{X} -t_{0{,}05 (2) \mathrm{(n-1)}}$ en $\overline{X} + t_{0{,}05 (2) \mathrm{(n-1)}}$ invalt. Oftewel, het 95%-BHI is: \[ \left[\overline{X} -t_{0{,}05 (2) \mathrm{(n-1)}}\, \mathrm{SE}_{\overline{X}}, \overline{X} + t_{0{,}05 (2) \mathrm{(n-1)}}\, \mathrm{SE}_{\overline{X}} \right]. \] Dit wordt ook vaak geschreven als \[ \mu_X = \overline{X} \pm t_{0{,}05 (2) \mathrm{(n-1)}}\, \mathrm{SE}_{\overline{X}}. \]

Laten we dit eindresultaat maar weer eens samenvatten:

95%-BHI voor de situatie dat $\sigma_{X}$ onbekend is

We gaan uit van het Normale Model, waarin variabele $X$ normaal verdeeld is in de populatie en de steekproef bestaat uit $n$ eenheden die onafhankelijk uit die normale verdeling zijn getrokken.

Het populatiegemiddelde $\mu_X$ en de standaarddeviatie $\sigma_{X}$ van de populatie zijn ons onbekend.

Dan wordt het 95%-BHI van het populatiegemiddelde $\mu_X$ gegeven door:

\[ \mu_X = \overline{X} \pm t_{0{,}05 (2) \mathrm{(n-1)}}\, \mathrm{SE}_{\overline{X}}, \tag{24.9}\] waarin \[ \mathrm{SE}_{\overline{X}} = \frac{s_X}{\sqrt{n}}. \]

Vergelijk Vergelijking 24.9 eens met Vergelijking 24.4, de formule die we hadden afgeleid voor de (onwaarschijnlijke) situatie dat $\sigma_{X}$ vooraf bekend is. De formules lijken erg op elkaar, maar er zijn twee verschillen.

Ten eerste, waar in Vergelijking 24.4 $\sigma_{\overline{X}}$ staat, staat in Vergelijking 24.9 de standaardfout $\mathrm{SE}_{\overline{X}}$. Ten tweede wordt in Vergelijking 24.9 de kritieke waarde van de $t$-verdeling gebruikt in plaats van de kritieke waarde van de standaard-normale verdeling 1,96. Omdat de kritieke $t$-waarde altijd groter is dan 1,96 zorgt dit voor bredere betrouwbaarheidsintervallen. Die bredere betrouwbaarheidsintervallen reflecteren het feit dat we zonder kennis van $\sigma_{X}$ ook onzekerder zijn over de waarde van $\mu_X$.

Maar als het aantal waarnemingen niet al te klein is, is dat effect klein. In Figuur 23.13 zagen we al dat de kritieke $t$-waarde al bij 10 vrijheidsgraden dicht bij 2 ligt. Als grove schatting is $\mu_X = \overline{X} \pm 2 \mathrm{SE}_{\overline{X}}$ dus meestal een redelijke benadering van het 95%- BHI

Om concreet te maken hoe je een 95%-BHI kunt berekenen, werken we nu een eenvoudig voorbeeld uit.

Voorbeeld 24.1 (95%-BHI voor de opbrengst van aardappelplanten)

StSP6A is een transgene lijn van aardappelplanten. Onderzoekers willen weten wat het effect van de mutatie is op de aardappelopbrengst. Daarom planten ze 5 planten in een kas en laten die 100 dagen groeien onder gecontroleerde condities. Daarna worden de aardappels gerooid en wordt van iedere plant het totale gewicht aan aardappels bepaald.

De resultaten zijn:

Plant	Tubergewicht (g)
1	112
2	125
3	119
4	89
5	102

We berekenen het 95%-BHI van het gemiddelde tubergewicht. We volgen de volgende stappen:

De steekproefgrootte is $n = 5$.
Het steekproefgemiddelde (Vergelijking 20.1) van de aardappelen is $\overline{x} = 109{,}4\textrm{\,g}$.
De standaarddeviatie (Vergelijking 20.6) van de steekproef is $s_X = 14{,}26\textrm{\,g}$.
De standaardfout (Vergelijking 24.6) is $\mathrm{SE}_{\overline{X}} = s_X / \sqrt{n} = 6{,}38\textrm{\,g}$.
Het aantal vrijheidsgraden is $\mathrm{df}= n - 1 = 4$.
Voor de kritieke waarde $t_{0{,}05 (2) \mathrm{4}}$ raadplegen we R (Paragraaf 23.10.3):
```
alpha <- 0.05
qt(1 - alpha/2, df = 4)
```
```
[1] 2.776445
```
De ondergrens van het 95%-BHI is $\overline{X} - t_{0{,}05 (2) \mathrm{4}} \mathrm{SE}_{\overline{X}} = 91{,}7\textrm{\,g}$.
De bovengrens van het 95%-BHI is $\overline{X} + t_{0{,}05 (2) \mathrm{4}} \mathrm{SE}_{\overline{X}} = 127{,}1\textrm{\,g}$.

Het betrouwbaarheidsinterval is dus $[91{,}7\textrm{\,g}; 127{,}1\textrm{\,g}]$, oftewel $(109{,}4 \pm 17{,}7)\textrm{\,g}$.

Oefening 24.5 (BHI van de gemiddelde spanwijdte van buizerds (alweer???))

Neem aan, net als in Oefening 24.4, dat de spanwijdte van buizerds in de populatie normaal verdeeld is. We weten het populatiegemiddelde $\mu_X$ niet. Deze keer weten we ook de standaarddeviatie van de populatie $\sigma_X$ niet.

Je vangt 9 buizerds (eenvoudig aselect) en meet hun spanwijdte. De meetwaarden zijn hetzelfde als in Oefening 24.4.

Volg de stappen 1 t/m 8 van Voorbeeld 24.1 om het BHI van de gemiddelde spanwijdte te berekenen. Je mag hiervoor een scriptje in R gebruiken.

Een BHI berekenen met R

In Voorbeeld 24.1 bestond de dataset uit maar vijf waarnemingen. Daarom was het goed mogelijk om het BHI met de hand uit te rekenen. Maar zelfs dan moesten we van R gebruik maken om de kritieke waarde op te zoeken. In de praktijk is het dus veel gemakkelijker en betrouwbaarder om de hele berekening in R uit te voeren.

In Oefening 24.5 heb je de berekening zelf ingeprogrammeerd. Maar het kan gemakkelijker: er bestaat standaard in R een functie die de hele berekening voor je uitvoert.

Het volgende blokje code laat zien hoe gemakkelijk dat is:

# creëer een vector met de gewichten (in g)
gewichten <- c(112, 125, 119, 89, 102)

# voer de functie t.test uit en sla de resultaten op
uitvoer_t.test <- t.test(gewichten)

# selecteer uit de resultaten het 95%-BHI
uitvoer_t.test$conf.int

[1]  91.69594 127.10406
attr(,"conf.level")
[1] 0.95

Je maakt gebruik van de functie t.test() die we in de volgende hoofdstukken nog vaker gaan terugzien. Deze functie berekent allerlei statistieken, waaronder het BHI. We slaan eerst de uitvoer van t.test() op in een variabele. In de code hierboven heet die variabele uitvoer_t.test. Vervolgens vragen we uit deze uitvoer het betrouwbaarheidsinterval (conf.int) op. Als het goed is, komt de uitvoer precies overeen met de getallen die in Voorbeeld 24.1 berekend hadden. Klopt dat?

Oefening 24.6 (Weer 9 buizerds)

Controleer je antwoord op Oefening 24.5 met behulp van de functie t.test().

24.8 Het Normale Model is ook geschikt voor herhaalde experimenten met meetfouten

In dit hoofdstuk hebben we het uitgebreid gehad over schattingen van populatieparameters op basis van een steekproef. In ander onderzoek in de biologie maken we geen gebruik van steekproeven maar proberen we een bepaalde grootheid te meten. In Hoofdstuk 21 hebben we besproken dat zulke metingen vaak worden geplaagd door meetfouten. Als de meetmethode goed gekalibreerd is, is de toevallige fout het grootste probleem.

Om de toevallige fout te verminderen wordt meestal dezelfde meting meerdere keren gerepliceerd (herhaald) en de uitkomsten gemiddeld. Als we een heel erg groot aantal replicaties zouden kunnen uitvoeren zou de toevallige fout uitmiddelen en het gemiddelde van de meetwaarden dus de échte waarde van grootheid aangeven. Maar in de praktijk moeten we het doen met een beperkt aantal replicaties.

In andere woorden, iedere meting is een kansexperiment, het meetresultaat is een kansvariabele, en we zijn geïnteresseerd in de verwachtingswaarde van die kansvariabele, $\mu_X$, die we zouden vinden als we de meting oneindig vaak zouden herhalen. Ons doel is om die parameter te schatten op basis van een eindig aantal experimenten.

Om dat te doen, hebben we weer een statistisch model nodig. Vaak worden daarbij de volgende aannames gemaakt:

Definitie: het Normale Model voor herhaalde metingen

De uitkomsten van de metingen zijn normaal verdeeld. Resultaten van herhaalde metingen hebben vaak een klokvormige verdeling. Dit model gaat ervan uit dat die verdeling benaderd kan worden door een normale verdeling.
Iedere meting $X_i$ is onafhankelijk uit die normale verdeling getrokken. Dit is het geval als de verschillende replicaties daadwerkelijk onafhankelijk van elkaar zijn uitgevoerd. Dat is niet altijd het geval. We gaan het hier uitgebreid over hebben in Paragraaf 24.9.1.

Als je deze modelaannames vergelijkt met het Normale Model voor steekproeven (Paragraaf 24.3) dan zie je dat die twee modellen precies overeenkomen! In beide gevallen behandelen we de gegevens als uitkomsten van een kansexperiment. In het ene geval bestaat dat kansexperiment uit het aselect trekken van een individu uit een populatie, in het andere geval uit het uitvoeren van een meting met een toevallige fout. Maar in beide gevallen veronderstelt het model dat de waarnemingen onafhankelijk getrokken zijn uit een normale verdeling.

Conclusie:

Het Normale Model voor steekproeven en het Normale Model voor herhaalde metingen komen op hetzelfde neer. Alle resultaten uit dit hoofdstuk, inclusief de formules voor BHI’s, kunnen dus net zo goed gebruikt worden voor gegevens die het resultaat zijn van een steekproef als voor gegevens die het resultaat zijn van herhaalde metingen.

Oefening 24.7 (Het meten van de smelttemperatuur van een DNA-sequentie)

Je meet de smelttemperatuur van een specifieke DNA-sequentie met behulp van UV-absorptie-metingen. Je repliceert de meting tien keer in onafhankelijke experimenten.

De resultaten zijn in de volgende vector gegeven (in $\degree$C):

T_m <- c(60.7, 58.3, 60.0, 62.7, 61.7, 60.5, 57.5, 59.1, 59.3, 55.9)

Ga ervan uit dat de toevallige fout normaal verdeeld is.

Bereken in R het 95%-betrouwbaarheidsinterval voor de smelttemeratuur.

24.9 Afwijkingen van het Normale Model

Bij de afleiding van de formule voor het BHI (Vergelijking 24.9) zijn we steeds uitgegaan van het Normale Model. In de praktijk zijn de aannames van het Normale Model zelden exact waar. Met name wijkt de verdeling van de variabele in de populatie vaak af van de normale verdeling. Kunnen we de formules die we geleerd hebben dan niet gebruiken? Daar gaan we in deze paragraaf op in.

Omdat het Normale Model niet altijd van toepassing is is het essentieel om, voordat je een BHI uitrekent met Vergelijking 24.9, eerst zorgvuldig na te gaan of dat wel verstandig is. We zullen het daarom eerst hebben over manieren waarop je afwijkingen van het Normale Model kunt herkennen. Daarna laten we zien dat zulke afwijkingen niet altijd een probleem zijn; de formule voor het BHI is behoorlijk “robuust” voor afwijkingen.

Is het Normale Model bij voorbaat wel plausibel?

Een eerste stap voordat je een BHI uitrekent is om te bedenken of de aannames van het Normale Model bij voorbaat wel plausibel zijn voor de toepassing.

Een belangrijke aanname is dat de variabele in de populatie normaal verdeeld is. Voor sommige variabelen is dat bij voorbaat niet aannemelijk. Dat kan verschillende redenen hebben; hier zijn er een paar:

Het is bekend dat de verdeling van de variabele scheef is. Bijvoorbeeld, het is algemeen bekend dat de verdeling van inkomen in Westerse landen heel erg rechts-scheef is. Hieronder is als voorbeeld in Figuur 24.6 de verdeling van persoonlijk inkomen in Nederland weergegeven, op basis van gegevens over 2022 van het Centraal Bureau voor de Statistiek.¹ Wanneer je een BHI wilt uitrekenen van zo’n variabele heb je dus bij voorbaat redenen om op je hoede te zijn.

Code

library(readr)
VerdelingPersoonlijkInkomen_CBS_2022 <- read_csv("data/VerdelingPersoonlijkInkomen-CBS-2022.csv", 
    locale = locale(decimal_mark = ","))

ggplot(VerdelingPersoonlijkInkomen_CBS_2022, 
       aes(y = Totaal, x = `persoonlijk inkomen (1000 euro)`)) +
  geom_line(color = "#e8601d", size = 1) +
  labs(x = "Persoonlijk inkomen (x 1000 euro)", y = "Aantal personen (x 1000)") +
  theme_minimal()

Figuur 24.6: Verdeling van persoonlijk inkomen in Nederland, 2022

De populatie bestaat uit meerdere deelpopulaties met vermoedelijk sterk verschillende gemiddelden. In dat geval zou de verdeling van de variabele multimodaal kunnen zijn of anderszins kunnen afwijken van een normale verdeling.
De variabele kan alleen waarden aannemen in een bepaald interval. Een normaal-verdeelde variabele kan iedere mogelijke waarde hebben, van $-\infty$ tot $\infty$. Vaak geldt dat niet voor biologische variabelen. Concentraties kunnen bijvoorbeeld nooit negatief zijn, en in veel toepassingen kunnen percentages alleen van 0 tot 100% lopen. Een normale verdeling kan dan nog steeds een goede benadering zijn, maar alleen als de bulk van de verdeling ver verwijderd is van de grenzen van het toelaatbare interval.
De variabele is discreet. Discrete variabelen kunnen bij benadering normaal verdeeld zijn als ze gemiddeld grote waarden hebben. Maar als de discrete waarden klein zijn is dat duidelijk niet het geval. Bijvoorbeeld, de vink (Fringilla Coelebs Africana) legt typisch 3 of 4 eieren; zie Figuur 24.7. Een normale verdeling is voor aantal eieren dus zeker geen model.

Code

dataVink <- data.frame(
  aantal_eieren <- c(rep(1, 3), rep(2, 2), rep(3, 15), rep(4, 14))
)

ggplot(dataVink, aes(aantal_eieren)) +
  geom_histogram(
    binwidth = 1, 
    color = "black", 
    fill = "#e8601d"
    ) +
  labs(x = "Aantal eieren per nest", y = "Frequentie") +
  theme_minimal()

Figuur 24.7: Aantal eieren in 34 nesten van de vink *Fringilla Coelebs Africana*².

De tweede aanname van het Normale Model was dat iedere waarde onafhankelijk uit de normale verdeling is getrokken. Ook over die aanname moet je goed nadenken. Op basis van de opzet van de steekproef of van het experiment heb je soms goede redenen om daaraan te twijfelen. Voorbeelden:

Pseudo-replicatie. In het kort verwijst pseudo-replicatie naar metingen of experimenten die meerdere keren zijn uitgevoerd onder dezelfde omstandigheden, zoals in één batch, op dezelfde dag, of met dezelfde materialen. Eventuele toevallige invloeden van die omstandigheden zijn daardoor bij al die pseudo-replicaties gelijk en dus zijn de replicaties niet onafhankelijk. De uitkomsten van zulke pseudo-replicaties zullen daardoor dichter bij elkaar liggen dan bij replicaties die daadwerkelijk onafhankelijk zijn uitgevoerd. Dat leidt onterecht tot een smaller BHI.
Tijdseries. Metingen over de tijd van dezelfde variabele zullen vaak niet onafhankelijk zijn. Een voorbeeld is de lichaamstemperatuur van één en dezelfde persoon gemeten op opeenvolgende dagen. Wie vandaag koorts heeft, heeft dat vaak morgen ook; metingen op opeenvolgende dagen zijn dus niet onafhankelijk. Je kunt bij het berekenen van een tijdsgemiddelde dus niet zomaar uitgaan van een model waarbij opeenvolgende metingen onafhankelijk zijn.
Steekproeven die bestaan uit groepjes. Een andere situatie is een steekproef die bestaat uit samenhangende groepjes. Bijvoorbeeld, stel dat een reisorganisatie alle reizigers om hun mening vraagt, maar de reizigers vaak reizen in paren of als gezin. Dan is er een gevaar dat de meningen van reisgezellen niet onafhankelijk zijn.

Kortom, het is belangrijk om goed na te denken over de aard van de variabele en de manier waarop de gegevens tot stand zijn gekomen om zo te bepalen of het Normale Model van toepassing lijkt.

Geven de gegevens aanleiding tot twijfel?

Vervolgens is het belangrijk om de gegevens zelf te bekijken om zo te bepalen of die aanleiding geven tot twijfel over het Normale Model.

Als de populatie een verdeling heeft die sterk afwijkt van een normale verdeling, dan verwacht je dat ook te zien in een steekproef uit die populatie (tenzij die steekproef heel klein is). Voordat je een BHI berekent of welke statistische analyse dan ook uitvoert moet je dus altijd eerst je gegevens goed bestuderen.

In Hoofdstuk 20 heb je geleerd hoe je de verdeling van een variabele kunt beschrijven en visualiseren. Voordat je een BHI berekent pas je dus altijd die methoden toe. Dat wil zeggen:

Bereken kengetallen.
- Bekijk het aantal waarnemingen $n$.
- Bereken maten van ligging en spreiding en bekijk die kritisch. Als de mediaan bijvoorbeeld sterk afwijkt van het gemiddelde, dan is dat een indicatie van een scheve verdeling of bijzondere uitbijters.
Visualiseer de gegevens.

Maak een histogram en/of een boxplot en bekijk die goed. Ziet de histogram er symmetrisch en klokvormig uit? Zijn de snorharen van de boxplot ongeveer symmetrisch, ligt de mediaan ongeveer in het midden van de box, en zijn er opvallende uitbijters? Als de verdeling er vreemd of scheef uitziet is het Normale Model mogelijk niet van toepassing.
Bestudeer een Normale Kwantiel-Kwantiel-plot.

Er bestaat een grafiek die speciaal ontworpen is om afwijkingen van de normale verdeling beter zichtbaar te maken. Deze grafiek heet de Normale Kwantiel-Kwantiel-plot (Normal QQ plot). We zullen die plot in het de volgende paragraaf even apart bespreken.

Door altijd de bovenstaande beschrijvende statistieken en figuren te bestuderen krijg je een goed beeld van de toepasbaarheid van het Normale Model.

Het Normale Kwantiel-Kwantiel-plot (normal QQ-plot)

Een Kwantiel-Kwantiel-plot (Quantile-Quantile plot, QQ-plot) is een grafiek waarmee je op het oog kunt controleren of een dataset een bepaalde verdeling volgt, zoals een normale verdeling.

Stel, je hebt in een experiment 101 waarnemingen verzameld. We sorteren de waarden van klein naar groot. De mediaan is dan de 51e waarde.

We willen controleren of deze gegevens uit een normale verdeling getrokken zijn. Als dat zo is, dan verwachten we dat de mediaan van de waarnemingen in de buurt ligt van het gemiddelde van de waarnemingen, want de normale verdeling is symmetrisch. Het eerste kwartiel (Q1) van de waarnemingen verwachten we ongeveer 0,67 standaarddeviaties onder het gemiddelde, want daar ligt voor de normale verdeling het eerste kwartiel.

Op deze manier kunnen we de kwantielen van de gegevens vergelijken met de kwantielen die we theoretisch hadden verwacht als de gegevens uit een normale verdeling getrokken waren.

Dat is in een notendop wat een QQ-plot doet. Voor iedere waarneming wordt op de $y$-as het kwantiel geplot dat bij die waarneming hoort, en op de $x$-as het kwantiel dat je had verwacht als de gegevens normaal verdeeld waren. Hier zijn twee voorbeelden:

(a) QQ-plot van een steekproef uit een normale verdeling

Als de gegevens uit een normale verdeling getrokken zijn komen de kwantielen van de waarnemingen grofweg overeen met de theoretische kwantielen. In dat geval vallen de datapunten ongeveer op een rechte lijn, zoals in Figuur 24.8 (a). Systematische afwijkingen van die rechte lijn geven aan dat de gegevens afwijkingen vertonen van een normale verdeling, zoals in Figuur 24.8 (b). Zo kun je beoordelen of de gegevens er “normaal” uitzien.

Het tekenen van een QQ-plot is met behulp van R erg gemakkelijk. Als gegevens een numerieke vector is, dan werkt het zo:

qqnorm(gegevens)  # normale kwantielen-kwantielen-plot (normal QQ-plot)
qqline(gegevens)  # voeg er een lijn aan toe die de theoretische voorspelling aangeeft

Oefening 24.8 (Normale QQ-plots en histogrammen) Hieronder vind je een figuur met vijf histogrammen en vijf normale QQ-plots van variabelen uit de NHANES dataset. Kun jij uitvinden welke QQ-plot bij welke histogram hoort?

De centrale limietstelling en de robuustheid van het Normale Model

Als de verdeling van de variabele in de populatie niet normaal is, gaat de afleiding van het BHI (Paragraaf 24.7) niet precies op. In deze paragraaf zullen we uitleggen dat de formule voor het BHI dan vaak toch nog een heel goede benadering is.

Om dat te laten zien, voeren we weer een aantal simulaties uit. In Figuur 24.9 (a) is de verdeling van lichaamsgewicht (Weight) uit de dataset NHANES weergegeven. We stellen ons voor dat deze verdeling de populatie voorstelt van personen uit de VS, waarvan wij het gemiddelde niet weten. Om iets te weten te komen over het gemiddelde van die populatie voeren we een steekproef uit. Figuur 24.9 (b) laat een kleine steekproef zien van $n = 3$ personen, en Figuur 24.9 (c) een iets grotere steekproef van $n = 30$. Paarse lijntjes geven het steekproefgemiddelde aan.

Wat zegt zo’n steekproefgemiddelde over het gemiddelde van de populatie? In Paragraaf 24.4.1 hebben we dat onderzocht door ons voor te stellen dat we heel vaak een steekproef trekken en telkens het steekproefgemiddelde uitrekenen. De verdeling van die steekproefgemiddelden vertelt ons hoe ver een steekproefgemiddelde typisch van het populatiegemiddelde afligt. Voor een variabele die in de populatie normaal verdeeld is, hadden we gezien dat het steekproefgemiddelde dat óók is. Dat feit hebben we gebruikt in de afleiding van de formule voor het BHI, Vergelijking 24.9.

De verdeling van lichaamsgewicht is echter bimodaal en bovendien scheef. De aannames van het Normale Model zijn dus niet van toepassing. Hoe zou de steekproevenverdeling van het gemiddelde er in dit geval uitzien?

Om dat te demonstreren nemen we heel vaak een steekproef van gewichten uit de NHANES dataset, berekenen steeds het gemiddelde, en tekenen een histogram van de steekproefgemiddelden. Het resultaat voor steekproeven met $n = 3$ is weergegeven in Figuur 24.9 (d), het resultaat voor $n = 30$ in Figuur 24.9 (e).

Het resultaat is heel verrassend: hoewel de populatie niet normaal verdeeld is, ziet de verdeling van de steekproefgemiddelden er wel behoorlijk normaal uit! Om dat te illustreren is in de figuren in het paars een normale verdeling toegevoegd.

Code

# Controleer en laad benodigde pakketten
if (!("ggplot2" %in% .packages())) {
  library(ggplot2)
}
if (!("NHANES" %in% .packages())) {
  library(NHANES)
}
if (!("patchwork" %in% .packages())) {
  library(patchwork)
}

set.seed(125)  # Voor reproduceerbaarheid

# Functie om een histogram te maken
maak_histogram <- function(
  data, xvariabele, x_label, binbreedte = 5, 
  x_bereik = c(0, 200), kleur = "#e8601d", alpha = 0.7, 
  voeg_normaal_toe = FALSE
) {
  # Basis histogram
  plot <- ggplot(data, aes_string(x = xvariabele)) +
    geom_histogram(
      fill = kleur, alpha = alpha, binwidth = binbreedte, 
      boundary = 0, color = "black"
    ) +
    labs(x = x_label, y = "Frequentie") +
    scale_x_continuous(
      limits = x_bereik,
      breaks = seq(x_bereik[1], x_bereik[2], by = 50)
    ) +
    theme_minimal()
  
  # Voeg normale verdeling toe indien gespecificeerd
  if (voeg_normaal_toe) {
    # Bereken parameters van de normale verdeling
    gemiddelde <- mean(data[[xvariabele]], na.rm = TRUE)
    standaarddev <- sd(data[[xvariabele]], na.rm = TRUE)
    schaal <- nrow(data) * binbreedte
    
    # Voeg normale verdeling toe aan de plot
    plot <- plot +
      stat_function(
        fun = function(x) {
          schaal * dnorm(x, mean = gemiddelde, sd = standaarddev)
        },
        color = "DarkOrchid", size = 1
      )
  }
  
  return(plot)
}

# Functie om een verticale lijn toe te voegen
voeg_lijn_toe <- function(plot, x_waarde, 
    lijnkleur = "DarkOrchid", linetype = "dashed", size = 1) {
  plot +
    geom_vline(
      aes(xintercept = x_waarde), color = lijnkleur,
      linetype = linetype, size = size
    )
}

# Functie om steekproef-histogram te maken
maak_steekproef_histogram <- function(data, xvariabele, 
    steekproefgrootte, x_label, binbreedte = 5, 
    x_bereik = c(0, 200), kleur = "#e8601d", alpha = 0.7) {
  steekproef <- sample(
    data[[xvariabele]], size = steekproefgrootte, replace = TRUE
  )
  steekproefgemiddelde <- mean(steekproef)
  
  plot <- maak_histogram(
    data.frame(steekproef), "steekproef", x_label,
    binbreedte, x_bereik, kleur, alpha
  )
  
  voeg_lijn_toe(plot, steekproefgemiddelde, lijnkleur = "DarkOrchid")
}

# Functie om histogram van steekproefgemiddelden te maken
maak_gemiddelden_histogram <- function(data, xvariabele, 
    steekproefgrootte, n_steekproeven, x_label, 
    binbreedte = 5, x_bereik = c(0, 200), kleur = "#e8601d", 
    alpha = 0.7, voeg_normaal_toe = FALSE) {
  steekproefgemiddelden <- replicate(
    n_steekproeven,
    mean(
      sample(
        data[[xvariabele]], size = steekproefgrootte, replace = FALSE
      )
    )
  )
  
  maak_histogram(
    data.frame(steekproefgemiddelden), "steekproefgemiddelden",
    x_label, binbreedte, x_bereik, kleur, alpha, voeg_normaal_toe
  )
}

# Parameters
steekproefgrootte <- 3
n_steekproeven <- 10000
x_label_boven_midden <- "Lichaamsgewicht (kg)"
x_label_onder <- "Steekproefgemiddelde Lichaamsgewicht (kg)"
binbreedte <- 5
x_bereik <- c(0, 200)

# Data; filter regels met missende gegevens
data <- NHANES[ !is.na( NHANES$Weight ), ]

# Maak de drie histogrammen
boven_plot <- maak_histogram(
  data, "Weight", x_label_boven_midden, binbreedte, x_bereik
) 
midden_plot <- maak_steekproef_histogram(
  data, "Weight", steekproefgrootte, x_label_boven_midden, 
  binbreedte, x_bereik
)
onder_plot <- maak_gemiddelden_histogram(
  data, "Weight", steekproefgrootte, n_steekproeven, 
  x_label_onder, binbreedte, x_bereik,
  voeg_normaal_toe = TRUE
)

# Parameters
steekproefgrootte <- 30  # Grootte van de steekproef
binbreedte_onder <- 2

# Maak twee histogrammen
midden_plot_b <- maak_steekproef_histogram(data, "Weight", 
    steekproefgrootte, x_label_boven_midden, binbreedte, x_bereik
)
onder_plot_b <- maak_gemiddelden_histogram(data, "Weight", 
    steekproefgrootte, n_steekproeven, x_label_onder, 
    binbreedte = binbreedte_onder, 
    x_bereik,
  voeg_normaal_toe = TRUE
)

boven_plot
midden_plot
midden_plot_b
onder_plot
onder_plot_b

(a) Histogram van lichaamsgewicht (`Weight`) in de NHANES dataset.

Waar komen die normale verdelingen vandaan? Intuïtief kun je je misschien voorstellen dat de bimodaliteit van de populatie niet terug te zien is in de gemiddelden van de steekproeven. Als de steekproef groot genoeg is, zal de steekproef typisch zowel kinderen bevatten (uit de linkerpiek van de verdeling) als volwassenen (uit de rechterpiek van de verdeling). Het gemiddelde komt daardoor meestal in het midden terecht.

Wiskundig kan daadwerkelijk worden aangetoond dat de verdeling van de steekproefgemiddelden meer op een normale verdeling gaat lijken naarmate de steekproef groter wordt; dat geldt ongeacht de verdeling van de populatie (behalve in heel pathologische gevallen waarover wij ons niet druk hoeven te maken). Dit volgt uit de centrale limietstelling, die we ook al kort tegenkwamen in Paragraaf 23.7.2.

Een consequentie van deze stelling is dat de afleiding van de formule voor het BHI (Paragraaf 24.7) bij goede benadering ook opgaat als de verdeling van de variabele in de populatie niet normaal is, mits $n$ groot genoeg is. Precies hoe groot $n$ moet zijn hangt af van de mate waarin de verdeling van de populatie afwijkt van de normale verdeling. Als de populatieverdeling heel erg scheef is en/of extreme uitbijters kent, dan is een grotere $n$ nodig. Als vuistregel kun je aannemen dat bij steekproeven van $n = 30$ of meer de formule voor het BHI ook van toepassing is als de populatie niet normaal verdeeld is. In Oefening 24.11 hieronder zul je zien dat zelfs bij behoorlijk vreemde verdelingen het BHI al bij veel kleinere $n$ betrouwbaar is.

Wat als het statistische model echt niet van toepassing is?

We hebben hierboven gezien dat het Normale Model vaak correcte resultaten geeft zelfs als de variabele niet daadwerkelijk normaal verdeeld is. Maar wat nou als de verdeling van de variabele sterk van de normale verdeling afwijkt en het aantal waarnemingen beperkt is?

Er zijn verschillende opties. In deze cursus zullen we die niet in detail kunnen behandelen, maar het is belangrijk om te weten dat er oplossingen bestaan die je in een vervolgcursus zou kunnen leren. Hier zijn een paar voorbeelden:

Andere statistische modellen: Als de normale verdeling geen goed model is voor de toepassing, kunnen we een ander statistisch model formuleren dat beter bij de toepassing past. In sommige gevallen is een model dat uitgaat van de Poisson-verdeling geschikt, in andere gevallen is wellicht een Binomiale verdeling van toepassing (zie Paragraaf 23.11). Voor zulke aangepaste modellen gelden andere formules voor het BHI, die we hier niet zullen behandelen.
Transformeren van de variabele: Als een variabele $X$ rechts-scheef is, dan kan het zijn dat $\log{X}$ of $\sqrt{X}$ wél ongeveer normaal verdeeld is. Als $X$ een grote waarde heeft, dan is $\log{X}$ een stuk kleiner, en $\sqrt{X}$ ook. Uitbijters in de rechter-staart van de verdeling van $X$ worden daardoor minder extreem, waardoor de staart als het ware wordt ingetrokken. Een functie zoals de log toepassen op alle waarnemingen heet het transformeren van de variabele. Als de getransformeerde variabele wél ongeveer normaal verdeeld is, kun je een het Normale Model dus wel op die variabele toepassen. Dat levert je een BHI voor de getransformeerde variabele.
Computer-intensieve methoden: Als het statistische model te ingewikkeld is, dan kan het moeilijk zijn om een formule af te leiden voor het BHI. Maar in die gevallen kun je vaak wel een BHI schatten door middel van simulaties. Voorbeelden van zulke computer-intensieve methoden zijn bootstrapping, likelihood-methoden, en methoden gebaseerd op Monte-Carlo-simulaties.

24.10 Samenvatting

Betrouwbaarheidsintervallen algemeen

Puntschatting: Een waarde die berekend is op basis van een steekproef of een reeks metingen, en als schatting dient voor een populatieparameter.
95%-betrouwbaarheidsinterval: Een interval om de puntschatter dat zo berekend is dat het in 95% van de gevallen de populatieparameter bevat. Geeft de onzekerheid in de schatting aan.

Het Normale Model

Van een steekproef van $n$ eenheden is de variabele $X$ gemeten.

Het Normale Model neemt dan aan dat:

De variabele $X$ is normaal verdeeld in de populatie.
Iedere waarde $X_i$ van de steekproef is onafhankelijk uit deze normale verdeling getrokken.

Hetzelde model kan ook gebruikt worden voor $n$ herhaalde metingen die door toevallige meetfouten van elkaar verschillen:

De uitkomsten van de metingen zijn normaal verdeeld.
Iedere meting $X_i$ is onafhankelijk uit die normale verdeling getrokken.

Notatie

Populatieparameters: Eigenschappen van de verdeling van de populatie waaruit de steekproef wordt getrokken. Parameters worden in de statistiek meestal aangeduid met een Griekse letter.

Voorbeelden: gemiddelde $\mu_X$ en standaarddeviatie $\sigma_X$ van de populatie.
Eigenschappen van steekproeven: Eigenschappen van een steekproef worden meestal aangeduid met een gewone (Romaanse) letter. Het gemiddelde van een steekproef wordt aangeduid met $\overline{X}$, de standaarddeviatie van de steekproef met $s_X$.

95%-betrouwbaarheidsintervallen voor het populatiegemiddelde

We hebben variabele $X$ gemeten voor een steekproef van grootte $n$ uit een populatie. We gebruiken die steekproef om het populatiegemiddelde $\mu_X$ te schatten.

Of:

We hebben een serie van $n$ metingen ($X_i$) van een grootheid en gebruiken die metingen om de verwachtingswaarde $\mu_X$ te schatten.
We maken gebruik van het Normale Model.

Dan is het 95%-betrouwbaarheidsinterval:

als de standaarddeviatie van de populatie $\sigma$ bekend is (zeer ongebruikelijk!):

\[ \mu_X = \overline{X} \pm 1,96 \sigma_{\overline{X}}, \]

met

\[ \sigma_{\overline{X}} = \frac{\sigma}{\sqrt{n}}.\]
als de standaarddeviatie van de populatie $\sigma$ onbekend is:

\[ \mu_X = \overline{X} \pm t_{0{,}05 (2) \mathrm{(n-1)}} \mathrm{SE}_{\overline{X}}. \]

Hier is

\[ \mathrm{SE}_{\overline{X}} = \frac{s_X}{\sqrt{n}} \]

de standaardfout van het gemiddelde.

Beoordelen of het Normale Model van toepassing is

Bedenk: Verwacht ik dat de Normale verdeling deze variabele redelijk kan beschrijven?

Voorbeelden van gevaren:

De verdeling is
- scheef
- multimodaal
- beperkt tot een interval
- discreet
Bedenk: Is het te verwachten dat de waarnemingen (bijna) onafhankelijk zijn?

Voorbeelden van gevaren:
- pseudoreplicatie
- tijdsseries
- de steekproef bestaat uit groepjes

Bestudeer de gegevens om afwijkingen van de normale verdeling te ontdekken

Bereken kengetallen en beoordeel ze.
Visualiseer de gegevens en beoordeel ze.
Maak een QQ-plot en beoordeel die.

Robuustheid van het Normale Model

Als $n$ groot is, is de steekproevenverdeling van het gemiddelde ongeveer normaal verdeeld, ook als de verdeling van de variabele in de populatie dat niet is. (Centrale limietstelling) Daardoor is het Normale Model vaak ook bruikbaar als de populatie niet normaal verdeeld is.

Vuistregel: als $n >30$ werken de formules voor het BHI ook goed voor populaties die niet normaal verdeeld zijn.

24.11 Terminologie

Tabel 24.1: Woordenlijst Hoofdstuk 24 .

Nederlands	Engels	Betekenis
betrouwbaarheidsinterval (BHI)	confidence interval (CI)	Een interval om de puntschatter dat zo berekend is dat het in 95% van de gevallen de populatieparameter bevat. Geeft de onzekerheid in de schatting aan.
centrale limietstelling	central limit theorem	Wiskundige stelling die zegt dat het gemiddelde van een groot aantal kansvariabelen convergeert naar een normale verdeling.
repliceren	replicate	herhalen (van een meting of experiment)
pseudoreplicatie	pseudoreplication	Het herhalen van een metingen of experiment op een manier waardoor de metingen niet onafhankelijk zijn.
puntschatting	point estimate	Een schatting zonder betrouwbaarheidsinterval.
standaardfout van het gemiddelde	Standard Error of the Mean (SEM)	Schatting van de standaarddeviatie van het steekproefgemiddelde.
statistische modellen	statistical models	Modellen van de kansprocessen die verantwoordelijk zijn voor de onzekerheid in de gegevens.
steekproevenverdeling van het gemiddelde	sampling distribution of the mean	Het nemen van een steekproef is een kansproces; het gemiddelde van de steekproef een kansvariabele. De steekproevenverdeling van het gemiddelde is de kansverdeling van steekproefgemiddelden.

24.12 Opgaven

Oefening 24.9 (BHI van een steekproef tomaten)

Van zestien aselect gekozen tomaten uit een kas is het gewicht bepaald. Het gemiddelde gewicht in de steekproef is gelijk aan 44,1 g; de standaardafwijking van de steekproef is gelijk aan 8,6 g. Bepaal het 95%-betrouwbaarheidsinterval voor het gemiddelde gewicht $\mu_X$ van dergelijke tomaten. (Je mag uit gaan van het Normale Model.)

Oefening 24.10 (Een schatting voor de gemiddelde rusthartslag.)

In Paragraaf 24.9 hebben we gezien hoe je kunt onderzoeken of een steekproef van een variabele afwijkingen laat zien van de normale verdeling. Als voorbeeld ga je de variabele Pulse (hartslagfrequentie) onderzoeken uit de NHANES-database.

Om te beginnen: zou je bij voorbaat verwachten dat een normale verdeling een redelijk model is voor de verdeling van hartslagfrequenties?

Zijn er redenen om daaraan te twijfelen?
Laad eerst het pakket NHANES in zoals je dat al vaak hebt gedaan.
Voor sommige personen in de dataset is geen hartslagfrequentie ingevoerd. In sommige regels staat daarom NA (not available). Selecteer alle gegevens die wel beschikbaar zijn:
```
# Verwijder missende waarden van de variabele 'Pulse'
pulse_data <- na.omit(NHANES$Pulse)
```
Hoeveel waarnemingen blijven er over?
Visualiseer de verdeling met een histogram en een boxplot. Zie je afwijkingen van een normale verdeling?
Bestudeer kengetallen voor de gegevens.
- Zijn de maten van locatie zoals verwacht?
- Vergelijk de mediaan met het gemiddelde; wat zegt je dat?
- Zijn de maten van spreiding (standaarddeviatie, IKA) ongeveer wat je had verwacht?
Voeg een Normale QQ-plot toe en beoordeel die.
Gebruik de functie t.test() om het BHI te berekenen.
Kun je het berekende BHI met een gerust hart gebruiken?

Oefening 24.11 (Robuustheid van het Normale Model controleren)

In Paragraaf 24.9.4 heb je geleerd over de Centrale Limietstelling. Die stelling zegt dat bij grote steekproeven de steekproevenverdeling van het gemiddelde bij benadering normaal verdeeld is, zelfs als de verdeling van populatie dat niet is. We gaan dat in deze opgave met een voorbeeld controleren.

We zagen dat de variabele lichaamsgewicht (Weight) in de dataset NHANES bimodaal was (Figuur 24.9 (a)). Het idee is om de personen in deze dataset te beschouwen als een populatie, en uit deze populatie een steekproef van $n = 30$ personen te nemen. Op basis van deze steekproef berekenen we een BHI. Als de berekening van het BHI robuust is (ondanks de bimodale verdeling), ligt het ware gemiddelde van de populatie 95% van de keren dat we dit doen binnen het BHI.

Maak een nieuw script in R. Laad daarin de library NHANES.
Voor sommige personen in de dataset mist het gewicht en is NA ingevuld. Sla de gewichten op in een vector gewichten, maar verwijder de NA’s, op dezelfde manier als in Oefening 24.10.

We gaan de gewichten in de dataset nu gebruiken als een populatie waaruit we een kleinere steekproef nemen.

Voeg een regel toe die het gemiddelde van deze “populatie” berekent.

Voeg een regel toe waarmee je een aselecte steekproef van 30 gewichten selecteert uit de populatie:

# replace = FALSE zorgt ervoor dat we niet twee keer dezelfde persoon
# in de steekproef kunnen opnemen
steekproefgrootte <- 30
steekproef <- sample(
    gewichten, size = steekproefgrootte, replace = FALSE
  )
steekproef

Bereken van die steekproef het BHI, zoals je dat hierboven geleerd hebt, met de functie t.test().

Ligt het ware populatiegemiddelde in het berekende BHI?

Als de formule van het Normale Model werkt, zou bij 95% van de keren dat je een steekproef en een BHI berekent, het ware populatiegemiddelde in het interval moeten liggen.

Voer de code een paar keer uit om te zien of dat ongeveer kan kloppen.

Om het BHI goed te controleren, willen we de code heel vaak uitvoeren en tellen hoe vaak het populatiegemiddelde in het interval ligt. Het script hieronder doet dat. Kopieer het, bestudeer het, en voer het uit op je laptop. Lijkt de formule voor het BHI goed te werken?

Code

library("NHANES")
data(NHANES)

# PARAMETERS
aantal_steekproeven <- 50000   # Aantal steekproeven
steekproefgrootte <- 30        # Grootte van elke steekproef

# POPULATIEGEGEVENS 
gewichten <- na.omit(NHANES$Weight)      # filter NA's weg
populatiegemiddelde <- mean(gewichten)   # bereken populatiegemiddelde

# SIMULATIE VOORBEREIDEN

# Functie om het BHI voor een steekproef te berekenen
# en te controleren of het populatiegemiddelde 
# binnen het betrouwbaarheidsinterval valt.
controleer_bhi <- function( steekproef, populatiegemiddelde) {
  resultaten.t.test <- t.test(steekproef)
  bhi <- resultaten.t.test$conf.int
  in_interval <- populatiegemiddelde >= bhi[1] && populatiegemiddelde <= bhi[2]

  return(in_interval)
}

# Logical vector om de resultaten in op te slaan
resultaten <- logical(aantal_steekproeven)

# SIMULEER STEEKPROEVEN
for (i in 1:aantal_steekproeven) {
  steekproef <- sample(
    gewichten, size = steekproefgrootte, replace = FALSE
  )
  resultaten[i] <- controleer_bhi(steekproef, populatiegemiddelde)
}

# Bereken en print percentage
print("Percentage steekproeven waarbij populatiegemiddelde binnen het BHI valt:")
print(mean(resultaten)*100)

Bij kleine steekproeven verwachten we dat het BHI zou kunnen afwijken. Voer het script opnieuw uit, maar dan met steekproeven van grootte 4. Wat concludeer je?

Bron: CBS, gedownload op 17 januari 2025.↩︎
Ramdani, K., M. Kouidri, M. L. Ouakid, and M. Houhamdi. “Breeding Biology of the Chaffinch Fringilla Coelebs Africana in the El Kala National Park (North East Algeria).” Arxius de Miscel·lània Zoològica, July 7, 2019, 109–21. Breeding biology.↩︎

24.1 Betrouwbaarheidsintervallen

24.2 Het interpreteren van onzekere gegevens vereist een statistisch model

24.3 Het Normale Model voor eenvoudige aselecte steekproeven uit grote populaties

24.4 De kansverdeling van het steekproefgemiddelde

Simulaties van het Normale Model

24.5 De standaarddeviatie van het steekproefgemiddelde.

24.6 95%-BHI van het gemiddelde als de standaarddeviatie van de populatie bekend is

24.7 BHI als de variantie van de populatie niet bekend is

De standaardfout van het steekproefgemiddelde

Transformeren naar \(t\)

Afleiding van het 95%-BHI

Een BHI berekenen met R

24.8 Het Normale Model is ook geschikt voor herhaalde experimenten met meetfouten

24.9 Afwijkingen van het Normale Model

Is het Normale Model bij voorbaat wel plausibel?

Geven de gegevens aanleiding tot twijfel?

Het Normale Kwantiel-Kwantiel-plot (normal QQ-plot)

De centrale limietstelling en de robuustheid van het Normale Model

Wat als het statistische model echt niet van toepassing is?

24.10 Samenvatting

Betrouwbaarheidsintervallen algemeen

Het Normale Model

Notatie

95%-betrouwbaarheidsintervallen voor het populatiegemiddelde

Beoordelen of het Normale Model van toepassing is

Bestudeer de gegevens om afwijkingen van de normale verdeling te ontdekken

Robuustheid van het Normale Model

24.11 Terminologie

24.12 Opgaven