---
title: "Categorische data analyse - full"  
output:
    html_document:
      code_download: true    
      theme: cosmo
      toc: true
      toc_float: true
      highlight: tango
      number_sections: true
---
# Koffie en zwangerschap

In de jaren 80 werd op basis van experimenten met ratten in de VS besloten om aan zwangere vrouwen af te raden koffie te drinken. Koffie tijdens de zwangerschap zou een verhoogd risiko op te laag geboortegewicht voor de baby met zich meebrengen. In 1982 rapporteerden Lin et al. in de New England Journal of Medicine een studie die rechtstreeks het effect van koffie op het geboortegewicht bij mensen onderzocht. Voor zo’n 12.197 borelingen die lukraak gekozen werden, werd genoteerd of het geboortegewicht te laag was of niet, en voor de moeder hoeveel koffie ze dronk.

De dataset `birthWeight` bestaat uit twee variabelen: Koffie en Gewicht. Elke variabele is een factor variabele, waarbij:

`Koffie`: 0 = gemiddeld 0 koppen koffie per dag, 1 = gemiddeld 1 tot 3 koppen koffie per dag, 2 = gemiddeld 4 of meer koppen koffie per dag.

`Gewicht`: 0 = Geen te laag geboortegewicht, 1 = Te laag geboortegewicht

```{r}
library(ggplot2)
library(dplyr)
```

## Lees de data in

```{r}
birthWeight <- data.frame(
Koffie = as.factor(rep(c(0,0,1,1,2,2), c(507,6410,388,4309,62,521))),
Gewicht = as.factor(rep(c(1,0,1,0,1,0), c(507,6410,388,4309,62,521)))
)
```

## Data Exploratie

```{r}
plot(birthWeight)
```

Alternatieve visualisatie
```{r}
frequencyTable <- as.data.frame(table(birthWeight))
# Stacked barplot with multiple groups
ggplot(data=frequencyTable, aes(x=Koffie, y=Freq, fill=Gewicht)) +
  geom_bar(stat="identity") +
  theme_bw()
```

### Relatief risico

Voor het berekenen van het relatief risico en de odds ratio, stellen we de gegeven data voor aan de hand van de volgende compacte en overzichtelijk vorm:

```{r}
tabel1 <- table(birthWeight$Gewicht, birthWeight$Koffie, dnn=c("Gewicht", "Koffie")) #dnn=benamingen in de frequentietabel
rownames(tabel1) <- c("Geen te laag geboortegewicht", "Te laag geboortegewicht")
colnames(tabel1) <-c("0 koppen", "1-3 koppen", "4 of meer koppen")
tabel1
```

Handig bij de berekening van relatief risico kan het toevoegen van rij- en kolomtotalen zijn. 
```{r}
tabel2 <- addmargins(tabel1)
```

#### Relatief risico voor koffiedrinkers van 1-3 koppen t.o.v. de referentieklasse

Het risico of 'de kans op te laag geboortegewicht' voor koffiedrinkers van 1-3 koppen per dag wordt
genoteerd met $\pi_{K13}$ en voor niet-koffiedrinkers met $\pi_{\bar{K}}$. In die notatie verwijst de index K13 naar ‘1-3 koppen koffie per dag’ en $\bar{K}$ naar 'geen koffie', niet blootgesteld aan koffie. Respectievelijk worden die kansen op basis van de gegevens geschat als 

- $\hat{\pi}_{K13}$ = `r tabel2[2,2]/tabel2[3,2]` = 8.3% voor koffiedrinkers 
- $\hat{\pi}_{\bar{K}}$ = `r tabel2[2,1]/tabel2[3,1]` = 7.3% en voor niet-koffiedrinkers. 

Het relatief risico in de doelpopulatie is $\pi_{K13}/\pi_{\bar{K}}$. Een schatting voor het
relatief risico (op te laag geboortegewicht) van koffiedrinkers van 1 tot 3 koppen per dag ten opzichte van
niet-koffiedrinkers is dus $\hat{\pi}_{K13}/\hat{\pi}_{\bar{K}}$ = 1.132.


