Gå till index

Analys med R

0% färdig
0/0 Steps
  1. Analys och forskning med R och Posit (Rstudio)
  2. Grunderna i R och Rstudio
    7 Ämnen
  3. Importera, exportera, spara och ladda data
    5 Ämnen
  4. Strängar och regular expressions (regex)
    1 Ämne
  5. Bearbetning av data med dplyr
    12 Ämnen
  6. Visualisera och presentera
    14 Ämnen
  7. Explorerande analyser
    6 Ämnen
  8. Prediktionsmodeller
    12 Ämnen
  9. Klassisk regressionsanalys
    8 Ämnen
  10. Machine learning (ML) och Artificiell Intelligens (AI)
    9 Ämnen
  11. Prediktionsmodeller: Tidymodels
  12. Hypotestester
    1 Ämne
Avsnitt 9, Ämne 4
Startad

Logistisk regression

Avsnitt Progress
0% färdig

Logistisk regression

Synonym: logit regression, logit model

Med linjär regression kan vi undersöka kontinuerliga utfall (\(Y\)), exempelvis kroppsvikt, blodtryck, inkomst eller pris. Linjär regression bör dock inte användas för att studerar binära utfall, dvs utfall med två kategorier. Exempel på binära utfall är sådana som klassificeras som ja/nej, 0/1, sjuk/frisk, bra/dålig, nöjd/missnöjd, lite/mycket, etc.

Det finns flera anledningar till att linjär regression inte bör användas för att studera varken binära utfall eller kategoriska utfall med fler än 2 kategorier. Den första anledningen är konceptuell och kan demonstreras med ett kategoriskt utfall med tre nivåer. Ponera att vi studerar sambandet mellan prediktorn \(X\) och en diagnos \(Y\), vilken kan ha en av följande kategorier: stroke, epilepsi eller Parkinsons sjukdom. För att använda linjär regression kodar vi dessa diagnoser med siffrorna 1 till 3. Vi definierar vårt utfall som följer:

\(
Y = \begin{cases} 1 = {stroke} \\ 2 = {epilepsi} \\ 3 = {parkinsons}
\end{cases}
\)

Det finns flera problem med att definiera utfallet på detta vis. Först och främst hävdar detta att det finns en inneboende ordning mellan diagnoserna, vilket det inte gör i detta fall. Dessutom innebär det att vi hävdar att skillnaden mellan stroke och epilepsi (numerisk skillnad på 1) är lika stor som skillnaden mellan epilepsi och Parkinsons sjukdom (numerisk skillnad på 1), vilket är en irrationell jämförelse.

Om syftet med regressionsmodellen är att studera hur diverse prediktorer (\(X\)) relaterar till risken för dessa diagnoser (Y) så leder den linjära modellen till orimliga resultat. För att illustrera detta skapar vi en linjär regression med data från PimaIndiansDiabetes, som finns i paketet mlbench. Paketet kan installeras med kommandot install.packages("mlbench"). I PimaIndiansDiabetes är variabeln diabetes kodad som neg (negative) eller pos (positive), vilket vi kodar om till 0 respektive 1 för att kunna göra linjär regression. I just detta fall skapar vi den linjära regressionen direkt i ggplot() (vanligtvis används funktionen lm()). Med ggplot() ritar vi ett spridningsdiagram med glucose på x-axeln och diabetes på y-axeln, därefter använder vi geom_smooth() för att skapa en linjär regression i data.

R
library(mlbench)
data("PimaIndiansDiabetes")

# Kodar om diabetes till 0 / 1
PimaIndiansDiabetes$diabetes <- ifelse(PimaIndiansDiabetes$diabetes=='pos', 1, 0)

# Linjär regressionsmodell
PimaIndiansDiabetes %>%
  ggplot(aes(glucose, diabetes)) +
  geom_point(alpha = 0.5, size=3) +
  geom_smooth(method = "lm") +
  labs(title = "Linjär regression med binärt utfall", 
       x = "Glucose",
       y = "Sannolikhet för diabetes")

Ovanstående graf är problematisk. Först och främst predicerar modellen negativa probabiliteter vid glucose <75. En probabilitet kan aldrig vara negativ. Dessutom skulle modellen predicera en probabilitet >1 vid höga värden (>200) på glucose. En probabilitet kan aldrig överstiga 1.0. Det är således tydligt att linjär regression är olämpligt för kategoriska utfall.

