Multipel regression (regressionsanalys)

Inom medicinsk forskning är det mycket vanligt att man samlar in information om flera variabler för varje individ/cell/enhet. I en klinisk studie kan detta 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 utmärkt verktyg. Regression tillåter oss göra följande:

• Studera relationen mellan en beroende variabel (utfallsmått, outcome) och en eller flera oberoende variabler (prediktorer).
• Studera hur stark relationen är mellan varje enskild prediktor och den beroende variabeln.
• Skapa en regressionmodell som kan användas för prediktion (dvs för att förutsäga värdet på den beroende variabeln). Exempelvis kan det vara av intresse att prediktera risken för cancer utifrån nivån på en cancermarkör i blodet.

Den fortsatta diskussionen kräver att du är bekant med regressionsanalys som diskuterades i föregående föreläsning: Regression och korrelation.

Målsättning med multipel regression

Som nämnt ovan kan regression användas för att (1) studera relationen mellan variabler, (2) studera sambandens styrka och (3) utföra prediktioner. Inom den kliniska forskningen 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. Man bör därför skilja på regressionsmodeller som syftar till prediktion samt modeller som syftar till effekt-estimering. Ett förtydligande följer.

Prediktion – Ponera att vi har samlat in data för 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 boende. Om vårt syfte är att kunna förutsäga vilka personer som kommer utveckla lungcancer, så är syftet med våra regressionsmodeller att göra prediktioner. När man konstruerar sådana prediktionsmodeller så står utmaningen i följande:

1. Välja rätt typ av regressionsmodell (exempel på alternativ: logistisk regression, linjär regression, Cox regression)
2. Välja vilka variabler som skall inkluderas i modellen. Detta är tämligen komplicerat och diskuteras i avsnittet Variabelselektion – Att välja variabler.
3. Utvärdera och förbättra modellens precision.

Effekt-estimering – Ponera att vi har samlat in data för 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 boende. Om vårt syfte är att studera relationen mellan (exempelvis) rökning och risken för lungcancer, så är syftet med modellen att estimera effekt av en prediktor, i detta fall rökning. Då är det önskvärt att sambandet mellan rökning och lungancer 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 lungcancer. Det vore därför önskvärt att ta höjd för – eller justera som det oftare kallas – socioekonomisk status. Detta görs helt enkelt genom att socioekonomisk status inkluderas som en prediktor i modellen. När man bygger regressionsmodeller är man alltså intresserad av hur en (alternativt flera) prediktorer relaterar till utfallsmåttet och som regel behöver man justera för faktorer som kan störa sambandet. När man bygger sådana modeller kan man – utifrån vetenskap och beprövad erfarenhet – själv avgöra vilka prediktorer som behöver finnas i modellen för att den skall vara meningsfull.

Även om dessa exempel baserades på ett utfallsmått som är dikotomt (att insjukna i lungcancer, ja/nej) så kommer den efterföljande diskussionen handla om utfallsmått som är kontinuerliga, dvs mätvärden eller andra numeriska värden. Regression med kontinuerliga utfallsmått kallas linjär regression. I senare avsnitt kommer andra typer av regression diskuteras.

Simpel linjär regression: ett utfallsmått och en prediktor.

Den vanliga regressionsekvationen ser ut som följer:

Y_i = β₀ + β₁X_i + e_i

Y_i är den beroende variabeln (utfallsmåttet, outcome).

β₀ är interceptet (dvs där regressionslinjen skär Y-axeln).

β₁ är koefficienten för prediktorn X_i, vilket innebär att det är värdet för linjens lutning.

e_i är modellens felprecision (residual error). Ingen modell är perfekt; det kommer alltid finnas ett visst fel i prediktionerna och detta felet är e_i.

Låt oss titta på dessa fyra parametrar i Figur 1.

