Inverse of the probability of censoring weights (IPCW) -harjoitus R:llä - monimuuttujamalli, johon ei koodata eri aikapisteissä mitattavaa, vaihtuvaa binaarista kovariaattia

Koodin tarkoitus

Luodaan esimerkkidata ja tehdään siitä survival-analyysi “inverse of the probability of censoring weights (IPCW)” -painotuksella.

Lähteet:

Grafféon ym. artikkeli (PMID: 31442762).
ipcwswitch-paketin helppidokumentti.

Tässä harjoituksessa toistetaan samat jutut kuin aiemmin julkaisemassani esimerkissä, seuraavin eroin:

ei analysoida seurannan aikana eri aikapisteissä mitattua kovariaattia
- sitä ei voi jättää ihan tyhjäksi, siitä lisää jäljempänä
kovariaateiksi vain baseline-muuttujia siis, mutta niitä muutama enemmän kuin aiemmassa esimerkissä

Huom! Nämä koodit ovat lähinnä omaa harjoitteluani varten, joten jos harjoittelet niillä, teet sen omalla vastuullasi.

Harjoittelussa tarvittavat datasetit

Latasin näitä harjoituksia varten seuraavat paketit sessioon. kableExtra-paketin latasin vain tämän web-sivun tekemiseen, sitä ei tarvitse koodia kokeillessa ladata. Käytän aina pacmania pakettien latailuun, joten se kannattaa ihan ensiksi asentaa, jos se ei jo ole R:ään asennettuna.

library(pacman)
p_load(dplyr, dlookr, ipcwswitch, ipw, ggplot2, survival, survminer, GGally, tidyr, kableExtra)

Esimerkkikoodi

Yllä mainitussa artikkelissa oli vain minimaalinen harjoitusdatasetti. Otetaan siitä mallia, mutta muokataan sen muotoa siten, että otetaan eri aikapisteissä seurannassa analysoitava kovariaatti (käytännössä) pois ja laitetaan muutama baseline-kovariaatti datasettiin lisää.

Huom! ipcwswitchin funktiot eivät toimi, jos seurannan aikana mitattava kovariaatti ei ole ollenkaan mukana mallissa. Siihen on olemassa ratkaisu, josta lisää jäljempänä.

Muuttujat:

id is the patient’s identifier,
randt is the date of the randomization visit,
lastdt is the date of latest news,
status is equal to 1 if the patient dies on lastdt (and 0 otherwise),
age is the patient’s age,
bas_chol is the patients cholesterol at baseline,
bas_sbp is the patients systolic blood pressure at baseline,
arm is the patient’s randomized arm,
swtrtdt is the date when the patient initiates the other arm treatment (NA if does not happen).
ps1, ps2, ps3 are the values of a covariate measured at dates randt, dt2, dt3, respectively

Kopioin manuaalisesti edellä mainitun artikkelin koodin tribbleksi, sen jälkeen lisäsin randomilla dataa, muuttuja kerrallaan, kunnes n oli 400. Esimerkkikoodia ei tarvitse sen enempää opiskella - riittää, että seuraavat rivit kopioi dataansa:

example_df <- tibble::tribble(
~id,~randt,~lastdt,~status,~age,~arm,~swtrtdt,
1L, "12.1.2018", "2.3.2018",    1L,  20L,  "A",  "1.3.2018",
2L, "4.11.2017", "15.12.2017",  1L,  50L,  "B",  NA,
3L, "20.5.2017", "4.1.2018",    0L,  40L,  "A",  NA
  )


# konvertoidaan pvm-sarakkeet oikeasti date-muotoisiksi:
example_df$randt <-   as.Date(example_df$randt, "%d.%m.%Y") 
example_df$lastdt <-  as.Date(example_df$lastdt, "%d.%m.%Y") 
example_df$swtrtdt <- as.Date(example_df$swtrtdt, "%d.%m.%Y") 

# konvertoidaan arm-sarake factor-muotoon:
example_df$arm <- as.factor(example_df$arm) 

# Luodaan lisää random-dataa.
#
# Ensin lisää tutkittavia (= id-numeroita) siten, että n=400:
example_df <- example_df %>% add_row(id = 4:400)

# "Uusille" tutkittaville pvm, jolloin randomoitu mukaan tutkimukseen:
example_df$randt[4:400] <- example_df$randt[1] + sample(1:30, 397, replace = T)

