16 Koko koodi
Alla on näkyvissä tässä harjoituksessa esitellyt koodit yhtenä putkena, jolloin ne kaikki voidaan ajaa kerralla.
# Load packages
library(pacman)
p_load(dplyr, dlookr, ipcwswitch, ipw, ggplot2, survival, survminer)
# Display toydata
knitr::kable(ipcwswitch::toydata)
# Esimerkkikoodin laajennus ad n=400
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 ))
# Yritetään antaa kaikille samat arvot eri aikapisteiden ps-kovariaattiin ihan testiksi:
toydata$ps1 <- 1
toydata$ps2 <- 1
toydata$ps3 <- 1
# ipswswitch nikottelee tibblejen kanssa,
# joten on muutettava datasettimme data.frame-muotoon:
toydata <- as.data.frame(toydata)
# Silmäillään nyt uudelleen toydataa
knitr::kable(head(toydata))
# Näytetään toydatan struktuuri
str(toydata)
# Survival-analyysin valmistelu: ensin *timesTokeep*-komennolla käsittely
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]]
# Seuraavaksi *long format* muotoon muuttaminen
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]])
# Näytetään 10 ensimmäistä riviä
knitr::kable(head(toy.long))
# Päivämäärät numeeriseen muotoon, jolloin tstart-muuttuja pisteessä 0
# 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
# 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 sensurointia ja sen jälkeen
knitr::kable(head(toy.long))
knitr::kable(head(toy.long2))
# Kaikkien event-aikojen yhdistäminen yhteen ja samaan dataan
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"])
# 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:
knitr::kable(head(toy.rep))
# 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')
# Tarkastellaan taas dataa
knitr::kable(head(toy.rep))
# Painotusten distribuution plottaaminen
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
# Truncated stabilized weights:
fit.stab.w <- coxph(Surv(tstart, tstop, event) ~ age + arm + cluster(id),
data = toy.rep, weights = toy.rep$weights.trunc)
# Summary mallista
summary(fit.stab.w)
# Survival-käyrän plottaaminen - tutkimuksen koko populaatio
fit_all <- survfit(fit.stab.w)
ggsurvplot(fit_all, data=toydata, color = "#2E9FDF",
ggtheme = theme_minimal())
# Survival-käyrän plottaaminen erikseen arm A:sta ja B:stä
# 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
fit2 <- survfit(fit.stab.w, newdata = apudf)
ggsurvplot(fit2,
data=toydata,
conf.int = TRUE,
legend.labs=c("Arm = A", "Arm = B"),
ggtheme = theme_minimal())