1. Inleiding: Waarom Statistiek?

Author

Lieven Clement

Published

Invalid Date

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

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 geen geurhinder veroorzaken.
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: aselecte steekproef $\rightarrow$ representatief voor populatie!
- de personen worden at random verdeeld over twee behandelingsgroepen:
  - placebo (enkel antibiotica)
  - transplantie (antibiotica en microbiële transplantatie).
  - Randomisatie (at random toewijzen van proefpersonen aan de behandeling) 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 gemiddeld 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/sbc/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/sbc/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 <-.

We pipen (via |>) het ap dataframe naar de group_by functie om de data te groeperen per treatment trt: group_by(trt).
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.

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 aesthetics aes functie. aes(x=trt,y=mean)
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).
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 dan 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) +
  theme_minimal()

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

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

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

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

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

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") +
  theme_minimal()

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 ($p << 0.001$, 95% BI [12.9%, 31.6%]).

2.3 Random steekproeftrekking (Random sampling)

Wat zijn de consequenties van het gebruik van een random steekproeftrekking?
Random steekproeftrekking 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 (herhaalde steekproeftrekking/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>        
# ℹ 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>, UrineVol2 <int>, …

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.

We pipen de dataset naar de function filter om de data te filteren volgens leeftijd.
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)") + 
  theme_minimal()

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 data in kleine steekproeven kan variëren van steekproef tot steekproef.

We filteren op leeftijd en verwijderen ontbrekenden gegevens (NA, Not Available).
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

Beide verdelingen:

lijken een normale verdeling te volgen
met dezelfde variantie

$\rightarrow$ ze verschillen dus enkel in gemiddelde.

Onderzoeksvraag: “Is er een verschil in lengte tussen mannen en vrouwen?”
kunnen vertalen naar 1 karakteristiek: het gemiddeld verschil in lengte

HeightSum$mean[1]-HeightSum$mean[2]

[1] -13.81445

Gemiddeld zijn volwassen Amerikaanse vrouwen 13.8 cm kleiner dan mannen.

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)") + 
  theme_minimal()

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

HeightSumExp1$mean[1]-HeightSumExp1$mean[2]

[1] -19.16

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

Boxplot is beter:

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

HeightSumExp2$mean[1]-HeightSumExp2$mean[2]

[1] -4.2

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)") +
  theme_minimal()

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

HeightSumExp3$mean[1]-HeightSumExp3$mean[2]

[1] 4.82

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)") + 
  theme_minimal()

3.4.1 Herhaal het experiment heel veel keer

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

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

3.4.2 Effect van vergroten van steekproef

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

# 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)
}

res2 <- data.frame(
  verschil=colMeans(femSamps) - colMeans(malSamps))

res2$ng <- 25
res$ng <- 5

rbind(res,res2) |>
  ggplot(aes(y=verschil)) +
  geom_boxplot() +
  ylab("Gemiddeld Verschil (cm)") +
  xlab("") +
  theme_minimal() +
  facet_wrap(~ng)

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.

$\rightarrow$ Bijgevolg zijn onze conclusies ook onzeker en kunnen deze wijzigen van steekproef tot steekproef.

$\rightarrow$ Met statistiek controleren we de kans op het trekken van foute conclusies.

We zagen ook dat we het werkelijke verschil meer precies kunnen schatten als we de steekproefgrootte vergroten.

$\rightarrow$ Design van experiment zal dus bepalen hoe makkelijk we het werkelijk verschil op kunnen pikken met de steekproef.

4 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?

4.1 NFIP Study

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

4.1.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?

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

4.3 Salk Study

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

4.3.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 cases en controles vergelijkbaar zijn, incidentie van respectievelijk 284 en 706 per miljoen.
De polio incidentie voor kinderen die geen toestemming geven blijft vergelijkbaar 453 en 463 per miljoen respectievelijk in de NFIP- en Salk-studie.

5 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
- random steekproeftrekking 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