# Uusille tutkittaville pvm, jolloin viimeksi potilaasta kuultiin,
# muotoa randomointi-pvm + 150 + 1-30:
example_df <- example_df %>%
  mutate(lastdt = case_when(is.na(lastdt) ~ example_df$randt + 150
                            + sample(1:300, 400, replace = T),
                            !is.na(lastdt) ~ lastdt))

# Uusille tutkittaville arm-arvot (A vs. B):
example_df$arm[4:400] <- sample(c("A","B"), 397, replace = TRUE)

# Uusille tutkittaville random-iät:
example_df$age[4:400] <- sample(20:50, 397, replace = T)

# Seuraavaksi uusille tutkittaville pvm,
# jolloin mahdollisesti siirtyivät toiseen tutkimus-armiin
# Laitetaan kuitenkin näitä siirtyjiä vain armiin "A",
# jotta näemme, toimiiko painotus.
# Ensin random-apumuuttuja: jos saa arvon 3,
# tutkittava vaihtaa ryhmää kesken tutkimuksen:
example_df$apumuuttuja <- sample(1:3, 400, replace = T)

example_df <- example_df %>%
  mutate(swtrtdt = case_when(apumuuttuja == 3 & arm =="A" ~ example_df$lastdt
                             - sample(1:3, 400, replace = T)))

# Kuvitellaan asetelma, jossa armissa "A" myös henkilöitä kuolee enemmän.
example_df[4:400,] <- example_df[4:400,] %>%
                      mutate(status = case_when(apumuuttuja == 1 & arm =="A" ~ 1,
                                                apumuuttuja == 2 & arm =="A" ~ 1,
                                                apumuuttuja == 3 & arm =="A" ~ 0,
                    
                                                apumuuttuja == 1 & arm =="B" ~ 1,
                                                apumuuttuja == 2 & arm =="B" ~ 0,
                                                apumuuttuja == 3 & arm =="B" ~ 0  ))


# Luodaan randomit syst. verenpaineet tutkittaville:
example_df$bas_sbp <- rnorm(400, 120, 7) 



# Luodaan randomit kolesteroliarvot tutkittaville
example_df$bas_chol <- rnorm(400, 180, 25) * 0.02586


# Siivotaan pois apumuuttuja
example_df <- example_df %>% select (-apumuuttuja)

Ongelmallinen eri aikapisteissä mitattava ja mallin toimimiselle pakollinen kovariaatti - what to do?

ipcwswitch ei toimi (antaa kaikille weightin 1.0), mikäli ajan mittaan vaihteleva kovariaatti ei ole mukana jollakin muuttuvalla arvolla.

Havaintoja:

Jos kyseistä kovariaattia ei ole mallissa, koko ajo keskeytyy.
Jos laittaa kaikille saman arvon eri aikapisteisiin (tai vaikka yhteen ainoaan aikapisteeseen ja muihin NA), ipcwswitch antaa kaikille tutkittaville painoarvon 1.0
eli tutkittaville on saatava ero aikaiseksi kys. kovariaatin suhteen

Oma ratkaisuni on, että mikäli tutkimuksessa ei oikeasti ole eri aikapisteissä mitattavaa binaarista kovariaattia, jokin baselinen binaarinen kovariaatti koodataan aikapisteessä nolla kyseiseksi mallista vielä puuttuvaksi kovariaatiksi. Täten ipcwswitch saadaan ajettua aidolla datalla läpi, siten että jokainen tutkittava saa häneltä oikeasti mitatun, todellisen binaarisen muuttuja-arvon (baselinessä). Joka tapauksessa ensimmäinen tällainen mittaus ipcwswitch:n dokumentaation mallin mukaisesti pitäisikin olla nollapisteessä (“randomointihetki”), joten jokin baselinen binaarinen muuttuja mielestäni soveltuu tähän tehtävään.

Koodataan seuraavassa “chunkissa” binaarinen muuttuja tupakointi kyseiseksi muuttujaksi.

# Alla seurannassa useassa eri aikapisteessä mitattava muuttuja,
# joka on pakko tavalla tai toisella olla mallissa, muuten malli ei toimi.

