Propensity score: teori, metoder och exempel i R (del 2)

Propensity score: teori, metoder och exempel i R

Beräkning av propensity score

Den verkliga sannolikheten för en person att få behandling kan man aldrig veta. Därför är propensity score alltid ett beräknat värde som anger den uppskattade sannolikheten att ha fått behandling. Som nämnts i ovanstående stycke, kan sannolikhet att ha fått behandling beräknas utifrån observerade variabler (patientkaraktäristika). I praktiken beräknar man oftast propensity score med hjälp av en logistisk regression. I den logistiska regressionen använder man behandlingsstatus (exempelvis 1 = behandlad och 0 = obehandlad) som utfallsmått, och prediktorerna (exponeringsvariablerna) utgörs av patientkaraktäristika (till exempel kön, ålder, genetik, matvanor, etnicitet, sociala faktorer etc.).

Den logistiska regressionen resulterar i en regressionskoefficient som kan användas för att beräkna sannolikheten att ha fått behandling – Med andra ord: man beräknar regressionskoefficienten för varje variabel med hjälp av en logistisk regression vilket sedan används för att beräkna propensity score (d.v.s. sannolikheten för behandling). Det finns även andra metoder för beräkning av propensity score, men logistisk regression torde vara en av de vanligaste.

Detta är Del 2 i vår serie om propensity score. Glöm inte läsa Del 1 – Propensity Score

Formeln för en logistisk regression kan skrivas på följande sätt:

ln[odds(Y = 1)] = β0 + β1X1 + β2X2 + βnXn

 

Till vänster om likhetstecknet får vi sannolikheten för utfallet. I detta fall blir det den beräknade sannolikheten för behandling.

 

Till höger om likhetstecknet anges β0 som representerar interceptet, samt β1-n som representerar regressionskoefficienterna. X1-n betecknar variablerna som representerar olika patientkaraktäristika (t.ex. kön, ålder, inkomst, etc.)

 

Den logistiska regressionen gör att man kan använda de beräknade regressionskoefficienterna (β) för varje variabel. Sannolikheten för behandling beror då på regressionskoefficienterna (β) samt variablerna som betecknar patientkaraktäristika (X1-n).

 

Eftersom vi vet värdet för varje variabel (X) för varje person, och även regressionskoefficienterna (β) kan vi beräkna sannolikheten för varje enskild person att få behandling.

 

 

Syftet med propensity score är att korrigera för ojämlikheter i patientkaraktäristika mellan behandlingsgrupper. Skillnader gällande patientkaraktäristika mellan behandlingsgrupper hindrar en från att dra korrekta slutsatser om behandlingseffekt. Om man däremot har möjlighet att bortse från eventuella skillnader i patientkaraktäristika kan behandlingseffekten beräknas på ett mer tillförlitligt sätt.

 

För att repetera föregående avsnitt:

Om en behandlad och en obehandlad person har samma propensity score, bör det i teorin enbart vara slumpen som skiljer de åt. Det är på så vis som man med hjälp av en propensity score kan få en observationsstudie att efterlikna en randomiserad studie.

 

Variabelselektion till propensity score

För att beräkna propensity score genomförs ofta en logistisk regression. Som nämnts tidigare utgörs utfallet av behandlingsstatus (d.v.s. behandlad eller obehandlad) och prediktorerna (exponeringsvariablerna) utgörs av patientkaraktäristika vid baseline. Detta innebär att man använder sig av så kallade baseline-variabler för att beräkna sannolikheten för att ha fått behandling.

Många menar att samtliga uppmätta baseline-variabler bör inkluderas i den logistiska regressionen för beräkning av propensity score. Det skall dock tilläggas att det inte finns någon säker konsensus gällande frågan om exakt vilka variabler som skall inkluderas. De flesta är nog överens om att variabler (patientkaraktäristika) som påverkar utfallet bör inkluderas. I praktiken kan det dock vara svårt att säkert avgöra exakt vilka variabler som möjligen kan påverka utfallet. Av det skälet menar många att man bör inkludera samtliga uppmätta baseline-variabler. Variabelselektionen vid propensity score är således en betydligt mindre rigorös process jämfört med traditionella metoder för justering där det oftast bör ske med omsorg.

 

