9. Categorische Data Analyse

1 Inleiding

• Vorige hoofdstukken:

  • continue uitkomst
    
  • categorische of continue predictor.

• Nu besluitvorming voor categorische uitkomst.

• Focus op associatie tussen categorische uitkomst en een categorische predictor.

• Gebruik van kruistabellen om associatie voor te stellen.

---

2 Toetsen voor een proportie

2.1 Saksen-studie

• Vrij gesloten populatie (weinig immigratie en emigratie)
• Wat is de kans dat een ongeboren kind mannelijk is?

boys <- 3175
n <- 6155

• Op 6155 ongeboren kinderen werden 3175 jongens geobserveerd.
• Verschil in kans dat ongeboren kind jongen of meisje is?

---

• Gegevens voorstellen als uitkomsten van een numerieke toevalsveranderlijke \(X\)

  • \(X=1\) voor jongens en
  • \(X=0\) voor meisjes.

• Merk op: telprobleem omdat de uitkomst een telling (nl. het aantal jongens)

• Formeel: populatie van ongeboren kinderen waarin elk individu gekenmerkt wordt door een 0 of een 1.

• Uitkomst variabele is binair.

---

2.2 Bernoulli verdeling

• Binaire data \(\longrightarrow\) Bernoulli verdeling: \[\begin{eqnarray*} X_i &\sim& B(\pi) \text{ met}\\ B(\pi)&=&\pi^{X_i}(1-\pi)^{(1-X_i)}, \end{eqnarray*}\]

• Distributie met 1 model parameter \(\pi\)

  • Verwachte waarde van \(X_i\): \(\text{E}[X_i]=\pi,\)
  • Proportie van ongeboren jongens (d.i. kinderen met een 1) in de populatie.
  • Bijgevolg is \(\pi\) kans dat lukraak getrokken individu een jongen is (een observatie die 1 oplevert).

• Variantie van Bernoulli data bepaald door kans \(\pi\). \[\text{Var}[X_i]=\pi (1-\pi).\]

---

Grafische weergave van enkele Bernoulli kansverdelingen

---

• In Saksenstudie worden lukraak 6155 observaties getrokken uit populatie.
• We schatten \(\pi\) als het steekproefgemiddelde : \[\hat \pi = \bar X = \frac{\sum\limits_{i=1}^n X_i}{n},\]

pi <- boys/n
pi

## [1] 0.5158408

In ons voorbeeld is \(\bar x =\) 3175 / 6155 = 51.6%.

---

2.3 Asymptotisch Betrouwbaarheidsinterval

• Saksenstudie heel grote steekproef: 6155 onafhankelijke observaties.

• CLT kan worden toegepast op gemiddelde: data volgt Bernoulli verdeling, maar gemiddelde o.b.v. onafhankelijke en identiek verdeelde observaties in heel grote steekproef volgt approximatief een normaal verdeling.
  

• \(\text{E}[\bar X] = \pi\)

• \(\text{Var}[\bar X] = \frac{\text{Var}[X]}{n} = \frac{\pi(1-\pi)}{n}\)

• \(\bar x\) = 0.516

se <- sqrt((pi*(1-pi))/n)

• SE = 0.0064\[ \text{SE} = \hat\sigma_{\bar x} = \sqrt{\frac{\bar x(1-\bar x)}{n}} \]

• BI: [0.503, 0.528] \[ [\bar x - z_{\alpha/2} \text{SE}_{\bar x}, \bar x + z_{\alpha/2} \text{SE}_{\bar x}]\]

• 0.5 valt niet binnen het 95% BI.

• De kans op een jongen is significant hoger dan 50% op het 5% significantie-niveau.

2.4 Asymptotische Test

• H\(_0\): \(\pi = 0.5\) vs H\(_1\): \(\pi \neq 0.5\)

• Merk op dat voor een Bernoulli verdeling de variantie onder H\(_0\) ook gekend is: \(\pi_0 (1-\pi_0)\).

• Dus onder \(H_0\) is standaard error op x:

\[ \text{SE}_{0, \bar x}=\sqrt{\frac{\pi_0 (1-\pi_0)}{n}} \]

• Statistiek \[ z = \frac{\bar x - \pi_0}{\text{SE}_{0, \bar x}}\]

• Volgt onder de nulhypothese dat er evenveel kans is op een jongen of een meisje (\(\pi_0=0.5\)) asymptotisch een normaal verdeling: CLT

• p-waarde? Merk op dat het een tweezijdige test is!

pi0 <- 0.5
se0 <- sqrt(pi0*(1-pi0)/n)
z <- (pi-pi0)/se0
z

## [1] 2.485539

pval <- pnorm(abs(z),lower=FALSE) *2
pval

## [1] 0.01293554

Er is significant meer kans dat een ongeboren kind mannelijk dan vrouwelijk is (p=0.013). De kans dat een ongeboren kind mannelijk is bedraagt 0.516 (95% BI [0.503, 0.528]).

2.5 Exacte test: Binomiale test

• Statistische toets voor \[H_0: \pi=1/2 \text{ versus } H_1: \pi\neq 1/2,\]

• Daarvoor moeten we verdeling van de

  • \(X\) en \(\bar X\)
  • of van de som \(S=n\bar X\) kennen.

• Als alle observaties onafhankelijk zijn en eenzelfde verdeling \(B(\pi)\) volgen dan volgt \(S\) een binomiaal verdeling.

---

• Stel \(H_0:\pi=1/2\) is waar (voorkomen van jongens en meisjes in populatie even waarschijnlijk)
• Lukrake trekking van één individu uit de populatie, kans op een jongen \[P(X=1)=\pi=1/2.\]
• Twee kinderen onafhankelijk van elkaar (en de populatie \(\approx \infty\)):
  • Kans \(\pi=1/2\) op jongen voor zowel eerste als tweede kind (onafhankelijk van elkaar)
  • Uitkomsten \((x_1, x_2)\) voor beide kinderen hebben dan 4 mogelijke waarden: \((0,0),(0,1),(1,0)\text{ en }(1,1).\)
  • Deze komen elk voor met kans \(1/4 = 1/2 \times 1/2\).