# Tässä esimerkissä valitaan kyseiseksi binaariseksi kovariaatiksi tupakointi.
# Baselinessä mitattu (haastateltu) arvo tupakoinnille saa arvon 0 tai 1.
# Tässä esimerkkidatassa kyseinen arvo tulee randomina.

example_df$ps1 <- sample(0:1, 400, replace = T)


# Muut kyseisen muuttujan (tupakointi) arvot ja aikapisteet koodataan NA:ksi.

example_df$ps2 <- NA
example_df$ps3 <- NA
example_df$dt2 <- NA
example_df$dt3 <- NA

ipcwswitch hyväksyy vain data.framen datan muodoksi

Havaitsin, että ipcwswitch ilmoittaa virheen ja keskeyttää ajon, jos data tarjotaan sille tibblenä.

Niinpä on tärkeä varmistua siitä, että datan muoto on oikea, ks. alla oleva chunk.

# ipswswitch nikottelee tibblejen kanssa,
# joten on muutettava datasettimme data.frame-muotoon:
example_df <- as.data.frame(example_df)

Pienemmän samplen ottaminen datastamme

Jos tahtoisimme pienemmän samplen tästä datasta, tekisimme sen seuraavilla käskyillä. Koska datamme on verraten pieni (n=400), emme sitä tällä kertaa tee.

#index <- sample(nrow(example_df), 100)

#example_df <- example_df[index, ]

Silmäillään esimerkkidataamme

Silmäillään hieman esimerkkidataa, varmistetaan, että kaikki kunnossa (vain muutama esimerkkirivi).

id	randt	lastdt	status	age	arm	swtrtdt	bas_sbp	bas_chol	ps1	ps2	ps3	dt2	dt3
1	2018-01-12	2018-03-02	1	20	A	2018-02-27	126.3828	5.998195	1	NA	NA	NA	NA
2	2017-11-04	2017-12-15	1	50	B	NA	120.0812	4.132616	0	NA	NA	NA	NA
3	2017-05-20	2018-01-04	0	40	A	NA	108.1852	4.795587	1	NA	NA	NA	NA
4	2018-01-21	2018-09-16	1	31	A	NA	114.4214	4.553998	0	NA	NA	NA	NA
5	2018-01-24	2018-10-23	0	46	B	NA	121.0980	4.707438	1	NA	NA	NA	NA
6	2018-01-13	2018-08-03	1	26	B	NA	124.2205	5.194685	0	NA	NA	NA	NA

Laajennetun esimerkkidatan struktuurin varmennus

On parasta katsoa myös esimerkkidatan struktuuri. ipswswitch nikottelee heti vastaan, jos esim. kvalitatiiviset muuttujat eivät ole tyyppiä factor. Kaikki näyttää alla olevassa listauksessa olevan kunnossa.

str(example_df)

## 'data.frame':    400 obs. of  14 variables:
##  $ id      : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ randt   : Date, format: "2018-01-12" "2017-11-04" ...
##  $ lastdt  : Date, format: "2018-03-02" "2017-12-15" ...
##  $ status  : int  1 1 0 1 0 1 1 0 0 0 ...
##  $ age     : int  20 50 40 31 46 26 37 39 48 37 ...
##  $ arm     : Factor w/ 2 levels "A","B": 1 2 1 1 2 2 1 1 1 2 ...
##  $ swtrtdt : Date, format: "2018-02-27" NA ...
##  $ bas_sbp : num  126 120 108 114 121 ...
##  $ bas_chol: num  6 4.13 4.8 4.55 4.71 ...
##  $ ps1     : int  1 0 1 0 1 0 1 0 1 1 ...
##  $ ps2     : logi  NA NA NA NA NA NA ...
##  $ ps3     : logi  NA NA NA NA NA NA ...
##  $ dt2     : logi  NA NA NA NA NA NA ...
##  $ dt3     : logi  NA NA NA NA NA NA ...

Survival-analyysin valmistelu: ensin timesTokeep-komennolla käsittely

Seuraavaksi ajetaan laajennettu toydata monen eri esivaiheen kautta, jotta sille voidaan lopuksi laskea painotukset ja tehdä aivan lopuksi analyysit.

kept.t <- timesTokeep(example_df, id       = "id",
                      tstart   = "randt",
                      tstop    = "lastdt",
                      mes.cov  = list(c("ps1", "ps2", "ps3")),
                      time.cov = list(c("randt", "dt2", "dt3")))