Den matematiska lösningen på ovanstående problem är att skapa en regressionsmodell som alltid producerar en prediktion mellan 0 och 1, oavsett värdet på prediktorn \(X\). För att göra detta måste utfallet hanteras som en sannolikhet (\(p\), probabilitet) vilken alltid är mellan 0 och 1. Detta kan åstadkommas med funktionen logit eller log odds \(𝑔(𝑝)\):

\(g(p) = \log\left(\frac{p}{1-p}\right)\text{ där }~ 0<p<1\)

I ovanstående formel är \(\frac{p}{1-p}\) ekvivalent med odds för utfallet. För att förstå den fortsatta diskussionen är begreppet odds viktigt.

Odds

Odds är centralt för logistisk regression. I en logistisk regression är utfallet binärt (0/1). Vi kan betrakta detta som en situation där en händelse kan utebli (0) eller inträffa (1). Den logistiska modellen beräknar odds för att händelsen skall inträffa. Odds är sannolikheten för att något skall inträffa dividerat med sannolikheten för att det inte skall inträffa. Nedan visas hur odds för ett utfall definierat som ja/nej, där \(p\) är sannolikheten.

\(odds = \frac{p(ja)}{p(nej)}\)

Detta är ekvivalent med:

\(odds = \frac{p}{1-p}\)

Exempel 1: Vad är oddsen för att få en jämn siffra (2, 4 eller 6) vid ett kast med tärning?

Svar: En tärning har sex siffror (1, 2, 3, 4, 5, 6). Sannolikheten för att få 2, 4 eller 6 är 3/6. Sannolikheten för att inte få 2, 4 eller 6 är också 3/6. Odds beräknas som följer:

\(odds = \frac{3/6}{3/6} = \frac{3}{6} \cdot \frac{6}{3} = \frac{18}{18} = 1\)

En odds på 1.0 innebär att det är lika stor chans att utfallet inträffar som att utfallet inte inträffar.

Exempel 2: Vad är oddsen för att få någon av siffrorna 1, 2, 3 eller 4 vid ett kast med tärning?

Svar: Sannolikheten för att få 1, 2, 3 eller 4 är 4/6 och sannolikheten för något annat resultat är 2/6. Odds beräknas som:

\(odds = \frac{4/6}{2/6} = \frac{4}{6} \cdot \frac{6}{2} = \frac{24}{12} = 2\)

Oddsen är således 2.0, vilket innebär att det är dubbelt så hög sannolikhet för att händelsen skall inträffa, jämfört med att den inte skall inträffa.

Vi återgår nu till diskussionen om logit.

I nedanstående figur visualiseras \(p\) mot logit (log odds). Figuren visar att probabiliteter kan översättas till vilken odds som helst med hjälp av logit.

Vi kan beräkna probabiliteter utifrån log odds med hjälp av den logistiska funktionen (logistic function):

\(g^{-1}(x) = \frac{e^{x}}{1+e^{x}}\)

Denna formelns funktion kan visualiseras som framgår av nedanstående figur. Observera att sannolikheten (y-axeln) är mellan 0 och 1, oavsett värdet på prediktorn \(X\).

Genom att använda logit för att transformera probabiliteten får vi ett värde som kan prediceras med en linjär ekvation, som följer:

\(g(p) = \log\big(\frac{p}{1-p}\big) = \beta_0 + \beta_1x\)

Logit funktionen gör det alltså möjligt att länka \(p\) till \(\beta_0 + \beta_1x\). Logit kan därför betraktas som en länkfunktion (link function). I logistisk regression är alltså logaritmen av odds en linjär funktion av prediktorerna.

Om denna funktionen används på samma data som ovan (prediktion av diabetes), så får regressionsmodellen följande funktion:

R
PimaIndiansDiabetes %>%
  ggplot(aes(glucose, diabetes)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "glm", method.args = list(family = "binomial")) +
  labs(title = "Logistisk regression med binärt utfall",
       x = "Glucose",
       y = "Sannolikhet för diabetes")

Denna modellen kommer alltid producera sannolikheter mellan 0 och 1.

Tolkning av logistisk regression

Regressionskoefficienten (\(\beta\)) representerar hur logaritmen av odds för händelsen ändras när prediktorn \(X\) ökar med 1 enhet. Detta påminner om linjär regression, där \(\beta\) representerar hur \(Y\) ändras när prediktorn \(X\)ökar med 1 enhet. Vid logistisk regression är det alltså logaritmen av odds för händelsen som ändras.

Odds ratio (oddskvot)

