11 Logistisk regression

Eksempel: prostata uge 8.

Den Logistiske regressionsmodel går også under betegnelsen 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})\)


Logistisk regression laver vi ved hjælp af proc logistic. Syntaksen er meget lig den vi kender for den generelle lineære model i proc glm.

For prostata-data er outcome-variablen gennemtraengning01, som er kodet 0/1. Vi laver som sædvanligt en model-linje, hvor outcome står på venstre side, forklarende variable på højre. Her skal vi dog angive om vi ønsker at modellere sandsynligheden for at outcome er 0 eller 1. Det gør vi ved at tilføje (event='0') eller (event=1) inden lighedstegnet. Glemmer vi det, benytter SAS pr default (event='0').

Den forklarende variabel involvering er kodet 0/1, dvs det er en dummy-variabel og vi kan derfor vælge om vi vil sætte den ind som kvantitativ eller kvalitativ.

  • Som kvantitativ behøver vi blot:
proc logistic data=prostata;
   model gennemtraengning01(event='1') = involvering;
run;
  • Som kvalitativ skal vi tilføje en class-linje. Her tilføjes param=glm efter en slash, hvorved vi får output på en form, som ligner den vi kender fra proc glm:
proc logistic data=prostata;
   class involvering (ref='0') / param=glm;
   model gennemtraengning01(event='1') = involvering;
run;

Estimater på log-odds skala og p-værdier aflæses i tabellen med titel Analysis of Maximum Likelihood Estimates, odds-ratio’er aflæses i tabellen med titel Odds Ratio Estimates

11.0.1 Rapportering af interaktioner

Interaktioner håndteres som sædvanligt ved at inkludere de to interaktionsvariable med en * imellem. Ønsker vi at kvantificere en interaktion, dvs beregne OR for den ene variabel for fastholdt værdi af den anden, kan vi imidlertid ikke gøre som vi plejer i proc glm med at udelade hovedeffekten af den ene variabel i modellinjen. I stedet skal vi tilføje en oddsratio-linje:

proc logistic data=prostata;
   class involvering (ref='0') alder65 / param=glm;
   model gennemtraengning01(event='1') = involvering alder65 involvering*alder65;
   oddsratio involvering / at (alder65="Under" "Over");
run;

Her angiver vi altså, at vi ønsker at bestemme OR for involvering for patienter “Under” hhv “Over” 65 år. Parentesen kan ikke udelades.

11.1 Modelkontrol

11.1.1 Hosmer-Lemeshow Goodness-Of-Fit test

Hosmer-Lemeshow testet udføres ved at tilføje lackfit i modellinjen:

proc logistic data=prostata;
   model gennemtraengning01(event='1') = psa / lackfit;
run;

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

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

11.1.2 Residualplot

Her plotter vi Pearson-residualerne mod kvantitative kovariater for at vurdere, om der er problemer med lineariteten

\[ res_i \ = \ \frac{y_i-\hat p_i}{\sqrt{\hat p_i\cdot(1-\hat p_i)}} \] Se i kapitel 11.2 nedenfor, hvordan vi får SAS til at beregne residualerne og lave plottene.

11.2 Beregning af prædiktioner og residualer

Vi kan få SAS til at bestemme prædiktioner og residualer for hvert individ ved at benytte et output-statement.

I den logistiske regressionsmodel kan vi prædiktere på to mulige 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)).

Vi kan få SAS til at beregne residualer og prædiktioner på log-odds såvel som sandsynlighedsskala ved at benytte en output-linje. Først angives navnet på det datasæt, som skal indeholde prædiktionerne, i et out= statement. Dernæst angiver vi hvilke beregninger vi ønsker foretaget, og hvad de skal hedde. Se evt SAS’ egen hjemmeside for hvilke muligheder der er her.
Her bestemmer vi prædiktionerne på sandsynlighedsskala (predicted=) og log-odds (xbeta=). Pearson residualerne fås ved stdreschi:

proc logistic data=prostata;
   model gennemtraengning01(event='1') = psa;
   output out=predictions predictions=pred xbeta=log_odds stdreschi=pearson_res;
run;

Dermed indeholder det nye datasæt med navn predictions en kopi af datasættet prostata, med to variable tilføjet: pred og log_odds. Tag gerne et kig på data, inden du plotter residualer, prædiktioner eller log-odds mod kovariaten (her har vi benyttet psa som forklarende variabel):

proc print data=predictions (obs=15);
run;

Residualplot:

proc sgplot data=predictions;
    loess x=psa y=pearson_res;
run;

Plot af prædiktioner på log-odds-skala (punkter og linje):

proc sgplot data=predictions;
    scatter x=psa y=log_odds;
    reg x=psa y=log_odds;
run;

Plot af prædiktioner på sandsynlighedsskala:

proc sgplot data=predictions;
    scatter x=psa y=gennemtraengning01;
    scatter x=psa y=pred;
run;

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.

11.2.1 Cook’s afstand og influence

I output-linjen beskrevet i kapitel 11.2 tilføjes dfbeta=_all_ hvorved vi får tilføjet Dfbeta for alle koefficienter i modellen.

Cooks distance virker ikke i nyeste version af SAS (i hvert fald ikke hos mig …) tidligere kunne man i output linjen også angive cooksd. Uden Cook må vi nøjes med at se på Dfbetas.

proc logistic data=prostata;
   model gennemtraengning01(event='1') = psa;
   output out=diagnostics dfbeta=_all_;
run;

I datasættet her navngivet diagnostics får vi tilføjet to variable til vores prostata-datasæt med navne DFBETA_Intercept og DFBETA_psa. Disse kan nu plottes mod f.eks. et id-nummer (hvis vi havde det) eller PSA-værdierne. Vælg gerne forskellig farve afhængigt af outcome, så det er nemt at se, hvilke observationer, der passer dårligt med modellen:

proc sgplot data=diagnostics;
    scatter x=psa y=dfbeta_psa / group=gennemtraengning01;
run;

11.3 ROC-kurver, AUC

Vi kan få tegnet en ROC kurve ved at tilføje plots(only)=roc i proc-linjen:

proc logistic data=prostata plots(only)=roc;
        model gennemtraengning01(event='1') = psa;
run;

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

‘Closest Top Left’ bestemmer tærskelværdien ud fra det punkt på kurven, som ligger tættest på øverste venstre hjørne.

Se SAS programmet hørende til slides uge 8 hvordan man finder disse størrelser, det kræver en lille smule programmering (du finder programmet i kapitel 14 her i e-bogen).