Esimerkiksi kolmannen tutkittavan päivämäärät näyttävät seuraavalta:

kept.t[[1]][[3]]

## [1] "2017-05-20" "2018-01-04"

Seuraavaksi long format muotoon muuttaminen

Long format -muoto aikaansaadaan wideToLongTDC-komennolla:

example_df.long <- wideToLongTDC(data = example_df, id = "id",
                          tstart = "randt", tstop = "lastdt",
                          event = "status",
                          bas.cov = c("age", "arm", "swtrtdt", "bas_sbp", "bas_chol"),
                          mes.cov = list(TDconf = c("ps1", "ps2", "ps3")),
                          time.cov = list(c("randt", "dt2", "dt3")),
                          times = kept.t[[1]])

Tämän jälkeen data näyttää seuraavalta. Huomaa, että muuttujaan TDconf tallentuu kovariaatin “ps” kaikki toistetut arvot, aggregaattina. Tässä esimerkissä toki tämä kovariaatti “ps” sai vain NA-arvoja, lukuun ottamatta ensimmäistä arvoa, joka kaikilla 1 (muuten koko funktio menettää toimintakykynsä ja kaikki saavat weightiksi 1).

id	tstart	tstop	event	age	arm	swtrtdt	bas_sbp	bas_chol	TDconf
1	2018-01-12	2018-03-02	1	20	A	2018-02-27	126.3828	5.998195	1
2	2017-11-04	2017-12-15	1	50	B	NA	120.0812	4.132616	0
3	2017-05-20	2018-01-04	0	40	A	NA	108.1852	4.795587	1
4	2018-01-21	2018-09-16	1	31	A	NA	114.4214	4.553998	0
5	2018-01-24	2018-10-23	0	46	B	NA	121.0980	4.707438	1
6	2018-01-13	2018-08-03	1	26	B	NA	124.2205	5.194685	0

Päivämäärät numeeriseen muotoon, jolloin tstart-muuttuja pisteessä 0

Tämä tehdään copy-pastettamalla seuraavat rivit koodiin.

# Put dates in numeric format with tstart at 0
example_df.long$tstart <- as.numeric(example_df.long$tstart)
example_df.long$tstop <- as.numeric(example_df.long$tstop)
example_df.long$swtrtdt <- as.numeric(example_df.long$swtrtdt)
tabi <- split(example_df.long, example_df.long$id)
L.tabi <- length(tabi)
tablist <- lapply(1:L.tabi, function(i){
  refstart <- tabi[[i]]$tstart[1]
  tabi[[i]]$tstart <- tabi[[i]]$tstart - refstart
  tabi[[i]]$tstop <- tabi[[i]]$tstop - refstart
  tabi[[i]]$swtrtdt <- tabi[[i]]$swtrtdt - refstart
  return(tabi[[i]])
})
example_df.long <- do.call( rbind, tablist )

Ryhmää vaihtavien tutkittavien sensorointi

Tämä tapahtuu copy-pastettamalla seuraava koodinpätkä mukaan:

# Patients are censored when initiating the other arm treatment, that is, at time swtrtdt
example_df.long2 <- cens.ipw(example_df.long, id = "id", tstart = "tstart", tstop = "tstop",
                      event = "event", arm = "arm",
                      realtrt = FALSE, censTime ="swtrtdt")

Datasetti ennen sensurointia ja sen jälkeen

Seuraavassa esitettynä muutamia rivejä datasetistä, ennen (ensimmäinen taulukko) ja jälkeen (jälkimmäinen taulukko) yllä mainitun sensuroinnin.

id	tstop	event	age	arm	swtrtdt	bas_sbp	bas_chol	TDconf
1	49	1	20	A	46	126.3828	5.998195	1
2	41	1	50	B	NA	120.0812	4.132616	0
3	229	0	40	A	NA	108.1852	4.795587	1
4	238	1	31	A	NA	114.4214	4.553998	0
5	272	0	46	B	NA	121.0980	4.707438	1
6	202	1	26	B	NA	124.2205	5.194685	0