I Figur 1-A redovisas relationen mellan två variabler (X, Y). Som framgår i Figur 1-A så finns ett positivt samband mellan X och Y, så till vida att när X ökar så ökar Y. Således kan vi omedelbart med en visuell bedömning konstatera att det finns ett samband mellan X och Y och sannolikt kan X användas för att prediktera Y.  För att på bästa sätt beskriva sambandet mellan X och Y så konstrueras en matematisk formel (samma formel som ovan):

Y = β₀ + β₁X₁ + e

Denna formel har också förklarats i föregående kapitel (Regression) samt i Figur 1. Formeln beskriver sambandet mellan den beroende variabeln Y och prediktorn X. Syftet är att studera relationen mellan X och Y, samt använda X för att förutsäga (prediktera) värdet på Y. Detta görs genom att konstruera en matematisk formel som definierar en rak linje genom datapunkterna. Denna linjen bör vara så bra som möjligt på att prediktera så många Y-värden som möjligt. Linjen skär Y-axeln (β₀), har en viss lutning för prediktorn X (β₁) och dessutom har den en viss felprecision (e) eftersom den inte kan pricka datapunkterna perfekt. Utmaningen ligger i att hitta värden för β₀, β₁ och e. Detta görs oftast genom metoden least-squares som är beskrivet i Figur 1-B. Kort sagt så innebär least-squares att linjen tippas fram och tillbaka till dess att den sammanlagda summan av alla error (se Figur 1-B) minimeras. När summan av alla error är minimal så erhåller man den linjen som bäst beskriver data.

Multipel regression

Detta exempel beskriver alltså hur en prediktor (X) relateras till den beroende variabeln (Y). Det är vanligt att vi önskar inkludera ytterligare prediktorer i regressionsmodellen och då talar vi om multipel regression. Med andra ord är multipel regression helt enkelt en modell med minst 2 prediktorer. Den matematiska formeln för multipel regression med två prediktorer ser ut som följer:

Y = β₀ + β₁X₁ + β₂X₂ + e

Detta samband kan inte längre presenteras med en rak regressionslinje. Sambandet mellan Y och prediktorerna X1 och X2 beskrivs nu istället med ett tredimensionellt regressionsplan. Figur 2 visar ett sådant regressionsplan.

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. Det är således önskvärt att erhålla ett objektivt mått på huruvida sambandet mellan X och Y är signifikant. Detta faller faktiskt inom ramen för hypotestestning och nollhypotesen. I dessa 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 Figur 3.

Signifikanstestning med t-test (t-statistic)

För att klargöra om det finns ett statistiskt signifikant samband mellan X och Y måste vi testa nollhypotesen. Som alltid är nollhypotesen att det inte finns något samband. Vi måste hitta en probabilitet (sannolikhet) 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 X är en statistiskt signifikant prediktor (alternativt att sambandet mellan X och Y är statistiskt signifikant).

Nollhypotesen testas med ett t-test. Värdet som beräknas kallas t-statistic, och det beräknas genom att koefficienten (β) 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 installerar och aktiverar sjPlot som underlättar utskrift av tabeller
install.packages("sjPlot")
library(sjPlot)
# Vi installerar och aktiverar mlbench som innehåller data från Pimaindianer (dessa indianer har mycket hög risk för diabetes)
install.packages("mlbench")
library(mlbench)
# med funktionen "data" anropar vi datafilen med Pimaindianer
data(PimaIndiansDiabetes)
# Vi inspekterar filen med View och summary för att få en överblick
View(PimaIndiansDiabetes)
summary(PimaIndiansDiabetes)
# Vi är intresserade av hur BMI, ålder, insulin och blodtryck relaterar till glukos och gör därför en scatter plot matrix med funktionen "scatterplotMatrix" från paketet "car"
install.packages("car")
library(car)
scatterplotMatrix(~glucose+mass+age+insulin+pressure, data=PimaIndiansDiabetes)
# Nu gör vi en linjär regression (multipel regression) där utfallsmåttet är glukos och vi har 4 prediktorer. Funktionen för linjär regression är "lm" som står för "linear model". Modellen specificeras genom att den beroende variabeln (glukos) placeras till vänster om "~" och prediktorerna till höger om "~". Därefter anges datafilen.
modell <- lm(glucose ~ mass + age + insulin + pressure, data=PimaIndiansDiabetes)
# Skriv ut en snygg tabell med resultaten
sjt.lm(modell)
# alternativt kan vi använda funktionen "summary"
summary(modell)
# Granska om modellen tillfredsställer antagandena för linjär regression.
Detta görs genom funktionen "plot". Du måste trycka på "Enter" för att bläddra igenom alla graferna som skall bedömas.
plot(modell)

