Regression

Inom all forskning och analys är det mycket vanligt att man samlar in information om flera variabler för att undersöka sambanden mellan dessa. Det kan innebära att man samlar in socioekonomiska data (inkomst, utbildning, etnicitet, etc), antropometriska data (vikt, längd, BMI, etc), sjukdomshistoria (förekomst av diabetes, reumatism, lungsjukdom, cancer, etc), biomarkörer (kolesterol, CRP, cancermarkörer etc) och andra mätvärden (blodtryck, ålder etc). Om man vill undersöka sambandet mellan dessa variabler så är regression ett kraftfullt verktyg. Regression tillåter oss göra följande:

• Studera relationen (associationen) mellan en beroende variabel (\(Y\)) och en eller flera oberoende variabler (prediktorer, \(X\)).
• Studera hur stark relationen (associationen) är mellan varje enskild prediktor (\(X\)) och den beroende variabeln (\(Y\)).
• Skapa en prediktionsmodell som kan användas för att predicera, vilket innebär att vi baserat på information om \(X\) kan förutsäga \(Y\).

I en regressionsanalys prediceras en variabel (\(Y\)) av minst en prediktor (\(X\)). Rent terminologiskt sägs \(Y\) vara beroende av \(X\); \(Y\) är därför den beroende variabel och \(X\) är oberoende variabel.

Prediktions vs. effektestimering

Som nämnt ovan kan regression användas för att studera relationen mellan variabler, studera sambandens styrka och/eller utföra prediktioner. I praxis innebär det att vi antingen studerar hur en eller flera prediktorer är associerade med utfallsmåttet, eller så skapar vi en prediktionsmodell för att kunna förutsäga ett värde eller en händelse (\(Y\)). Man bör särskilja regressionsmodeller som syftar till prediktion (dvs göra prediktioner) från modeller som syftar till effektestimering (dvs studera sambandens styrka och utseende).

Prediktion – Ponera att vi har samlat in data på 10000 personer och följt dessa under flera år för att studera insjuknandet i cancer. För varje person har vi samlat in uppgifter om bland annat ålder, kön, yrke, socioekonomi, rökning, sjukdomar, läkemedelsanvändning och boendesituation. Om syftet är att kunna förutsäga vilka personer som kommer utveckla lungcancer, så är syftet med regressionsanalysen att göra prediktioner.

Effektestimering – Ponera att vi har samlat in data på 10000 personer och följt dessa under flera år för att studera insjuknandet i lungcancer. För varje person har vi samlat in uppgifter om bland annat ålder, kön, yrke, socioekonomi, rökning, sjukdomar, läkemedelsanvändning och boendesituation. Om vårt syfte är att studera relationen mellan (exempelvis) rökning och risken för cancer så är syftet med regressionsanalysen att estimera effekten av en prediktor (rökning). Då är det viktigt att sambandet mellan rökning och cancer inte störs av andra faktorer. Det kan exempelvis vara så att rökare tenderar att ha lägre socioekonomisk status, som också påverkar risken för cancer. Det vore därför önskvärt att ta höjd för – med andra ord justera för – socioekonomisk status. Detta kan göras genom att inkludera socioekonomisk status som en prediktor i modellen. När man konstruerar regressionsmodeller för effektestimering är det viktigt att beakta confounders (se kapitel Kausal Inferens) av denna typen.

Simpel linjär regression: en beroende variabel och en prediktor

Detta är en kort repetition av kapitlet Linjär Regression.

Den vanliga regressionsekvationen ser ut som följer:

\(Y = b_1 X_1 + b_0 + \epsilon_i\)

• \(b_0\) är interceptet, vilket är punkten där regressionslinjen skär Y-axeln (se figur ovan).
• \(b_1\) är koefficienten för mass, vilket innebär att det är värdet på regressionslinjens lutning.
• \(\epsilon_i\) är skillnaden mellan modellens predicerade värde på \(Y\) och det observerade värdet på \(Y\). Detta är alltså modellens felprecision (error).

Följande figur visar sambandet mellan prediktorn mass (kroppsvikt) och utfallet glucose (blodsocker).

Ovanstående regressionslinje visar ett positivt samband mellan \(X\) och \(Y\); om \(X\) ökar så ökar \(Y\). Det visuella intrycket av ovanstående figur är således att \(X\) kan användas för att predicera \(Y\). Vi kan nu skriva om ekvationen för linjär regression med variablerna mass och glucose:

\(glucose_i = b_1 mass_i + b_0 + \epsilon_i\)