• Toevalsveranderlijke \(S\) die som van uitkomsten weergeeft kan volgende waarden aannemen:

\((x_1,x_2)\)	\(s\)	\(P(S = s)\)
(0,0)	0	1/4
(0,1), (1,0)	1	1/2
(1,1)	2	1/4

---

2.5.1 Algemeen: \(n\) onafhankelijke observaties

• Kans \(\pi\) op “succes” (uitkomst 1) voor elke observatie

• Totaal aantal successen \(S\) (som van alle 1-en) kan \(n+1\) mogelijke waarden hebben \[S=k\text{, met }k=0,\ldots,n\]

• Verdeling van \(S\)? \[\begin{equation} P(S=k) = \left ( \begin{array}{c} n \\ k \\ \end{array} \right ) \pi^k (1-\pi)^{n-k} \end{equation}\]

• \(1-\pi\): kans op mislukking in 1 enkele trekking (uitkomst met 0 genoteerd) en

• binomiaalcoëfficient \[\begin{equation*} \left ( \begin{array}{c} n \\ k \\ \end{array} \right ) = \frac{n \times (n-1) \times ...\times (n-k+1) }{ k!} = \frac{ n!}{ k!(n-k)! } \end{equation*}\]

• In R kan je de kansen van binomiale verdeling voor elke \(S=k\) opvragen met dbinom(k,n,p)

---

2.5.2 Binomiale Verdeling

De toevalsveranderlijke \(S\) volgt:

• Binomiaal verdeelde toevalsveranderlijke met bijhorende Binomiale kansverdeling