id	tstop	event	age	arm	swtrtdt	bas_sbp	bas_chol	TDconf	cens
1	46	0	20	A	46	126.3828	5.998195	1	1
2	41	1	50	B	NA	120.0812	4.132616	0	0
3	229	0	40	A	NA	108.1852	4.795587	1	0
4	238	1	31	A	NA	114.4214	4.553998	0	0
5	272	0	46	B	NA	121.0980	4.707438	1	0
6	202	1	26	B	NA	124.2205	5.194685	0	0

Kaikkien event-aikojen yhdistäminen yhteen ja samaan dataan

Painotukset pitää saada kaikille event- ja treatment-sensoroinnin ajoille, joten seuraavat rivit täytyy ajaa.

rep.times1 <- unique(c(example_df.long2$tstop[example_df.long2$cens==1 & example_df.long2$arm == "A"],
                       example_df.long2$tstop[example_df.long2$event==1]))
rep.times2 <- unique(c(example_df.long2$tstop[example_df.long2$cens==1 & example_df.long2$arm == "B"],
                       example_df.long2$tstop[example_df.long2$event==1]))

# to put times in same order as arms levels
# (Huom. tätä riviä ei ollut referoidussa artikkelissa,
# mutta R-dokumentaation mukaan se kannattaa ajaa.) 
levels(example_df.long2[, "arm"])

## [1] "A" "B"

# Now, we can replicate the rows
example.rep <- replicRows(example_df.long2, tstart = "tstart", tstop = "tstop",
                      event = "event", cens = "cens",
                      times1 = rep.times1, times2 = rep.times2,
                      arm = "arm")

Lopputuloksena seuraavan näköinen datasetti (vain muutama esimerkkirivi):

id	age	arm	swtrtdt	bas_sbp	bas_chol	TDconf	cens	tstart	tstop
1	20	A	46	126.3828	5.998195	1	0	0	41
1	20	A	46	126.3828	5.998195	1	1	41	46
3	40	A	NA	108.1852	4.795587	1	0	0	41
3	40	A	NA	108.1852	4.795587	1	0	41	46
3	40	A	NA	108.1852	4.795587	1	0	46	152
3	40	A	NA	108.1852	4.795587	1	0	152	153

Vihdoinkin painotusten laskeminen

Painotukset lasketaan komennolla ipcw:

example.rep <- ipcw(example.rep, id = "id", tstart = tstart, tstop = tstop,
                cens = cens, arm = "arm", bas.cov = c("age", "bas_sbp", "bas_chol"),
                conf = c("TDconf"), trunc = 0.01, type = 'kaplan-meier')

Tämä analyysi heittää warningseja runsaasti, mutta ajoin oleellisesti ihan samat analyysit R:n helppidokumentissa käytetyllä esimerkkidatasetillä (data(“SHIdat”)) ja samassa dokumentissa listatulla koodilla läpi, ja aivan samantyyppiset warningit näyttivät silläkin tulevan.

Saadaan aikaiseksi seuraavalta näyttävä datasetti (vain muutama esimerkkiriviä):

id	age	arm	swtrtdt	bas_sbp	bas_chol	TDconf	cens	weights	tstart	tstop	weights.trunc
1	20	A	46	126.3828	5.998195	1	0	1.000000	0	41	1.000000
1	20	A	46	126.3828	5.998195	1	1	1.000000	41	46	1.000000
3	40	A	NA	108.1852	4.795587	1	0	1.000000	0	41	1.000000
3	40	A	NA	108.1852	4.795587	1	0	1.000000	41	46	1.000000
3	40	A	NA	108.1852	4.795587	1	0	0.999901	46	152	0.999901
3	40	A	NA	108.1852	4.795587	1	0	0.999901	152	153	0.999901

Painotusten distribuution plottaaminen

Luodut painotukset pystyy tämän jälkeen plottaamaan seuraavasti. (Allekirjoittaneelle on epäselvää, kuinka yleistä on, että näitä julkaistaan esim. supplemental-osuudessa.)

Un-truncated stabilized weights:

ipwplot(weights = example.rep$weights,
        timevar = example.rep$tstart,
        binwidth = 50,
        logscale = T,
        xlab = "Time since enrolment (days)",
        ylab = "Logarithm of the stabilized\n un-truncated weights")

Truncated stabilized weights:

ipwplot(weights = example.rep$weights,
        timevar = example.rep$tstart,
        binwidth = 50,
        logscale = T,
        xlab = "Time since enrolment (days)",
        ylab = "Logarithm of the stabilized\n un-truncated weights")