Många menar att samtliga uppmätta baseline-variabler bör inkluderas i den logistiska regressionen för beräkning av propensity score.

 

Praktisk tillämpning av propensity score

Med hjälp av propensity score kan man ta hänsyn till skillnader i patientkaraktäristika mellan behandlade och obehandlade och på så sätt minska bias. Följande är exempel på vanliga metoder för att minska bias med hjälp av propensity score:1) Matchning utifrån propensity score, 2) stratifiering utifrån propensity score, 3) användning av propensity score som variabel i en multipel regressionsmodell.

 

Matchning utifrån propensity score

Matchning utifrån propensity score innebär att man matchar par (behandlade och obehandlade personer) som har liknande propensity score. Med andra ord, man skapar par med så likvärdiga propensity score som möjligt. Målsättningen är således att en behandlad person matchas till en obehandlad person med samma, eller ungefär samma, propensity score. Andra metoder, än just matchning, kan tillämpas med propensity score (var god se nedan) för att minska bias. En del menar dock att just matchning utifrån propensity score är det mest tillförlitliga sättet.

Två vanliga metoder kan tillämpas för att matcha behandlade och obehandlade med liknande propensity score. Dessa två metoder är nearest neighbor matching och nearest neighbor matching med restriktion, så kallad specifik caliper distance. Nearest neighbor matchning innebär att en behandlad person med ett visst propensity score matchas till den obehandlade personen vars propensity score är mest lik. Om det finns flera obehandlade personer med samma propensity score väljs en utav dessa slumpmässigt. Matchning med nearest neighbour med restriktion med en viss caliper distance är väldigt likt detta, men man lägger till ytterligare ett krav (det är detta krav som restriktionen syftar till). Kravet (restriktionen) som man lägger till är att skillnaden i propensity score mellan behandlade och obehandlade par får inte vara större än ett visst förutbestämt värde. Restriktionen är alltså ett förutbestämt värde som kallas för caliper distance. Detta innebär att man, från en grupp med obehandlade personer som har en propensity score inom ett visst intervall (caliper distance), väljer ut den som har det mest lika värdet som den behandlade personen som den skall matchas till. Om det inte finns obehandlade personer som har propensity score inom det förutbestämda intervallet för en given behandlad person innebär det att den behandlade personen inte kommer att matchas och exkluderas därmed från analyserna. Ett högt caliper-värde innebär att fler matchningar kommer att kunna genomföras, men innebär samtidigt att olikheterna mellan behandlade och obehandlade blir större. Omvänt gäller att ett lågt caliper-värde ger medför att behandlade och obehandlade är mer lika varandra, men samtidigt att färre matchningar kan genomföras. Gällande det exakta värdet för restriktionen, det vill säga caliper-värdet, finns inga bestämda regler. Ofta rekommenderas dock att man använder ett caliper-värde som är 0.2 standardavvikelser av den naturliga logaritmen av det estimerade propensity score.

Ofta används en metod som kallas 1:1 matchning vilket innebär att en behandlad och en obehandlad person utgör ett matchat par. Man kan dock även matcha en obehandlad person till flera behandlade (på engelska kallas detta för matching  with replacement). En annan metod för matchning utifrån propensity score är det som kallas greedy matching. Detta innebär att man slumpmässigt väljer en behandlad person och sedan matchar till den obehandlade personen vars propensity score är mest lik den slumpmässigt utvalda, behandlade, personen. Därefter upprepas detta fram tills det att alla behandlade personer har matchats, eller tills det att man inte längre finner någon obehandlad person att para ihop med. Matchning kan också göras med hjälp av en metod som kallas för optimal matching. Optimal matching innebär att varje par, bestående av en behandlad och obehandlad person, väljs ut (matchas) på ett sätt som ger minsta möjliga skillnad i propensity score.