Proporties t.o.v. rij- en kolomtotalen kunnen ook berekend worden als volgt: 
```{r}
prop.table(tabel1,1) #proportie tov rijtotaal
prop.table(tabel1,2) #proportie tov kolomtotaal
```


Het risico op een boreling met laag geboortegewicht is dus 13.2% hoger voor moeders die 1 tot 3 koppen koffie per dag drinken dan voor niet-koffiedrinkende moeders. 

#### Relatief risico voor koffiedrinkers van 4 of meer koppen t.o.v. de referentieklasse

Het risico of 'de kans op te laag geboortegewicht' voor koffiedrinkers van 4 of meer koppen per dag wordt
genoteerd met $\pi_{K4}$ en voor niet-koffiedrinkers met $\pi_{\bar{K}}$. In die notatie verwijst de index K4 naar ‘4 of meer koppen koffie per dag’ en $\bar{K}$ naar 'geen koffie', niet blootgesteld aan koffie. Respectievelijk worden die kansen op basis van de gegevens geschat als 

- $\hat{\pi}_{K4}$ = `r tabel2[2,3]/tabel2[3,3]` = 10.6% voor koffiedrinkers 
- $\hat{\pi}_{\bar{K}}$ = `r tabel2[2,1]/tabel2[3,1]` = 7.3% en voor niet-koffiedrinkers. 

Het relatief risico in de doelpopulatie is $\pi_{K4}/\pi_{\bar{K}}$. Een schatting voor het
relatief risico (op te laag geboortegewicht) van koffiedrinkers van 4 of meer koppen per dag ten opzichte van
niet-koffiedrinkers is dus $\hat{\pi}_{K4}/\hat{\pi}_{\bar{K}}$ = 1.451.

Het risico op een boreling met laag geboortegewicht is dus 45.1% hoger voor moeders die 4 of meer koppen koffie per dag drinken dan voor niet-koffiedrinkende moeders. 


### Odds ratio

Een andere manier om over associatie te spreken tussen 2 categorische veranderlijken is via de odds ratio.

Odds ratio op laag geboortegewicht voor vrouwen die 1 tot 3 koppen koffie per dag drinken versus vrouwen die geen koffie drinken, wordt als volgt berekend:

\[OR = \frac{\text{odds}_{K13}}{\text{odds}_{\bar{K}}} = \frac{\frac{\pi_{K13}}{1-\pi_{K13}}}{\frac{\pi_{\bar{K}}}{1-\pi_{\bar{K}}}}\]

- Odds op laag geboortegewicht voor koffiedrinkers van 1-3 koppen koffie per dag:

\[\frac{\hat{\pi}_{K13}}{1-\hat{\pi}_{K13}} = \frac{388/4697}{4309/4697} = 0.090\]

- Odds op laag geboortegewicht voor niet-koffiedrinkers:

\[\frac{\hat{\pi}_{\bar{K}}}{1-\hat{\pi}_K} =\frac{507/6917}{6410/6917} = \frac{507}{6410} = 0.079\]

Daarmee vinden we dat de Odds ratio op laag geboortegewicht tussen koffiedrinkers van 1-3 koppen koffie per dag en niet-koffiedrinkers gelijk is aan $\frac{0.090}{0.079}$ = 1.139.

De odds op een boreling met laag geboortegewicht is 13.9 % hoger bij moeders die gemiddeld 1-3 koppen koffie per drinken dan bij niet-koffiedrinkende moeders (Koffiedrinkende moeders die 1-3 koppen koffie drinken per dag hebben dus een hogere kans op een boreling met te laag geboortegewicht dan niet-koffiedrinkende moeders).

\textbf{Opmerking} 

Op analoge manier kan de odds ratio berekend worden op:

- laag geboortegewicht voor vrouwen die 4 of meer koppen koffie per dag drinken versus vrouwen die geen koffie per dag drinken

- laag geboortegewicht voor vrouwen die 4 of meer koppen koffie per dag drinken versus vrouwen die 1-3 koppen koffie per dag drinken