--- title: "1. Inleiding: Waarom Statistiek?" author: "Lieven Clement" date: "statOmics, Ghent University (https://statomics.github.io)" format: html: toc: true number-sections: true theme: cosmo highlight-style: tango code-tools: true embed-resources: true --- # Introductie ```{r include=FALSE} knitr::opts_chunk$set(include = TRUE, comment = NA, echo = TRUE, message = FALSE, warning = FALSE) library(tidyverse) ``` ![](../figures/ways_to_capture_b2b_customer_data_5bacd2dbe1a09.png) ![](../figures/wpGraph.jpeg) ![](../figures/mvpc.jpeg) - 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. *** # Case study: oksel microbiome ![](https://images.vrt.be/width1280/2018/06/27/4e72098d-7a0d-11e8-abcc-02b7b76bf47f.jpg){width=70%} - <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 geen geurhinder veroorzaken. - 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) *** ## Proefopzet (experimental design) ```{r out.width='80%',fig.asp=.8, fig.align='center',echo=FALSE} if ("pi"%in%ls()) rm("pi") kopvoeter<-function(x,y,angle=0,l=.2,cex.dot=.5,pch=19,col="black") { angle=angle/180*pi points(x,y,cex=cex.dot,pch=pch,col=col) lines(c(x,x+l*cos(-pi/2+angle)),c(y,y+l*sin(-pi/2+angle)),col=col) lines(c(x+l/2*cos(-pi/2+angle),x+l/2*cos(-pi/2+angle)+l/4*cos(angle)),c(y+l/2*sin(-pi/2+angle),y+l/2*sin(-pi/2+angle)+l/4*sin(angle)),col=col) lines(c(x+l/2*cos(-pi/2+angle),x+l/2*cos(-pi/2+angle)+l/4*cos(pi+angle)),c(y+l/2*sin(-pi/2+angle),y+l/2*sin(-pi/2+angle)+l/4*sin(pi+angle)),col=col) lines(c(x+l*cos(-pi/2+angle),x+l*cos(-pi/2+angle)+l/2*cos(-pi/2+pi/4+angle)),c(y+l*sin(-pi/2+angle),y+l*sin(-pi/2+angle)+l/2*sin(-pi/2+pi/4+angle)),col=col) lines(c(x+l*cos(-pi/2+angle),x+l*cos(-pi/2+angle)+l/2*cos(-pi/2-pi/4+angle)),c(y+l*sin(-pi/2+angle),y+l*sin(-pi/2+angle)+l/2*sin(-pi/2-pi/4+angle)),col=col) } par(mar=c(0,0,0,0),mai=c(0,0,0,0)) plot(0,0,xlab="",ylab="",xlim=c(0,10),ylim=c(0,10),col=0,xaxt="none",yaxt="none",axes=FALSE) rect(0,6,10,10,border="red",lwd=2) text(.5,8,"population",srt=90,col="red",cex=2) symbols (3, 8, circles=1.5, col="red",add=TRUE,fg="red",inches=FALSE,lwd=2) set.seed(330) grid=seq(0,1.3,.01) for (i in 1:50) { angle1=runif(n=1,min=0,max=360) angle2=runif(n=1,min=0,max=360) radius=sample(grid,prob=grid^2*pi/sum(grid^2*pi),size=1) kopvoeter(3+radius*cos(angle1/180*pi),8+radius*sin(angle1/180*pi),angle=angle2) } text(7.5,8,"Microbiome in population",col="red",cex=1.2) ``` *** ```{r out.width='80%',fig.asp=.8, fig.align='center',echo=FALSE} if ("pi"%in%ls()) rm("pi") kopvoeter<-function(x,y,angle=0,l=.2,cex.dot=.5,pch=19,col="black") { angle=angle/180*pi points(x,y,cex=cex.dot,pch=pch,col=col) lines(c(x,x+l*cos(-pi/2+angle)),c(y,y+l*sin(-pi/2+angle)),col=col) lines(c(x+l/2*cos(-pi/2+angle),x+l/2*cos(-pi/2+angle)+l/4*cos(angle)),c(y+l/2*sin(-pi/2+angle),y+l/2*sin(-pi/2+angle)+l/4*sin(angle)),col=col) lines(c(x+l/2*cos(-pi/2+angle),x+l/2*cos(-pi/2+angle)+l/4*cos(pi+angle)),c(y+l/2*sin(-pi/2+angle),y+l/2*sin(-pi/2+angle)+l/4*sin(pi+angle)),col=col) lines(c(x+l*cos(-pi/2+angle),x+l*cos(-pi/2+angle)+l/2*cos(-pi/2+pi/4+angle)),c(y+l*sin(-pi/2+angle),y+l*sin(-pi/2+angle)+l/2*sin(-pi/2+pi/4+angle)),col=col) lines(c(x+l*cos(-pi/2+angle),x+l*cos(-pi/2+angle)+l/2*cos(-pi/2-pi/4+angle)),c(y+l*sin(-pi/2+angle),y+l*sin(-pi/2+angle)+l/2*sin(-pi/2-pi/4+angle)),col=col) } par(mar=c(0,0,0,0),mai=c(0,0,0,0)) plot(0,0,xlab="",ylab="",xlim=c(0,10),ylim=c(0,10),col=0,xaxt="none",yaxt="none",axes=FALSE) rect(0,6,10,10,border="red",lwd=2) text(.5,8,"population",srt=90,col="red",cex=2) symbols (3, 8, circles=1.5, col="red",add=TRUE,fg="red",inches=FALSE,lwd=2) set.seed(330) grid=seq(0,1.3,.01) for (i in 1:50) { angle1=runif(n=1,min=0,max=360) angle2=runif(n=1,min=0,max=360) radius=sample(grid,prob=grid^2*pi/sum(grid^2*pi),size=1) kopvoeter(3+radius*cos(angle1/180*pi),8+radius*sin(angle1/180*pi),angle=angle2) } text(7.5,8,"Microbiome in population",col="red",cex=1.2) rect(0,0,10,4,border="blue",lwd=2) text(.5,2,"sample",srt=90,col="blue",cex=2) symbols (3, 2, circles=1.5, col="red",add=TRUE,fg="blue",inches=FALSE,lwd=2) for (i in 0:1) for (j in 0:4) { kopvoeter(2.1+j*(3.9-2.1)/4,1.1+i/2,col="black") } for (i in 2:3) for (j in 0:4) { kopvoeter(2.1+j*(3.9-2.1)/4,1.6+i/2,col="black") } text(7.5,2,"Microbiome in sample",col="blue",cex=1.2) arrows(3,5.9,3,4.1,col="black",lwd=3) text(1.5,5,"EXP. DESIGN (1)",col="black",cex=1.2) ``` *** - Experiment: - 20 personen worden at random uit de populatie getrokken van personen met een okselgeur: aselecte steekproef $\rightarrow$ representatief voor populatie! - de personen worden at random verdeeld over twee behandelingsgroepen: - placebo (enkel antibiotica) - transplantie (antibiotica en microbiële transplantatie). - Randomisatie (at random toewijzen van proefpersonen aan de behandeling) 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://journals.plos.org/plosone/article/figure/image?size=large&download=&id=10.1371/journal.pone.0070538.g001){width=70%} <https://doi.org/10.1371/journal.pone.0070538> *** Vertaal onderzoeksvraag naar iets wat we kunnen quantificeren: Is er gemiddeld een verschil in relatieve abundantie van *Staphylococcus spp.* in het microbiome van de placebo groep en de transplantatie *** ```{r out.width='80%',fig.asp=.8, fig.align='center',echo=FALSE} if ("pi"%in%ls()) rm("pi") kopvoeter<-function(x,y,angle=0,l=.2,cex.dot=.5,pch=19,col="black") { angle=angle/180*pi points(x,y,cex=cex.dot,pch=pch,col=col) lines(c(x,x+l*cos(-pi/2+angle)),c(y,y+l*sin(-pi/2+angle)),col=col) lines(c(x+l/2*cos(-pi/2+angle),x+l/2*cos(-pi/2+angle)+l/4*cos(angle)),c(y+l/2*sin(-pi/2+angle),y+l/2*sin(-pi/2+angle)+l/4*sin(angle)),col=col) lines(c(x+l/2*cos(-pi/2+angle),x+l/2*cos(-pi/2+angle)+l/4*cos(pi+angle)),c(y+l/2*sin(-pi/2+angle),y+l/2*sin(-pi/2+angle)+l/4*sin(pi+angle)),col=col) lines(c(x+l*cos(-pi/2+angle),x+l*cos(-pi/2+angle)+l/2*cos(-pi/2+pi/4+angle)),c(y+l*sin(-pi/2+angle),y+l*sin(-pi/2+angle)+l/2*sin(-pi/2+pi/4+angle)),col=col) lines(c(x+l*cos(-pi/2+angle),x+l*cos(-pi/2+angle)+l/2*cos(-pi/2-pi/4+angle)),c(y+l*sin(-pi/2+angle),y+l*sin(-pi/2+angle)+l/2*sin(-pi/2-pi/4+angle)),col=col) } par(mar=c(0,0,0,0),mai=c(0,0,0,0)) plot(0,0,xlab="",ylab="",xlim=c(0,10),ylim=c(0,10),col=0,xaxt="none",yaxt="none",axes=FALSE) rect(0,6,10,10,border="red",lwd=2) text(.5,8,"population",srt=90,col="red",cex=2) symbols (3, 8, circles=1.5, col="red",add=TRUE,fg="red",inches=FALSE,lwd=2) set.seed(330) grid=seq(0,1.3,.01) for (i in 1:50) { angle1=runif(n=1,min=0,max=360) angle2=runif(n=1,min=0,max=360) radius=sample(grid,prob=grid^2*pi/sum(grid^2*pi),size=1) kopvoeter(3+radius*cos(angle1/180*pi),8+radius*sin(angle1/180*pi),angle=angle2) } text(7.5,8,"Microbiome in population",col="red",cex=1.2) rect(0,0,10,4,border="blue",lwd=2) text(.5,2,"sample",srt=90,col="blue",cex=2) symbols (3, 2, circles=1.5, col="red",add=TRUE,fg="blue",inches=FALSE,lwd=2) for (i in 0:1) for (j in 0:4) { kopvoeter(2.1+j*(3.9-2.1)/4,1.1+i/2,col="purple") } for (i in 2:3) for (j in 0:4) { kopvoeter(2.1+j*(3.9-2.1)/4,1.6+i/2,col="orange") } text(7.5,2,"Microbiome in sample",col="blue",cex=1.2) arrows(3,5.9,3,4.1,col="black",lwd=3) text(1.5,5,"EXP. DESIGN (1)",col="black",cex=1.2) text(7.5,.5,"DATA EXPLORATION &\nDESCRIPTIVE STATISTICS (2)",col="black",cex=1.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. *** ### 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. ```{r} read_lines("https://raw.githubusercontent.com/statOmics/sbc/master/data/armpit.csv") ``` - 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. ```{r} ap <- read_csv("https://raw.githubusercontent.com/statOmics/sbc/master/data/armpit.csv") ap ``` *** ### 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. ```{r} 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 ```{r} apRelSum ``` We kunnen ook een tabel in de webpagina of het pdf bestand integreren via het commando kable van het knitr pakket: ```{r} knitr::kable(apRelSum) ``` *** ### 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. #### 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 aesthetics `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 dan deze van de bar. ```{r} apRelSum |> ggplot(aes(x=trt,y=mean)) + geom_bar(stat="identity") + geom_errorbar(aes(ymin=mean-sd,ymax=mean+sd),width=.2) + theme_minimal() ``` - Is deze plot informatief?? *** #### 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()` ```{r} ap |> ggplot(aes(x=trt,y=rel)) + geom_boxplot() + theme_minimal() ``` *** Merk op dat we de plot ook op kunnen slaan als een object. ```{r} apBoxplot <- ap |> ggplot(aes(x=trt,y=rel)) + geom_boxplot() + theme_minimal() ``` De plot wordt dan niet gemaakt. Om de plot weer te geven kunnen we het object aanroepen ```{r} 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. ```{r} ap |> ggplot(aes(x=trt,y=rel)) + geom_boxplot(outlier.shape=NA) + geom_point(position="jitter") + theme_minimal() ``` 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. ```{r} apBoxplot + geom_point(position="jitter") ``` We hadden het resultaat ook opnieuw op kunnen slaan voor later hergebruik. ```{r} 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? *** ```{r out.width='80%',fig.asp=.8, fig.align='center',echo=FALSE} if ("pi"%in%ls()) rm("pi") kopvoeter<-function(x,y,angle=0,l=.2,cex.dot=.5,pch=19,col="black") { angle=angle/180*pi points(x,y,cex=cex.dot,pch=pch,col=col) lines(c(x,x+l*cos(-pi/2+angle)),c(y,y+l*sin(-pi/2+angle)),col=col) lines(c(x+l/2*cos(-pi/2+angle),x+l/2*cos(-pi/2+angle)+l/4*cos(angle)),c(y+l/2*sin(-pi/2+angle),y+l/2*sin(-pi/2+angle)+l/4*sin(angle)),col=col) lines(c(x+l/2*cos(-pi/2+angle),x+l/2*cos(-pi/2+angle)+l/4*cos(pi+angle)),c(y+l/2*sin(-pi/2+angle),y+l/2*sin(-pi/2+angle)+l/4*sin(pi+angle)),col=col) lines(c(x+l*cos(-pi/2+angle),x+l*cos(-pi/2+angle)+l/2*cos(-pi/2+pi/4+angle)),c(y+l*sin(-pi/2+angle),y+l*sin(-pi/2+angle)+l/2*sin(-pi/2+pi/4+angle)),col=col) lines(c(x+l*cos(-pi/2+angle),x+l*cos(-pi/2+angle)+l/2*cos(-pi/2-pi/4+angle)),c(y+l*sin(-pi/2+angle),y+l*sin(-pi/2+angle)+l/2*sin(-pi/2-pi/4+angle)),col=col) } par(mar=c(0,0,0,0),mai=c(0,0,0,0)) plot(0,0,xlab="",ylab="",xlim=c(0,10),ylim=c(0,10),col=0,xaxt="none",yaxt="none",axes=FALSE) rect(0,6,10,10,border="red",lwd=2) text(.5,8,"population",srt=90,col="red",cex=2) symbols (3, 8, circles=1.5, col="red",add=TRUE,fg="red",inches=FALSE,lwd=2) set.seed(330) grid=seq(0,1.3,.01) for (i in 1:50) { angle1=runif(n=1,min=0,max=360) angle2=runif(n=1,min=0,max=360) radius=sample(grid,prob=grid^2*pi/sum(grid^2*pi),size=1) kopvoeter(3+radius*cos(angle1/180*pi),8+radius*sin(angle1/180*pi),angle=angle2) } text(7.5,8,"Microbiome in population",col="red",cex=1.2) rect(0,0,10,4,border="blue",lwd=2) text(.5,2,"sample",srt=90,col="blue",cex=2) symbols (3, 2, circles=1.5, col="red",add=TRUE,fg="blue",inches=FALSE,lwd=2) for (i in 0:1) for (j in 0:4) { kopvoeter(2.1+j*(3.9-2.1)/4,1.1+i/2,col="purple") } for (i in 2:3) for (j in 0:4) { kopvoeter(2.1+j*(3.9-2.1)/4,1.6+i/2,col="orange") } text(7.5,2,"Microbiome in sample",col="blue",cex=1.2) arrows(3,5.9,3,4.1,col="black",lwd=3) text(1.5,5,"EXP. DESIGN (1)",col="black",cex=1.2) text(7.5,.5,"DATA EXPLORATION &\nDESCRIPTIVE STATISTICS (2)",col="black",cex=1.2) arrows(7,4.1,7,5.9,col="black",lwd=3) text(8.5,5,"ESTIMATION &\nINFERENCE (3)",col="black",cex=1.2) ``` *** - 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. ```{r} t.test(rel~trt,data=ap) ``` Conclusie: Gemiddeld is de relatieve abundantie van Staphylococcus in het microbiome van personen met een zweetgeur `r format(apRelSum$mean[2]-apRelSum$mean[1],digits=3)`% hoger na de transplantie dan na de placebo behandeling ($p << 0.001$, 95% BI [`r paste0(round(sort(abs(t.test(rel~trt,data=ap)$conf.int)),1),"%")`]). *** ## Random steekproeftrekking (Random sampling) - Wat zijn de consequenties van het gebruik van een random steekproeftrekking? - Random steekproeftrekking 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? *** # 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 (`herhaalde steekproeftrekking`/`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` *** ```{r} library(NHANES) head(NHANES) glimpse(NHANES) ``` *** ## 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. ```{r} NHANES|> filter(Age >= 18 & !is.na(Height)) |> ggplot(aes(x = Height))+ geom_histogram() + facet_grid(Gender ~ .) + xlab("Lengte (cm)") + theme_minimal() ``` *** 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 data 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. ```{r} 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). ```{r} HeightSum <- nhanesSub |> group_by(Gender) |> summarize_at("Height", list(mean = mean, sd = sd) ) knitr::kable( HeightSum |> mutate_if(is.numeric, round, digits=1) ) ``` Beide verdelingen: - lijken een normale verdeling te volgen - met dezelfde variantie $\rightarrow$ ze verschillen dus enkel in gemiddelde. - Onderzoeksvraag: "Is er een verschil in lengte tussen mannen en vrouwen?" - kunnen vertalen naar 1 karakteristiek: het gemiddeld verschil in lengte ```{r} HeightSum$mean[1]-HeightSum$mean[2] ``` Gemiddeld zijn volwassen Amerikaanse vrouwen 13.8 cm kleiner dan mannen. *** ## 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. *** ```{r} 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 ``` *** Data Exploratie ```{r} samp1 |> ggplot(aes(x=Height)) + geom_histogram() + facet_grid(Gender~.) + xlab("Lengte (cm)") + theme_minimal() HeightSumExp1 <- samp1 |> group_by(Gender) |> summarize_at("Height", list(mean = mean, sd = sd) ) HeightSumExp1 HeightSumExp1$mean[1]-HeightSumExp1$mean[2] ``` Histogram is niet zinvol als we maar zo weinig datapunten hebben. *** Boxplot is beter: ```{r echo=FALSE} samp1 |> ggplot(aes(x = Gender,y = Height)) + geom_boxplot(outlier.shape = NA) + geom_point(position = "jitter") + geom_point( aes(x = 1, y = HeightSumExp1$mean[1]), size = 3, pch = 17, color="darkred") + geom_point( aes(x = 2, y = HeightSumExp1$mean[2]), size = 3, pch = 17, color = "darkred") + ylab("Height (cm)") + theme_minimal() ``` ## Herhaal het experiment Als we het experiment herhalen selecteren we andere mensen en verkrijgen we andere resultaten. ```{r} 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 HeightSumExp2$mean[1]-HeightSumExp2$mean[2] 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)") + theme_minimal() ``` ## Herhaal het experiment opnieuw ```{r} 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 HeightSumExp3$mean[1]-HeightSumExp3$mean[2] 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)") + theme_minimal() ``` ### Herhaal het experiment heel veel keer ```{r} 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)) res |> ggplot(aes(y=verschil)) + geom_boxplot() + ylab("Gemiddeld Verschil (cm)") + xlab("") + theme_minimal() ``` ### Effect van vergroten van steekproef ```{r} set.seed(15142) # Aantal simulaties en steekproefgrootte per groep nSim <- 10000 nSamp <- 25 # 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) } res2 <- data.frame( verschil=colMeans(femSamps) - colMeans(malSamps)) res2$ng <- 25 res$ng <- 5 rbind(res,res2) |> ggplot(aes(y=verschil)) + geom_boxplot() + ylab("Gemiddeld Verschil (cm)") + xlab("") + theme_minimal() + facet_wrap(~ng) ``` ## 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. $\rightarrow$ Bijgevolg zijn onze conclusies ook onzeker en kunnen deze wijzigen van steekproef tot steekproef. $\rightarrow$ Met statistiek controleren we de kans op het trekken van foute conclusies. - We zagen ook dat we het werkelijke verschil meer precies kunnen schatten als we de steekproefgrootte vergroten. $\rightarrow$ Design van experiment zal dus bepalen hoe makkelijk we het werkelijk verschil op kunnen pikken met de steekproef. *** # 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? *** ## NFIP Study ### 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. ### Data ```{r} 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) ``` Vergelijk de polio incidentie? *** ```{r} nfip <- nfip |> mutate(incidencePM = round(nfip$polio/nfip$total*1e6,0)) knitr::kable(nfip) ``` Wat kunnen we concluderen? *** ## Confounding ```{r,echo=FALSE, fig.align = "center",out.width = '50%'} plot(c(0,0,1),c(-2,2,0),pch=c("S","V","P"),xaxt="none",yaxt="none",axes=FALSE,xlab="",ylab="",cex=4,ylim=c(-2.2,2.2)) arrows(x0=0.1,x1=.9,y0=1.8,y1=0.1,lwd=4) arrows(x0=0.1,x1=.9,y0=-1.8,y1=-0.2,lwd=4) arrows(x0=0,x1=0,y0=-1.4,y1=1.4,lwd=4) ``` - 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. *** ## Salk Study ### 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 *** ### Data ```{r} 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) ``` - We observeren een veel groter effect nu cases en controles vergelijkbaar zijn, incidentie van respectievelijk `r salk$incidencePM[1]` en `r salk$incidencePM[2]` per miljoen. - De polio incidentie voor kinderen die geen toestemming geven blijft vergelijkbaar `r nfip$incidencePM[3]` en `r salk$incidencePM[3]` per miljoen respectievelijk in de NFIP- en Salk-studie. *** # 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 - random steekproeftrekking 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