Inverse of the probability of censoring weights (IPCW) -harjoitus R:llä

Koodin tarkoitus

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

Käytännössä tehdään sama prosessi mutta isommalla datalla kuin mitä tehtiin artikkelissa Grafféon ja kumppaneiden artikkelissa (PMID: 31442762).

Artikkelin koodissa oli joitakin pieniä virheitä, ja koodia piti paikata samaisten henkilöiden tekemästä R:n helppidokumentista.

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. 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)

Esimerkkikoodi

Yllä mainitussa artikkelissa oli vain minimaalinen harjoitusdatasetti, joka näkyy alla. Sen muoto/muuttujat ovat seuraavanlaiset (lähde: Grafféon ym. (PMID: 31442762):

This Toy Example Data is given in wide format. It contains the following variables:

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,
ps1, ps2, ps3 are the values of a covariate measured at dates randt, dt2, dt3, respectively
arm is the patient’s randomized arm,
swtrtdt is the date when the patient initiates the other arm treatment (NA if does not happen).

id	randt	lastdt	status	age	ps1	ps2	ps3	dt2	dt3	arm	swtrtdt
1	2018-01-12	2018-03-02	1	20	0	0	0	2018-02-02	2018-03-01	A	2018-03-01
2	2017-11-04	2017-12-15	1	50	1	NA	2	NA	2017-12-12	B	NA
3	2017-05-20	2018-01-04	0	40	0	0	1	2017-08-02	2018-01-02	A	NA

Esimerkkikoodin laajennus ad n=400

Kopioin manuaalisesti 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:

toydata <- tibble::tribble(
~id,~randt,~lastdt,~status,~age,~ps1,~ps2,~ps3,~dt2,~dt3,~arm,~swtrtdt,
1L, "12.1.2018", "2.3.2018", 1L, 20L, 0L, 0L, 0L, "2.2.2018", "1.3.2018", "A", "1.3.2018",
2L, "4.11.2017", "15.12.2017", 1L, 50L, 1L, 1L, 2L, "1.12.2017", "12.12.2017", "B", NA,
3L, "20.5.2017", "4.1.2018", 0L, 40L, 0L, 0L, 1L, "2.8.2017", "2.1.2018", "A", NA
  )


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

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

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

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

# Uusille tutkittaville ensimmäinen "labrakontrollipvm",
# muotoa randomointi-pvm + 60 + 1-30:
toydata <- toydata %>%
  mutate(dt2 = case_when(is.na(dt2) ~ toydata$randt + 60
                         + sample(1:30, 400, replace = T),
                         !is.na(dt2) ~ dt2))

# Uusille tutkittaville toinen "labrakontrollipvm",
# muotoa randomointi-pvm + 100 + 1-30:
toydata <- toydata %>%
  mutate(dt3 = case_when(is.na(dt3) ~ toydata$randt + 100
                         + sample(1:30, 400, replace = T),
                         !is.na(dt3) ~ dt3))

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

# Uusille tutkittaville puuttuvat "labratulokset" (kovariaatit):
toydata$ps1[4:400] <- sample(0:2, 397, replace = T)
toydata$ps2[4:400] <- sample(0:2, 397, replace = T)
toydata$ps3[4:400] <- sample(0:2, 397, replace = T)

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

# Uusille tutkittaville random-iät:
toydata$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:
toydata$apumuuttuja <- sample(1:3, 400, replace = T)

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

# Kuvitellaan asetelma, jossa armissa "A" myös henkilöitä kuolee enemmän.
toydata[4:400,] <- toydata[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  ))


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

Laajennetun esimerkkidatan silmäily

Silmäillään hieman esimerkkidataa, varmistetaan, että kaikki kunnossa. Katsotaan head-komennolla 10 ensimmäistä riviä:

id	randt	lastdt	status	age	ps1	ps2	ps3	dt2	dt3	arm	swtrtdt	apumuuttuja
1	2018-01-12	2018-03-02	1	20	0	0	0	2018-02-02	2018-03-01	A	NA	2
2	2017-11-04	2017-12-15	1	50	1	1	2	2017-12-01	2017-12-12	B	NA	3
3	2017-05-20	2018-01-04	0	40	0	0	1	2017-08-02	2018-01-02	A	NA	1
4	2018-01-26	2019-02-20	1	26	2	2	1	2018-04-16	2018-05-16	A	NA	1
5	2018-01-18	2018-10-19	0	44	1	0	1	2018-04-15	2018-05-22	A	2018-10-18	3
6	2018-02-05	2018-11-19	0	20	2	2	1	2018-04-09	2018-05-25	A	2018-11-16	3

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(toydata)

