Logistische regressie

In de statistische literatuur bestaat ook een raamwerk voor het modelleren van binaire data (vb. kanker vs geen kanker): logistische regressie-modellen. Net zoals we bij het regressiemodel voor Normale continue gegevens hebben gezien in Hoofdstuk 6, 7 en 10, laten logistische regressie modellen toe om binaire gegevens te modelleren a.d.h.v. continue en/of dummy variabelen.

De modellen veronderstellen dat de 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)\\\\ \log \frac{\pi_i}{1-\pi_i}&=&\beta_0 + \beta_1X_{i1} + \ldots + \beta_p X_{ip} \end{array}\right. \end{equation}\]

Categorische predictor

We illustreren logistische regressie aan de hand van het borstkanker voorbeeld waar we bestuderen of de BRCA 1 variant geassocieerd is met het krijgen van borstkanker (zie Sectie 9.3.2).

Net zoals in de anova context, kunnen we een factor in het regressieraamwerk introduceren door gebruik te maken van dummy variabelen. We zullen hierbij 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 variant} \\ 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.

Het model wordt als volgt in R gefit:

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!

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 kanker en de genetische variant.

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

Als er een significante associatie was geweest, dan hadden we post-hoc tests kunnen uitvoeren om te evalueren welke odds ratio’s verschillend zijn. Voor het BRCA1 voorbeeld zouden we uiteraard geen post-hoc testen uitvoeren omdat de globale $\chi^2$-test voor het testen van associatie niet significant is. We illustreren de post-hoc analyse toch zodat jullie over de code beschikken voor het uitvoeren van de analyses.

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.3245
## 95% family-wise confidence level
##  
## 
## Linear Hypotheses:
##                        Estimate lwr      upr     
## leu/leu - pro/pro == 0  0.21197 -0.22770  0.65164
## pro/leu - pro/pro == 0  0.13802 -0.13099  0.40704
## pro/leu - leu/leu == 0 -0.07395 -0.51378  0.36589

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:

OR <- exp(posthocBI$confint)
OR

##                    Estimate       lwr      upr
## leu/leu - pro/pro 1.2361111 0.7963615 1.918690
## pro/leu - pro/pro 1.1480000 0.8772237 1.502358
## pro/leu - leu/leu 0.9287191 0.5982270 1.441792
## attr(,"conf.level")
## [1] 0.95
## attr(,"calpha")
## [1] 2.324504

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 13: 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. Ze is dus enkel correct voor grote steekproeven.