1. Inleiding: Waarom Statistiek?

1 Introductie

• We leven in het big data era
• Gegevens over locatie, surfgedrag, aankoopgedrag, sociale media
• Wetenschap: expressie meten van duizenden genen, eiwitten,… voor elk subject
• Chemische proces controle: groot aantal sensoren die continu een proces meten zodat het proces kan worden bijgestuurd
• Data driven journalism
• …

Statistiek is de wetenschap om te leren uit empirische gegevens

Statistische geletterdheid is cruciaal om resultaten en publicaties kritisch te kunnen interpreteren.

---

2 Case study: oksel microbiome

• https://www.vrt.be/vrtnws/nl/2018/10/22/gezocht-mensen-met-penetrante-lijfgeur-om-probiotische-deodor/
• https://youtu.be/9RIFyqLXdVw

---

• Okselgeur wordt niet veroorzaakt door het zweet zelf. De geur is afkomstig van specificiek micro-organismen van de groep Corynebacterium spp. die het zweet metaboliseren. Staphylococcus spp. zijn een andere groep bacteriën die ook abundant zijn in het microbiome van de oksel die zweet metaboliseren naar metabolieten die niet stinken.

• De CMET-groep aan UGENT doet onderzoek naar het transplanteren van de microbiële gemeenschap, het microbiome, om mensen van geurende oksels af te helpen.

• Voorgestelde therapie

  1. Verwijder het oksel microbiome met antibiotica
  2. Beïnvloed het oksel microbiome door microbiome te transplanteren van een individue die geen geurende oksels heeft (https://youtu.be/9RIFyqLXdVw)

---

2.1 Proefopzet (experimental design)

---

---

• Experiment:

  • 20 personen worden at random uit de populatie getrokken van personen met een okselgeur: steekproef representatief voor populatie!

  • de personen worden at random verdeeld over twee behandelingsgroepen:

    • placebo (enkel antibiotica)
    • transplantie (antibiotica en microbiële transplantatie).
    • Randomisatie is belangrijk om ervoor te zorgen dat de groepen vergelijkbaar zijn.

  • Het microbiome wordt bemonsterd 6 weken na de behandeling.

  • The relative abundantie van Staphylococcus spp. op Corynebacterium spp. + Staphylococcus spp. in het microbiome wordt gemeten via DGGE (Denaturing Gradient Gel Electrophoresis).

---

DGGE

https://doi.org/10.1371/journal.pone.0070538

---

Vertaal onderzoeksvraag naar iets wat we kunnen quantificeren: Is er een verschil in relatieve abundantie van Staphylococcus spp. in het microbiome van de placebo groep en de transplantatie

---

---

2.2 Data Exploratie en beschrijvende statistiek

• Data exploratie is heel belangrijk om inzicht te krijgen in de data en is een essentiële eerste stap om te leren uit data.
• Het wordt vaak ondergewaardeerd of over het hoofd gezien.

---

2.2.1 Importeer de data

• Data in deze cursus wordt verwerkt via het statistisch software pakket R.

• Dit pakket laat toe om te leren uit data.

• In deze cursus gaan we data eerst exploreren om inzicht te verwerven in de gegevens om die vervolgens statistisch te verwerken.

• Vooraleer we hiermee van start kunnen gaan moeten we de data eerst importeren in R.

• Via het volgende commando kunnen we enkele regels van een data bestand inlezen om de structuur van het data bestand te weten te komen.

read_lines("https://raw.githubusercontent.com/statOmics/sbc20/master/data/armpit.csv")

 [1] "trt,rel"                       "placebo,54.99207606973059"    
 [3] "placebo,31.84466019417476"     "placebo,41.09947643979057"    
 [5] "placebo,59.52063914780293"     "placebo,63.573407202216075"   
 [7] "placebo,41.48648648648649"     "placebo,30.44041450777202"    
 [9] "placebo,42.95676429567643"     "placebo,41.7391304347826"     
[11] "placebo,33.896515311510036"    "transplant,57.218124341412015"
[13] "transplant,72.50900360144058"  "transplant,61.89258312020461" 
[15] "transplant,56.690140845070424" "transplant,76"                
[17] "transplant,71.7357910906298"   "transplant,57.757296466973884"
[19] "transplant,65.1219512195122"   "transplant,67.53424657534246" 
[21] "transplant,77.55359394703657" 

• Gegevens in het bestand zijn door comma’s gescheiden.
• Elke rij bevat de gegevens voor 1 proefpersoon
• Verschillende variabelen worden gemeten per persoon en zijn van elkaar gescheiden door een comma. Het bestand is csv formaat: “comma separated values”.
• We kunnen bestanden met dit formaat inlezen R via het commando read_csv.
• We slaan de data op in R in het object met naam ap. Hiervoor gebruiken we de <- operator.
• We geven de data tabel terug door het object aan te roepen door zijn naam te typen.

ap <- read_csv("https://raw.githubusercontent.com/statOmics/sbc20/master/data/armpit.csv")
ap

# A tibble: 20 × 2
   trt          rel
   <chr>      <dbl>
 1 placebo     55.0
 2 placebo     31.8
 3 placebo     41.1
 4 placebo     59.5
 5 placebo     63.6
 6 placebo     41.5
 7 placebo     30.4
 8 placebo     43.0
 9 placebo     41.7
10 placebo     33.9
11 transplant  57.2
12 transplant  72.5
13 transplant  61.9
14 transplant  56.7
15 transplant  76  
16 transplant  71.7
17 transplant  57.8
18 transplant  65.1
19 transplant  67.5
20 transplant  77.6

---

2.2.2 Beschrijvende statistiek

• In artikels en de media worden resultaten uit een steekproef vaak gerapporteerd a.d.h.v. gemiddelde en de standaardafwijking

• We vatten de data eerst samen. We berekenen het gemiddelde en de standaard deviatie (een maat voor de spreiding, zie volgende hoofdstukken). We slaan het resultaat hiervan op in het object apRelSum via apRelSum <-.

1. We pipen (via %>%) het ap dataframe naar de group_by functie om de data te groeperen per treatment trt: group_by(trt).

2. We pipen het resultaat naar de summarize_at function om de “rel” variable samen te vatten en berekenen hierbij het gemiddelde en standaardafwijking. Omdat we de data eerst hebben gegroepeerd zullen we het gemiddelde en de standaard deviatie berekenen per groep.

apRelSum <- ap %>%
  group_by(trt) %>%
  summarize_at("rel",
               list(mean=mean,
                    sd=sd
                    )
                )

We tonen vervolgens het resultaat door het object apRelSum aan te roepen

apRelSum

# A tibble: 2 × 3
  trt         mean    sd
  <chr>      <dbl> <dbl>
1 placebo     44.2 11.5 
2 transplant  66.4  7.88

We kunnen ook een tabel in de webpagina of het pdf bestand integreren via het commando kable van het knitr pakket:

knitr::kable(apRelSum)

trt	mean	sd
placebo	44.15496	11.543251
transplant	66.40127	7.880175

---

2.2.3 Grafieken

We maken in deze cursus gebruik van het pakket ggplot2 om grafieken te maken.
Met de ggplot2 bibliotheek kunnen we gemakkelijk grafieken opbouwen in lagen (layers). Hierdoor leest de code veel makkelijker.

2.2.3.1 barplot

Bar plots worden heel veel gebruikt in artikels om resultaten weer te geven.

1. We pipen de samengevatte data naar de functie ggplot. Dat is de basis van elke ggplot. We selecteren de variabele met de behandeling trt als x variabele en de variabele met naam mean als y-variabele voor de plot. We doen dit steeds via de aestetics aes functie. aes(x=trt,y=mean)

2. We maken een barplot door een laag toe te voegen via de geom_bar function. De statistiek is stat="identity" omdat de hoogte van de bar gelijk is aan de waarde voor y (hier het gemiddelde voor de relatieve abundantie).

3. We voegen foutenvlaggen toe om de onzekerheid op het gemiddelde weer te geven. We doen dit via de geom_errorbar functie en specifiëren het minimum en maximum van de error bar. Het width argument wordt gebruikt om de breedte van de error bar smaller te maken dat deze van de bar.

apRelSum %>%
  ggplot(aes(x=trt,y=mean)) +
  geom_bar(stat="identity") +
  geom_errorbar(aes(ymin=mean-sd,ymax=mean+sd),width=.2)

• Is deze plot informatief??

---

2.2.3.2 boxplots

Barplots zijn geen goede grafieken:

• Ze zijn niet informatief

• Ze visulariseren een samenvatting van de data door twee punten: two point summary: gemiddelde en de standaard deviatie hierop en kunnen beter in een tabel worden opgenomen.

• Ze laten niet toe om anomalieën in de data zoals meetfouten op te sporen.

• Ze verbruiken veel ruimte, b.v. van nul tot de kleinste waarde voor de relatieve abundantie) waar geen data in ligt.