• parameters

  • \(n\) (d.i. het aantal trekkingen of, equivalent, de maximale uitkomstwaarde)
    
  • \(\pi\) (de kans op een `succes’ bij elke trekking).

• Kans berekenen dat \(k\) gebeurtenissen zich voordoen op \(n\) onafhankelijke experimenten waarbij kans op 1 zo’n gebeurtenis per experiment, \(\pi\) bedraagt.

• Voor analyse van gegevens die slechts 2 mogelijke waarden kunnen aannemen.

• Bijvoorbeeld: al dan niet besmet met HIV, wild type van een gen vs een mutant,…

• Gebruik: Proporties of risico’s op een gebeurtenis van een bepaald type vergelijken tussen verschillende groepen.

---

2.5.3 Een grafische weergave van enkele Binomiale kansverdelingen.

Binomiale verdelingen.

---

Toetsstatistieken voor \[H_0:\pi=1/2\text{ versus }H_1:\pi\neq 1/2\]

• \(\bar X-1/2\) of, equivalent,
• \(\Delta=n(\bar X-\pi_0)=S-s_0\).
• Verdeling van deze laatste toetsstatistiek volgt rechtstreeks uit de Binomiale verdeling:
  • We observeren \(s=\) 3175 en dus \(\delta=s-s_0=\) 3175 \(-\) 6155 \(\times 0.5=\) 97.5.
  • In veronderstelling dat jongens en meisjes even waarschijnlijk zijn (d.i. onder de nulhypothese \(H_0:\pi=1/2\)), bekomen we de bijhorende tweezijdige p-waarde: \[p=\text{P}_0\left[S-s_0\geq \vert \delta\vert \right] + \text{P}_0\left[S-s_0\leq - \vert \delta\vert \right].\]
• Merk op dat we dit kunnen herschrijven in termen van S. \[p=\text{P}_0\left[S\geq s_0+ \vert \delta\vert \right] + \text{P}_0\left[S \leq s_0 - \vert \delta\vert \right].\]

---

• Voor ons voorbeeld kunnen we deze kansen als volgt berekenen:

\[\begin{eqnarray*} \text{P}_0\left[S\geq s_0+ \vert \delta\vert \right] &=& P(S \geq 6155 \times 0.5 + \vert 3175 - 6155 \times 0.5\vert ) \\&=& P(S \geq 3175)\\ &= &P(S= 3175) + P(S=3176) + ... + P(S=6155)\\ & =& 0.0067\\\\ \text{P}_0\left[S \leq s_0 - \vert \delta\vert \right] &=& P(S \leq 6155 \times 0.5 - \vert 3175- 6155 \times 0.5\vert) \\&=& P(S \leq 2980)\\ &= &P(S=0) + ... + P(S=2980) \\ &=&0.0067 \end{eqnarray*}\]

• Binomiale distributie is symmetrisch als \(\pi=1/2\): \[\text{P}_0\left[S\geq s_0+ \vert \delta\vert \right] = \text{P}_0\left[S \leq s_0 - \vert \delta\vert \right]\]
• Dat is niet langer het geval wanneer \(\pi\) afwijkt van 0.5.

---

pi0 <- 0.5; s0 <- pi0 *n
delta <- abs(boys- s0)
delta

## [1] 97.5

sUp <- s0 + delta
sDown <- s0 -delta
c(sDown,sUp)

## [1] 2980 3175

pUp <- 1-pbinom(sUp-1,n,pi0)
pDown <- pbinom(sDown,n,pi0)
p <- pUp+pDown
c(pUp,pDown, p)

## [1] 0.006699883 0.006699883 0.013399766

• Als \(\pi= 1/2\), kans om door toeval minstens \(\delta=\) 97.5 jongens meer of minder te observeren dan het gemiddelde onder \(H_0: s_0=\) 3077.5 , slechts 1.34% is: de \(p\)-waarde van de binomiale test.
  
• Kleine kans om een dergelijk groot aantal jongens te observeren als in realiteit jongens en meisjes even veel voorkomen in de populatie.
• Drukt uit dat de onderstelling dat jongens en meisjes even veel voorkomen, weinig gesteund wordt door de data.

---

---

De test kan worden uitgevoerd a.d.h.v. de binom.test functie in R.

binom.test(x=boys,n=n,p=pi0)

## 
##  Exact binomial test
## 
## data:  boys and n
## number of successes = 3175, number of trials = 6155, p-value = 0.0134
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.5032696 0.5283969
## sample estimates:
## probability of success 
##              0.5158408

Op het 5% significantie-niveau besluiten we dat er gemiddeld meer kans is dat een ongeboren kind mannelijk dan vrouwelijk is.

• De Binomiale test heeft ook een exact BI weer op een proportie.

• Het exacte BI is te verkiezen boven het BI dat gebaseerd is op de CLT.

• Voor Saksen-studie ligt BI o.b.v. CLT heel dicht bij exacte BI: grote steekproef (\(n=\) 6155).

• Merk op dat het testen voor een proportie kan gezien worden als het equivalent van een one-sample t-test voor binaire data.

---

2.6 Conclusie

• Voor de Saksen populatie besluiten we op het 5% significantieniveau dat er meer kans is dat een ongeboren kind mannelijk dan vrouwelijk is (\(p=\) 0.013). De kans dat een ongeboren kind mannelijk is, bedraagt 51.6% (95% BI [50.3, 52.8]%).

---

3 Toets voor associatie tussen 2 kwalitatieve variabelen

3.1 Gepaarde gegevens

• 2 keer zelfde individu meten
• bijvoorbeeld, vóór en na blootstelling aan de experimentele stof
• telkens de categorische uitkomst te observeren.
• Hier enkel: gepaarde binaire uitkomsten
• Statistische analyse moet rekening houden met de paring.

---

3.1.1 Voorbeeld: partnerkeuze van seksueel mature vrouwelijke Campbelli dwerghamster (Rogovin et al. 2017)

---

• Na 3 minuten, scheidingswand weg

• aggressief vs niet-agressief mannentje

• Elk vrouwtje onderging tweemaal de test: na verblijf in

  • vijandige omgeving (hoge populatie, weinig voedsel, veel concurrentie)
  • vriendelijkere omgeving

Kruistabel van partnerkeuze bij dwerghamster.

	vriendelijk-agressief	vriendelijk-niet-agressief	totaal
vijandig-agressief	3 (e)	17 (f)	20
vijandig-niet-agressief	1 (g)	13 (h)	14
totaal	4	30	34

---

• \(\pi_2=P[\text{agressief mannetje } \vert \text{ verblijf vijandige omgeving}]\)

• \(\hat \pi_2=(e+f)/n\), waarbij \(n=e+f+g+h\).

• \(\pi_1=P[\text{agressief mannetje } \vert \text{ verblijf vriendelijke omgeving}]\)

• \(\hat \pi_1=(e+g)/n\)

• Absoluut riscoverschil (ARV) \[\begin{equation*} \widehat{\text{ARV}}=\hat\pi_2-\hat\pi_1=\frac{e+f}{n}-\frac{e+g}{n}=\frac{f-g}{n} \end{equation*}\]

• Enkel beïnvloed door aantallen discordante paren \(f\) en \(g\)

---

• Standaard error op ARV \[\begin{equation*} \text{SE}_{\widehat{\text{ARV}}}=\frac{1}{n}\sqrt{f+g-\frac{(f-g)^2}{n}} \end{equation*}\]

• Als er voldoende gegevens zijn, kan men een \((1-\alpha)100\%\) BI op ARV \[\left[\widehat{\text{ARV}}-z_{\alpha/2}\text{SE}_{\widehat{\text{ARV}}},\widehat{\text{ARV}}-z_{\alpha/2}\text{SE}_{\widehat{\text{ARV}}}\right]\] of \[\left[\frac{f-g}{n}-\frac{z_{\alpha/2}}{n}\sqrt{f+g-\frac{(f-g)^2}{n}},\frac{f-g}{n}+\frac{z_{\alpha/2}}{n}\sqrt{f+g-\frac{(f-g)^2}{n}}\right] \]

---

hamster <- matrix(c(3,17,1,13),ncol=2,byrow=TRUE)
rownames(hamster) <- c("vijandig-agressief", "vijandig-niet-agressief")
colnames(hamster) <- c("vriendelijk-agressief","vriendelijk-niet-agressief")

f <- hamster[1,2]
g <- hamster[2,1]
n <- sum(hamster)
riskdiff <- (f-g)/n
riskdiff

## [1] 0.4705882

se <- sqrt(f+g-(f-g)^2/n)/n
se

## [1] 0.09517144

bi <- riskdiff + c(-1,1)*qnorm(0.975)*se
bi

## [1] 0.2840556 0.6571208

---

\[\begin{equation*} \widehat{\text{ARV}}=\frac{17-1}{34}=0.471 \end{equation*}\] of 47.1%.

• De standaard error \[\begin{equation*} \text{SE}_{\widehat{\text{ARV}}}=\frac{1}{34}\sqrt{17+1-\frac{(17-1)^2}{34}}=0.0952 \end{equation*}\]

• Een 95% betrouwbaarheidsinterval voor het absolute risicoverschil op de keuze van een agressief mannetje tussen een verblijf in een vijandige en vriendelijke omgeving is bijgevolg \[\begin{equation*} \left[0.471-1.96\times 0.0952,0.471+1.96\times 0.0952\right]=[0.284,0.658] \end{equation*}\]

• We hebben besluiten dus dat absolute risico verschil in het interval [28.4,65.8]% ligt en weten dat dergelijke intervallen voor grote steekproeven het werkelijke absolute risco verschil met een kans van 95% bevatten.

• Op basis van het BI besluiten we op het 5% niveau dat het meer waarschijnlijk is dat vrouwelijke hamsters een aggressieve partner kiezen nadat ze in een vijandige omgeving hebben verbleven dan wanneer ze in een vriendelijke omgeving verbleven. De kans op de keuze van een agressief mannetje is 47% hoger nadat ze in een vijandige omgeving heeft verbleven (95% BI 28.4,65.8]%)

---

3.1.2 McNemar test

	vriendelijk-agressief	vriendelijk-niet-agressief	totaal
vijandig-agressief	3 (e)	17 (f)	20
vijandig-niet-agressief	1 (g)	13 (h)	14
totaal	4	30	34

• Toetsen of de risico’s verschillen tussen de vijandige en vriendelijke omgeving.
• Enkel de discordante paren leveren hier informatie over.
• \(f>g\) indicatie tegen \(H_0\): partnerkeuze niet geassocieerd met omgeving
• Kans evalueren dat in een lukraak discordant paar, vrouwtje na verblijf in een vijandige omgeving kiest voor het agressieve mannetje.
• Deze kans wordt geschat als \[\frac{f}{f+g}\]

---

\[\begin{eqnarray*} \text{E}\left[f/(f+g)\right]&\stackrel{H_0}{=}& 0.5\\ f & \stackrel{H_0}{\sim}& \text{Binom}(n=f+g,\pi=0.5)\\ \text{SE}_{\frac{f}{f+g}} & \stackrel{H_0}{=}& \sqrt{(f+g)\times 0.5\times 0.5}=\frac{\sqrt{f+g}}{2} \end{eqnarray*}\]

---

• Asymptotisch one-sample z-test (o.b.v. normale verdeling)

\[\begin{equation*} z <- \frac{f-(f+g)/2}{\sqrt{f+g}/2}=\frac{f-g}{\sqrt{f+g}} \end{equation*}\]

• De Normale benadering is goed als \[f \times g/(f+g) \geq 5\]

• In kleine steekproeven is het meer aangewezen om een continuïteitscorrectie te gebruiken d.m.v. de toetsingsgrootheid \[\begin{equation*} \frac{|f-g|-1}{\sqrt{f+g}} \end{equation*}\]

De Mc Nemar test analogon van de gepaarde t-test voor binaire, kwalitatieve i.p.v. continue variabelen.

---

In R mcnemar.test functie

mcnemar.test(hamster)

## 
##  McNemar's Chi-squared test with continuity correction
## 
## data:  hamster
## McNemar's chi-squared = 12.5, df = 1, p-value = 0.000407

• We verwerpen bijgevolg de nulhypothese op het 5% significantieniveau en

• Besluiten dat de parternkeuze extreem significant geassocieerd is met de omgeving.

• We zien dat hier eveneens de continuïteitscorrectie werd uitgevoerd en dat we exact dezelfde p-waarde bekomen.

---

• Normale benadering van deze toetstatistiek niet ideaal is omdat \(f \times g/(f+g) < 5\).

• Steeds beter: binomiale test

binom.test(x = f,
  n = f + g,
  p = 0.5)

## 
##  Exact binomial test
## 
## data:  f and f + g
## number of successes = 17, number of trials = 18, p-value = 0.000145
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.7270564 0.9985944
## sample estimates:
## probability of success 
##              0.9444444

---

3.1.3 Conclusie

• Op basis van de exacte test besluiten we eveneens dat de parternkeuze extreem significant geassocieerd is met de omgeving (\(p<0.001\)).
• De kans op de keuze van een agressief mannetje ligt 47.1% hoger als een dwerghamster vrouwtje zich in een vijandige omgeving bevindt dan wanneer ze zich in een vriendelijke omgeving bevindt (95% BI [28.4, 65.7]%).

---

3.2 Ongepaarde gegevens

3.2.1 Genetische associatie studie

• Genetische associatiestudie polymorfismen in het BRCA1 gen geassocieerd is met borstkanker?
• Retrospectieve case-controle studie met 800 borstkankercases en 572 controles

library(tidyverse)

## ── Attaching packages ─────────────────────────────────────────────── tidyverse 1.3.0 ──

## ✔ ggplot2 3.3.2     ✔ purrr   0.3.4
## ✔ tibble  3.0.4     ✔ dplyr   1.0.2
## ✔ tidyr   1.1.2     ✔ stringr 1.4.0
## ✔ readr   1.4.0     ✔ forcats 0.5.0

## ── Conflicts ────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()

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

## 
## ── Column specification ────────────────────────────────────────────────────────────────
## cols(
##   cancer = col_character(),
##   variant = col_character(),
##   variant2 = col_character()
## )

head(brca)

## # A tibble: 6 x 3
##   cancer  variant variant2
##   <chr>   <chr>   <chr>   
## 1 control pro/pro other   
## 2 control pro/pro other   
## 3 control pro/pro other   
## 4 control pro/pro other   
## 5 control pro/pro other   
## 6 control pro/pro other

summary(brca)

##     cancer            variant            variant2        
##  Length:1372        Length:1372        Length:1372       
##  Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character

---

Genotype	Controles	Cases	Totaal
Pro/Pro	266 (a)	342 (d)	608 (a+d)
Pro/Leu	250 (b)	369 (e)	619 (b+e)
Leu/Leu	56 (c)	89 (f)	145 (c+f)
Totaal	572 (a+b+c)	800 (d+e+f)	1372 (n)

• In case-controle studies kiest men een vast aantal cases en controles en spoort men voor hen op welke blootstellingen ze in het verleden ondervonden hebben.
• Dergelijke studies noemt men ook retrospectief
• Onmogelijk om het risico’s and riscoverschillen op borstkanker te schatten: proportie van cases en controles weerspiegelt populatie niet!

---

Genotype	Controles	Cases	Totaal
Pro/Pro	266 (a)	342 (d)	608 (a+d)
Pro/Leu	250 (b)	369 (e)	619 (b+e)
Leu/Leu	56 (c)	89 (f)	145 (c+f)
Totaal	572 (a+b+c)	800 (d+e+f)	1372 (n)

• Wel mogelijk om kans te schatten om allel Leu/Leu
  • cases: \(\pi_1=f/(d+e+f)=89/800=11.1\%\)
  • controles:\(\pi_0=c/(a+b+c)=56/572=9.8\%\)
• Relatief risico op blootstelling voor cases versus controles is bijgevolg \(11.1/9.8=1.14\).
• Vrouwen met borstkanker hebben dus 14% meer kans om de allelcombinatie Leu/Leu te hebben op het BRCA1 gen dan vrouwen zonder borstkanker.
• Dit suggereert een associatie, maar drukt iet uit hoeveel hoger het risico op borstkanker is voor vrouwen met de allelcombinatie Leu/Leu dan voor andere vrouwen
• Andere risicomaat?

---

\[\begin{equation*} Odds=\frac{p}{1-p} \end{equation*}\] waarbij \(p\) de kans is op die gebeurtenis.

Transformatie van het risico, met volgende eigenschappen:

• de odds neemt waarden aan tussen nul en oneindig.

• de odds is gelijk aan 1 als en slechts als de kans zelf gelijk is aan 1/2.

• de odds neemt toe als de kans toeneemt.

• populair bij gokkers: hoeveel waarschijnlijker het is om te winnen dan om te verliezen

---

Genotype	Controles	Cases	Totaal
Pro/Pro	266 (a)	342 (d)	608 (a+d)
Pro/Leu	250 (b)	369 (e)	619 (b+e)
Leu/Leu	56 (c)	89 (f)	145 (c+f)
Totaal	572 (a+b+c)	800 (d+e+f)	1372 (n)

Odds op allel Leu/Leu

• Cases: \(\mbox{odds}_1=\frac{ f/(d+e+f)}{(d+e)/(d+e+f)}=f/(d+e)=89/711=0.125\). Vrouwen met borstkanker hebben ongeveer 8 keer meer kans om de allelcombinatie Leu/Leu niet te hebben op het BRCA1 gen dan om het wel te hebben.

• Controles: \(\mbox{odds}_2=c/(a+b)=56/516=0.109\).

• Associatie tussen blootstelling en uitkomst: \[ OR_{Leu/Leu}=\frac{\mbox{odds}_T}{\mbox{odds}_C}= \frac{f/(d+e)}{c/(a+b)}=\frac{f/(d+e)}{c/(a+b)}=1.15 \]

---

Genotype	Controles	Cases	Totaal
Pro/Pro	266 (a)	342 (d)	608 (a+d)
Pro/Leu	250 (b)	369 (e)	619 (b+e)
Leu/Leu	56 (c)	89 (f)	145 (c+f)
Totaal	572 (a+b+c)	800 (d+e+f)	1372 (n)

• Was de bovenstaande studie echter een volledig lukrake steekproef geweest (waarbij het aantal cases en controles niet per design werden vastgelegd),
• dan konden we daar ook de odds ratio op borstkanker berekenen voor mensen met versus zonder het allel Leu/leu. \[\begin{equation*} OR_{case}=\frac{ \frac{ f}{c}}{ \frac{(d+e)}{(a+b)}} = \frac{f(a+b)}{c(d+e)}=OR_{Leu/Leu}=1.15, \end{equation*}\]
• OR is een symmetrische maat! OR op borstkanker kan wel worden geschat!
• De odds op borstkanker is bijgevolg 15% hoger bij vrouwen met die specifieke allelcombinatie.

---

• Is verschil groot genoeg zodat we het effect die we in de steekproef zien kunnen veralgemenen naar de populatie toe.

• Hiertoe zullen we de kruistabel eerst herschrijven tot een 2x2 tabel

Genotype	Controles	Cases	Totaal
andere	516 (a)	711 (c)	1227 (a+c)
Leu/Leu	56 (b)	89 (d)	145 (b+d)
Totaal	572 (a+b)	800 (c+d)	1372 (n)

---

3.3 De Pearson Chi-kwadraat test voor ongepaarde gegevens

• Testen voor associatie tussen de categorische blootstelling (bvb. variant, X) en de categorische uitkomst (bvb. ziekte, Y). \[H_0: \text{Er is geen associatie tussen } X \text{ en } Y \text{ vs } H_1: X \text{ en } Y \text{ zijn geassocieerd}\]

• Beschouw de rijtotalen \(n_\text{andere}=a+c\), \(n_\text{leu,leu}=b+d\) enerzijds en

• de kolomtotalen \(n_\text{contr}=a+b\) en \(n_\text{case}=c+d\) anderzijds.

• Zij verstrekken informatie over de marginale verdeling van de blootstelling (bvb. variant, X) en de uitkomst (bvb. ziekte, Y), maar niet over de associatie tussen die veranderlijken.

• Onder \(H_0\) zijn \(X\) en \(Y\) onafhankelijk zijn en verwacht men een proportie \((b+d)/n\) van \(a+b\) controles met een Leu/Leu variant, of dat \((a+b)(b+d)/n\) een Leu/Leu variant hebben

• Analoog kan men verwachte aantal \(E_{ij}\) berekenen dat onder de nulhypothese in elke cel van de \(2\times 2\) tabel zou liggen.

---

• \(E_{11}\) = het verwachte aantal onder \(H_0\) in de (1,1)-cel = 145 \(\times\) 800/1372 = 84.55 ;

• \(E_{12}\) = het verwachte aantal onder \(H_0\) in de (1,2)-cel = 145 \(\times\) 572/1372 = 60.45 ;

• \(E_{21}\) = het verwachte aantal onder \(H_0\) in de (2,1)-cel = 1227 \(\times\) 800/1372 = 715.5 ;

• \(E_{22}\) = het verwachte aantal onder \(H_0\) in de (2,2)-cel = 1227 \(\times\) 572/1372 = 511.5 ;

Toetsstatistiek: \[\begin{eqnarray*} X^2 &=& \frac{\left (|O_{11} - E_{11}| - .5 \right)^2 }{ E_{11}} + \frac{ \left ( |O_{12} - E_{12}| - .5 \right)^2 }{E_{12} }+ \\ &&\quad\quad \frac{ \left ( |O_{21} - E_{21}| - .5 \right)^2 }{E_{21}}+ \frac{ \left ( |O_{22} - E_{22}| - .5 \right)^2 }{E_{22} }\\ X^2 &\stackrel{H_0}{\longrightarrow}& \chi^2(df=1) \end{eqnarray*}\]

---

---

• Toetsingsgrootheid groot \(\longrightarrow\) indicatie voor afwijking van \(H_0\).
• De p-waarde voor een 2-zijdige toets is kans om grotere waarde voor toetsingsgrootheid te observeren dan geobserveerde waarde \(x^2\) als nulhypothese waar is \(\longrightarrow\) kans dat \(\chi^2_1\)-verdeelde toevalsveranderlijke waarden groter dan \(x^2\) aanneemt.

---

expected <- matrix(0,nrow=2,ncol=2)
for (i in 1:2)
    for (j in 1:2)
        expected[i,j] <-
            sum(brcaTab2[i,])*sum(brcaTab2[,j])/sum(brcaTab2)
expected

##          [,1]     [,2]
## [1,]  84.5481  60.4519
## [2,] 715.4519 511.5481

x2 <- sum((abs(brcaTab2-expected) - .5)^2/expected)
1-pchisq(x2,1)

## [1] 0.481519

---

• Omdat observaties \(O_{ij}\) discrete getallen zijn, kan toetsingsgrootheid \(X^2\) slechts discrete waarden aannemen en kan continue verdeling zoals \(\chi^2_1\)-verdeling slechts een benadering zijn.
• continuïteitscorrectie: 0.5 aftrekken van verwachte aantallen in elke cel \(\longrightarrow\) om beter aan te laten sluiten bij \(\chi^2_1\)-verdeling
• Pearson Chi-kwadraat toets met Yates correctie.
• Wanneer de correctie niet gebruikt wordt (d.w.z. wanneer de getallen `0.5’ in de uitdrukking voor \(X^2\) door 0 vervangen worden), dan spreekt men van de Pearson Chi-kwadraat toets.