Om vi önskar analysera sambandet mellan mass och glucose justerat för ålder så krävs multipel regression. Det innebär att multipel regression används för att studera sambandet mellan minst två prediktorer och ett utfall.

Multipel regression

Multipel regression är en modell med multipla, minst 2, prediktorer. Den matematiska formeln för multipel regression med två prediktorer ser ut som följer:

\(Y = \beta_1 X_1 + \beta_2 X_2 + \beta_0 + \epsilon_i\)

Detta samband kan inte längre presenteras med en rak regressionslinje. Sambandet mellan \(Y\) och prediktorerna \(X_1\) och \(X_2\) utgör nu ett tredimensionellt regressionsplan, som exemplifieras i ovanstående figur.

Statistisk signifikans: är sambandet mellan X och Y statistiskt signifikant?

En central fråga är om sambandet mellan \(X\) och \(Y\) är statistiskt signifikant. Denna fråga kan delvis besvaras genom en visuell bedömning av det grafiska sambandet. Visuella bedömningar är dock subjektiva och det är dessutom svårt att bedöma tredimensionella regressionsplan (eller regressionsplan med ännu fler dimensioner [variabler]). Det behövs ett formellt mått på huruvida sambandet mellan \(X\) och \(Y\) är signifikant. Detta kan göras inom ramen för hypotestestning, dvs genom att etablera en nollhypotes. I detta sammanhang är nollhypotesen att det inte finns något samband mellan \(X\) och \(Y\), vilket innebär att lutningen (koefficienten) för prediktorn \(X\) är lika med 0. Om lutningen är 0, så finns inget samband mellan \(X\) och \(Y\). Ett sådant samband illustreras i nedanstående figur.

Signifikanstestning med t-test (t-statistic)

För att undersöka om det finns ett statistiskt signifikant samband mellan \(X\) och \(Y\) måste vi testa nollhypotesen. Nollhypotesen är att koefficienten för prediktorn inte har någon lutning, vilket innebär att värdet på \(Y\) inte ändras om X ökas eller minskas. Uppgiften blir att beräkna probabiliteten (sannolikheten) för att nollhypotesen stämmer. Om den sannolikheten är <5% (p<0.05) så förkastar vi nollhypotesen och säger därmed att det finns ett samband, dvs \(X\) är en statistiskt signifikant prediktor.

Nollhypotesen kan testas med ett t-test. Värdet som beräknas kallas t-statistic och det beräknas genom att koefficienten (\(\beta\)) divideras med dess standard error. Det aktuella t-värdet kan sedan översättas till ett p-värde. Alla statistikprogram, inklusive R, gör dessa beräkningar automatiskt. Det innebär att när man skapar en regressionsmodell så erhåller man per automatik t-tester och p-värden för samtliga prediktorer i modellen.

Härnäst skall vi skapa en regressionsmodell och tolka resultaten i detalj.

Övningsexempel med multipel regression

Här följer två övningar med multipel regression och tonvikten ligger på tolkning av resultaten. Övningarna genomförs i R. Du kan kopiera koden nedan till R och göra samma körningar.

Vi kommer använda flera paket för att underlätta och effektivisera arbetet. Multipel linjär regression finns i base R (standardinstallationen av R). Funktionen heter lm(). Därutöver behöver vi paketen sjPlot, broom och tidyverse för att extrahera och presentera resultaten. Vi kommer använda PimaIndiansDiabetes data (från paketet mlbench), som innehåller information om Pimaindianer, vilka har hög risk för diabetes. För att installera och aktivera paketen, samt data skrives följande kommandon:

BLOCKERAD KOD, BILD ELLER TEXT


Du måste logga in för att komma åt allt material.

pregnant glucose  pressure triceps insulin mass  pedigree age diabetes
6	 148	  72	   35	   0	   33.6	 0.627	  50  pos
1	 85	  66	   29	   0	   26.6	 0.351	  31  neg
8	 183	  64	   0	   0	   23.3	 0.672	  32  pos
1	 89	  66	   23	   94	   28.1	 0.167	  21  neg
0	 137	  40	   35	   168	   43.1	 2.288	  33  pos

I nästa steg skapar vi en regressionsmodell. Detta görs med funktionen lm() (linear model), som används när den beroende variabeln är kontinuerlig. I detta fall är den beroende variabeln (utfallet, Y) blodsocker (glucose).

Med funktionen lm() används tecknet ~ för att specificera modellen. Till vänster om ~ placeras den beroende variabeln (glucose) och till höger specificeras alla prediktorer separerade med +. Resultatet sparar vi i objektet my_model. Vi kör sedan objektet för att se resultatet.