Sannolikhet (probabilitet) är ett svårgreppat koncept. Att tolka logaritmerade sannolikheter är ännu svårare. Därför transformerar vi regressionskoefficienten genom att exponentiera den och därmed skapa en oddskvot (odds ratio, OR). Odds ratio är betydligt enklare att tolka. Exempel följer:

  • För prediktorn \(X\) erhåller vi en regressionskoefficient på 0.5, vilket innebär att när prediktorn \(X\) ökar med 1 enhet så ökar logaritmen av odds för att händelsen skall inträffa med 0.5.
  • Vi exponentierar koefficienten: \(e^{0.5}=1.65\)
  • Detta innebär att när prediktorn \(X\) ökar med 1 enhet så ökar odds för händelsen med 65%.

Man kan alltså omvandla en oddskvot till procent! Detta görs med formeln: \(100 \cdot (OR–1)\).

Exempel 1:
Odds ratio = 1.10
Beräkning: 100•(1.10 – 1) = 10%.
Tolkning: För varje enhet prediktorn X ökar så ökar odds för händelsen med 10%.
Exempel 2:
Odds ratio = 0.90
Beräkning: 100•(0.9 – 1) = –10%.
Tolkning: För varje enhet prediktorn X ökar så minskar odds för händelsen med 10%.

En odds ratio (oddskvot) på 1.0 innebär att prediktorn \(X\) inte har ett samband med \(Y\) eftersom en ökning/minskning av \(X\) inte påverkar odds för händelsen \(Y\). Om odds ratio är större än 1.0 så ökar odds för händelsen och om odds ratio är mindre än 1.0 så minskar odds för händelsen. Odds ratio på 1.50 innebär 50% ökad odds och odds ratio 0.5 innebär 50% minskad odds.

När man använder odds ratio så är procent-ändringen i odds samma varje gång prediktorn \(X\) ökar med 1 enhet. Detta beror på att den logistiska regressionen förutsätter att förhållandet mellan \(X\) och log(odds) är linjärt.

Exempel på logistisk regression

Vi skapar en logistisk regression med data från PimaIndiansDiabetes, som finns i paketet mlbench. Installera paketet med kommandot install.packages("mlbench"). I data är variabeln diabetes kodad som pos (positive) eller neg (negative). Den logistiska regressionen skapas med hjälp av funktionen glm() som finns i grundinstallationen av R. GLM står för generalized linear model, vilken estimerar linjära modeller med hjälp av en länkfunktion (link function). I logistisk regression är log odds (logit) länkfunktionen. Modellen specificeras som följer:

  • diabetes="pos" ~ glucose + age innebär att det är kategorin "pos" vi vill predicera med hjälp av variablerna glucose och age.
  • family = binomial(link = logit) specificerar vilken länkfunktion GLM skall använda.
R
library(mlbench)
data("PimaIndiansDiabetes")

min_modell <- glm(diabetes=="pos" ~ glucose + age,
                  family = binomial(link = logit),
                  data=PimaIndiansDiabetes)

Resultatet av modellen med base R funktionen summary():

R
summary(min_modell)
Resultat
Call:
glm(formula = diabetes == "pos" ~ glucose + age, family = binomial(link = logit), 
    data = PimaIndiansDiabetes)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.3367  -0.7775  -0.5087   0.8367   3.1630  

Coefficients:
             Estimate Std. Error z value             Pr(>|z|)    
(Intercept) -5.912449   0.462620  -12.78 < 0.0000000000000002 ***
glucose      0.035644   0.003290   10.83 < 0.0000000000000002 ***
age          0.024778   0.007374    3.36             0.000778 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 993.48  on 767  degrees of freedom
Residual deviance: 797.36  on 765  degrees of freedom
AIC: 803.36

Number of Fisher Scoring iterations: 4

Under rubriken Coefficients presenteras regressionskoefficienterna. Båda koefficienterna (se Estimate för glucose och age) är positiva (>0). Dessutom är P-värdena <0.005, vilket innebär att koefficienterna är statistiskt signifikanta.

I en vetenskaplig publikation är man oftast intresserad av att presentera odds ratio (OR) med 95% konfidensintervall. För att göra detta använder vi funktionen tidy från paketet broom. Erinra att koefficienten är log odds för Y, vilket innebär att koefficienten måste exponentieras för att erhålla odds ratio. Med kommandot nedan begär vi både exponentiering och 95% konfidensintervall:

R
library(broom)
tidy(min_modell, exponentiate=T, conf.int=T)
Resultat
# A tibble: 3 × 7
  term        estimate std.error statistic  p.value conf.low conf.high
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
1 (Intercept)  0.00271   0.463      -12.8  2.11e-37  0.00106   0.00653
2 glucose      1.04      0.00329     10.8  2.37e-27  1.03      1.04   
3 age          1.03      0.00737      3.36 7.78e- 4  1.01      1.04  

Resultatet ovan tolkas som följer:

  • Estimate är odds ratio för händelse Y.
    • När glucose ökar med 1 enhet så ökar risken för diabetes med 4%.
    • När age ökar med 1 enhet så ökar risken för diabetes med 3%.
  • conf.low och conf.high inderikerar ett 95% konfidensintervall. Med 95% säkerhet är odds ratio mellan 1.03 och 1.04 för glucose, vilket innebär att 1 enhet ökning i glucose ger med 95% säkerhet 3% till 4% högre sannolikhet för diabetes.

Modeldiagnostik kan erhållas effektivt med funktionen check_model() från paketet easystats:

R
library(easystats)
check_model(min_modell)

Predicera med logistisk regression

Vi kan göra prediktioner med regressionsmodellen med hjälp av funktionen predict(). Här följer två exempel där vi predicerar på en enda person genom att skapa en ny dataframe med personens värde på glucose och age:

R
# Definiera en ny observation och kalla objektet "new_patient"
new_patient = data.frame(age=50,
                         glucose=150)

# Predicera sannolikheten för diabetes
predict(min_modell, newdata=new_patient, type="response")
Resultat
0.00422564

Ovanstående siffra är sannolikheten (\(p\)) för att personen har diabetes.

ROC kurva för modellen

Att utvärdera modellens prediktiva förmåga görs som regel med hjälp av en ROC-kurva. Här nedan skapar vi tre varianter av modellen och jämför de med hjälp av ROC. Observera dock att en jämförelse av olika modeller bör helst göras med någon typ av resampling, förslagsvis korsvalidering (se kapitel Utvärdera och Jämföra Modeller).

R
# En modell med glucose och age
min_modell_1 <- glm(diabetes=="pos" ~ glucose + age,
                  family = binomial(link = logit),
                  data=PimaIndiansDiabetes)

# En modell med glucose, age, insulin och pressure
min_modell_2 <- glm(diabetes=="pos" ~ glucose + age + insulin + pressure,
                  family = binomial(link = logit),
                  data=PimaIndiansDiabetes)

# En modell med samtliga tillgängliga prediktorer i data
# Punkt (.) är kortkommando för "samtliga prediktorer"
min_modell_3 <- glm(diabetes=="pos" ~ .,
                  family = binomial(link = logit),
                  data=PimaIndiansDiabetes)

library(easystats)
plot(performance_roc(min_modell_1, min_modell_2, min_modell_3))

Multipel logistisk regression

Som framgår i exemplen ovan kan man medbringa multipla prediktorer i den logistiska modellen och den accepterar både kategoriska och kontinuerliga prediktorer. Om modellen inbegriper flera prediktorer så tolkas koefficienten för varje enskild prediktor som dennes "effekt" när övriga prediktorer hålls konstanta. "Effekten" av prediktorn \(X\) är således den som observeras när justering görs för övriga prediktorer i modellen.

Maximum Likelihood Estimation (MLE): Beräkning av regressionskoefficienter i logistisk regression

Erinra att den linjära modellen har en enkel formel för att hitta modellens parametrar (intercept och regressionskoefficienter). Ekvationen är den regressionslinje som bäst beskriver data. Detta görs genom att minimera summan av alla errors (se Linjär Regression). För logistisk regression finns ingen enkel ekvation att lösa. Istället måste mjukvaran hitta parametrarna med hjälp av optimeringsfunktioner (optimizer functions). Den vanligaste funktionen för att göra detta är MLE (maximum likelihood estimation).

Maximum likelihood estimation (MLE, maximimetoden) används generellt för att beräkna parametrarna i en sannolikhetsfördelning. Metoden väljer de parametrar som maximerar sannolikheten för att erhålla de observerade utfallen. Metodern kräver att en probabilitetsdistribution för \(X\) används, vilket är Bernoullifördelning i fallet med logistisk regression. Givet denna fördelning kan en likelihood funktion beräknas för observerad i data. Funktionen kan sedan optimeras för att hitta parametrarna som bäst beskriver data. Den exakta beräkningen av dessa parametrar ligger utanför ramen för denna kurs.