---

In R kan je deze toetsen uitvoeren door de optie op TRUE of FALSE te zetten:

chisq.test(brcaTab2)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  brcaTab2
## X-squared = 0.49542, df = 1, p-value = 0.4815

chisq.test(brcaTab2, correct=FALSE)

## 
##  Pearson's Chi-squared test
## 
## data:  brcaTab2
## X-squared = 0.62871, df = 1, p-value = 0.4278

---

• Zelfs met continuïteitscorrectie is \(\chi^2_1\) benadering slechts verantwoord als in geen enkele van de cellen het verwachte aantal onder \(H_0\) kleiner is dan 5.
• Wanneer de \(\chi^2\)-benadering niet verantwoord is, kan men een Fisher’s exact test uitvoeren.
• De nulhypothese van deze test is eveneens dat \(X\) en \(Y\) onafhankelijk zijn, en de alternatieve hypothese dat \(X\) en \(Y\) afhankelijk zijn.

fisher.test(brcaTab2)

## 
##  Fisher's Exact Test for Count Data
## 
## data:  brcaTab2
## p-value = 0.4764
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.7998798 1.6738449
## sample estimates:
## odds ratio 
##   1.153279

---

3.3.1 Uitbreiding naar categorische variabelen met meerdere niveaus

• \(\chi^2\)-toets kan ook als minstens 1 van de discrete variabelen \(X\) en \(Y\) meer dan 2 mogelijke waarden aanneemt