Lopultakin Coxin MSM-mallin kokoaminen painotusten kera

Nyt kun painotukset on vihdoin viimein saatu laskettua, ne saadaan otettua Coxin MSM-malliin (marginal structural models) mukaan. Käytetään coxph-funktiota paketista survival, ja painotukset määritetään funktion argumenteissa.

Huomaa seuraavat 2 asiaa:

tässä käytettiin painoina truncated stabilized weights
tupakointia (muuttuja ps1) ei kirjoiteta tähän malliin vaan sen vaikutus on jo weights-määreessä “sisällä” (näin coxph-malli kasattiin ipcwswitch:n dokumentaatiossa).

fit.stab.w <- coxph(Surv(tstart, tstop, event) ~ age + bas_sbp + bas_chol + arm + cluster(id),
                    data = example.rep, weights = example.rep$weights.trunc)

Mallin tulos saadaan laajemmin esiin komennolla summary:

summary(fit.stab.w)

## Call:
## coxph(formula = Surv(tstart, tstop, event) ~ age + bas_sbp + 
##     bas_chol + arm, data = example.rep, weights = example.rep$weights.trunc, 
##     cluster = id)
## 
##   n= 29828, number of events= 187 
## 
##               coef exp(coef)  se(coef) robust se      z Pr(>|z|)    
## age       0.007858  1.007889  0.008087  0.007367  1.067 0.286120    
## bas_sbp   0.001190  1.001191  0.011054  0.012233  0.097 0.922513    
## bas_chol  0.151413  1.163477  0.128099  0.126491  1.197 0.231297    
## armB     -0.585240  0.556972  0.153143  0.154096 -3.798 0.000146 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##          exp(coef) exp(-coef) lower .95 upper .95
## age          1.008     0.9922    0.9934    1.0225
## bas_sbp      1.001     0.9988    0.9775    1.0255
## bas_chol     1.163     0.8595    0.9080    1.4908
## armB         0.557     1.7954    0.4118    0.7534
## 
## Concordance= 0.609  (se = 0.023 )
## Likelihood ratio test= 19.17  on 4 df,   p=7e-04
## Wald test            = 18.28  on 4 df,   p=0.001
## Score (logrank) test = 19.16  on 4 df,   p=7e-04,   Robust = 18.91  p=8e-04
## 
##   (Note: the likelihood ratio and score tests assume independence of
##      observations within a cluster, the Wald and robust score tests do not).

Vertailun vuoksi sama Coxin MSM-malli ilman painotuksia

fit.stab.w.weightless  <- coxph(Surv(tstart, tstop, event) ~ age + bas_sbp
                          + bas_chol + arm + cluster(id),
                          data = example.rep)

summary(fit.stab.w.weightless)

## Call:
## coxph(formula = Surv(tstart, tstop, event) ~ age + bas_sbp + 
##     bas_chol + arm, data = example.rep, cluster = id)
## 
##   n= 29828, number of events= 187 
## 
##               coef exp(coef)  se(coef) robust se      z Pr(>|z|)    
## age       0.007853  1.007884  0.008089  0.007369  1.066 0.286603    
## bas_sbp   0.001164  1.001165  0.011055  0.012218  0.095 0.924082    
## bas_chol  0.151848  1.163983  0.128083  0.126462  1.201 0.229851    
## armB     -0.585385  0.556891  0.153154  0.154101 -3.799 0.000145 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##          exp(coef) exp(-coef) lower .95 upper .95
## age         1.0079     0.9922    0.9934    1.0225
## bas_sbp     1.0012     0.9988    0.9775    1.0254
## bas_chol    1.1640     0.8591    0.9085    1.4914
## armB        0.5569     1.7957    0.4117    0.7533
## 
## Concordance= 0.609  (se = 0.023 )
## Likelihood ratio test= 19.18  on 4 df,   p=7e-04
## Wald test            = 18.29  on 4 df,   p=0.001
## Score (logrank) test = 19.17  on 4 df,   p=7e-04,   Robust = 18.92  p=8e-04
## 
##   (Note: the likelihood ratio and score tests assume independence of
##      observations within a cluster, the Wald and robust score tests do not).

Survival-käyrän plottaaminen - tutkimuksen koko populaatio