Matchning utifrån propensity score medför att matchade par är väldigt lika varandra med avseende på patientkaraktäristika. Med hjälp av denna metod kan man beräkna behandlingseffekten på ett mer tillförlitligt sätt än om det inte hade gjorts. Det är alltså ett sätt att kontrollera för bias, d.v.s. ojämlikheter i patientkaraktäristika. Effekten av en behandling kan då beräknas genom att direkt jämföra behandlade med obehandlade. På så vis är en propensity score matchade studie tämligen lik en randomiserad studie. Gällande kontinuerliga variabler kan man jämföra medelvärde och median och för kategoriska variabler andel/proportioner. Gällande signifikanstester rekommenderar man ofta att använda parat t-test (för kontinuerliga variabler) och McNemars test (för kategoriska variabler). Skälet till detta är att många anser att matchade par inte utgör oberoende observationer (både t-test och chi-två-test har som villkor att observationerna är oberoende).

Om antalet personer i studien är liten kan det trots en propensity score matchning kvarstå skillnader mellan behandlade och obehandlade. I sådana fall är det möjligt att använda sig av regressionsanalyser på det matchade datasetet för att korrigera ytterligare. Man kan alltså kombinera propensity score matchning med regressionsanalys och kovariatjustering för ökad precision.

 

Stratifiering utifrån propensity score

Vid stratifiering delas behandlade och obehandlade personer upp i grupper (strata), där alla personer inom ett visst stratum har samma eller liknande propensity score. En vanlig metod är att dela upp behandlade och obehandlade i fem grupper utifrån propensity score kvintiler (femtedelar). Detta medför att inom varje grupp (stratum) kan man direkt jämföra behandlingseffekten eftersom behandlade och obehandlade har liknande propensity score och därmed torde patientkaraktäristika vara jämt fördelade.

Efter att ha stratifierat utifrån propensity score kan man beräkna en gemensam riskdifferens för behandlade och obehandlade i samtliga strata. Man bör då vikta varje enskild strata utifrån antalet personer som de innehåller.

 

Kovariatjustering med propensity score

En tredje metod som kan användas för att minska bias är kovariatjustering med hjälp av propensity score. Detta innebär att man beräknar propensity score för varje behandlad och obehandlad person och sedan inkluderar värdet i en multipel regression. Följaktligen inkluderar man i regressionsmodellen en behandlingsvariabel (exponering), en utfallsvariabel och propensity score. Resultatet blir att man uppskattar effekten av en viss exponering på utfallet, justerat för propensity score. Detta innebär alltså att man, med hjälp av propensity score, tar hänsyn till skillnader i patientkaraktäristika mellan behandlade och obehandlade. Med andra ord, behandlingseffekten uppskattas oberoende av uppmätta skillnader i mellan behandlade och obehandlade personer. Efter att ha beräknat propensity score kan olika typer av regressionsmodeller tillämpas. Den exakta typen av regressionsmodell får givetvis avgöras utifrån utfallet som man önskar att studera.

 

Kontrollera propensity score har beräknats korrekt

Som nämnts tidigare kan vi aldrig med säkerhet veta den sanna propensity score – det vill säga den faktiska sannolikheten för en person att få behandling eller ej. Vid beräkning av propensity sker således enbart en uppskattning av det sanna värdet utifrån uppmätta variabler (patientkaraktäristika). Av det skälet är det viktigt att bedöma att beräkningen (specificeringen) av propensity score har gjorts korrekt.

Skillnader i patientkaraktäristika mellan behandlade och obehandlade, som har samma eller närliggande värden utav propensity score, kan tala för beräkningen av själva propensity score inte har utförts korrekt. Behandlade personer med ett visst propensity score bör vara tämligen lika obehandlade personer med samma propensity score. Om så inte är fallet, kan det indikera problem med uträkningen av propensity score.