Het is beter om de data zo ruw mogelijk voor te stellen zodat we inzicht krijgen in de verdeling van de gegevens. Hiervoor zullen we ondermeer boxplots gebruiken.

---

We maken nu een boxplot voor de ap data

1. We pipen het ap dataframe naar ggplot
2. We selecteren de data voor de plot via ggplot(aes(x=trt,y=rel))
3. We voegen laag toe voor de boxplot dmv de functie geom_boxplot()

ap %>%  
  ggplot(aes(x=trt,y=rel)) +
  geom_boxplot()

---

Merk op dat we de plot ook op kunnen slaan als een object.

apBoxplot <- ap %>%  
  ggplot(aes(x=trt,y=rel)) +
  geom_boxplot()

De plot wordt dan niet gemaakt.

Om de plot weer te geven kunnen we het object aanroepen

apBoxplot

---

• Merk op dat we niet zoveel gegevens hebben. Enkel 10 per groep.

• Het is altijd beter om de data zo ruw mogelijk te tonen!

Omdat er niet zoveel gegevens zijn kunnen we de data toevoegen aan de plot zonder dat die te druk wordt.

• Merk op dat we het argument outlier.shape op NA (not available) zetten outlier.shape=NA in the geom_boxplot functie omdat we anders outliers twee keer weer zullen geven. Eerst via de boxplot laag en daarna omdat we een laag met alle ruwe data toevoegen aan de plot.
• We geven de ruwe data weer via de geom_point(position="jitter") functie. We gebruiken hierbij het argument position=‘jitter’ zodat we wat random ruis toevoegen aan de x-cordinaat zodat de gegevens elkaar niet overlappen.

ap %>%  
  ggplot(aes(x=trt,y=rel)) +
  geom_boxplot(outlier.shape=NA) +
  geom_point(position="jitter")

Dit is een informatieve plot!

---

Aangezien we de plot ook hebben opgeslagen konden we dit ook doen door de plot op te roepen en nog een laag toe te voegen.

apBoxplot +
  geom_point(position="jitter")

We hadden het resultaat ook opnieuw op kunnen slaan voor later hergebruik.

apBoxplot <- apBoxplot +
  geom_point(position="jitter")

---

• We zagen duidelijk een effect van de transplantatie op de relatieve abundantie van Staphylococcus.

• Is dat effect nu groot genoeg om te kunnen concluderen dat de behandeling werkt?

---

---

• Inductie: Door middel van statistische besluitvorming (inference) kunnen we uitspraken doen over de populatie op basis van een steekproef.

• De prijs die we hiervoor betalen is onzekerheid!

• We kunnen op basis van een steekproef nooit absoluut zeker zijn van onze conclusies.

---

• Met data kunnen we niet bewijzen dat een behandeling werkt.

• Falsificatie principe van Popper: Data kunnen enkel een hypothese of een theorie ontkrachten.

• Met statistiek kunnen we dus niet aantonen dat de behandeling werkt.

• Statistiek zal ons wel toelaten om het omgekeerde te falsifiëren: als we veronderstellen dat er geen effect van de behandeling, spreekt de data in de steekproef dit tegen?