BLOCKERAD KOD, BILD ELLER TEXT


Du måste logga in för att komma åt allt material.

Call:
lm(formula = glucose ~ mass + age + insulin + pressure, data = pima)

Coefficients:
(Intercept)         mass          age      insulin     pressure  
   68.89351      0.58340      0.72475      0.08665      0.03373  

Ovan ser vi modellens specifikation under rubriken Call. Under rubriken Coefficients ser vi koefficienterna. Det var således en tämligen sparsmakad presentation av modellen. Låt oss tillämpa funktionen summary() på objektet istället.

BLOCKERAD KOD, BILD ELLER TEXT


Du måste logga in för att komma åt allt material.

Call:
lm(formula = glucose ~ mass + age + insulin + pressure, data = pima)

Residuals:
    Min      1Q  Median      3Q     Max 
-123.65  -17.78   -2.50   16.23   86.57 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 68.893507   5.491146  12.546  < 2e-16 ***
mass         0.583397   0.138432   4.214 2.81e-05 ***
age          0.724746   0.090398   8.017 4.06e-15 ***
insulin      0.086647   0.009137   9.483  < 2e-16 ***
pressure     0.033732   0.057169   0.590    0.555    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 28.52 on 763 degrees of freedom
Multiple R-squared:  0.2087,	Adjusted R-squared:  0.2045 
F-statistic:  50.3 on 4 and 763 DF,  p-value: < 2.2e-16

Detta extraherar mer information om modellen. För varje prediktor presenteras nu koefficienten (Estimate), standard error (Std. Error), t statistic (t value) och p-värde (Pr(>|t|)). Resultaten tolkas som följer:

• (Intercept): Vi är inte intresserade av interceptet. Det är ytterst sällan som interceptet är av intresse. Interceptet indikerar vad Y-värdet (glucose) är om alla X-värden är 0 och detta är inte biologiskt meningsfullt (t ex är det ingen människa som väger 0 kilo). Interceptet redovisas ändå alltid.
• Prediktorerna:
  • mass: Om mass ökar 1 enhet, så ökar glucose 0.583397 enheter.
  • age: Om age ökar 1 enhet, så ökar glucose 0.724746 enheter.
  • Likaledes är övriga prediktorer också positivt associerade med glucose, eftersom alla koefficienterna är positiva.
  • Alla prediktorer, förutom pressure, har signifikanta p-värden (p<0.05). Det betyder att sambandet mellan pressure och glucose inte är statistiskt signifikant.
  • Association för varje prediktor är justerad för alla andra prediktorer. Det betyder att varje koefficient visar hur prediktorn påverkar glucose, om övriga prediktorer hålls konstanta.
• R² och Adjusted R²: R² är modellens förklaringskapacitet. I detta fall är R² 0.2087, vilket innbär att modellen förklarar 20.87% av all variation i glucose.

Det finns ytterligare funktioner och paket som kan extrahera och presentera modellens resultat. Vi börjar med paketet sjPlot, som har flera användbara funktioner för ändamålet. Med tab_model() presenteras modellens resultat i en tabell:

BLOCKERAD KOD, BILD ELLER TEXT


Du måste logga in för att komma åt allt material.

I resultaten ovan presenteras även ett konfidensintervall som indikerar inom vilket intervall den sanna koefficienten finns med 95% sannolikhet. Med 95% sannolikhet är koefficienten för mass mellan 0.31 och 0.86. Notera att detta intervall inte inkluderar 0 (noll), vilket innebär att vi med 95% säkerhet kan säga att en ökning av mass ger en ökning av glukos. Detta stämmer med P-värdet som är <0.001. Således är sambandet mellan mass och glukos statistiskt signifikant. Denna effekten av mass är justerad för effekten av de övriga prediktorerna (kovariaterna) i modellen.

Koefficienten för mass indikerar hur stor effekt 1 enhet ökning av mass har på glucose, när alla andra prediktorer (kovariater) hålls konstanta. Samma princip gäller för övriga prediktorer.

Notera raden Observations i tabellen. Den redovisar hur många observationer som inkluderades i regressionsmodellen. Kom ihåg att om en observationer saknar uppgift (dvs har "missing") på någon av prediktorerna som ingår i modellen, så exkluderas den observationer från hela modellen. Det innebär att om en observationer exempelvis saknar uppgift på mass, men har värden för alla övriga variabler, så exkluderas den personen från modellen och vi tappar en observation.

Med funktionen plot_model() får vi en grafisk presentation av resultaten. Funktionen har en lång rad argument som du kan använda för att justera grafen.