## 'data.frame':    400 obs. of  13 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 0 0 1 0 0 ...
##  $ age        : int  20 50 40 26 44 20 50 27 42 40 ...
##  $ ps1        : int  0 1 0 2 1 2 2 2 0 0 ...
##  $ ps2        : int  0 1 0 2 0 2 1 0 0 0 ...
##  $ ps3        : int  0 2 1 1 1 1 1 1 1 1 ...
##  $ dt2        : Date, format: "2018-02-02" "2017-12-01" ...
##  $ dt3        : Date, format: "2018-03-01" "2017-12-12" ...
##  $ arm        : Factor w/ 2 levels "A","B": 1 2 1 1 1 1 2 2 2 2 ...
##  $ swtrtdt    : Date, format: NA NA ...
##  $ apumuuttuja: int  2 3 1 1 3 3 2 1 3 2 ...

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(toydata, id       = "id",
                      tstart   = "randt",
                      tstop    = "lastdt",
                      mes.cov  = list(c("ps1", "ps2", "ps3")),
                      time.cov = list(c("randt", "dt2", "dt3")))

Esimerkiksi tutkittavan id=3 päivämäärät näyttävät seuraavalta:

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

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

Seuraavaksi long format muotoon muuttaminen

Long format -muoto aikaansaadaan wideToLongTDC-komennolla:

toy.long <- wideToLongTDC(data = toydata, id = "id",
                          tstart = "randt", tstop = "lastdt",
                          event = "status",
                          bas.cov = c("age", "arm", "swtrtdt"),
                          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 (tässä näytetään 10 ensimmäistä riviä). Huomaa, että muuttujaan TDconf tallentuu kaikkien kovariaattien ps1, ps2 ja p3 toistetut arvot, aggregaattina.

id	tstart	tstop	event	age	arm	swtrtdt	TDconf
1	2018-01-12	2018-03-02	1	20	A	NA	0
2	2017-11-04	2017-12-12	0	50	B	NA	1
2	2017-12-12	2017-12-15	1	50	B	NA	2
3	2017-05-20	2018-01-02	0	40	A	NA	0
3	2018-01-02	2018-01-04	0	40	A	NA	1
4	2018-01-26	2018-05-16	0	26	A	NA	2

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
toy.long$tstart <- as.numeric(toy.long$tstart)
toy.long$tstop <- as.numeric(toy.long$tstop)
toy.long$swtrtdt <- as.numeric(toy.long$swtrtdt)
tabi <- split(toy.long, toy.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]])
})
toy.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
toy.long2 <- cens.ipw(toy.long, id = "id", tstart = "tstart", tstop = "tstop",
                      event = "event", arm = "arm",
                      realtrt = FALSE, censTime ="swtrtdt")

Datasetti ennen sensurointi ja sen jälkeen

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

id	tstart	tstop	event	age	arm	swtrtdt	TDconf
1	0	49	1	20	A	NA	0
2	0	38	0	50	B	NA	1
2	38	41	1	50	B	NA	2
3	0	227	0	40	A	NA	0
3	227	229	0	40	A	NA	1
4	0	110	0	26	A	NA	2

id	tstart	tstop	event	age	arm	swtrtdt	TDconf
1	0	49	1	20	A	NA	0
2	0	38	0	50	B	NA	1
2	38	41	1	50	B	NA	2
3	0	227	0	40	A	NA	0
3	227	229	0	40	A	NA	1
4	0	110	0	26	A	NA	2

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(toy.long2$tstop[toy.long2$cens==1 & toy.long2$arm == "A"],
                       toy.long2$tstop[toy.long2$event==1]))