• Opnieuw: nulhypothese \(H_0: X\) en \(Y\) zijn onafhankelijk (niet-geassocieerd), ten opzichte van het tweezijdig alternatief \(H_A: X\) en \(Y\) zijn niet onafhankelijk (geassocieerd).

• Als de variabele voorgesteld op de rijen \(r\) mogelijke uitkomsten heeft en die op de kolommen \(c\) mogelijke uitkomsten, dan noemt men de kruistabel die \(X\) tegenover \(Y\) uitzet, een \(r \times c\) tabel.

• Zoals voorheen vergelijkt men het aantal geobserveerde waarden in cel \((i,j)\), \(O_{ij}\) genoteerd, met het aantal verwachte waarden onder de nulhypothese, \(E_{ij}\) -Opnieuw is \(E_{ij}\) product van het \(i\)-de rijtotaal met het \(j\)-de kolomtotaal gedeeld door het algemene totaal.

\[\begin{equation*} X^2 = \sum_{ij} \frac{\left (O_{ij} - E_{ij}\right)^2 }{ E_{ij}} \end{equation*}\]

---

• Men kan aantonen dat ze een Chi-kwadraat verdeling volgt met \((r-1) \times (c-1)\) vrijheidsgraden als de nulhypothese waar is.
• De continuïteitscorrectie wordt meestal niet gebruikt bij meer dan 2 rijen of kolommen.
• Pearson \(\chi^2\) test is analogon van de one-way variantie-analyse voor kwalitatieve i.p.v. continue variabelen.