BLOCKERAD KOD, BILD ELLER TEXT


Du måste logga in för att komma åt allt material.

Linjäritet (Linearity)

Användning av linjär regression kräver att varje prediktor är linjärt associerad med Y. Figuren nedan visar en linjär regressionslinje genom (A) data med icke-linjärt samband, (B) data utan samband och (C) data med linjärt samband. I fall (A) misslyckas den linjära modellen med att fånga upp ett existerande samband. I fall (B) finns inget nämnvärt samband att fånga upp.

Att logaritmera X hjälper inte för att fånga sambandet i fall (A), även om logaritmering kan göra icke-linjära samband mer linjära. För att fånga sambandet i fall (A) krävs en mer flexibel linje. Metoder för att hantera icke-linjaritet har diskuterats i kapitlet Hantering av Kontinuerliga Variabler.

Modelspecifikation: att specificera modellen korrekt

Om en regressionsmodell är felaktigt specificerad så är koefficienterna inte tillförlitliga. Att utelämna viktiga kovariater, interaktioner mellan variabler eller att ignorera icke-linjära samband är vanliga orsaker till att modellen är felaktigt specificerad. Detta kan leda till att effekten av en eller flera prediktorer överskattas, underskattas eller feltolkas. En interaktion innebär att två (eller flera) prediktorer interagerar med varandra, så till vida att effekten av ena prediktorn beror på den andra prediktor.

Interaktioner

Ponera att läkemedel X används för att behandla huvudvärk. Läkemedlet fungerar enbart på kvinnor eftersom substansen aktiveras av ett enzym som finns i bröstkörtlar. Effekten av läkemedlet beror således på kön. Det finns alltså en interaktion mellan läkemedlets effekt och kön. Att inte beakta denna interaktionen leder till en felaktig modell. Den felaktiga modellen specificeras som följer.

\(Y = \beta_0 + \beta_{Läkemedel} X_{Läkemedel} + \beta_{Kön} X_{Kön} + \epsilon_i\)

Modellen är felaktig eftersom den behandlar kön och läkemedel som två separata parametrar (prediktorer). Den korrekta modellen specificeras med interaktionen mellan läkemedel och kön. Interaktionen innebär att man multiplicerar läkemedel med kön (rent matematiskt innebär det att kombinationen av varje nivå av kön och varje nivå av läkemedel introduceras som en prediktor):

\(Y = \beta_0 + \beta_{Läkemedel} X_{Läkemedel} • \beta_{Kön} X_{Kön} + \epsilon_i\)

För att avgöra om interaktionen är signifikant så bedömer man P-värdet för interaktionstermen. All statistisk mjukvara (inklusive R) redovisar P-värde (och 95% konfidensintervall) för interaktionstermer. Om interaktionstermen är signifikant (p < 0.05) så är det bevis för att det finns en interaktion mellan prediktorerna.

Multicollinearity (Multikolinjaritet)

Syftet med multipel regression är att kvantifiera den oberoende effekten av varje enskild prediktor. Det finns dock situationer då en prediktor inte kan göras sig oberoende av en eller flera andra prediktorer. Ponera att vi undersöker hur BMI (body mass indeX) påverkar blodsocker. Vi skapar en regression där blodsocker är utfallsmått och BMI samt ålder är prediktorer:

\(Y_{blodsocker} = \beta_0 + \beta_{BMI}X_{BMI} + \beta_{ålder}X_{ålder} + \epsilon_i\)

Låt oss säga att effekten (koefficienten) för BMI är 0.5, vilket innebär att för varje enhet BMI ökar så stiger blodsockret med 0.5 enheter. Denna effekten av BMI är oberoende av ålder, eftersom vi justerat för ålder. Man kan också säga att detta är effekten av BMI när vi håller ålder konstant. Om vi skulle introducera ytterligare en kovariat, nämligen vikt, så kan multicollinearity uppstå. BMI och vikt är nämligen mycket starkt relaterade (BMI är en funktion av vikt och längd), så till vida att när BMI ökar så tenderar även kroppsvikten öka. Man kan då inte hävda att koefficiencenten för BMI representerar effekten av BMI när vikt hålls konstant, eftersom vikt med största sannolikhet också ökar när BMI ökar. Detta är ett exempel på multikolinjaritet, dvs när två eller flera prediktorer är starkt korrelerade med varandra.

Denna typ av korrelation mellan prediktorer försämrar koefficienternas precision. Det innebär också att vi använder mer information än nödvändigt. Om informationen som finns i BMI även finns i vikt, och vice versa, så räcker det med att en av dessa prediktorer inkluderas i modellen.

