Eén benadering is het berekenen van de normale benadering van de binomiale verdeling (Wald-interval), zoals beschreven in het OpenIntro-boek. Dit is eigenlijk hetzelfde als de ‘handmatige’ berekening, maar dan in R.
Als voorbeeld gebruiken we een enquête onder kiezers waarin 23% aangeeft de Conservatieven te steunen. We schatten een 95% betrouwbaarheidsinterval voor deze schatting:
# Input definiëren (verander deze waarden)
p_hat <- 0.23
n <- 1000
confidence <- 0.95
# Voer berekeningen uit (je hoeft hier niets te veranderen)
se <- sqrt(p_hat * (1 - p_hat)/ n) # Standaardfout berekenen
z_star <- qnorm((1 + confidence) / 2) # Z*-waarde berekenen
lower <- p_hat - z_star * se # Ondergrens betrouwbaarheidsinterval
upper <- p_hat + z_star * se # Bovengrens betrouwbaarheidsinterval
c(lower, upper)[1] 0.203917 0.256083Wanneer je alleen beknopte gegevens hebt, d.w.z. informatie over de steekproefgrootte en de steekproefproportie, raden we aan prop.test te gebruiken voor het berekenen van een betrouwbaarheidsinterval voor een enkele proportie.1 De functie prop.test maakt deel uit van het package stats, dat een van de weinige packages is die R automatisch laadt wanneer het start.
# Input definiëren (verander deze waarden)
p_hat <- 0.23
n <- 1000
# Voer berekeningen uit (pas zo nodig alleen conf.level aan)
prop.test(x = p_hat * n,
n = n,
conf.level = 0.95)prop.testDit is de functie die een test uitvoert voor een proportie en het bijbehorende betrouwbaarheidsinterval.
x = p_hat * nDe eerste opgegeven waarde moet het aantal successen zijn. In ons voorbeeld is dat het aantal mensen dat de Conservatieven steunt. We kunnen dit berekenen als de proportie in de steekproef ($ $) keer de steekproefomvang.
n = nDe tweede opgegeven waarde is het aantal proeven (lees: het aantal waarnemingen in onze dataset). In ons geval: het aantal respondenten in de dataset.
conf.level = 0.95Dit bepaalt het betrouwbaarheidsniveau. De standaardwaarde is 0.95, wat overeenkomt met een betrouwbaarheidsinterval van 95%.
1-sample proportions test with continuity correction
data: p_hat * n out of n, null probability 0.5
X-squared = 290.52, df = 1, p-value < 2.2e-16
alternative hypothesis: true p is not equal to 0.5
95 percent confidence interval:
0.2045014 0.2576046
sample estimates:
p
0.23 Er is een kleine afwijking t.o.v. de handmatige berekening, omdat R automatisch een continuïteitscorrectie toepast en een Wilson score interval gebruikt (een iets geavanceerdere versie van het betrouwbaarheidsinterval).
Als je een dataframe hebt met een variabele die voor elk geval het succes (of niet) weergeeft, kunnen we de volgende procedure gebruiken. In ons voorbeeld hebben we een variabele die de stemintentie van een respondent registreert, die ofwel Conservative ofwel Other party is:
# Voor dit voorbeeld maken we een dataset van 1000 respondenten waarvan er 230 een voorkeur voor de Conservatieven hebben.
example_data <- data.frame(party_choice = factor(c(rep("Conservatives", 230),
rep("Other party", 770))))
table(example_data$party_choice)
Conservatives Other party
230 770 Als we dergelijke gegevens hebben, kunnen we direct het betrouwbaarheidsinterval berekenen door prop.test uit te voeren voor de tabel:
prop.test(table(example_data$party_choice),
conf.level = 0.95)table(example_data$party_choice)Het eerste argument is een tabel met onze variabele, die we met deze code kunnen aanmaken. Als je je eigen gegevens gebruikt, vervang dan example_data door de naam van je dataframe en party_choice door de naam van de variabele. 2
Noot: De variabele/tabel mag slechts twee categorieën bevatten. R zal het betrouwbaarheidsinterval berekenen voor de eerste categorie in de tabel (in ons geval: Conservatieven). 3
conf.level = 0.95Dit bepaalt het betrouwbaarheidsniveau. De standaardwaarde is 0.95, wat overeenkomt met een betrouwbaarheidsinterval van 95%.
1-sample proportions test with continuity correction
data: table(example_data$party_choice), null probability 0.5
X-squared = 290.52, df = 1, p-value < 2.2e-16
alternative hypothesis: true p is not equal to 0.5
95 percent confidence interval:
0.2045014 0.2576046
sample estimates:
p
0.23 De resultaten zijn identiek aan de berekening op basis van de samengevatte gegevens.
Het 95% betrouwbaarheidsinterval voor de steun van kiezers voor de Conservatieve partij loopt van 20,5% tot 25,8% (N = 1000).
In deze peiling (N = 1000) was de steun van kiezers voor de Conservatieve partij gelijk aan 23%, 95% CI [20,5%; 25,8%].
Neem dus naast een omschrijving in je rapportage op:
De onder- en bovengrens van het interval.
Om welk betrouwbaarheidsinterval het gaat: 90%, 95%, 99% (of nog anders)
Zeker bij peilingen is het gebruikelijk en aanbevelingswaardig om de steekproefomvang te vermelden.
Wanneer je alleen informatie hebt over de steekproefgrootte en de proportie, raden we aan prop.test te gebruiken voor het uitvoeren van een hypothesetoets.4
# Input definiëren (verander deze waarden)
p_hat <- 0.23
n <- 1000
# Voer berekeningen uit
prop.test(x = p_hat * n,
n = n,
p = 0.25,
alternative = "two.sided")prop.testDit is de functie die een test uitvoert voor een proportie en het bijbehorende betrouwbaarheidsinterval.
x = p_hat * nDe eerste opgegeven waarde moet het aantal successen zijn. In ons voorbeeld is dat het aantal mensen dat de Conservatieven steunt. We kunnen dit berekenen als de proportie in de steekproef (\(hat{p}\)) keer de steekproefomvang.
n = nDe tweede opgegeven waarde is het aantal proeven (lees: het aantal waarnemingen in onze dataset). In ons geval: het aantal respondenten in de dataset.
p = 0.25Dit specificeert de waarde onder de nulhypothese. In dit voorbeeld is dit 0.25 of 25%, dus testen we de nulhypothese \(mu_0 = 0.25\).
alternative = "two.sided"Bepaalt of we een tweezijdige of een eenzijdige test willen gebruiken. Opties zijn “tweezijdig” (standaard), “less” (als \(H_1: \mu < p\)) of “greater” (als \(H_1: \mu > p\)).
1-sample proportions test with continuity correction
data: p_hat * n out of n, null probability 0.25
X-squared = 2.028, df = 1, p-value = 0.1544
alternative hypothesis: true p is not equal to 0.25
95 percent confidence interval:
0.2045014 0.2576046
sample estimates:
p
0.23 Als je een dataframe hebt met een variabele die voor elk geval het succes (of niet) weergeeft, kunnen we de volgende procedure gebruiken. In ons voorbeeld hebben we een variabele die de stemintentie van een respondent registreert, die ofwel Conservative ofwel Other Party is (voor het voorbereiden van de data zie hierboven Gegevens in een dataframe).
prop.test(table(example_data$party_choice),
p = 0.25,
alternative = "two.sided")table(example_data$party_choice)Het eerste argument is een tabel met onze variabele, die we met deze code kunnen aanmaken. Als je je eigen gegevens gebruikt, vervang dan example_data door de naam van je dataframe en party_choice door de naam van de variabele.
Noot: De variabele/tabel mag slechts twee categorieën bevatten. R voert de test uit dat de proportie voor de eerste categorie (in ons geval ‘Conservatieven’) gelijk is aan de hypothetische waarde (p, zie hieronder).
p = 0.25Dit specificeert de waarde onder de nulhypothese. In dit voorbeeld is dit 0.25 of 25%, dus testen we de nulhypothese \(mu_0 = 0.25\).
alternative = "two.sided"Bepaalt of we een tweezijdige of een eenzijdige test willen gebruiken. Opties zijn “tweezijdig” (standaard), “less” (als \(H_1: \mu < p\)) of “greater” (als \(H_1: \mu > p\)).
1-sample proportions test with continuity correction
data: table(example_data$party_choice), null probability 0.25
X-squared = 2.028, df = 1, p-value = 0.1544
alternative hypothesis: true p is not equal to 0.25
95 percent confidence interval:
0.2045014 0.2576046
sample estimates:
p
0.23 of, equivalent (waarbij we de z-waarde rapporteren, die gelijk is aan de wortel van bovengenoemde \(\chi^2\)-waarde):
Neem dus naast een inhoudelijke omschrijving in je rapportage op:
De twee percentages/proporties die worden vergeleken.
De z-waarde of Chikwadraatwaarde. Als je de laatste gebruikt ook het aantal vrijheidsgraden tussen haken, dat is voor dit soort berekeningen altijd 1.
p = p-waarde. Maar: schrijf nooit \(p = 0,000\). Want de p-waarde is nooit precies nul, maar heel klein. Het is dan beter om te zeggen \(p < 0,001\). Als de hypothesetoets handmatig uitvoert schrijf je p < α-waarde, bijvoorbeeld \(p < 0,05\).
Hierboven hebben we een betrouwbaarheidsinterval berekend voor één proportie. In deze paragraaf behandelen we twee (onafhankelijke) proporties. De onderstaande berekeningen gelden voor het geval dat de twee proporties afkomstig zijn uit verschillende steekproeven of deelsteekproeven, bijvoorbeeld de steun voor een partij in twee afzonderlijke opiniepeilingen (deze groepen zijn onafhankelijk omdat de samenstelling van de steekproef doorgaans verschilt van steekproef tot steekproef) of het links-rechtsstandpunt van respondenten die in de hoofdstad wonen en respondenten die daar niet wonen (deze groepen zijn onafhankelijk omdat je er maar tot één kunt behoren).
We gebruiken de CPR dataset uit het openintro package (zie p. 218). De dataset bevat gegevens over patiënten die willekeurig werden verdeeld in een behandelingsgroep waar ze een bloedverdunner kregen of de controlegroep waar ze geen bloedverdunner kregen. De uitkomstvariabele was of de patiënten ten minste 24 uur overleefden.
library(tidyverse)
library(openintro)
library(flextable) data(cpr)
# Factorniveaus van variabelen aanpassen om ervoor te zorgen dat de volgorde hetzelfde is als in de tabel in het tekstboek
cpr <- cpr |>
mutate(group = factor(group, levels = c("control", "treatment"), ordered = TRUE),
outcome = factor(outcome, levels = c("survived", "died"), ordered = TRUE))
table_example <- proc_freq(x = cpr,
row = "group",
col = "outcome",
include.row_percent = FALSE,
include.column_percent = FALSE,
include.table_percent = FALSE)
table_examplegroup | outcome | ||
|---|---|---|---|
survived | died | Total | |
control | 11 | 39 | 50 |
treatment | 14 | 26 | 40 |
Total | 25 | 65 | 90 |
Merk op dat we gewoonlijk de onafhankelijke variabele in de kolommen zetten. Het boek volgt deze conventie echter niet.
Om een 90 % betrouwbaarheidsinterval van het verschil voor de overlevingspercentages (\(p_{1}\) en \(p_{2}\)) met de hand te berekenen, kunnen we schrijven:
# Definieer invoer (pas deze aan voor je eigen berekeningen)
p_hat_1 <- 14 / 40 # Treatment group
p_hat_2 <- 11 / 50 # Control group
n_1 <- 40
n_2 <- 50
confidence = 0.9
# Voer berekeningen uit (je hoeft hier niets te veranderen)
p_hat <- p_hat_1 - p_hat_2
se <- sqrt((p_hat_1*(1-p_hat_1)/n_1) + (p_hat_2*(1-p_hat_2)/n_2))
z_star <- qnorm((1 + confidence) / 2) # Bereken z*-value
lower <- p_hat - z_star * se # Ondergrens betrouwbaarheidsinterval
upper <- p_hat + z_star * se # Bovengrens betrouwbaarheidsinterval
c(lower, upper)[1] -0.02707706 0.28707706Merk op dat we de volgorde van het boek volgen (p1 is de behandelingsgroep en p2 is de controlegroep).
Wanneer je alleen samenvattende gegevens hebt, d.w.z. informatie over de steekproefgrootte en de steekproefproportie, raden we aan prop.test te gebruiken om een betrouwbaarheidsinterval te berekenen voor het verschil tussen twee proporties in twee steekproeven. De functie prop.test maakt deel uit van het package stats, dat een van de weinige packages is die R automatisch laadt bij het starten.
prop.test(x = c(14, 11),
n = c(40, 50),
conf.level = 0.90,
correct=FALSE)prop.testDit is de functie die hypothesetoetsen voor porporties uitvoert en bijbehorende betrouwbaarheidsintervallen berekent.
x = c(14, 11)Hier geven we het aantal “successen” voor elke groep aan (in dit geval het aantal overlevenden in elke groep).
n = c(40, 50)De tweede opgegeven waarde moet het aantal casussen zijn voor elk van de twee steekproeven. In ons geval: het aantal deelnemers in beide groepen van de studie.
conf.level = 0.90Dit bepaalt het betrouwbaarheidsniveau. De standaardwaarde is 0.95, wat overeenkomt met een betrouwbaarheidsinterval van 95%. Wij stellen het interval in op 0.90.
correct = FALSEDit zet de continuïteitscorrectie op FALSE. Dit zou dezelfde resultaten moeten opleveren als de handmatige berekeningen hierboven.
2-sample test for equality of proportions without continuity correction
data: c(14, 11) out of c(40, 50)
X-squared = 1.872, df = 1, p-value = 0.1712
alternative hypothesis: two.sided
90 percent confidence interval:
-0.02707706 0.28707706
sample estimates:
prop 1 prop 2
0.35 0.22 of
Neem dus naast een omschrijving in je rapportage op:
De onder- en bovengrens van het interval.
Om welk betrouwbaarheidsinterval het gaat: 90%, 95%, 99% (of nog anders)
Zeker bij peilingen is het gebruikelijk en aanbevelingswaardig om de steekproefomvang te vermelden.
We gebruiken de mammogram dataset uit het openintro package. De dataset bevat gegevens van een experiment waarbij 89.835 vrouwen werden gerandomiseerd om al dan niet een mammogram te ondergaan. De gemeten respons was of zij binnen 25 jaar aan borstkanker waren overleden.
library(tidyverse)
library(openintro)
library(flextable)
data(mammogram)
# Factorniveaus van variabelen aanpassen om ervoor te zorgen dat de volgorde hetzelfde is als in de tabel in het tekstboek
mammogram <- mammogram |>
mutate(breast_cancer_death = factor(breast_cancer_death, levels = c("yes", "no"), ordered=TRUE)) |>
mutate(treatment = factor(treatment, levels = c("mammogram", "control"), ordered = TRUE))
table_example <- proc_freq(x = mammogram,
row = "treatment",
col = "breast_cancer_death",
include.row_percent = FALSE,
include.column_percent = FALSE,
include.table_percent = FALSE)
table_exampletreatment | breast_cancer_death | ||
|---|---|---|---|
yes | no | Total | |
mammogram | 500 | 44,425 | 44,925 |
control | 505 | 44,405 | 44,910 |
Total | 1,005 | 88,830 | 89,835 |
Deze tabel is dezelfde als tabel 6.2 op blz. 219 van Openintro Statistics.
Om de hypothese of er een verschil was in borstkankersterfte in de twee groepen met de hand te toetsen kunnen we schrijven:
# Definieer invoer (pas deze aan voor je eigen berekeningen)
p_hat_1 <- 500 / (500 + 44425)
p_hat_2 <- 505 / (505 + 44405)
n_1 <- 500 + 44425 # Totaal aantal deelnemers in de eerste groep
n_2 <- 505 + 44405 # Totaal aantal deelnemers in de tweede groep
null_value <- 0
# Voer berekeningen uit (je hoeft hier niets te veranderen)
p_hat_pooled <- (p_hat_1 * n_1 + p_hat_2 * n_2) / (n_1 + n_2)
point_est <- p_hat_1 - p_hat_2
se <- sqrt((p_hat_pooled*(1-p_hat_pooled)/n_1)+(p_hat_pooled*(1-p_hat_pooled)/n_2))
z <- (point_est - null_value)/se
pnorm(z)[1] 0.434892De oppervlakte onder de onderste staart is 0,4349 (klein verschil met het boek door afronding). De p-waarde, gevormd door de beide staarten samen, is \(0,434892 * 2 = 0,869784\). Omdat deze waarde \(p > 0,05\), verwerpen we de nulhypothese niet.
We kunnen dit als volgt visualiseren:
library(visualize)
visualize.norm(stat = c(-z, z), section = "tails")
Als alternatief kunnen we de R-functie prop.test() gebruiken om de significantie van het verschil tussen de twee groepen te testen. De functie prop.test maakt deel uit van het package stats, dat een van de weinige packages is die R automatisch laadt bij het opstarten.
We hebben prop.test al gebruikt voor een hypothesetoets bij één proportie, maar we kunnen deze functie ook gebruiken voor een toets met twee proporties:
result <- prop.test(x = c(500, 505),
n = c(44925, 44910),
alternative = "two.sided",
correct = FALSE)
resultprop.test(Dit is de functie die een test voor een proportie uitvoert.
x = c(500, 505),Hier geven wij het aantal gevallen aan dat zich in onze twee groepen bevindt. In dit voorbeeld: 500 patiënten die een mammografie kregen en stierven en 505 patiënten die geen mammografie kregen en stierven. Omdat dit twee waarden zijn, combineren we ze met c().
n = c(44925, 44910),De tweede opgegeven waarde moet het totale aantal mensen per studiegroep zijn. Omdat dit twee waarden zijn, combineren we ze met c().
alternative = "two.sided"Bepaal of we een tweezijdige of een eenzijdige test willen gebruiken. Opties zijn “two.sided” (standaard), “less” (als \(H_1: \mu < p\)) of “greater” (als \(H_1: \mu > p\)).
correct = FALSE)Hiermee geven we aan dat we het betrouwbaarheidsinterval niet willen berekenen met een ‘continuïteitscorrectie’. Standaard staat deze waarde op “TRUE”.
2-sample test for equality of proportions without continuity correction
data: c(500, 505) out of c(44925, 44910)
X-squared = 0.026874, df = 1, p-value = 0.8698
alternative hypothesis: two.sided
95 percent confidence interval:
-0.001490590 0.001260488
sample estimates:
prop 1 prop 2
0.01112966 0.01124471 We kunnen bovenstaande analyse ook uitvoeren op de gegevens uit het dataframe.
prop.test(table(mammogram$treatment, mammogram$breast_cancer_death),
alternative = "two.sided",
correct = FALSE)
2-sample test for equality of proportions without continuity correction
data: table(mammogram$treatment, mammogram$breast_cancer_death)
X-squared = 0.026874, df = 1, p-value = 0.8698
alternative hypothesis: two.sided
95 percent confidence interval:
-0.001490590 0.001260488
sample estimates:
prop 1 prop 2
0.01112966 0.01124471 prop.test(table(mammogram$treatment, mammogram$breast_cancer_death),Dit is de functie die een test voor een proportie uitvoert. We selecteren de onafhankelijke (mammogram$treatment) en afhankelijke variabele (mammogram$breast_cancer_death). Let op: de volgorde van de variabelen is van belang voor het betrouwbaarheidsinterval: vermeld dus eerst de onafhankelijke en dan de afhankelijke variabele.
alternative = "two.sided"Bepaal of we een tweezijdige of een eenzijdige test willen gebruiken. Opties zijn “two.sided” (standaard), “less” (als \(H_1: \mu < p\)) of “greater” (als \(H_1: \mu > p\)).
correct = FALSE)Hiermee geven we aan dat we het betrouwbaarheidsinterval niet willen berekenen met een ‘continuïteitscorrectie’. Standaard staat deze waarde op “TRUE”.
De uitkomst is exact hetzelfde als hierboven onder ‘samengevatte gegevens’.
Neem dus naast een inhoudelijke omschrijving in je rapportage op:
De twee percentages/proporties die worden vergeleken.
De z-waarde of Chikwadraatwaarde. Als je de laatste gebruikt ook het aantal vrijheidsgraden tussen haken, dat is voor dit soort berekeningen altijd 1. Aangezien R de chikwadraatwaarde rapporteert, is deze het meest logisch om op te nemen.
p = p-waarde. Maar: schrijf nooit \(p = 0,000\). Want de p-waarde is nooit precies nul, maar heel klein. Het is dan beter om te zeggen \(p < 0,001\). Als de hypothesetoets handmatig uitvoert schrijf je p < α-waarde, bijvoorbeeld \(p < 0,05\).
Deze gebruikt een iets ingewikkelder formule voor de berekening van het betrouwbaarheidsinterval (Wilson score-interval), waarin ook een “continuïteitscorrectie” is opgenomen (die rekening houdt met het feit dat we een continue normale verdeling gebruiken om een discreet fenomeen te benaderen, namelijk het aantal successen).↩︎
Als je echt van pipes houdt, kun je ook gebruiken: example data |> pull(party_choice) |> table().↩︎
Als je het betrouwbaarheidsinterval voor de tweede categorie wilt berekenen, kun je bijvoorbeeld relevel gebruiken om de categorieën van de eerste variabele opnieuw te rangschikken:
prop.test(table(relevel(example_data$party_choice, ref = "Other party")),
conf.level = 0.95)
1-sample proportions test with continuity correction
data: table(relevel(example_data$party_choice, ref = "Other party")), null probability 0.5
X-squared = 290.52, df = 1, p-value < 2.2e-16
alternative hypothesis: true p is not equal to 0.5
95 percent confidence interval:
0.7423954 0.7954986
sample estimates:
p
0.77 Dit gebruikt een iets ingewikkelder formule voor de berekening van het betrouwbaarheidsinterval (Wilson score-interval), waarin ook een “continuïteitscorrectie” is opgenomen (die rekening houdt met het feit dat we een continue normale verdeling gebruiken om een discreet fenomeen te benaderen, d.w.z. het aantal successen).↩︎