Tolkning av resultatet (koefficienterna)

Tabellen visar resultat från ovanstående regressionmodell.

		Estimate	95% CI	p
(Intercept)		68.89	58.11 – 79.67	<0.001
mass		0.58	0.31 – 0.86	<0.001
age		0.72	0.55 – 0.90	<0.001
insulin		0.09	0.07 – 0.10	<0.001
pressure		0.03	-0.08 – 0.15	0.555
Observations		768
R2 / adj. R2		.209 / .205

Den första kolumnen presenterar interceptet, prediktorerna (mass, age, insulin, pressure), antalet observationer och R². För interceptet och prediktorerna har vi fått uppgifter om “Estimate”, “95% CI” och “P” (P-värde). Notera att om Du använder summary()-funktionen så ser utskriften annorlunda ut men innehållet är detsamma. Härnäst tolkar vi vardera parameter som presenterats för oss.

• Intercept: Vi är inte intresserade av interceptet. Faktum är att det är ytterst sällan som interceptet är av intresse. Interceptet indikerar endast vad Y-värdet (glukos) ä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.
• För mass är estimate 0.58. Detta är koefficienten (“lutningen”) för mass. Det innebär att när mass ökar med 1 enhet, så ökar glukos med 0.58 enheter. Koefficienten redovisar alltså hur mycket 1 enhet ökning av prediktorn (mass) påverkar den beroende variabeln (glukos). Vi får även ett 95% CI för mass och det indikerar var den sanna koefficienten finns med 95% sannolikhet. Så med 95% sannolikhet är koefficienten mellan 0.31 och 0.86. Notera att detta intervall inte inkluderar 0.0, 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 som finns i modellen. Man kan också uttrycka det så här: koefficienten för mass indikerar hur stor effekt 1 enhets ökning av mass har på glukos, när alla andra variabler hålls konstanta. Dessa principer gäller för alla prediktorer!
• Age och insulin tolkas på samma sätt som mass.
• Pressure har en koefficient på 0.03, vilket innebär att när pressure stiger med 1 enhet så stiger glukos med 0.03 enheter. Effekten av pressure tycks således vara mindre än effekten av övriga prediktorer. Av konfidensintervallet framgår också att det inte med säkerhet kan sägas att pressure faktiskt har en effekt. Konfidensintervallet sträcker sig nämligen från -0.08 till 0.15 vilket inkluderar 0.0 och då är inte koefficienten statistiskt signifikant. Detta stämmer med P-värdet som är 0.555 (alltså större än 0.05 som är gränsen för signifikans).
• Observations redovisar endast hur många personer som inkluderades i regressionen. Kom ihåg att om en person saknar uppgift (saknat mätvärde) på någon av prediktorerna som ingår i modellen, så exkluderas den personen från hela modellen och du tappar då en observation. Det innebär att om en person (exempelvis) saknar uppgift om insulin, så exkluderas den personen från modellen och du tappar en observation.
• R² är modellens förklaringskapacitet. I detta fall är R² 0.209, vilket innbär att modellen förklarar 20.9% av all variation i glukos. Detta kan tyckas vara väldigt lite men det är faktiskt helt OK för en regressionsmodell.