• Met statistiek kunnen we berekenen hoe waarschijnlijk het is om in een random steekproef (wanneer we het experiment dus opnieuw uit zouden voeren) een gemiddeld verschil in relatieve abundantie te zien tussen placebo en transplantatiegroep dat minstens zo groot is als in onze steekproef als de behandeling geen effect zou hebben.

• Die kans wordt een p-waarde genoemd.

• Als p heel klein is, dan is het heel onwaarschijnlijk om een dergelijk effect te observeren in een steekproef door toeval.

• p wordt meestal vergeleken met 5%. Als er geen effect is van de behandeling dan tolereren 5% vals positieve conclusies.

• Om de kans p te berekenen is het nodig om de data te modelleren met statistische modellen.

---

In latere hoofdstukken zullen we zien dat we t-test kunnen gebruiken om hetgeen we observeren in de microbiome dataset te veralgemenen naar de populatie.

t.test(rel~trt,data=ap)

    Welch Two Sample t-test

data:  rel by trt
t = -5.0334, df = 15.892, p-value = 0.0001249
alternative hypothesis: true difference in means between group placebo and group transplant is not equal to 0
95 percent confidence interval:
 -31.62100 -12.87163
sample estimates:
   mean in group placebo mean in group transplant 
                44.15496                 66.40127 

Conclusie: Gemiddeld is de relatieve abundantie van Staphylococcus in het microbiome van personen met een zweetgeur 22.2% hoger na de transplantie dan na de placebo behandeling.

---

2.3 Randomisatie

• Wat wordt zijn de consequenties van het gebruik van een steekproef en van randomisatie?

• Randomisatie is sterk gerelateerd met het concept van de populatie en scope van de studie.

• De scope van de studie moet goed worden omschreven voor de start van het experiment.

• Omdat de statistische analyse valide zou zijn is het noodzakelijk dat de subjecten volledig random worden getrokken uit de populatie naar waar we onze conclusies wensen te veralgemenen.

• Volledig random trekken van de populatie impliceert dat:

  • alle subjecten van de populatie evenveel kans hebben om in de steekproef te worden opgenomen
  • de selectie van een subject onafhankelijk is van de andere subjecten in de steekproef.

• De steekproef is dan representatief voor de populatie, maar is nog steeds random.

• Wat betekent dit?

---

3 Case study: Lengte van mannen en vrouwen - Variabiliteit van steekproef tot steekproef

• Om te begrijpen dat een steekproef random is zouden we hetzelfde experiment veel keer moeten kunnen herhalen (repeated sampling).

• Dan zouden we inzicht kunnen krijgen hoe de gegevens veranderen van steekproef tot steekproef.

• Om dit te illustreren zullen we gebruik maken van een hele grote studie.

• Uit die studie zullen we dan herhaaldelijk kleine steekproeven trekken om te begrijpen hoe de gegevens en statistieken veranderen van steekproef tot steekproef. Of om met andere woorden na te gaan wat de variabiliteit is tussen steekproeven.

---

National Health And Nutrition Examination Study (NHANES)

• Sinds 1960 worden elk jaar mensen van alle leeftijden geïnterviewd bij hen thuis.
• Er maakt ook een gezondheidsonderzoek deel uit van de study die in een mobiel onderzoekscentrum wordt afgenomen.
• We zullen deze grote studie gebruiken om at random personen te selecteren van de Amerikaanse populatie.
• Dat zal inzicht geven in hoe de gegevens en resultaten van een analyse zullen variëren van steekproef tot steekproef.
• De data van deze studie is terug te vinden in het R pakket NHANES

---

library(NHANES)
head(NHANES)

# A tibble: 6 × 76
     ID SurveyYr Gender   Age AgeDecade AgeMonths Race1 Race3 Education   
  <int> <fct>    <fct>  <int> <fct>         <int> <fct> <fct> <fct>       
1 51624 2009_10  male      34 " 30-39"        409 White <NA>  High School 
2 51624 2009_10  male      34 " 30-39"        409 White <NA>  High School 
3 51624 2009_10  male      34 " 30-39"        409 White <NA>  High School 
4 51625 2009_10  male       4 " 0-9"           49 Other <NA>  <NA>        
5 51630 2009_10  female    49 " 40-49"        596 White <NA>  Some College
6 51638 2009_10  male       9 " 0-9"          115 White <NA>  <NA>        
# … with 67 more variables: MaritalStatus <fct>, HHIncome <fct>,
#   HHIncomeMid <int>, Poverty <dbl>, HomeRooms <int>, HomeOwn <fct>,
#   Work <fct>, Weight <dbl>, Length <dbl>, HeadCirc <dbl>, Height <dbl>,
#   BMI <dbl>, BMICatUnder20yrs <fct>, BMI_WHO <fct>, Pulse <int>,
#   BPSysAve <int>, BPDiaAve <int>, BPSys1 <int>, BPDia1 <int>, BPSys2 <int>,
#   BPDia2 <int>, BPSys3 <int>, BPDia3 <int>, Testosterone <dbl>,
#   DirectChol <dbl>, TotChol <dbl>, UrineVol1 <int>, UrineFlow1 <dbl>, …

glimpse(NHANES)