---

brcaTab <- table(brca$variant,brca$cancer)
chisq.test(brcaTab)

## 
##  Pearson's Chi-squared test
## 
## data:  brcaTab
## X-squared = 2.0551, df = 2, p-value = 0.3579

• Om te onderzoeken of het BRCA1 gen geassocieerd is met borstkanker, berekenen we de Pearson chi-kwadraat toets voor de case-controle studie.
• De toetsingsgrootheid bedraagt nu 2.055 en volgt een Chi-kwadraat verdeling met 2 vrijheidsgraden. De kans dat zo’n \(\chi^2\)- verdeelde toevalsveranderlijke extremer is dan 2.055, bedraagt 36%.
• Op het 5% significantieniveau kunnen we dus niet besluiten dat het BRCA1 gen geassocieerd is met borstkanker.

---

4 Logistische regressie

• Raamwerk voor het modelleren van binaire data (vb. kanker vs geen kanker): logistische regressie-modellen.

• Binaire gegevens modelleren a.d.h.v. continue en/of dummy variabelen.

• De modellen veronderstellen dat observaties voor subject \(i=1,\ldots,n\) onafhankelijk zijn en een Bernoulli verdeling volgen.

• Het logaritme van de odds wordt dan gemodelleerd d.m.v. een lineair model, ook wel lineaire predictor genoemd: \[\begin{equation} \left\{ \begin{array}{ccl} Y_i&\sim&B(\pi_i)\\\\ E[Y_i] &=& \pi_i\\\\ \log \frac{\pi_i}{1-\pi_i}&=&\beta_0 + \beta_1X_{i1} + \ldots + \beta_p X_{ip} \end{array}\right. \end{equation}\]

---

4.1 Categorische predictor

• Borstkanker voorbeeld: is BRCA 1 variant geassocieerd is met het krijgen van borstkanker.

• Net zoals in de anova context, factor in het regressieraamwerk d.m.v. dummy variabelen.

• 1 dummy variable minder nodig hebben dan er groepen zijn.

• Voor het BRCA 1 voorbeeld zijn dus twee dummy variabelen nodig en kunnen we de data dus modelleren met onderstaande lineaire predictor:

\[\begin{eqnarray*} \log \frac{\pi_i}{1-\pi_i} &=& \beta_0+\beta_1 x_{i1} +\beta_2 x_{i2} \end{eqnarray*}\]

---

• Waarbij de predictoren dummy-variabelen zijn: \[x_{i1} = \left\{ \begin{array}{ll} 1 & \text{ als subject $i$ heterozygoot is, Pro/Leu variant} \\ 0 & \text{ als subject $i$ homozygoot is, (Pro/Pro of Leu/Leu variant)} \end{array}\right. .\] \[x_{i2} = \left\{ \begin{array}{ll} 1 & \text{ als subject $i$ homozygoot is in de Leucine mutatie: Leu/Leu } \\ 0 & \text{ als subject $i$ niet homozygoot is in de Leu/Leu variant} \end{array}\right. .\]

• Homozygositeit in het wild type allel Pro/Pro wordt voor dit model de referentiegroep.

---

We fitten het model in R.

• Merk op dat we functie as.factor gebruiken om de cancer variabele om te zetten naar een factor.
• en dat we de functie relevel gebruiken om de controles als referentie groep de definiëren: eerste niveau van de factor.

brca <- brca %>%
  mutate(
    cancer = cancer %>%
      as.factor %>%
      relevel("control"),
    variant = variant %>%
      as.factor %>%
      relevel("pro/pro")
  )

brcaLogit <- glm(
  cancer ~ variant,
  data = brca,
  family = binomial)

summary(brcaLogit)

## 
## Call:
## glm(formula = cancer ~ variant, family = binomial, data = brca)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -1.379  -1.286   1.017   1.017   1.073  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)   
## (Intercept)     0.25131    0.08175   3.074  0.00211 **
## variantleu/leu  0.21197    0.18915   1.121  0.26243   
## variantpro/leu  0.13802    0.11573   1.193  0.23302   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1863.9  on 1371  degrees of freedom
## Residual deviance: 1861.9  on 1369  degrees of freedom
## AIC: 1867.9
## 
## Number of Fisher Scoring iterations: 4

Het intercept is de log-odds op kanker in de referentieklasse (Pro/Pro) en de hellingstermen zijn log odds ratio’s tussen de behandeling en de referentieklasse: \[\begin{eqnarray*} \log \text{ODDS}_\text{Pro/Pro}&=&\beta_0\\\\ \log \text{ODDS}_\text{Pro/Leu}&=&\beta_0+\beta_1\\\\ \log \text{ODDS}_\text{Leu/Leu}&=&\beta_0+\beta_2\\\\ \log \frac{\text{ODDS}_\text{Pro/Leu}}{\text{ODDS}_\text{Pro/Pro}}&=&\log \text{ODDS}_\text{Pro/Leu}-\log ODDS_{Pro/Pro}\\ &=&\beta_0+\beta_1-\beta_0=\beta_1\\\\ \log \frac{\text{ODDS}_\text{Leu/Leu}}{\text{ODDS}_\text{Pro/Pro}}&=&\beta_2 \end{eqnarray*}\]

• De analyse laat dus toe om de resultaten onmiddellijk te interpreteren in termen van Odds’es en Odds-ratio’s!

---

anova(brcaLogit, test = "Chisq")

## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: cancer
## 
## Terms added sequentially (first to last)
## 
## 
##         Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL                     1371     1863.9         
## variant  2   2.0562      1369     1861.9   0.3577

De \(\chi^2\)-test op het logistische regressiemodel geeft eveneens aan dat er geen significante associatie is tussen de uitkomst (voorkomen van kanker) en de factor ( de genetische variant van het BRCA gen) (\(p=\) 0.358). De p-waarde is bijna equivalent aan de p-waarde van de \(\chi^2\)-test uit de vorige sectie.

---

• Significante associatie? Post-hoc tests om te evalueren welke odds ratio’s verschillend zijn.
• Voor het BRCA1 voorbeeld zouden we uiteraard geen post-hoc testen
• Toch illustratie zodat jullie over de code beschikken

---

suppressPackageStartupMessages({library(multcomp)})
posthoc <- glht(brcaLogit, linfct = mcp(variant = "Tukey"))
posthocTests <- summary(posthoc)
posthocTests

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: glm(formula = cancer ~ variant, family = binomial, data = brca)
## 
## Linear Hypotheses:
##                        Estimate Std. Error z value Pr(>|z|)
## leu/leu - pro/pro == 0  0.21197    0.18915   1.121    0.493
## pro/leu - pro/pro == 0  0.13802    0.11573   1.193    0.449
## pro/leu - leu/leu == 0 -0.07395    0.18922  -0.391    0.917
## (Adjusted p values reported -- single-step method)

---

posthocBI <- confint(posthoc)
posthocBI

## 
##   Simultaneous Confidence Intervals
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: glm(formula = cancer ~ variant, family = binomial, data = brca)
## 
## Quantile = 2.3232
## 95% family-wise confidence level
##  
## 
## Linear Hypotheses:
##                        Estimate lwr      upr     
## leu/leu - pro/pro == 0  0.21197 -0.22745  0.65139
## pro/leu - pro/pro == 0  0.13802 -0.13084  0.40688
## pro/leu - leu/leu == 0 -0.07395 -0.51353  0.36564

• Door middel van de confint functie worden BI’s verkregen op de log-odds ratios die gecorrigeerd zijn voor multiple testing.

---

• BI’s kunnen als volgt worden teruggetransformeerd naar odds ratios:

OR <- exp(posthocBI$confint)
OR

##                    Estimate       lwr      upr
## leu/leu - pro/pro 1.2361111 0.7965606 1.918210
## pro/leu - pro/pro 1.1480000 0.8773578 1.502128
## pro/leu - leu/leu 0.9287191 0.5983766 1.441432
## attr(,"conf.level")
## [1] 0.95
## attr(,"calpha")
## [1] 2.323183

• De odds ratios die worden bekomen met het logistisch regressiemodel zijn exact gelijk aan de odds ratios die we zouden bekomen op basis van Tabel:

• vb. \(\text{OR}_\text{Leu/Leu-Pro/Pro}=89\times 266/(56\times 342)=\) 1.236.

• Merk op dat de statistische besluitvorming bij logistische modellen beroep doet op asymptotische theorie.

---

4.2 Continue predictor

• Toxicologisch effect van koolstofdisulfide (CS\(_2\)) op kevers.

• De centrale onderzoeksvraag is of de concentratie van CS\(_2\) een effect heeft op de mortaliteit (i.e. kans op sterven) van de kevers?

Design

• 32 onafhankelijk experimenten
• Telkens 1 kever blootgesteld aan één van 8 concentraties (mg/l) van CS\(_2\) voor een gegeven periode.
• De uitkomst van het experiment is: de kever sterft (\(y=1\)) of de kever overleeft (\(y=0\)).

beetles <- read_csv("https://raw.githubusercontent.com/statomics/sbc20/master/data/beetles.csv")

## 
## ── Column specification ────────────────────────────────────────────────────────────────
## cols(
##   dose = col_double(),
##   status = col_double()
## )

head(beetles)

## # A tibble: 6 x 2
##    dose status
##   <dbl>  <dbl>
## 1  169.      1
## 2  169.      0
## 3  169.      0
## 4  169.      0
## 5  172.      1
## 6  172.      0

table(beetles$dose, beetles$status)

##         
##          0 1
##   169.07 3 1
##   172.42 3 1
##   175.52 3 1
##   178.42 2 2
##   181.13 1 3
##   183.69 0 4
##   186.1  0 4
##   188.39 0 4

---

Logistisch regressiemodel dat log odds modelleert in functie van dosis \(x_i\): \[\log \frac{\pi_i}{1-\pi_i}=\beta_0+\beta_1 \times x_i.\]

beetleModel<-glm(
  status~dose,
  data = beetles,
  family = binomial)

summary(beetleModel)

## 
## Call:
## glm(formula = status ~ dose, family = binomial, data = beetles)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.7943  -0.7136   0.2825   0.5177   2.1670  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept) -53.1928    18.0046  -2.954  0.00313 **
## dose          0.3013     0.1014   2.972  0.00296 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 42.340  on 31  degrees of freedom
## Residual deviance: 26.796  on 30  degrees of freedom
## AIC: 30.796
## 
## Number of Fisher Scoring iterations: 5

---

• Intercept: log odds op mortaliteit wanneer er geen \(\text{CS}_2\) gas wordt toegediend.

• Erg lage odds op sterfte (\(\pi/(1-\pi)=\exp(-53.2)\)) en dus op een kans die nagenoeg nul is.

• Merk op: heel sterke extrapolatie: minimum dosis in de dataset 169.07 mg/l.

• Geschatte odds ratio voor dosis effect op de mortaliteitskans is \(\exp(0.3013)=1.35\).

• Dus bij een toename van de dosis CS\(_2\) met 1 mg/l, is de odds ratio voor de mortaliteit \(1.35\).

---

• We besluiten dat dit effect heel significant is (\(p=\) 0.003).
• Een toename in de CS\(_2\) dosis doet de kans op sterven toenemen.

beetlesTab<-table(beetles) %>% data.frame

data.frame(
    grid = seq(
      min(beetles$dose),
      max(beetles$dose),.1)
  ) %>%
  mutate(piHat = predict(beetleModel,
          newdata = data.frame(dose=grid),
          type = "response")
  ) %>%
  ggplot(aes(grid, piHat))+
  geom_line() +
  xlab("dose") +
  ylab("probability (dead)") +
  geom_text(
    aes(x = dose %>%
          as.character %>%
          as.double,
        y = status %>%
          as.character %>%
          as.double,
        label = Freq),
    beetlesTab %>%
      filter(status==0)
  ) +
  geom_text(
    aes(x = dose %>%
          as.character %>%
          as.double,
        y = status %>%
          as.character %>%
          as.double,
       label=Freq),
   beetlesTab %>%
    filter(status==1)
  )

---

Home