Linjär regression innebär alltså att man studerar ett utfallsmått som är kontinuerligt. Linjär regression kan var antingen simpel (då finns endast en prediktor) eller multipel (då finns minst två prediktorer). Genom att använda linjär regression kan vi studera magnituden och kvaliteten på samband mellan prediktorer och utfallsmåttet. Vi är också intresserade av att klargöra om sambandet mellan varje prediktor och utfallsmåttet är statistiskt signifikant och för detta använder vi P-värden och 95% konfidensintervall. Om prediktorn inte är statistiskt signifikant så är den som regel inte av betydelse för utfallsmåttet. Vid multipel regression är varje prediktors effekt justerad för de andra prediktorerna som ingår i modellen. Det innebär att om modellen inkluderar prediktorerna X₁, X₂ och X₃, så är den effekten som noteras för (exempelvis) X₁ justerad för effekten av X₂ och X₃.

Användning av linjär regression kräver att varje prediktor är linjärt associerad med Y. Figuren nedan visar en simpel regression som i ena fallet (A) visar ett linjärt samband mellan X och Y och i andra fallet (B) uppvisar ett icke-linjärt samband. I fall (B) kan linjär regression inte användas eftersom en rak linje inte kan fånga data.

Linjäritet (eng. linearity)

Linjär regression kräver alltså att X och Y skall vara linjärt associerade. Det innebär egentligen att koefficienten för X skall vara linjärt associerad med Y. I figuren ovan ses ett linjärt (A) och ett icke-linjärt (B) samband mellan X och Y. I (A) kan sambandet beskrivas med en rak (linjär) linje. I (B) är sambandet inte rakt utan en mer flexibel linje behövs för att fånga sambandet mellan X och Y. Denna typ av spridningsdiagram (scatterplot) kan rekommenderas för att bedöma linjäritet mellan en prediktor och ett utfallsmåt. Om sambandet mellan X och Y inte är linjärt så kan man transformera X och/eller Y. Ett exempel på transformation är att logaritmera X eller Y och därefter undersöka sambandet igen. Logaritmen kan ofta göra icke-linjära samband mer linjära. Om man inte önskar logaritmera så finns flertal metoder för att hantera icke-linjaritet men dessa diskuteras inte här.

Modelspecifikation: att specificera modellen korrekt

Om en regressionsmodell är felaktigt specificerad så är koefficienterna inte pålitliga. Att utelämna viktiga prediktorer 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 eller underskattas. En annan viktig anledning till att modeller felspecificeras är interaktioner. En interaktion innebär att två (eller flera) prediktorer interagerar med varandra, så till vida att effekten av den ena prediktorn beror på den andra prediktor. Exempel följer:

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

Y = β₀ + β•Läkemedel + β•Kön + e

Modellen är felaktig eftersom den behandlar kön och läkemedel som två separata prediktorer. Den korrekta modellen specificeras med interaktionen mellan läkemedel och kön. Interaktionen innebär oftast 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 = β₀ + (β•Läkemedel)•(β•Kön) + e

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

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 påverkar blodsocker. Vi skapar en regression där blodsocker är utfallsmått och BMI samt ålder är prediktorer. Låt oss säga att effekten (koefficienten) för BMI var 0.2, vilket innebär att för varje enhet BMI ökar så stiger blodsockret med 0.2 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 prediktor, 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 säga att “man håller vikt konstant”, eftersom vikt med största sannolikhet också ökar när BMI ökar.  Det innebär alltså att om två eller flera prediktorer är starkt korrelerade så kan multicollinearity uppstå.

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 endera används i modellen.

Risken för multicollinearity är hög om korrelationen 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, men det diskuteras inte här då det sällan är nödvändigt i praxis.

Kategoriska prediktorer

Hittills har diskussionen enbart avhandlat prediktorer som är kontinuerliga variabler. Självfallet kan man studera prediktorer som är kategoriska variabler, såsom kön, sjukdomstyp, celltyp, kalenderår, antibiotikum osv. Det är vanligt att en regression inkluderar en blandning av kontinuerliga och kategoriska prediktorer.

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