Rows: 10,000
Columns: 76
$ ID               <int> 51624, 51624, 51624, 51625, 51630, 51638, 51646, 5164…
$ SurveyYr         <fct> 2009_10, 2009_10, 2009_10, 2009_10, 2009_10, 2009_10,…
$ Gender           <fct> male, male, male, male, female, male, male, female, f…
$ Age              <int> 34, 34, 34, 4, 49, 9, 8, 45, 45, 45, 66, 58, 54, 10, …
$ AgeDecade        <fct>  30-39,  30-39,  30-39,  0-9,  40-49,  0-9,  0-9,  40…
$ AgeMonths        <int> 409, 409, 409, 49, 596, 115, 101, 541, 541, 541, 795,…
$ Race1            <fct> White, White, White, Other, White, White, White, Whit…
$ Race3            <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ Education        <fct> High School, High School, High School, NA, Some Colle…
$ MaritalStatus    <fct> Married, Married, Married, NA, LivePartner, NA, NA, M…
$ HHIncome         <fct> 25000-34999, 25000-34999, 25000-34999, 20000-24999, 3…
$ HHIncomeMid      <int> 30000, 30000, 30000, 22500, 40000, 87500, 60000, 8750…
$ Poverty          <dbl> 1.36, 1.36, 1.36, 1.07, 1.91, 1.84, 2.33, 5.00, 5.00,…
$ HomeRooms        <int> 6, 6, 6, 9, 5, 6, 7, 6, 6, 6, 5, 10, 6, 10, 10, 4, 3,…
$ HomeOwn          <fct> Own, Own, Own, Own, Rent, Rent, Own, Own, Own, Own, O…
$ Work             <fct> NotWorking, NotWorking, NotWorking, NA, NotWorking, N…
$ Weight           <dbl> 87.4, 87.4, 87.4, 17.0, 86.7, 29.8, 35.2, 75.7, 75.7,…
$ Length           <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ HeadCirc         <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ Height           <dbl> 164.7, 164.7, 164.7, 105.4, 168.4, 133.1, 130.6, 166.…
$ BMI              <dbl> 32.22, 32.22, 32.22, 15.30, 30.57, 16.82, 20.64, 27.2…
$ BMICatUnder20yrs <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ BMI_WHO          <fct> 30.0_plus, 30.0_plus, 30.0_plus, 12.0_18.5, 30.0_plus…
$ Pulse            <int> 70, 70, 70, NA, 86, 82, 72, 62, 62, 62, 60, 62, 76, 8…
$ BPSysAve         <int> 113, 113, 113, NA, 112, 86, 107, 118, 118, 118, 111, …
$ BPDiaAve         <int> 85, 85, 85, NA, 75, 47, 37, 64, 64, 64, 63, 74, 85, 6…
$ BPSys1           <int> 114, 114, 114, NA, 118, 84, 114, 106, 106, 106, 124, …
$ BPDia1           <int> 88, 88, 88, NA, 82, 50, 46, 62, 62, 62, 64, 76, 86, 6…
$ BPSys2           <int> 114, 114, 114, NA, 108, 84, 108, 118, 118, 118, 108, …
$ BPDia2           <int> 88, 88, 88, NA, 74, 50, 36, 68, 68, 68, 62, 72, 88, 6…
$ BPSys3           <int> 112, 112, 112, NA, 116, 88, 106, 118, 118, 118, 114, …
$ BPDia3           <int> 82, 82, 82, NA, 76, 44, 38, 60, 60, 60, 64, 76, 82, 7…
$ Testosterone     <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ DirectChol       <dbl> 1.29, 1.29, 1.29, NA, 1.16, 1.34, 1.55, 2.12, 2.12, 2…
$ TotChol          <dbl> 3.49, 3.49, 3.49, NA, 6.70, 4.86, 4.09, 5.82, 5.82, 5…
$ UrineVol1        <int> 352, 352, 352, NA, 77, 123, 238, 106, 106, 106, 113, …
$ UrineFlow1       <dbl> NA, NA, NA, NA, 0.094, 1.538, 1.322, 1.116, 1.116, 1.…
$ UrineVol2        <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ UrineFlow2       <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ Diabetes         <fct> No, No, No, No, No, No, No, No, No, No, No, No, No, N…
$ DiabetesAge      <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ HealthGen        <fct> Good, Good, Good, NA, Good, NA, NA, Vgood, Vgood, Vgo…
$ DaysPhysHlthBad  <int> 0, 0, 0, NA, 0, NA, NA, 0, 0, 0, 10, 0, 4, NA, NA, 0,…
$ DaysMentHlthBad  <int> 15, 15, 15, NA, 10, NA, NA, 3, 3, 3, 0, 0, 0, NA, NA,…
$ LittleInterest   <fct> Most, Most, Most, NA, Several, NA, NA, None, None, No…
$ Depressed        <fct> Several, Several, Several, NA, Several, NA, NA, None,…
$ nPregnancies     <int> NA, NA, NA, NA, 2, NA, NA, 1, 1, 1, NA, NA, NA, NA, N…
$ nBabies          <int> NA, NA, NA, NA, 2, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ Age1stBaby       <int> NA, NA, NA, NA, 27, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ SleepHrsNight    <int> 4, 4, 4, NA, 8, NA, NA, 8, 8, 8, 7, 5, 4, NA, 5, 7, N…
$ SleepTrouble     <fct> Yes, Yes, Yes, NA, Yes, NA, NA, No, No, No, No, No, Y…
$ PhysActive       <fct> No, No, No, NA, No, NA, NA, Yes, Yes, Yes, Yes, Yes, …
$ PhysActiveDays   <int> NA, NA, NA, NA, NA, NA, NA, 5, 5, 5, 7, 5, 1, NA, 2, …
$ TVHrsDay         <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ CompHrsDay       <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ TVHrsDayChild    <int> NA, NA, NA, 4, NA, 5, 1, NA, NA, NA, NA, NA, NA, 4, N…
$ CompHrsDayChild  <int> NA, NA, NA, 1, NA, 0, 6, NA, NA, NA, NA, NA, NA, 3, N…
$ Alcohol12PlusYr  <fct> Yes, Yes, Yes, NA, Yes, NA, NA, Yes, Yes, Yes, Yes, Y…
$ AlcoholDay       <int> NA, NA, NA, NA, 2, NA, NA, 3, 3, 3, 1, 2, 6, NA, NA, …
$ AlcoholYear      <int> 0, 0, 0, NA, 20, NA, NA, 52, 52, 52, 100, 104, 364, N…
$ SmokeNow         <fct> No, No, No, NA, Yes, NA, NA, NA, NA, NA, No, NA, NA, …
$ Smoke100         <fct> Yes, Yes, Yes, NA, Yes, NA, NA, No, No, No, Yes, No, …
$ Smoke100n        <fct> Smoker, Smoker, Smoker, NA, Smoker, NA, NA, Non-Smoke…
$ SmokeAge         <int> 18, 18, 18, NA, 38, NA, NA, NA, NA, NA, 13, NA, NA, N…
$ Marijuana        <fct> Yes, Yes, Yes, NA, Yes, NA, NA, Yes, Yes, Yes, NA, Ye…
$ AgeFirstMarij    <int> 17, 17, 17, NA, 18, NA, NA, 13, 13, 13, NA, 19, 15, N…
$ RegularMarij     <fct> No, No, No, NA, No, NA, NA, No, No, No, NA, Yes, Yes,…
$ AgeRegMarij      <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 20, 15, N…
$ HardDrugs        <fct> Yes, Yes, Yes, NA, Yes, NA, NA, No, No, No, No, Yes, …
$ SexEver          <fct> Yes, Yes, Yes, NA, Yes, NA, NA, Yes, Yes, Yes, Yes, Y…
$ SexAge           <int> 16, 16, 16, NA, 12, NA, NA, 13, 13, 13, 17, 22, 12, N…
$ SexNumPartnLife  <int> 8, 8, 8, NA, 10, NA, NA, 20, 20, 20, 15, 7, 100, NA, …
$ SexNumPartYear   <int> 1, 1, 1, NA, 1, NA, NA, 0, 0, 0, NA, 1, 1, NA, NA, 1,…
$ SameSex          <fct> No, No, No, NA, Yes, NA, NA, Yes, Yes, Yes, No, No, N…
$ SexOrientation   <fct> Heterosexual, Heterosexual, Heterosexual, NA, Heteros…
$ PregnantNow      <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…