Gällande sample som har matchats utifrån propensity score bör man kontrollera att medelvärden, medianvärden och proportioner är lika mellan behandlade och obehandlade par som har matchats. För att efterforska detta vidare kan man, för både kontinuerliga och kategoriska variabler, beräkna standardiserade skillnader. Beräkning av standardiserade skillnader innebär att skillnaden i värdet för en viss variabel, mellan en behandlad och obehandlad person, uttrycks i standardavvikelser. En standardiserad skillnad för en viss variabel som är lika med eller mindre än 0.1 anses indikera att variabeln är jämt fördelad. Detta innebär att om en viss variabel jämförs mellan en behandlad och en obehandlad person och den standardiserade skillnaden är mindre än 0.1 indikerar detta att man med hjälp av propensity score har på ett adekvat sätt lyckats balansera fördelningen jämt gällande just den variabeln. I tillägg till detta bör man även granska varians, skewness och kurtosis. Fördelningen bör också granskas grafiskt till exempel med hjälp av box plots och density plots.

Om väsentliga skillnader kvarstår mellan behandlade och obehandlade, med liknande propensity score, kan man beräkna om en ny propensity score till exempel genom att inkludera fler variabler och/eller eventuella interaktioner.

 

Propensity score i R: Praktisk vägledning

I R kan man genomföra propensity score direkt med hjälp av befintliga basfunktioner, eller med enkelhet, med hjälp av särskilda paket, till exempel nonrandom. I detta exempel går vi igenom båda metoder. Data som vi använder oss av i det här exemplet heter Right Heart Catheterization Dataset, och finns att ladda ner via denna länk (observera dock att man behöver inte ladda ner datan, utan det går även att hämta det direkt till R genom att skriva in hemsidans URL i R Studio, var god se nedan). Vi använder samma data för detta exempel som i kapitlet om logistisk regression.

Dataset Right Heart Catheterization innehåller uppgifter om patienter som vårdats på intensivvårdsavdelning och erhållit, en så kallad PA-kateter (På engelska kallas det för Right heart catheter, eller Swan-Ganz catheter. En PA-kateter kan användas, hos en svårt sjuk person, för att på ett exakt sätt mäta hjärtats pumpförmåga samt blodtryck-och blodflöde).

Datasetet Right Heart Catheterization innehåller 63 variabler  (patientkaraktäristika) på 5735 patienter. Den ursprungliga studien som använde dessa data publicerades av Connors et al i tidskriften JAMA 1996. Studien var en observationsstudie som undersökte överlevnad i relation till att ha fått en PA-kateter. Man fann att patienter som erhållit en PA-kateter hade ökad risk för död jämfört med patienter som inte erhållit en PA-kateter. Eftersom att studien var icke-randomiserad (observationsstudie) valde man att analysera data med hjälp av propensity score.

I det här exemplet är vi framförallt intresserade av två varibler: 1) prediktorvariabeln (exponering) och 2) utfallsvariabeln. Prediktorn (exponeringen) är huruvida patienten har fått en PA-kateter, och denna variabel är döpt till swang1 (efter den engelska förkortningen för Swan-Ganz catheter). Utfallsvariabeln är död, det vill säga huruvida patienten har avlidit eller ej, och är i datasetet döpt till death.

### Installera paket, hämta dataset och granska data

### Börja med att installera följande paket

install.packages("pROC")

install.packages("ggplot2")

install.packages("moonBook")

install.packages(”nonrandom”)

### Ladda paketen

library(pROC)

library(ggplot2)

library(moonBook)

library(nonrandom)

### Vi använder data som kan nås via denna länk. För att läsa in datan utan att ladda ner filen kan vi använda följande kod. Vi döper vår data till rhc