Risken för multicollinearity är hög om korrelationen (se kapitel Korrelation) mellan två prediktorer är 0.9 eller högre. Multicollinearity kan dock bli ett problem även vid lägre korrelationsgrader. Med variance inflation factor (VIF) kan man undersöka multicollinearity.

Test för multikollinearitet (multicolinearity)

Det finns flera olika mått för att testa multikollinearitet. Detta kan testas globalt för hela modellen och för enskilda prediktorer i modellen. Ett vanligt mått på multikollinearitet är VIF (variance inflation factor). VIF är ett mått på hur uppblåst en prediktors varians är till följd av multikollinearitet i modellen. Om VIF är större än 5 finns multikollinearitet. Om VIF är >10 bör variabeln exkluderas från modellen.

Vi skapar nu en modell med data från Pimaindianerna. Denna gången använder vi alla variabler i data för att prediera blodsocker (glucose). Vi aktiverar paketet mlbench för att ladda data med Pimaindianer.

BLOCKERAD KOD, BILD ELLER TEXT


Du måste logga in för att komma åt allt material.

Därefter installerar vi paketet performance och använder funktionen check_model() för att extrahera VIF grafiskt och i tabellerat format:

BLOCKERAD KOD, BILD ELLER TEXT


Du måste logga in för att komma åt allt material.

# Check for Multicollinearity

Low Correlation

     Term  VIF   VIF 95% CI Increased SE Tolerance Tolerance 95% CI
 pregnant 1.46 [1.34, 1.61]         1.21      0.69     [0.62, 0.75]
 pressure 1.19 [1.11, 1.31]         1.09      0.84     [0.76, 0.90]
  triceps 1.48 [1.36, 1.63]         1.21      0.68     [0.61, 0.74]
  insulin 1.27 [1.18, 1.40]         1.13      0.79     [0.71, 0.85]
     mass 1.36 [1.26, 1.50]         1.17      0.74     [0.67, 0.80]
 pedigree 1.08 [1.03, 1.22]         1.04      0.92     [0.82, 0.97]
      age 1.55 [1.42, 1.71]         1.24      0.65     [0.58, 0.71]
 diabetes 1.22 [1.14, 1.35]         1.11      0.82     [0.74, 0.88]

Som framgår ovan fanns ingen prediktor med VIF >5, så kolinjaritet är inget problem i vår modell.

Kategoriska prediktorer (factors)

Hittills har diskussionen enbart avhandlat prediktorer som är kontinuerliga variabler. Självfallet kan man studera prediktorer som är kategoriska variabler, såsom kön, behandling, stad, land, etc. Regressionsmodeller brukar innehålla både kontinuerliga och kategoriska prediktorer.

Låt oss säga att vi studerar hur ålder, vikt och födelseland påverkar blodtryck. Modellen ser ut som följer:

\(Y_{blodtryck} = \beta_0 + \beta_{ålder}X_{ålder} + \beta_{vikt}X_{vikt} + \beta_{land}X_{land} + \epsilon_i\)

Ålder och vikt är kontinuerliga variabler. Födelseland (land) är däremot en kategorisk variabel av vilken det finns följande kategorier: (1) Sverige, (2) Norden och (3) Europa. När vi inkluderar land i modellen så måste en av dessa kategorier användas som referenskategori och de övriga kategorierna jämförs med den referensen. Låt oss nu tolka resultaten från regressionsmodellen ovan:

Vi ser att det finns koefficienter (Estimate) för Norden och Europa, men på Sverige finns bara texten "ref". Dett förklaras av att kategorin Sverige är referenskategori, med vilken Norden och Europa jämförs. Så effekterna som presenteras för Norden och Europa är i jämförelse med Sverige. Då blir tolkningen som följer:

• För varje enhet BMI ökar så ökar blodtrycket med 0.91 enheter och detta är statistiskt signifikant.
• För varje år man blir äldre så ökar blodtrycket med 1.14 enheter och detta är statistiskt signifikant.
• Jämfört med Svenskar så har personer födda i Norden 3.0 enheter högre blodtryck. Detta är statistiskt signifikant (p<0.001) och det sanna värdet är mellan 2.0 till 4.1 enheter lägre blodtryck.
• Jämfört med Svenskar så har personer födda i Europa 2.1 enheter högre blodtryck. Detta är statistiskt signifikant (p<0.001) och det sanna värdet är mellan 1.5 till 2.7 enheter lägre blodtryck.