## Nul- en alternative hypothese

De Pearson Chi-kwadraat test wordt gebruikt om het verband na te gaan tussen koffiegebruik van de moeders en het geboortegewicht van hun baby. We toetsen de nulhypothese van ‘onafhankelijkheid’ of ‘geen associatie’ tussen 2 discrete (binaire) veranderlijken.

H0: Er is geen associatie tussen het koffiegebruik en het geboortegewicht

HA:  Er is een associatie tussen het koffiegebruik en het geboortegewicht

## Statistische test en assumpties

Alvorens de test uit te voeren zullen we de assumpties van de test nagaan. De assumpties van de Chi-kwadraat test zijn

- De gegevens dienen onafhankelijk van elkaar verzameld te zijn.
- Minstens 80% van de verwachte waarden onder de nulhypothese in de kruistabel dienen groter dan 5 te zijn.

We nemen de assumptie van onafhankelijkheid aan via het design. De verwachte waarden van de kruistabel onder de nulhypothese kan men bekomen via de `chisq.test` functie, die alsook de Chi-kwadraat test uitvoert.

```{r}
tabel1 <- table(birthWeight$Gewicht, birthWeight$Koffie, dnn=c("Gewicht", "Koffie")) #dnn=benamingen in de frequentietabel
rownames(tabel1) <- c("Geen te laag geboortegewicht", "Te laag geboortegewicht")
colnames(tabel1) <-c("0 koppen", "1-3 koppen", "4 of meer koppen")
tabel1
```

```{r}
test <- chisq.test(tabel1)

# verwachte waarden
test$expected
```

We zien dat alle verwachte waarden groter zijn dan 5 en gaan dus verder met de output van de test te bekijken.

```{r}
test
```

Merk op dat we de $\chi^2$-test op twee manieren kunnen uitvoeren. Enerzijds door een kruistabel mee te geven zoals hierboven, anderzijds door de dataset met de twee variabelen mee te geven aan de `chisq.test` functie.

```{r}
chisq.test(birthWeight$Gewicht, birthWeight$Koffie)
```

De p-waarde van de statistische test bedraagt 0.007. Er is dus voldoende bewijs om de nulhypothese van onafhankelijkheid te verwerpen op het 5% significantieniveau.

### Conclusie

Er is een heel significante associatie tussen het koffiegebruik van moeders en tussen het voorkomen van een laag geboortegewicht van hun baby's (p=`r round(chisq.test(birthWeight$Gewicht, birthWeight$Koffie)$p.value,3)`).

# Logistische regressie

Net zoals we t-testen kunnen veralgemenen naar een lineair model, kunnen we voor binaire data de Chi-kwadraat test veralgemenen naar een logistisch regressiemodel. Dit heeft verschillende voordelen, zoals een onmiddellijke interpretatie van de parameters in termen van odds ratio’s ten opzichte van het referentie niveau (weergegeven door het intercept).

## Model

```{r}
birthWeightLogit <- glm(
  Gewicht ~ Koffie,
  data = birthWeight,
  family = binomial)

summary(birthWeightLogit)
```

Het intercept is de log-odds op laag geboortegewicht in de referentieklasse "Koffie0" of dus in de klasse waarbij vrouwen geen koffie drinken. De hellingstermen zijn log odds ratio’s tussen de klassen van koffiedrinkers en de referentieklasse (niet-koffiedrinkers).
De analyse laat dus toe om de resultaten onmiddellijk te interpreteren in termen van odds’es en odds-ratio’s.

- De odds op een laag geboortegewicht is bij moeders die geen koffie gebruiken is bijvoorbeeld: 

$$
\widehat{\text{ODDS}}_{\bar K}=\exp(\hat{\beta}_0)  = `r round(exp(birthWeightLogit$coef[1]),3)`
$$

- De odds op een laag geboortegewicht is bij moeders die dagelijks 1-3 koppen koffie gebruiken is bijvoorbeeld: 

