9 Logistisk regression

Eksempel: prostata uge 8.

Logistisk regression laver vi vha glm() - Generalized Linear Model. Dette er ikke det samme som General Linear Model, idet sidstnævnte dækker over modeller for kvantitative outcomes hvor kovariaterne kan være både kvantitative og kvalitative (herunder lineær regression, ANOVA, ANCOVA, …). Generaliseret Lineær Model refererer til at vi specificerer modeller som vi plejer ved lineære udtryk (a la a+b*x), men at sammenhængen mellem den parameter vi er interesseret i (for logistisk regression odds eller sandsynlighed) og det lineære udtryk ikke er lineært. For logistisk regression er sammenhængen mellem sandsynligheden \(p_i\) for at outcome for person nr i antager værdien 1 (ud af de mulige udfald 0 og 1)

\[ \log (\frac{p_i}{1-p_i}) = a + b \cdot x_i \]

Som funktion af \(p_i\) er dette ikke lineært i \(p_i\) - men udtrykket på højre side er lineært som funktion af kovariaten \(x_i\). Derfor udtrykket generaliseret lineær model. Funktionen på venstre side af lighedstegnet kaldes link-funktionen, fordi den angiver et ‘link’ mellem det lineære udtryk og sandsyndsynligheden for at vores outcome antager værdien 1. Linkfunktionen her hørende til den logistiske regressionsmodel kaldes også for logit-funktionen: \(\mbox{logit}(p)=\log(\frac{p}{1-p})\)


Som for lm() specificerer vi modellen ved formel notation. På venstre side af en ~ sætter vi outcome / respons (0-1-variabel, f.eks. gennemtraengning01 for prostata-data), på højre side en eller flere forklarende variable adskilt af +’-er - helt efter samme principper som for lm(). Dertil et data-argument og et ekstra argument family='binomial' - sidste argument er ekstremt vigtigt: Glemmer vi det, får vi bare lavet en almindelig lineær regression (dvs helt samme resultater, som hvis vi benyttede lm()).

Dvs syntaks a la

glm1 <- glm( gennemtraengning01 ~ involvering + psa,
              data=prostate, family='binomial')

Responsvariablen gennemtræning01 skal være kodet 0/1, og herved modelleres odds (og dermed også sandsynligheden) for at gennemtraengning01 er 1, dvs odds for gennemtrængning (frem for odds for ikke at have gennemtrængning).

Vi får desværre ikke så meget information ud af glm() men skal selv trække det hele ud. Hertil benytter vi igen og igen funktionerne:

coef( glm1 )             # estimerede koefficienter
summary( glm1 )          # estimerede koefficienter + p-værdier
confint.default( glm1 )  # konfidensintervaller for koefficienterne

Med disse funktioner får vi koefficienterne frem på log-odds (logit)-skala.

Bemærk: Vi benytter her confint.default() til at bestemme CI, hvor vi for lm() benyttede confint(). Med confint.default() får vi bestemt Wald CI (dvs ud fra formlen +/-1.96*SE). Det sikrer os, at konfidensintervallerne stemmer overens med p-værdierne angivet i vores summary. Glemmer vi det, og kommer til at benytte confint(glm1) i stedet får vi stadig konfidensintervaller frem, de er dog beregnet ved en mere avanceret (og bedre …) metode. Men vi vil helst have at CI og p-værdi passer sammen (dumt med en p-værdi som angiver signifikans, CI som ikke angiver signifikans (dette forekommer typisk kun når vi har p-værdier tæt på 0.05)).

Vi får desværre ikke OR og CI for OR frem automatisk, men skal selv bestemme disse ved at benytte antilog (exp()) til de estimerede koefficienter og CI’er:

exp( coef(glm1) )
exp( confint.default(glm1) )

9.1 Modelkontrol

9.1.1 Hosmer-Lemeshow Goodness-Of-Fit test

Hosmer-Lemeshow testet er implementeret i funktionen hoslem.test() i pakken resourceSelection

install.packages("ResourceSelection") # kør kun denne linje 1 gang
library(ResourceSelection)

I testet sammenligner vi i et antal grupper (ofte 10 nogenlunde lige store grupper) det observerede antal med outcome=1 (gennemtraengning01=1) med de prædikterede værdier (dvs ud fra vores model, sandsynlighederne for at observere gennemtraengning01=1 (plejer vi at kalde for p)).

Det angiver vi således

hl <- hoslem.test( prostate$gennemtraengning01, fitted(glm1), g=10)
hl

hvor vi altså har angivet at vi ønsker at sammenligne vores outcome med prædiktionerne fra modellen glm1 (fitted() beregner prædiktionerne for hvert individ på sandsynlighedsskala). Ekstra argument g= specificerer antallet af grupper, her 10.

Med ovenstående får vi blot en p-værdi frem - er denne større end 0.05 er der ikke tegn på at modellen ikke passer, er den mindre end 0.05 bør vi undersøge modellen nærmere.