Survival-käyrän plottaamista ei käsitelty yllä referoidussa lähdeartikkelissa. Sain kuitenkin käyrää aikaiseksi kaikkien tutkittavien ryhmästä (sekä arm A että B) seuraavalla tavalla:

fit_all <- survfit(fit.stab.w)
ggsurvplot(fit_all, data=example_df, color = "#2E9FDF",
           ggtheme = theme_minimal())

## Warning: Now, to change color palette, use the argument palette= '#2E9FDF'
## instead of color = '#2E9FDF'

Survival-käyrän plottaaminen erikseen arm A:sta ja B:stä

Erityistä pohdintaa aiheutti tässä otsikossa esitetty kysymys. Ryhmien erottamiseksi toisistaan omiksi survival-käyrikseen täytyy käsitykseni mukaan tehdä “apu-dataframe”, johon kuvaillaan millintarkasti, mitä muuttujia tahdotaan plotata. Jos ja kun erotellaan ryhmät erikseen, pitää niiden levelien olla millintarkasti plotattavaa datasettiä vastaavat, ja kyseinen apumuuttuja täytyy muuttaa factor-muotoon. Jos joukossa on muita kovariaatteja, niiden meania käytetään survival-käyrän mallissa:

# Apudatasetin luominen 
apudf <- with(example_df,
              data.frame(arm = c("A", "B"), 
                         age = rep(mean(age, na.rm = TRUE), 2),
                         bas_sbp = rep(mean(bas_sbp, na.rm = TRUE), 2),
                         bas_chol = rep(mean(bas_chol, na.rm = TRUE), 2)))

apudf$arm <- as.factor(apudf$arm) 

apudf

##   arm     age  bas_sbp bas_chol
## 1   A 35.7575 119.8857 4.630074
## 2   B 35.7575 119.8857 4.630074

fit2 <- survfit(fit.stab.w, newdata = apudf)

ggsurvplot(fit2,
           data=example_df,
           conf.int = TRUE,
           legend.labs=c("Arm = A", "Arm = B"),
           #risk.table = T,
           ggtheme = theme_minimal())

Siistimmän survival-käyrän plottaaminen erikseen arm A:sta ja B:stä - alkuasetukset

Mielestäni ggsurvplot tuottaa rujonnäköisen survival-käyrän. Luu esittää tapoja, joilla käyrää saa kauniimmaksi. Ne eivät toimi tässä suoraan, koska survfit-komennolla käsiteltiin omassa esimerkissämme coxph:sta tullut fit.

Käytämme ensimmäiseksi Luun esittämää metodia, jolla ihan aluksi voidaan leikata (truncate) luomaamme fitiä ja siis kuvaajan oikeaa laitaa, jos esimerkiksi koko toinen kohortti on jo menehtynyt esimerkiksi aikapisteessä 480. Tämä ei tietenkään ole pakollista, mutta esitellään tämäkin tekniikka tässä.

# Tehdään summary fit2:sta ja tallennetaan s:ksi 

s <- summary(fit2, times = seq(0, 480, 1), extend = T)

Fitit pakkaavat olemaan list-tyyppisiä, joten niiden tarkastelu on hiukan erilaista kuin dataframejen. Mutta jos katsomme tekemäämme fitiä esim. komennolla View(fit2), voimme nähdä, että listan kohta time on muotoa double[x-luku] kun kohdassa surv taas näkyy double[x-luku x 2]. Tästä voimme päätellä, että kutakin yhteistä aikapistettä kohden kohdassa surv on kullekin stratalle omat informaationsa.

Yhtä hyvä ajatus on tarkastella listaa komennolla str. Tämä on hiukan sotkuisempaa mielestäni lukea kuin komennolla View aikaansaatu näkymä, mutta ihan samat tiedot kuin edellisessä kappaleessa kuvailtiin nähdään kohdissa time ja surv:

str(s)