---

3.1 Data exploratie

Onderzoeksvraag: hoe verschilt de lengte van volwassen mannen en vrouwen.

1. We pipen de dataset naar de function filter om de data te filteren volgens leeftijd.
   
2. We plotten de lengte metingen.
   • We selecteren de data data met het commando ggplot(aes(x=lengte))
   • We voegen een histogram toe met het commando geom_histogram()
   • We maken twee vertikale panels met het commando facet_grid(Gender~.)
   • We veranderen het label van de x-as met de xlab functie.

NHANES%>%  
  filter(Age >= 18 & !is.na(Height)) %>%
  ggplot(aes(x = Height))+
  geom_histogram() +
  facet_grid(Gender ~ .) +
  xlab("Lengte (cm)")

---

We zien dat de data nu min of meer symmetrisch verdeeld zijn en een klokvorm hebben.
Dat zal ons toe laten om de data verder samen te vatten door gebruik te maken van twee statistieken: het gemiddelde en de standaard deviatie wat een maat is voor de spreiding van de gegevens rond het gemiddelde.

---

We maken nu een subset van de data die we zullen gebruiken om aan te tonen hoe de variabiliteit in kleine steekproeven kan variëren van steekproef tot steekproef.

1. We filteren op leeftijd en verwijderen ontbrekenden gegevens (NA, Not Available).
2. We selecteren enkel het geslacht en Lengte zodat de dataset geen onnodige variabelen bevat.

nhanesSub <- NHANES %>%
  filter(Age >= 18 & !is.na(Height)) %>%
  select(c("Gender","Height"))

---

We berekenen het gemiddelde en de standaard deviatie voor de lengte voor mannen en vrouwen in de grote dataset. We groeperen de data hiervoor op basis van het geslacht (variable Gender).

HeightSum <- nhanesSub %>%
  group_by(Gender) %>%
  summarize_at("Height",
               list(mean = mean,
               sd = sd)
              )

knitr::kable(
  HeightSum %>%
  mutate_if(is.numeric, round, digits=1)
  )

Gender	mean	sd
female	162.1	7.3
male	175.9	7.5

---

3.2 Experiment

• Stel dat we geen toegang hebben tot de metingen van de NHANES studie.

• We zouden dan een experiment op moeten zetten om metingen bij mannen en vrouwen te doen.

• Veronderstel dat we budget hebben om metingen bij 5 mannen en 5 vrouwen te doen.

• We zouden dan 5 mannen en 5 vrouwen boven de 18 jaar at random selecteren uit de Amerikaanse populatie.

• We kunnen dit experiment simuleren door 5 vrouwen en 5 mannen at random te selecteren uit de NHANES studie.

---

set.seed(1023)
nSamp <- 5
fem <- nhanesSub %>%
  filter(Gender=="female") %>%
  sample_n(size=5)

mal <- nhanesSub %>%
  filter(Gender=="male") %>%
  sample_n(size=5)

samp1 <- rbind(fem,mal)

samp1

# A tibble: 10 × 2
   Gender Height
   <fct>   <dbl>
 1 female   164 
 2 female   160.
 3 female   159 
 4 female   154.
 5 female   156.
 6 male     170.
 7 male     183.
 8 male     183.
 9 male     185.
10 male     170.

---

Data Exploratie

samp1 %>%
  ggplot(aes(x=Height)) +
  geom_histogram() +
  facet_grid(Gender~.) +
  xlab("Lengte (cm)")

HeightSumExp1 <- samp1 %>%
  group_by(Gender) %>%
  summarize_at("Height",
               list(mean = mean,
                      sd = sd)
                  )
HeightSumExp1

# A tibble: 2 × 3
  Gender  mean    sd
  <fct>  <dbl> <dbl>
1 female  159.  3.76
2 male    178.  7.55

Histogram is niet zinvol als we maar zo weinig datapunten hebben.

---

Boxplot is beter:

---

We voeren hier ook een t-test uit.

t.test(Height~Gender,data=samp1)

    Welch Two Sample t-test

data:  Height by Gender
t = -5.0783, df = 5.8695, p-value = 0.00242
alternative hypothesis: true difference in means between group female and group male is not equal to 0
95 percent confidence interval:
 -28.441927  -9.878073
sample estimates:
mean in group female   mean in group male 
              158.72               177.88 

In het experiment zijn vrouwen gemiddeld 19.16 cm kleiner dan mannen. En als we een statistische test uitvoeren (zie hoofdstuk 5: Statistische besluitvorming) kunnen we besluiten dat dit verschil statistisch significant is.

---

3.3 Herhaal het experiment

Als we het experiment herhalen selecteren we andere mensen en verkrijgen we andere resultaten.

