Logistisk regression
Eduq
2017-07-24
Logistisk regression
Vi kommer göra en logistisk regression med data från Vanderbilt University. Från Vanderbilt laddar vi ner ett dataset med patienter som har hjärtsjukdom. Själva data och de ingående variablerna är inte av intresse, utan endast (1) metoden för att skapa en logistisk modell och (2) tolkning av resultaten.
# Ladda in (Direkt från Internet) dataset med hjärtpatienter
rhc <- read.csv("http://biostat.mc.vanderbilt.edu/wiki/pub/Main/DataSets/rhc.csv")
# Se datastruktur och variabler
str(rhc)## 'data.frame': 5735 obs. of 63 variables:
## $ X : int 1 2 3 4 5 6 7 8 9 10 ...
## $ cat1 : Factor w/ 9 levels "ARF","CHF","Cirrhosis",..: 6 9 8 1 9 6 8 1 8 1 ...
## $ cat2 : Factor w/ 6 levels "Cirrhosis","Colon Cancer",..: NA NA 6 NA NA NA NA 3 NA NA ...
## $ ca : Factor w/ 3 levels "Metastatic","No",..: 3 2 3 2 2 2 1 2 3 3 ...
## $ sadmdte : int 11142 11799 12083 11146 12035 12389 12381 11453 12426 11381 ...
## $ dschdte : int 11151 11844 12143 11183 12037 12396 12423 11487 12437 11400 ...
## $ dthdte : int NA 11844 NA 11183 12037 NA NA 11491 NA NA ...
## $ lstctdte: int 11382 11844 12400 11182 12036 12590 12616 11490 12560 11590 ...
## $ death : Factor w/ 2 levels "No","Yes": 1 2 1 2 2 1 1 2 1 1 ...
## $ cardiohx: int 0 1 0 0 0 0 0 0 0 0 ...
## $ chfhx : int 0 1 0 0 0 1 0 0 0 0 ...
## $ dementhx: int 0 0 0 0 0 0 0 0 0 0 ...
## $ psychhx : int 0 0 0 0 0 0 0 0 0 0 ...
## $ chrpulhx: int 1 0 0 0 0 1 0 0 0 0 ...
## $ renalhx : int 0 0 0 0 0 0 0 0 0 0 ...
## $ liverhx : int 0 0 0 0 0 0 0 0 0 0 ...
## $ gibledhx: int 0 0 0 0 0 0 0 0 0 0 ...
## $ malighx : int 1 0 1 0 0 0 1 0 0 1 ...
## $ immunhx : int 0 1 1 1 0 0 0 0 0 0 ...
## $ transhx : int 0 1 0 0 0 0 0 1 0 0 ...
## $ amihx : int 0 0 0 0 0 0 0 0 0 0 ...
## $ age : num 70.3 78.2 46.1 75.3 67.9 ...
## $ sex : Factor w/ 2 levels "Female","Male": 2 1 1 1 2 1 2 2 1 1 ...
## $ edu : num 12 12 14.07 9 9.95 ...
## $ surv2md1: num 0.641 0.755 0.317 0.441 0.437 ...
## $ das2d3pc: num 23.5 14.8 18.1 22.9 21.1 ...
## $ t3d30 : int 30 30 30 30 2 30 30 30 30 30 ...
## $ dth30 : Factor w/ 2 levels "No","Yes": 1 1 1 1 2 1 1 1 1 1 ...
## $ aps1 : int 46 50 82 48 72 38 29 25 47 48 ...
## $ scoma1 : int 0 0 0 0 41 0 26 100 0 0 ...
## $ meanbp1 : num 41 63 57 55 65 115 67 128 53 73 ...
## $ wblc1 : num 22.1 28.9 0.05 23.3 29.7 ...
## $ hrt1 : int 124 137 130 58 125 134 135 102 118 141 ...
## $ resp1 : num 10 38 40 26 27 36 10 34 30 40 ...
## $ temp1 : num 38.7 38.9 36.4 35.8 34.8 ...
## $ pafi1 : num 68 218 276 157 478 ...
## $ alb1 : num 3.5 2.6 3.5 3.5 3.5 ...
## $ hema1 : num 58 32.5 21.1 26.3 24 ...
## $ bili1 : num 1.01 0.7 1.01 0.4 1.01 ...
## $ crea1 : num 1.2 0.6 2.6 1.7 3.6 ...
## $ sod1 : int 145 137 146 117 126 138 136 136 136 146 ...
## $ pot1 : num 4 3.3 2.9 5.8 5.8 ...
## $ paco21 : num 40 34 16 30 17 68 45 26 40 30 ...
## $ ph1 : num 7.36 7.33 7.36 7.46 7.23 ...
## $ swang1 : Factor w/ 2 levels "No RHC","RHC": 1 2 2 1 2 1 1 1 1 2 ...
## $ wtkilo1 : num 64.7 45.7 0 54.6 78.4 ...
## $ dnr1 : Factor w/ 2 levels "No","Yes": 1 1 1 1 2 1 1 1 1 1 ...
## $ ninsclas: Factor w/ 6 levels "Medicaid","Medicare",..: 2 6 5 6 2 2 5 5 5 1 ...
## $ resp : Factor w/ 2 levels "No","Yes": 2 1 1 2 1 2 1 2 1 1 ...
## $ card : Factor w/ 2 levels "No","Yes": 2 1 2 1 2 1 1 1 1 1 ...
## $ neuro : Factor w/ 2 levels "No","Yes": 1 1 1 1 1 1 1 2 1 2 ...
## $ gastr : Factor w/ 2 levels "No","Yes": 1 1 1 1 1 1 1 1 1 2 ...
## $ renal : Factor w/ 2 levels "No","Yes": 1 1 1 1 1 1 2 1 1 1 ...
## $ meta : Factor w/ 2 levels "No","Yes": 1 1 1 1 1 1 1 1 1 1 ...
## $ hema : Factor w/ 2 levels "No","Yes": 1 1 1 1 1 1 1 1 2 1 ...
## $ seps : Factor w/ 2 levels "No","Yes": 1 2 1 1 1 1 1 2 1 1 ...
## $ trauma : Factor w/ 2 levels "No","Yes": 1 1 1 1 1 1 1 1 1 1 ...
## $ ortho : Factor w/ 2 levels "No","Yes": 1 1 1 1 1 1 1 1 1 1 ...
## $ adld3p : int 0 NA NA NA NA 0 NA NA NA NA ...
## $ urin1 : num NA 1437 599 NA 64 ...
## $ race : Factor w/ 3 levels "black","other",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ income : Factor w/ 4 levels "> $50k","$11-$25k",..: 4 4 3 2 4 4 3 3 4 4 ...
## $ ptid : int 5 7 9 10 11 12 13 14 16 17 ...# Skapa en baseline-tabell med funktionen mytable() från paketet moonBook
# Vi gör ett urval av variabler som vi är intresserade av och delar upp tabellen på kön.
# Glöm inte installera moonBook om du saknar paketet
# install.packages("moonBook")
library(moonBook)
mytable(sex ~ age + race + income + alb1, data=rhc)##
## Descriptive Statistics by 'sex'
## _______________________________________________
## Female Male p
## (N=2543) (N=3192)
## -----------------------------------------------
## age 62.1 ± 17.1 60.8 ± 16.3 0.006
## race 0.000
## - black 465 (18.3%) 455 (14.3%)
## - other 157 ( 6.2%) 198 ( 6.2%)
## - white 1921 (75.5%) 2539 (79.5%)
## income 0.000
## - > $50k 170 ( 6.7%) 281 ( 8.8%)
## - $11-$25k 479 (18.8%) 686 (21.5%)
## - $25-$50k 378 (14.9%) 515 (16.1%)
## - Under $11k 1516 (59.6%) 1710 (53.6%)
## alb1 3.1 ± 0.9 3.1 ± 0.7 0.822
## -----------------------------------------------# Kontrollera hur utfallsmåttet är kodat
levels(rhc$death)## [1] "No" "Yes"Det verkar som att “death” är kodat som “Yes” eller “No”. Detta måste specificeras i den logistiska regressionen!
# kontrollera kategorier för kategoriska variabler. Kom ihåg att första kategorin i utskriften kommer användas som referenskategori vid regressionen
levels(rhc$sex)## [1] "Female" "Male"levels(rhc$race)## [1] "black" "other" "white"levels(rhc$income)## [1] "> $50k" "$11-$25k" "$25-$50k" "Under $11k"# Skapa en logistisk regression där död är utfallsmåttet
# logistisk regression görs med funktionen glm(), där argumentet family skall vara "binomial"
# Vi anger att det är händelsen death = Yes som vi vill relatera till kön, ålder, ras och inkomst
modell <- glm(death=="Yes" ~ sex + age + race + income, data=rhc, family="binomial")
# Rita en graf (forest plot) med resultaten. Detta görs med funktionen ORplot från moonBook
ORplot(modell)
# Kontrollera modellens goodness-of-fit med ROC-kurva. Du behöver installera paketet "Deducer" för detta.
# install.packages("Deducer")
library(Deducer)## Loading required package: ggplot2## Loading required package: JGR## Loading required package: rJava## Loading required package: JavaGD## Loading required package: iplots## Note: On Mac OS X we strongly recommend using iplots from within JGR.
## Proceed at your own risk as iplots cannot resolve potential ev.loop deadlocks.
## 'Yes' is assumed for all dialogs as they cannot be shown without a deadlock,
## also ievent.wait() is disabled.
## More recent OS X version do not allow signle-threaded GUIs and will fail.##
## Please type JGR() to launch console. Platform specific launchers (.exe and .app) can also be obtained at http://www.rforge.net/JGR/files/.## Loading required package: car## Loading required package: MASS##
##
## Note Non-JGR console detected:
## Deducer is best used from within JGR (http://jgr.markushelbig.org/).
## To Bring up GUI dialogs, type deducer().
rocplot(modell)