## List of 16
##  $ n            : int 29828
##  $ time         : num [1:481] 0 1 2 3 4 5 6 7 8 9 ...
##  $ n.risk       : num [1:481] 400 400 400 400 400 400 400 400 400 400 ...
##  $ n.event      : num [1:481] 0 0 0 0 0 0 0 0 0 0 ...
##  $ n.censor     : num [1:481] 0 0 0 0 0 0 0 0 0 0 ...
##  $ surv         : num [1:481, 1:2] 1 1 1 1 1 1 1 1 1 1 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:2] "1" "2"
##  $ cumhaz       : num [1:481, 1:2] 0 0 0 0 0 0 0 0 0 0 ...
##  $ std.err      : num [1:481, 1:2] 0 0 0 0 0 0 0 0 0 0 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:2] "1" "2"
##  $ logse        : logi TRUE
##  $ lower        : num [1:481, 1:2] 1 1 1 1 1 1 1 1 1 1 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:2] "1" "2"
##  $ upper        : num [1:481, 1:2] 1 1 1 1 1 1 1 1 1 1 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:2] "1" "2"
##  $ conf.type    : chr "log"
##  $ conf.int     : num 0.95
##  $ call         : language survfit(formula = fit.stab.w, newdata = apudf)
##  $ table        : num [1:2, 1:9] 29828 29828 400 400 400 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:2] "1" "2"
##   .. ..$ : chr [1:9] "records" "n.max" "n.start" "events" ...
##  $ rmean.endtime: num 450
##  - attr(*, "class")= chr "summary.survfit"

Listasta saa kopioitua elementtejä “ulos” painamalla RStudiossa kyseisen listan View-näkymässä kunkin elementin kohdalla oikeassa laidassa olevaa nappia. Listan elementti hakasulkuineen siirtyy konsoliin, jolloin se voidaan edelleen siirtää R-skriptiin. Katso tarkemmin kuva alla.

Listan elementteihin niiden indekseillä viittaaminen RStudiossa

Näin toimimalla voimme luoda dataframen, jolle annamme nimeksi surv_a_and_b. “Stratat” eli surv-elementtien kolumnit eivät saa mitään järkevää nimeä, joten nimeämme ne manuaalisesti dplyr:n rename-komennolla muotoon “A” ja “B”:

# Otetaan listasta s irti oleelliset tiedot eli aika ja kummankin ryhmän survivalit

surv_a_and_b      <- as.data.frame(s[["surv"]])
surv_a_and_b$time <- s[["time"]]

surv_a_and_b <- surv_a_and_b %>% rename(A = 1, B = 2)

tail(surv_a_and_b)

##              A         B time
## 476 0.01880586 0.1093526  475
## 477 0.01880586 0.1093526  476
## 478 0.01880586 0.1093526  477
## 479 0.01880586 0.1093526  478
## 480 0.01880586 0.1093526  479
## 481 0.01880586 0.1093526  480

Siistimmät survival-käyrät ggplotilla - loppusilaus

Ennen kuin muotoilemme käyrät ggplotilla haluttuun muotoon, data kannattaa muuttaa long-muotoon, sillä värien ja legendien tekeminen on rutkasti hankalampaa muussa tapauksessa. Long-muoto saavutetaan kuitenkin helposti tidyr-paketin komennolla gather.

surv_a_and_b_long <- gather(surv_a_and_b, "strata", "surv", -time)

head(surv_a_and_b_long)

##   time strata surv
## 1    0      A    1
## 2    1      A    1
## 3    2      A    1
## 4    3      A    1
## 5    4      A    1
## 6    5      A    1

tail(surv_a_and_b_long)

##     time strata      surv
## 957  475      B 0.1093526
## 958  476      B 0.1093526
## 959  477      B 0.1093526
## 960  478      B 0.1093526
## 961  479      B 0.1093526
## 962  480      B 0.1093526

Nyt data on vihdoin muodossa, joka on helppo syöttää ggplot-komennolle. Värit ja viivatyypit määrittyvät automaattisesti, kun aes-kohdissa ohjataan color ja linetype vaihtelemaan muuttujan strata mukaan:

ggplot(surv_a_and_b_long, aes(x = time, y = surv, color = strata)) +
  geom_step(aes(linetype = strata), size = 1) +
  theme_classic() +
  labs(x = "Time",
       y = "Probability of survival")

Jos tahdotaan kokonaan mustavalkoinen kuvaaja, määritetään geom_step-komennolle color = "black" - huom. tämä on tehtävä aes-funktion (sulkumerkkien) ulkopuolella:

ggplot(surv_a_and_b_long, aes(x = time, y = surv, color = strata)) +
  geom_step(aes(linetype = strata), size = 1, color="black") +
  theme_classic() +
  labs(x = "Time",
       y = "Probability of survival")