set.seed(1024)
fem <- nhanesSub %>%
  filter(Gender=="female") %>%
  sample_n(size=5)

mal <- nhanesSub %>%
  filter(Gender=="male") %>%
  sample_n(size=5)

samp2 <- rbind(fem,mal)

HeightSumExp2 <- samp2 %>%
  group_by(Gender) %>%
  summarize_at("Height",
               list(mean=mean,
                    sd=sd)
              )
HeightSumExp2

# A tibble: 2 × 3
  Gender  mean    sd
  <fct>  <dbl> <dbl>
1 female  166.  9.89
2 male    170.  8.24

samp2 %>%
  ggplot(aes(x = Gender,y = Height)) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = "jitter") +
  geom_point(
    aes(x = 1, y = HeightSumExp2$mean[1]),
    size = 3,
    pch = 17,
    color="darkred") +
  geom_point(
    aes(x = 2, y = HeightSumExp2$mean[2]),
    size = 3,
    pch = 17,
    color = "darkred") +
  ylab("Height (cm)")

t.test(Height ~ Gender, data=samp2)

    Welch Two Sample t-test

data:  Height by Gender
t = -0.7295, df = 7.747, p-value = 0.4872
alternative hypothesis: true difference in means between group female and group male is not equal to 0
95 percent confidence interval:
 -17.552343   9.152343
sample estimates:
mean in group female   mean in group male 
              166.26               170.46 

In de nieuwe steekproef zijn vrouwen gemiddeld 4.2 cm kleiner dan mannen. En dit verschil is statistisch niet significant

---

3.4 Herhaal het experiment opnieuw

seed <- 88605
set.seed(seed)
fem <- nhanesSub %>%
  filter(Gender=="female") %>%
  sample_n(size=5)

mal <- nhanesSub %>%
  filter(Gender=="male") %>%
  sample_n(size=5)

samp3 <- rbind(fem,mal)

HeightSumExp3 <- samp3 %>%
  group_by(Gender) %>%
  summarize_at("Height",
               list(mean=mean,
                    sd=sd)
              )
HeightSumExp3

# A tibble: 2 × 3
  Gender  mean    sd
  <fct>  <dbl> <dbl>
1 female  173.  1.97
2 male    168.  2.84

samp3 %>%
  ggplot(aes(x = Gender,y = Height)) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = "jitter") +
  geom_point(
    aes(x = 1, y = HeightSumExp3$mean[1]),
    size = 3,
    pch = 17,
    color="darkred") +
  geom_point(
    aes(x = 2, y = HeightSumExp3$mean[2]),
    size = 3,
    pch = 17,
    color = "darkred") +
  ylab("Height (cm)")

t.test(Height ~ Gender, data=samp3)

    Welch Two Sample t-test

data:  Height by Gender
t = 3.1182, df = 7.136, p-value = 0.01648
alternative hypothesis: true difference in means between group female and group male is not equal to 0
95 percent confidence interval:
 1.178916 8.461084
sample estimates:
mean in group female   mean in group male 
              172.60               167.78 

In de nieuwe steekproef zijn vrouwen gemiddeld 4.82 cm groter dan mannen. En dit verschil is statistisch significant

---

3.5 Samenvatting

• We trokken at random andere proefpersonen in elke steekproef

• Hierdoor verschillen lengtemetingen van steekproef tot steekproef.

• Dus ook de geschatte gemiddeldes en standaard deviaties.

• Bijgevolg zijn onze conclusies ook onzeker en kunnen deze wijzigen van steekproef tot steekproef.

• Voor het lengte voorbeeld zijn steekproeven waarbij het effect tegengesteld is aan dat in de populatie en waarbij we besluiten dat het verschil significant echter zeldzaam.

\(\rightarrow\) Met statistiek controleren we de kans op het trekken foute conclusies.

3.6 Controle van fouten

Hoe controleert statistiek de kans op het trekken van foute conclusies?

• In onderstaande code trekken we 10000 herhaalde steekproeven van 5 vrouwen en 5 mannen uit de NHANES studie.

set.seed(15152)
# Aantal simulaties en steekproefgrootte per groep
nSim <- 10000
nSamp <- 5

# We filteren de data vooraf zodat we dit niet telkens opnieuw hoeven te doen
fem <- nhanesSub %>%
  filter(Gender == "female")

mal <- nhanesSub %>%
  filter(Gender == "male")

# Simulatie studie
# Om snelle functies te kunnen gebruiken nemen we eerst nSim steekproeven en berekenen we daarna alles.

femSamps <- malSamps <-matrix(NA, nrow=nSamp, ncol=nSim)
for (i in 1:nSim)
{
  femSamps[,i] <- sample(fem$Height, nSamp)
  malSamps[,i] <- sample(mal$Height, nSamp)
}

res <- data.frame(
  verschil=colMeans(femSamps) - colMeans(malSamps),
  Rfast::ttests(femSamps, malSamps)
  )

sum(res$pvalue < 0.05 & res$verschil < 0)

[1] 7234

sum(res$pvalue >= 0.05)

[1] 2766

sum(res$pvalue < 0.05 & res$verschil>0)

[1] 0

res %>%
  ggplot(aes(x=verschil,y=-log10(pvalue),color=pvalue < 0.05)) +
  geom_point() +
  xlab("Gemiddeld Verschil (cm)") +
  ylab("Statistische Significantie (-log10 p)")

res %>%
  ggplot(aes(y = verschil)) +
  geom_boxplot() +
  ylab("Gemiddeld Verschil (cm)")

  xlab("")

$x
[1] ""

attr(,"class")
[1] "labels"

Op basis van 10000 steekproeven van 5 mannen en 5 vrouwen zagen we dat in 7234 steekproeven vrouwen gemiddeld significant kleiner zijn dan mannen. In 2766 steekproeven besluiten we dat vrouwen en mannen gemiddeld niet significant verschillen in lengte. En in 0 besluiten we dat vrouwen gemiddeld significant groter zijn dan mannen.