$$
\widehat{\text{ODDS}}_{1-3}=\exp(\hat\beta_0+\hat\beta_1) = \exp(`r round(birthWeightLogit$coef[1],2)` + `r round(birthWeightLogit$coef[2],2)`) = \exp(`r round(sum(birthWeightLogit$coef[1:2]),2)`) = `r round(exp(sum(birthWeightLogit$coef[1:2])),3)`
$$

- De odds ratio op een laag geboortegewicht tussen moeders die dagelijks 1-3 koppen koffie gebruiken en moeders die geen koffie drinken is bijvoorbeeld:

$$
\log \widehat{\text{OR}} = \log \widehat{\text{ODDS}}_{13} - \log \widehat{\text{ODDS}}_{\bar K} = \hat\beta_0 + \hat\beta_1 -\hat\beta_0= \hat\beta_1
$$

$$
\widehat{\text{OR}} = \exp(\hat{\beta}_1) = `r round(exp(birthWeightLogit$coef[2]),3)`
$$

Merk op dat deze waarden identiek zijn aan wat we bekwamen in de data exploratie! 

## Inferentie

Net zoals bij een ANOVA analyse bij een continue response, kan de anova functie worden gebruikt voor het testen of er een associatie is tussen het voorkomen van laag geboortegewicht en het koffiegebruik.

```{r}
anova(birthWeightLogit, test = "Chisq")
```

De $\chi^2$-test in het logistische regressiemodel geeft eveneens aan dat er een significante associatie is tussen de uitkomst (voorkomen van laag geboortegewicht) en de factor (de hoeveelheid koffie die gedronken wordt) (p=0.009). De p-waarde is bijna equivalent aan de p-waarde van de $\chi^2$-test die we eerder uitvoerden (p-waarde = 0.007).

Aangezien er een heel significante associatie is, kunnen we post-hoc tests kunnen uitvoeren om te evalueren welke odds ratio’s verschillend zijn.

```{r}
library(multcomp)
posthoc <- glht(birthWeightLogit,
  linfct = mcp(Koffie = "Tukey"))

posthocTests <- summary(posthoc)
posthocTests

posthocBI <- confint(posthoc)
posthocBI
```

Door middel van de confint functie worden BI’s verkregen op de log-odds ratios die gecorrigeerd zijn voor multiple testing. Deze kunnen als volgt worden teruggetransformeerd naar odds ratios:

```{r}
OR <- exp(posthocBI$confint)
OR
```

## Conclusie

Er is een heel significantie associatie tussen het voorkomen van een laag geboortegewicht bij borelingen en het koffiegebruik van de moeders (p-waarde = `r round(anova(birthWeightLogit, test = "Chisq")$Pr[2],3)`). 
De odds op een boreling met laag geboortegewicht is anderhalf keer hoger bij moeders die 4 of meer koppen koffie drinken per dag dan bij moeders die geen koffie drinken (p-waarde = `r round(posthocTests$test$pvalues[2],3)`, 95% BI op OR [`r round(exp(posthocBI$confint[2,2:3]),3)`]). De overige odds ratio's zijn niet significant verschillend op het 5% significantieniveau.

Merk op dat we hier niet kunnen concluderen dat er een causaal verband is tussen het koffieverbruik van moeders en het voorkomen van het baren van kinderen met een laag geboortegewicht vermits het observationele studie betreft. De vrouwen uit de verschillende koffiegebruiksgroepen die deelnamen aan de studie verschillen mogelijks ook in andere factoren. Rookgedrag en alcoholgebruik zijn bijvoorbeeld typisch geassocieerd met het koffiegebruik en kunnen ook een nefaste invloed hebben op de ontwikkeling van de foetus. 

### Logistische regressie via de binomiale familie

Voor het uitvoeren van logistische regressie is onafhankelijkheid van de data vereist. 
Wat hier het geval is, de baby's werden immers allemaal lukraak gekozen. 

Dat laat ons ook toe om de binaire random variabele om te zetten tot binomiaal verdeelde random variabele voor elke groep door het aantal babies met laag geboortegewicht te berekenen. Deze variabele volgt dan een bionomiale verdeling met parameters $\pi_k$, de kans op een baby met ondergewicht en $n_k$ het aantal vrouwen per groep $k$. 