Advarsel: Hosmer-Lemeshow testet er følsomt over for, hvordan de 10 grupper defineres, og det afhænger af, hvilket software-program man bruger. Derudover kan p-værdien også ændres, hvis man i stedet interesserede sig for odds for den modsatte hændelse (outcome=0 i stedet for outcome=1). Det betyder at man skal tage testet med et gran salt.

9.1.2 Linearitet - residualer

Vi kan kontrollere lineariteten af de kvantitative variable ved at bestemme Pearson-residualerne for hvert individ:

\[ res_i \ = \ \frac{y_i-\hat p_i}{\sqrt{\hat p_i\cdot(1-\hat p_i)}} \]

Disse plottes derefter mod de kvantitative variable, gerne sammen med en loess-smoother (/glidende gennemsnit) så det er lettere at spotte, om der er problemer med lineariteten. F.eks. for prostata-data og vores model glm1, som benytter psa som kvantitativ forklarende variabel:

r <- resid( glm1, type = 'pearson') 
scatter.smooth( r ~ prostate$psa, lpars=list(col='blue',lwd=2))

Er der problemer med lineariteten kan vi vælge at benytte spline-funktioner eller transformere kovariaten (ofte med en logaritme…). Præcis som for den lineære model.

9.1.3 Cook’s afstand

Vi kan her gøre præcis som vi gjorde for Cook for den lineære regressionsmodel, se kapitel 5.7.1:

plot( glm1, which=4)
cooks.distance( glm1 )

9.1.4 Influence

Helt som for den linære model, se kapitel 5.7.2:

dfbetas( glm1 )

9.2 Prædiktion

Vi prædikterer som for den lineære model vha predict(). Når vi laver logistisk regression kan vi prædiktere på to skalaer:

  • logit / log-odds-skala: Her beregnes prædiktionerne ud fra modelformlen, e.g. a+b*psa. Dette kaldes også den linære prædiktor og angives ved at benytte ekstra argument type='link' (link-funktionen henviser her til logit-funktionen, som angiver sammenhængen mellem den lineære prædiktor og sandsynligheden (p))
  • sandsynlighedsskala: Her beregnes prædiktionerne ud fra formlen, exp(a+b*psa)/(1+exp(a+b*psa)). Angives ved at benytte estra argument type='response'.

Vi kan derfor benytte predict() på følgende to måder

predict( glm1, newdata=new, type='link')     # prædiktion log-odds-skala
predict( glm1, newdata=new, type='response') # prædiktion sandsynligheds-skala

Som tidligere skal datasættet new benyttet som argument til newdata indeholde de samme variable som angivet i modelformlen i glm1.

Bemærk at det her ikke giver mening at benytte argumentet interval='prediction' som vi tidligere benyttede for normalfordelte data, fordi vores outcome selvfølgelig ikke kan være normalfordelt, når det er født som et 0-1-outcome. Faktisk virker interval-argumentet slet ikke her, så vi kan ikke engang få beregnet konfidensintervaller for vores prædiktioner som tidligere, hvor vi benyttede interval='confidence'.

9.3 ROC-kurver, AUC

ROC-kurver kan tegnes vha pakken pROC.

install.packages('pROC')
library(pROC)

Derefter tegnes ROC-kurven vha roc()-funktionen. Modellens prædiktioner på sandsynlighedsskala gemmes i en variabel, her p, og vores 0/1-outcome og prædiktionen angives med formelnotation:

p <- predict( glm1, type='response')
roc( prostate$gennemtraengning01 ~ p, percent=F, plot=T )

Ekstra argument percent=F sikrer at skala’en på plottet går fra 0-1 (ellers bliver det 0-100%), plot=T skal angives for at få lavet plottet.

AUC incl CI kan bestemmes ved at gemme resultaterne fra roc()

ROC1 <- roc( prostate$gennemtraengning01 ~ p, percent=F, plot=T )

og dernæst benytte

auc(ROC1)
ci.auc(ROC1)

9.3.1 Optimale tærskelværdier

Youden’s index er bestemt ved det punkt på ROC-kurven, som har den maksimale lodrette afstand ned til diagonalen. Det findes ved funktionen coords (COORDinateS)

coords(ROC1, "best", ret="threshold", best.method="youden",
            transpose=T)
  • Første argument skal være et objekt fra roc()-funktionen, her ROC1
  • I andet argument angiver vi at vi ønsker den optimale ("best") tærskelværdi
  • ret står for RETurn og her fortæller vi, at vi ønsker tærskelværdien
  • metoden angives i best.method
  • transpose=T er for at undgå en overflødig warning

Ønskes i stedet ‘Closest Top Left’, hvor tærskelværdien er bestemt ud fra det punkt på kurven, som ligger tættest på øverste venstre hjørne, benyttes i stedet argument best.method='closest.topleft'