Blodtryck = β₀ + β_Ålder + β_BMI + β_{Blodtrycksmedicin} + e

Ålder och BMI är kontinuerliga variabler. Blodtrycksmedicin är däremot en kategorisk variabel av vilken det finns följande kategorier: (1) ACE-hämmare, (2) diuretika, (3) beta-blockerare och (4) AR-blockerare. När vi introducerar blodtrycksmedicin i modellen så kommer en av dess kategorier att användas som referenskategori och de övriga kategorierna kommer jämföras med den referensen. Vilken av kategorierna som är referensnivå beror på mjukvaran man använder. Om man kodat läkemedelsgrupp med siffror så brukar mjukvaran använda den lägsta siffran som referensnivå. All mjukvara tillåter att man ändrar vilken kategori som är referens nivå (exempel följer senare). Låt oss nu tolka resultaten från en regressionsmodell med en kategorisk variabel, nämligen blodtrycksmedicin. Utfallsmåttet är blodtryck och vi har även med BMI och ålder som prediktorer.

		Estimate	95% CI	P values
(Intercept)		83.68	73.47 – 90.90	<.001
BMI		0.91	0.76 – 1.05	<.001
Ålder		1.14	1.01 – 1.30	<.001
Blodtrycksmedicin
diuretika		-1.60	-2.00 – -1.19	<.001
beta-blockerare		2.11	-2.01 – -2.31	<.001
AR-blockerare		-3.11	-2.71 – -4.51	<.001
Observations		150
R2 / adj. R2		.837 / .832

I detta fall ser vi att det finns koefficienter för tre av blodtrycksmedicinerna, nämligen diuretika, beta-blockerare och AR-blockerare. Det finns ingen koefficient för ACE-hämmare och detta beror på att ACE-hämmare är vår referenskategori! Så de övriga tre (diuretika, beta-blocekrare, AR-hämmare) jämförs med ACE-hämmare. 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 ACE-hämmare, så har personer med diuretika 1.60 enheter lägre blodtryck. Direutika är alltså assoceirat med lägre blodtyck, jämfört med ACE-hämmare och detta är statistiskt signifikant.
• Jämfört med ACE-hämmare så har personer med betablockerare 2.11 enheter högre blodtryck. Beta-blockerare är alltså associerat med högre blodtryck, jämfört med ACE-hämmare och det är statistiskt signifikant.
• Jämfört med ACE-hämmare så har personer med AR-blockerare 3.11 enheter lägre blodtryck och resten av tolkningen följer samma metod som ovan.

Här följer ett exempel på regression med en kategorisk variabel i R. Denna gången använder vi ett dataset med blommor för enkelhetens skull. Om du vill se hela outputen så klickar du här.

# Vi installerar och aktiverar "datasets" som inkluderar "iris"
install.packages("datasets")
library(datasets)
# Laddar datasettet "iris"
data(iris)
# Vi granskar tabellen
View(iris)
# Vi undersöker vilka variabler som finns i iris. Detta görs med str()-funktionen. Denna funktionen visar att det finns 5 variabler och variabeln "Species" är en "Factor" med 3 kategorier. Det är denna variabeln som vi är intresserade av.
str(iris)
# Vi kan dubbelkolla att "Species" faktiskt är en kategorisk variabel (faktorvariabel) genom att använda class()-funktionen
class(iris$Species)
# med funktionen levels() får vi veta alla kategorier som Species har och den kategorin som skrivs ut först är referenskategorin!
levels(iris$Species)
# Då ser vi att "setosa" är referenskategorin för Species.
# Nu gör vi en linjär regression där vi studerar hur tre prediktorer påverkar Sepal.Length och vi är primärt intresserade av Species.
modell <- lm(Sepal.Length ~ Petal.Length + Petal.Width + Species, data=iris)
# Skriv ut en snygg tabell med resultaten. Resultaten redovisas nedan.
sjt.lm(modell)