• De steekproef die we toonden waaruit we zouden besluiten dat vrouwen significant groter zijn dan mannen is heel onwaarschijnlijk. Er moesten 88605 steekproeven worden getrokken om deze extreme steekproef te vinden.

Het feit dat we in veel steekproeven resultaten vinden die statistisch niet significant zijn komt omdat de statistische toets een te lage kracht heeft om het verschil te detecteren wanneer er maar 5 observaties zijn per groep.

---

3.6.1 Grotere steekproef?

Wat gebeurt er als we de steekproef verhogen naar 50 observaties per groep?

set.seed(11145)
# Aantal simulaties en steekproefgrootte per groep
nSim <- 10000
nSamp <- 50

# We filteren de data vooraf zodat we dit niet telkens opnieuw hoeven te doen
fem <- nhanesSub %>%
  filter(Gender == "female")

mal <- nhanesSub %>%
  filter(Gender == "male")

# Simulatie studie
# Om snelle functies te kunnen gebruiken nemen we eerst nSim steekproeven en berekenen we daarna alles.

femSamps <- malSamps <- matrix(NA, nrow = nSamp, ncol = nSim)
for (i in 1:nSim)
{
  femSamps[,i] <- sample(fem$Height, nSamp)
  malSamps[,i] <- sample(mal$Height, nSamp)
}

res <- data.frame(
  verschil = colMeans(femSamps) - colMeans(malSamps),
  Rfast::ttests(femSamps, malSamps)
  )

sum(res$pvalue < 0.05 & res$verschil < 0)

[1] 10000

sum(res$pvalue >= 0.05)

[1] 0

sum(res$pvalue < 0.05 & res$verschil > 0)

[1] 0

res %>%
  ggplot(aes(x=verschil,y=-log10(pvalue),color=pvalue<0.05)) +
  geom_point() +
  xlab("Gemiddeld Verschil (cm)") +
  ylab("Statistische Significantie (-log10 p)")

res %>%
  ggplot(aes(y=verschil)) +
  geom_boxplot() +
  ylab("Gemiddeld Verschil (cm)")

  xlab("")

$x
[1] ""

attr(,"class")
[1] "labels"

• We zien dus dat we de kans om een verschil te vinden als er in werkelijkheid een verschil is in de populatie kunnen beïnvloeden in de design fase: aan de hand van de steekproefgrootte.

• Hoe meer gegevens hoe makkelijker we het werkelijk verschil oppikken in de steekproef.

3.6.2 Controle van vals positieven

Wat gebeurt er als er geen verschil is tussen beide groepen?

• We moeten hiervoor experimenten simuleren waarbij de groepen gelijk zijn.

• Hiervoor zullen we twee groepen vergelijken waarvoor de lengte gemiddeld niet verschillend is.

• Dat kunnen we doen door een steekproef te trekken waarbij we voor beide groepen at random subjecten trekken uit de subset van vrouwen in de NHANES studie.

• We doen dit opnieuw voor een steekproef met 5 subjecten per groep

set.seed(13245)
# Aantal simulaties en steekproefgrootte per groep
nSim <- 10000
nSamp <- 5

# We filteren de data vooraf zodat we dit niet telkens opnieuw hoeven te doen
fem <- nhanesSub %>%
  filter(Gender == "female")

# Simulatie studie
# Om snelle functies te kunnen gebruiken nemen we eerst nSim steekproeven en berekenen we daarna alles.

femSamps <- femSamps2 <-matrix(NA, nrow=nSamp, ncol=nSim)
for (i in 1:nSim)
{
  femSamps[,i] <- sample(fem$Height, nSamp)
  femSamps2[,i] <- sample(fem$Height, nSamp)
}

res <- data.frame(
  verschil=colMeans(femSamps) - colMeans(femSamps2),
  Rfast::ttests(femSamps, femSamps2)
  )

sum(res$pvalue < 0.05 & res$verschil < 0)

[1] 213

sum(res$pvalue >= 0.05)

[1] 9558

sum(res$pvalue < 0.05 & res$verschil>0)

[1] 229

res %>%
  ggplot(aes(x=verschil,y=-log10(pvalue),color=pvalue < 0.05)) +
  geom_point() +
  xlab("Gemiddeld Verschil (cm)") +
  ylab("Statistische Significantie (-log10 p)")

res %>%
  ggplot(aes(y = verschil)) +
  geom_boxplot() +
  ylab("Gemiddeld Verschil (cm)")

  xlab("")

$x
[1] ""

attr(,"class")
[1] "labels"

Op basis van 10000 steekproeven zien we dat we in 442 steekproeven ten onrechte besluiten dat er een verschil is in gemiddelde lengte tussen twee groepen vrouwen.

Met de statistische analyse controleren we dus het aantal vals positieve resultaten correct op 5%.

Wat gebeurt er als we het aantal observaties verhogen?

We simuleren opnieuw experimenten met 50 subjecten per groep maar we trekken de subjecten opnieuw telkens uit de populatie van vrouwen.

set.seed(1345)
# Aantal simulaties en steekproefgrootte per groep
nSim <- 10000
nSamp <- 50

# We filteren de data vooraf zodat we dit niet telkens opnieuw hoeven te doen
fem <- nhanesSub %>%
  filter(Gender == "female")

# Simulatie studie
# Om snelle functies te kunnen gebruiken nemen we eerst nSim steekproeven en berekenen we daarna alles.

femSamps <- femSamps2 <-matrix(NA, nrow=nSamp, ncol=nSim)
for (i in 1:nSim)
{
  femSamps[,i] <- sample(fem$Height, nSamp)
  femSamps2[,i] <- sample(fem$Height, nSamp)
}

res <- data.frame(
  verschil=colMeans(femSamps) - colMeans(femSamps2),
  Rfast::ttests(femSamps, femSamps2)
  )

sum(res$pvalue < 0.05 & res$verschil < 0)

[1] 271

sum(res$pvalue >= 0.05)

[1] 9501

sum(res$pvalue < 0.05 & res$verschil>0)

[1] 228

res %>%
  ggplot(aes(x=verschil,y=-log10(pvalue),color=pvalue < 0.05)) +
  geom_point() +
  xlab("Gemiddeld Verschil (cm)") +
  ylab("Statistische Significantie (-log10 p)")