rhc <- read.csv("http://biostat.mc.vanderbilt.edu/wiki/pub/Main/DataSets/rhc.csv")

### För att få översikt, se alla variabelnamn och första sex patienter kan vi använda oss av nedanstående kod

View(rhc)

names(rhc)

head(rhc)

### Se antal personer som har avlidit (death) och antal personer som har fått PA-kateter (swang1).

summary(rhc$death)

summary(rhc$swang1)

### För att se antalet personer som har avlidit i förhållande till antalet personer som har fått PA-kateter kan man skriva nedanstående kod. I kolumnerna betecknar ”No” eller ”Yes” antalet personer som avlidit, och i raderna betecknas antalet personer som har fått PA-kateter med RHC respektive No RHC.

tab <- table(rhc$swang1, rhc$death)

tab

### Nedanstående visar förhållandet mellan andelen personer som avlidit bland de som har fått (eller inte fått) PA-kateter samt p-värde beräknat med Chi-två-test.

prop.table(tab, 1)

chisq.test(tab)

### Ett enklare och mer effektivt sätt att få fram samma resultat är att använda paktetet ”moonBook” med följande kod

mytable(swang1 ~ death, data=rhc)

### Vill man dessutom granska ytterligare variabler, till exempel ålder (age), kön (sex), etnicitet (race) och inkomst (income) lägger man till följande.

mytable(swang1 ~ death + age + sex + race + income, data=rhc)

# Gör en ojusterad logistisk regression på effekten av PA-kateter på risk för död.

ojusterad <- glm(formula = death ~ swang1,

family  = binomial(link = "logit"),

data = rhc)

# Använd nedanstående kod för att beräkna oddskvot (odds ratio) för den ojusterade effekten av PA-kateter på risken för död.

exp(coef(ojusterad))

## Första steget gällande beräkning av propensity score är att skapa en logistisk regressionsmodell som inkluderar alla variabler. Vi döper modellen till propmodell. Genom denna kod skapar vi en logistisk regression som, utifrån alla andra variabler, beräknar sannolikheten att ha fått en PA-kateter (swang1).

propmodell <- glm(formula = swang1 ~ age + sex + race + edu + income + ninsclas + cat1 + das2d3pc + dnr1 + ca + surv2md1 + aps1 + scoma1 + wtkilo1 + temp1 + meanbp1 + resp1 + hrt1 + pafi1 + paco21 + ph1 + wblc1 + hema1 + sod1 + pot1 + crea1 + bili1 + alb1 + resp + card + neuro + gastr + renal + meta + hema + seps + trauma + ortho + cardiohx + chfhx + dementhx + psychhx + chrpulhx + renalhx + liverhx + gibledhx + malighx + immunhx + transhx + amihx,

family  = binomial(link = "logit"),

data = rhc)

## Sedan använder vi funktionen predict för att beräkna propensity score för varje enskild individ utifrån ovanstående modell. Vi lägger till resultatet av propensity score till datasetet och kallar det för ps.

rhc$ps <- predict(propmodell, type = "response")

### Ett värde för propensity score bör nu finnas för varje enskild patient som kan ses i den sista kolumnen till höger, döpt till ps.

View(rhc)

### Vi kan även se medianvärde, samt minsta och högsta värde för ps för hela datasetet genom att skriva följande

summary(rhc$ps)

### Granska ROC (Receiver operator characteristic) kurva och AUC (area under curve) för att bedöma modellens prediktiva värde

rocPsModel <- roc(swang1 ~ ps, data = rhc)

plot(rocPsModel, legacy.axes = TRUE)

rocPsModel

# Visa distribution och överlappning av propensity score mellan behandlade och icke-behandlade. Den första raden ger dig histogram och den andra raden en density plot.

ggplot(rhc, aes(x=ps, fill=swang1)) +

geom_histogram(binwidth=.03, alpha=.5, position="identity")

ggplot(rhc, aes(x=ps, colour=swang1)) + geom_density()