Dat laat ons toe om de dataset op een meer efficiënte manier voor te stellen, namelijk als een tabel met daarin de aantallen per categorie. 

Een binomiaal verdeelde variabele kunnen we onmiddellijk in R modelleren met de `glm`-functie, we moeten dan het aantal successen (1, hier aantal baby's met ondergewicht) en het aantal failures (0, hier aantal baby's zonder ondergewicht) meegeven aan R. 

Wanneer we de output van deze fit vergelijken met de logistische regressie die we hierboven verkregen, zien we hetzelfde resultaat. 

Het voordeel van deze werkwijze is dat we de logistische regressie op een meer efficiënte manier hebben gecodeerd d.m.v. de binomiale familie.

```{r}
birthWeightSum <- data.frame(
Koffie = as.factor(c(0,1,2)),
laag = c(507, 388,62),
hoog = c(6410,4309,521)
)
birthWeightLogitSum<- glm(
cbind(laag,hoog) ~ Koffie,
data = birthWeightSum,
family = binomial)
summary(birthWeightLogitSum)
```


## Illustratie Beperkte vergelijking

Ter illustratie tonen we hier ook een analyse met een twee op twee tabel. 
We beperken hiervoor de oorspronkelijke vraagstelling tot het vergelijken van het geboortegewicht van baby's tussen 2 groepen moeders, namelijk zij die koffie drinken en zij die geen koffie drinken. We kunnen dit eenvoudig doen door het aggregeren van de groepen moeders die 1-3 koppen koffie per drinken en zij die 4 of meer koppen koffie per dag drinken.

```{r}
birthWeight$Koffie2 <- ifelse(birthWeight$Koffie == 0, 0,1)

tabel3 <- table(birthWeight$Gewicht, birthWeight$Koffie2)
rownames(tabel3) <- c("Geen te laag geboortegewicht", "Te laag geboortegewicht")
colnames(tabel3) <-c("0 koppen", "1 of meer koppen")
tabel3
```

Als voorwaarde zullen we opnieuw de verwachte waarden onder de nulhypothese bekijken.

```{r}
chisq.test(tabel3)$expected
```

Aangezien alle waarden groter zijn dan 5, is aan de voorwaarden voldaan, en kunnen we de test uitvoeren en interpreteren.

```{r}
chisq.test(tabel3)
```

Uit deze analyse kunnen we besluiten dat er een significante associatie is tussen het koffiegebruik van moeders  en het voorkomen van ondergewicht van hun babies (`r round(chisq.test(tabel3)$p.value, 3)`). 

Merk op dat de output van de test ditmaal vermeldt `Pearson's Chi-squared test **with Yates' continuity correction**`. Aangezien we hier met een zekere vorm van teldata werken, kan de bekomen teststatistiek van een Pearson Chi-kwadraat test slechts enkele waarden aannemen; de distributie van de teststatistiek is dus discreet, terwijl we deze onder de nulhypothese wel benaderen met een continue $\chi^2$ distributie. Deze benadering is gepaster voor grotere kruistabellen met veel data (vandaar de assumptie dat 80% van de waarden onder H0 groter dan 5 dienen te zijn). Voor de kleinst mogelijke kruistabellen, zoals deze 2×2 kruistabel, wordt de continuïteitscorrectie van Yates toegepast, omdat hierdoor deze benadering beter opgaat. `R` doet dit dus automatisch voor ons bij kleine kruistabellen, maar dit kan ook zelf aangepast worden via het argument `correct` in de `chisq.test` functie.

### Fisher exact test

Merk op dat voor 2×2 tabellen ook een exacte test mogelijk is, de Fisher exact test. Deze is echter conservatief (d.w.z. dat ze een kleinere kans heeft op een Type I fout dan vooropgesteld, en bijgevolg een grotere kans op Type II fout). Merk op dat men deze test dus enkel zal doen wanneer de $\chi^2$ benadering voor de verdeling van de toetsingsgrootheid niet verantwoord is, m.a.w. in het geval in geen enkele van de cellen het verwachte aantal onder $H_0$ kleiner is dan 5.