res %>%
  ggplot(aes(y = verschil)) +
  geom_boxplot() +
  ylab("Gemiddeld Verschil (cm)")

  xlab("")

$x
[1] ""

attr(,"class")
[1] "labels"

Op basis van 10000 steekproeven zien we dat we in 499 steekproeven ten onrechte besluiten dat er een verschil is in gemiddelde lengte tussen twee groepen vrouwen.

Met de statistische analyse controleren we dus ook bij het nemen van een grote steekproef het aantal vals positieve resultaten correct op 5% (Vals positief: Op basis van de steekproef besluiten dat er gemiddeld een verschil is in lengte tussen beide groepen terwijl er in werkelijkheid geen verschil is in de populatie.)

3.6.3 Conclusies

• De statistische analyse controleert steeds de kans op het nemen van een vals positieve conclusie.

• De statistische analysis controleert de kans op het nemen van een vals negatieve conclusies niet, maar in de design fase kunnen we deze kans beïnvloeden ondermeer d.m.v. de steekproefgrootte.

---

#Case study: Salk vaccin

• In 1916, brak de eerste grote polio epidemie uit in de USA.
• Begin de jaren 50 ontwikkelde John Salk een vaccin met belovende resultaten in het lab.
• In 1954, heeft de National Foundation for Infantile Paralysis (NFIP) een grote studie opgezet om de effectiviteit van het Salk vaccin na te gaan.
• Veronderstel dat de NFIP in 1954 een groot aantal kinderen zou hebben gevaccineerd, wat zouden ze dan kunnen besluiten als de polio incidentie in 1954 lager was dan in 1953?

---

3.7 NFIP Study

3.7.1 Design

• Grote simultane studie met gevaccineerde kinderen (cases) en ongevaccineerde kinderen (controles).
• In scholen van districten met hoge polio incidentie.
• Cases: kinderen van de tweede graad van het lager onderwijs waarvan de ouders toestemden met vaccinatie.
• Controles: kinderen van de eerste en derde graad.

3.7.2 Data

nfip <- tibble(
  group=c("cases","controls","noConcent"),
  grade=c("g2","g1g3","g2"),
  vaccin=c("yes","no","no"),
  total=c(221998,725173,123605),
  polio=c(54,391,56)
  ) %>%
  mutate(noPolio = total - polio)
knitr::kable(nfip)

group	grade	vaccin	total	polio	noPolio
cases	g2	yes	221998	54	221944
controls	g1g3	no	725173	391	724782
noConcent	g2	no	123605	56	123549

Vergelijk de polio incidentie?

---

nfip <- nfip %>%
  mutate(incidencePM = round(nfip$polio/nfip$total*1e6,0))
knitr::kable(nfip)

group	grade	vaccin	total	polio	noPolio	incidencePM
cases	g2	yes	221998	54	221944	243
controls	g1g3	no	725173	391	724782	539
noConcent	g2	no	123605	56	123549	453

Wat kunnen we concluderen?

---

3.8 Confounding

• We observeren een lagere polio (P) incidentie voor kinderen bij wie de ouders geen toestemming gaven dan in de controle groep.

• Toestemming voor vaccinatie (V) is geassocieerd met de socio-economische status (S).

• Kinderen van lagere socio-economische status zijn meer resistent tegen de ziekte.

• De groepen van cases en controles zijn niet vergelijkbaar

  • verschil in leeftijd
  • verschil in socio-economische status en
  • verschil in vatbaarheid voor de ziekte.

---

3.9 Salk Study

3.9.1 Design

Een nieuwe studie werd uitgevoerd: dubbel blinde gerandomiseerde studie.

• Kinderen worden at random toegewezen aan controle of case arm van het experiment nadat de ouders toestemden met vaccinatie.
• Controle: vaccinatie met placebo
• Treatment: vaccinatie met vaccin
• Double blinding:
  • ouders en kinderen weten niet of ze werden gevaccineerd of niet
  • medische staf en onderzoekers weten niet of het kind het vaccin of de placebo kreeg

---

3.9.2 Data

salk <- data.frame(
  group=c("cases","control","noConcent"),
  treatment=c("vaccine","placebo","none"),
  total=c(200745,201229, 338778),polio=c(57,142,157)
  ) %>%
  mutate(
    noPolio = total-polio,
    incidencePM = round(polio/total*1e6,0)
    )
knitr::kable(salk)

group	treatment	total	polio	noPolio	incidencePM
cases	vaccine	200745	57	200688	284
control	placebo	201229	142	201087	706
noConcent	none	338778	157	338621	463

• We observeren een veel groter effect nu dat cases en controles vergelijkbaar zijn, incidentie van respectievelijk 284 and 706 per miljoen.

• De polio incidentie voor kinderen die geen toestemming geven blijft vergelijkbaar 453 and 463 per miljoen respectievelijk in the NFIP and Salk study.

---

4 Rol van Statistiek

• We hebben gezien dat
  • het belangrijk is om de scope van de studie goed te specifiëren voor de start van het experiment
  • randomisatie nodig is om een representatieve steekproef te nemen
  • steekproef grootte is heel belangrijk
  • we moeten ons bewust zijn van Confounding
  • een goede controle is belangrijk

\(\rightarrow\) Goede proefopzet is cruciaal!

---

• We hebben ook geobserveerd dat er variabiliteit is in de populatie
• We kunnen maar een beperkte steekproef nemen uit de populatie

\(\rightarrow\) onzekerheid in de resultaten \(\rightarrow\) onzekerheid in de conclusies

---

• Statistiek is de wetenschap voor het
  1. verzamelen (experimenteel design),
  2. exploren (data exploration) en
  3. leren van data zodat we hetgeen we observeren in de steekproef zouden kunnen veralgemenen naar de populatie terwijl we de onzekerheid quantificeren, controleren en rapporteren (statistisch modelleren en statistische besluitvorming).
• Statistiek speelt daarom een heel belangrijke rol in zowat alle wetenschappen

Home