rep.times2 <- unique(c(toy.long2$tstop[toy.long2$cens==1 & toy.long2$arm == "B"],
                       toy.long2$tstop[toy.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(toy.long2[, "arm"])

## [1] "A" "B"

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

Lopputuloksena seuraavan näköinen datasetti:

id	age	arm	swtrtdt	tstart	tstop	event
1	20	A	NA	0	41	0
1	20	A	NA	41	49	1
3	40	A	NA	0	41	0
3	40	A	NA	41	49	0
3	40	A	NA	49	151	0
3	40	A	NA	151	153	0

Vihdoinkin painotusten laskeminen

Painotukset lasketaan komennolla ipcw:

toy.rep <- ipcw(toy.rep, id = "id", tstart = tstart, tstop = tstop,
                cens = cens, arm = "arm", bas.cov = c("age"),
                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:

id	age	arm	swtrtdt	weights	tstart	tstop	event	weights.trunc
1	20	A	NA	1.0000000	0	41	0	1.0000000
1	20	A	NA	1.0000000	41	49	1	1.0000000
3	40	A	NA	1.0000000	0	41	0	1.0000000
3	40	A	NA	1.0000000	41	49	0	1.0000000
3	40	A	NA	1.0000000	49	151	0	1.0000000
3	40	A	NA	0.9996188	151	153	0	0.9996188

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 = toy.rep$weights,
        timevar = toy.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 = toy.rep$weights,
        timevar = toy.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.

Truncated stabilized weights:

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

Mallin tulos saadaan laajemmin esiin komennolla summary:

summary(fit.stab.w)

## Call:
## coxph(formula = Surv(tstart, tstop, event) ~ age + arm, data = toy.rep, 
##     weights = toy.rep$weights.trunc, cluster = id)
## 
##   n= 35655, number of events= 213 
## 
##           coef exp(coef)  se(coef) robust se      z Pr(>|z|)    
## age  -0.013133  0.986952  0.008139  0.007935 -1.655 0.097880 .  
## armB -0.571294  0.564794  0.146832  0.147006 -3.886 0.000102 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##      exp(coef) exp(-coef) lower .95 upper .95
## age     0.9870      1.013    0.9717    1.0024
## armB    0.5648      1.771    0.4234    0.7534
## 
## Concordance= 0.584  (se = 0.023 )
## Likelihood ratio test= 16.68  on 2 df,   p=2e-04
## Wald test            = 16.35  on 2 df,   p=3e-04
## Score (logrank) test = 16.48  on 2 df,   p=3e-04,   Robust = 17.76  p=1e-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=toydata, 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(toydata,
              data.frame(arm = c("A", "B"), 
                         age = rep(mean(age, na.rm = TRUE), 2)))

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

apudf

##   arm    age
## 1   A 34.005
## 2   B 34.005

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

ggsurvplot(fit2,
           data=toydata,
           conf.int = TRUE,
           legend.labs=c("Arm = A", "Arm = B"),
           ggtheme = theme_minimal())

id	tstart	tstop	event	age	arm	swtrtdt	TDconf
1	0	49	1	20	A	NA	0
2	0	38	0	50	B	NA	1
2	38	41	1	50	B	NA	2
3	0	227	0	40	A	NA	0
3	227	229	0	40	A	NA	1
4	0	110	0	26	A	NA	2

id	tstart	tstop	event	age	arm	swtrtdt	TDconf
1	0	49	1	20	A	NA	0
2	0	38	0	50	B	NA	1
2	38	41	1	50	B	NA	2
3	0	227	0	40	A	NA	0
3	227	229	0	40	A	NA	1
4	0	110	0	26	A	NA	2

id	age	arm	swtrtdt	tstart	tstop	event
1	20	A	NA	0	41	0
1	20	A	NA	41	49	1
3	40	A	NA	0	41	0
3	40	A	NA	41	49	0
3	40	A	NA	49	151	0
3	40	A	NA	151	153	0

id	tstart	tstop	event	age	arm	swtrtdt	TDconf
1	0	49	1	20	A	NA	0
2	0	38	0	50	B	NA	1
2	38	41	1	50	B	NA	2
3	0	227	0	40	A	NA	0
3	227	229	0	40	A	NA	1
4	0	110	0	26	A	NA	2

id	tstart	tstop	event	age	arm	swtrtdt	TDconf
1	0	49	1	20	A	NA	0
2	0	38	0	50	B	NA	1
2	38	41	1	50	B	NA	2
3	0	227	0	40	A	NA	0
3	227	229	0	40	A	NA	1
4	0	110	0	26	A	NA	2

id	age	arm	swtrtdt	tstart	tstop	event
1	20	A	NA	0	41	0
1	20	A	NA	41	49	1
3	40	A	NA	0	41	0
3	40	A	NA	41	49	0
3	40	A	NA	49	151	0
3	40	A	NA	151	153	0

id	tstart	tstop	event	age	arm	swtrtdt	TDconf
1	0	49	1	20	A	NA	0
2	0	38	0	50	B	NA	1
2	38	41	1	50	B	NA	2
3	0	227	0	40	A	NA	0
3	227	229	0	40	A	NA	1
4	0	110	0	26	A	NA	2

id	tstart	tstop	event	age	arm	swtrtdt	TDconf
1	0	49	1	20	A	NA	0
2	0	38	0	50	B	NA	1
2	38	41	1	50	B	NA	2
3	0	227	0	40	A	NA	0
3	227	229	0	40	A	NA	1
4	0	110	0	26	A	NA	2

id	age	arm	swtrtdt	tstart	tstop	event
1	20	A	NA	0	41	0
1	20	A	NA	41	49	1
3	40	A	NA	0	41	0
3	40	A	NA	41	49	0
3	40	A	NA	49	151	0
3	40	A	NA	151	153	0