```{r}
fisher.test(tabel3)
```


# Plantenbloei in functie van diens orientatie

Niet alle planten bloeien jaarlijks, al dan niet door hun biologie, of door omgevingsvariabelen.
Inderdaad, indien de condities niet optimaal zijn, zal een plant soms niet bloeien, omdat alle kostbare energie nodig is om te overleven, en er geen resterende energie is voor reproductie.
Dit is belangrijk voor tuinplanten, waarvan de eigena(a)r(es) typisch graag heeft dat de planten jaarlijks bloeien.
Vele planten aarden beter indien die gericht staan naar het Zuiden.
Een nieuw plantenras wordt bestudeerd om te kijken of diens oriëntatie een invloed heeft op de bloeicapaciteit.
Hiervoor plaatst men $50$ volwassen planten eens een zomer aan de zuidkant van een gebouw, en eens een zomer aan de oostkant.
Voor elk jaar wordt genoteerd of de plant al dan niet bloeit.

```{r}
data <- matrix(c(20, 10,
                 8, 12),
               nrow=2, ncol=2, byrow=TRUE,
               dimnames=list(c("Zuidkant: bloeien", "Zuidkant: Niet bloeien"),
                             c("Oostkant: bloeien", "Oostkant: Niet bloeien")))
data
```

Merk op dat we hier met gepaarde data werken, elke plant werd namelijk tweemaal geobserveerd over de periode van twee jaar.
Hierdoor is de Chi-kwadraat test niet langer gepast.

De kans op bloeien met een Zuiderse oriëntatie bedraagt $(20+10)/50 = 60$%
De kans op bloeien met een Oosterse oriëntatie bedraagt $(20+8)/50 = 56$%

Het absoluut risicoverschil (ARV) bedraagt dus **60% - 56% = 4%**.
Dit kon men ook bekomen via de **discordante paren**, dit is, de planten die van bloeistatus veranderden afhankelijk van hun oriëntatie. Inderdaad, **10/50 - 8/50 = 2/50 = 4%**.

Men kan nagaan of de oriëntatie een significant effect heeft op de bloeistatus door een asymptotisch betrouwbaarheidsinterval op te stellen voor het ARV.

```{r}
f=10
g=8
n=50
arv <- (f-g)/n
print(paste0("Absoluut risicoverschil = ", arv))
seARV <- sqrt(f+g - (f-g)^2/n)/n
print(paste0("Standaardfout op het absoluut risicoverschil = ", seARV))

lower <- arv - qnorm(.975)*seARV
upper <- arv + qnorm(.975)*seARV
bi <- c(lower, upper)

cat("95% BI voor ARV = ", bi)
```

Een formele test kan uitgevoerd worden via de Binomiale distributie.
Indien de oriëntatie geen invloed heeft op de bloeistatus, dan zouden er evenveel planten zijn die beginnen bloeien als ze van Zuid naar Oost worden geplaatst, als dat er planten zijn die stoppen met bloeien als ze van Zuid naar Oost worden geplaatst.
Met andere woorden, in termen van de discordante paren, verwacht men dat $f=g$.
Men kan dit ook herschrijven als $f/(f+g) = 1/2$, wat men kan testen via een binomiale test via de `binom.test` functie.
De nul- en alternatieve hypothese zijn:

- $H_0$: De oriëntatie heeft geen invloed op de bloeistatus.
- $H_1$: De oriëntatie heeft invloed op de bloeistatus.

Merk op dat we de nul- en alternatieve hypothese hier causaal kunnen  formuleren aangezien we data gebruiken uit een experimentele studie. 

```{r}
binom.test(x=f, n=f+g, p=0.5)
```

## Conclusie

Er is onvoldoende bewijs in de studie om te kunnen besluiten dat de oriëntatie een significante invloed heeft op de bloeistatus op het 5% significantieniveau (p=`r round(binom.test(x=f, n=f+g, p=0.5)$p.value,2)`).

