# SVall.r
# written by JuG
# July 27 2020
#' Compute weights for all time points
#' @author JuG
#' @description modifed ipwtm function to compute the weights for all time points in the dataset (from Graffeo et al DOI: 10.1002/bimj.20160223)
#' @param exposure exposure variable
#' @param numerator numerator as in ipw package
#' @param denominator denominator as in ipw package
#' @param id patient id
#' @param tstart time interval start
#' @param timevar timeinterval stop
#' @param data data.frame
#' @param trunc truncation value
#' @details
#' @references Graffeo et al DOI: 10.1002/bimj.20160223
#' @examples
#'
#'
#' @return
#' @export
SVall <- function(exposure, numerator = NULL, denominator,
id, tstart, timevar, data, trunc= NULL){
#######################################################################################
# Filename : SVall.R
#
# Description : modifed ipwtm function
# computes the weights for all time points in the dataset,
# that is event times and time changes in the treatment of interest status
#
# Required package : survival
#
# Usage : SVall(exposure, numerator, denominator,
# id, tstart, timevar, data, trunc)
# Same usage as in ipwtm (of the ipw package) with
# family = survival & type = all
#
# Value : same as the ones of ipwtm
########################################################################################
tempcall <- match.call()
tempdat <- data.frame(id = data[ , as.character(tempcall$id)],
tstart = data[ , as.character(tempcall$tstart)],
timevar = data[ , as.character(tempcall$timevar)],
exposure = data[ , as.character(tempcall$exposure)])
row.data <- nrow(tempdat)
# if numerator = NULL (unstabilized weights), then the weight equals 1
if (is.null(tempcall$numerator))
tempdat$w.numerator <- 1
# Numerator
else {
# Andersen-Gill model
mod1 <- coxph(formula = eval(parse(text = paste("Surv(",
deparse(tempcall$tstart),
",",
deparse(tempcall$timevar),
",",
deparse(tempcall$exposure),
")",
deparse(tempcall$numerator),
"+ cluster (",
deparse(tempcall$id),
")"))),
data = data)
# Relative risks
tempdat$risk.numerator <- predict(mod1, type = "risk")
# Baseline hazard
tauxBasCum <- basehaz(mod1, centered = TRUE)
# all tstart and timevar
tempo.times <- sort(unique(c(tempdat$tstart[tempdat$tstart >= 0], tempdat$timevar)))
if( ! all(tempo.times %in% tauxBasCum$time)){
tauxBasCum <- data.frame(time = tempo.times,
hazard = approx(x = tauxBasCum$time, y = tauxBasCum$hazard,
xout = tempo.times,
method = "constant", rule = 2)$y)
}
# split by patient
tabi <- split(tempdat, tempdat$id)
L.tabi <- length(tabi)
tablist <- lapply(1:L.tabi, function(i){
lignes.tabi <- nrow(tabi[[i]])
# for each [t1;t2]: \Lambda(t2) - \Lambda(t1)
bashaz.numerator <- vector(length = lignes.tabi)
# when time = t1
taux.cum.start.numerator <- approx(x = tauxBasCum$time, y = tauxBasCum$hazard,
xout = tabi[[i]][ ,"tstart"])$y
# when time = t2
taux.cum.stop.numerator <- approx(x = tauxBasCum$time, y = tauxBasCum$hazard,
xout = tabi[[i]][ ,"timevar"])$y
bashaz.numerator[1] <- taux.cum.stop.numerator[1]
if(lignes.tabi > 1){
bashaz.numerator[2:lignes.tabi] <- (taux.cum.stop.numerator - taux.cum.start.numerator)[-1]
}
indices0 <- which(tabi[[i]]$exposure == 0) # non-exposed patients
indices1 <- which(tabi[[i]]$exposure == 1) # exposed patients
tabi[[i]]$p.numerator <- vector("numeric", lignes.tabi)
# when exposition = 0 => exp( - (\Lambda(t2) - \Lambda(t1)) * RR )
tabi[[i]]$p.numerator[indices0] <- exp( - bashaz.numerator[indices0] * tabi[[i]]$risk.numerator[indices0])
# when exposition = 1 => 1 - exp ( - (\Lambda(t2) - \Lambda(t1)) * RR )
tabi[[i]]$p.numerator[indices1] <- 1 - exp( - bashaz.numerator[indices1] * tabi[[i]]$risk.numerator[indices1])
# weights obtained by the cumulative product
tabi[[i]]$w.numerator <- cumprod(tabi[[i]]$p.numerator)
return(tabi[[i]])
} )
rm(tauxBasCum)
result.numerator <- do.call(rbind, tablist)
}
# Denominator
# Andersen-Gill model
mod2 <- coxph(formula = eval(parse(text = paste("Surv(",
deparse(tempcall$tstart),
",",
deparse(tempcall$timevar),
",",
deparse(tempcall$exposure),
")",
deparse(tempcall$denominator),
"+ cluster (",
deparse(tempcall$id),
")"))),
data = data)
# Relative risks
tempdat$risk.denominator <- predict(mod2, type = "risk")
# Baseline hazard
tauxBasCum <- basehaz(mod2, centered = TRUE)
# all tstart and timevar
tempo.times <- sort( unique( c( tempdat$tstart[tempdat$tstart >= 0] , tempdat$timevar ) ) )
tempo.times <- sort( unique( c( tempdat$tstart[ tempdat$tstart >= 0 ] , tempdat$timevar ) ) )
if( ! all( tempo.times %in% tauxBasCum$time ) ){
tauxBasCum <- data.frame( time = tempo.times,
hazard = approx( x = tauxBasCum$time, y = tauxBasCum$hazard,
xout = tempo.times,
method = "constant", rule = 2 )$y )
}
# split by patient
tabi <- split(tempdat, tempdat$id)
L.tabi <- length(tabi)
tablist <- lapply(1:L.tabi, function(i){
lignes.tabi <- nrow(tabi[[i]])
# for each [t1;t2]: \Lambda(t2) - \Lambda(t1)
bashaz.denominator <- vector(length = lignes.tabi)
# when time = t1
taux.cum.start.denominator <- approx(x = tauxBasCum$time, y = tauxBasCum$hazard,
xout = tabi[[i]][ , "tstart"])$y
# when time = t2
taux.cum.stop.denominator <- approx(x = tauxBasCum$time, y = tauxBasCum$hazard,
xout = tabi[[i]][ , "timevar"])$y
if(min(tabi[[i]]$tstart ) < 0){
bashaz.denominator[1] <- tauxBasCum[1, "hazard"]
if(lignes.tabi > 1){
bashaz.denominator[2:lignes.tabi] <- (taux.cum.stop.denominator - taux.cum.start.denominator)[-1]
}
}else{
bashaz.denominator[1] <- taux.cum.stop.denominator[1]
if(lignes.tabi > 1){
bashaz.denominator[2:lignes.tabi] <- (taux.cum.stop.denominator - taux.cum.start.denominator)[-1]
}
}
indices0 <- which(tabi[[i]]$exposure == 0) # non-exposed patients
indices1 <- which(tabi[[i]]$exposure == 1) # exposed patients
tabi[[i]]$p.denominator <- vector("numeric", lignes.tabi)
# when exposition = 0 => exp( - (\Lambda(t2) - \Lambda(t1)) * RR )
tabi[[i]]$p.denominator[indices0] <- exp( - bashaz.denominator[indices0] * tabi[[i]]$risk.denominator[indices0])
# when exposition = 1 => 1 - exp ( - (\Lambda(t2) - \Lambda(t1)) * RR )
tabi[[i]]$p.denominator[indices1] <- 1 - exp( - bashaz.denominator[indices1] * tabi[[i]]$risk.denominator[indices1])
# weights obtained by the cumulative product
tabi[[i]]$w.denominator <- cumprod(tabi[[i]]$p.denominator)
return( tabi[[i]] )
} )
rm(tauxBasCum)
result.denominator <- do.call( rbind, tablist )
# Final weights
if (is.null(tempcall$numerator)){
tempdat$ipw.weights <- tempdat$w.numerator / result.denominator$w.denominator
}else{
tempdat$w.numerator <- result.numerator$w.numerator
tempdat$w.denominator <- result.denominator$w.denominator
tempdat$ipw.weights <- result.numerator$w.numerator / result.denominator$w.denominator
}
# Truncated weights (optional)
if (!(is.null(tempcall$trunc))) {
tempdat$weights.trunc <- tempdat$ipw.weights
tempdat$weights.trunc[tempdat$ipw.weights <= quantile(tempdat$ipw.weights,
0 + trunc)] <- quantile(tempdat$ipw.weights, 0 +
trunc)
tempdat$weights.trunc[tempdat$ipw.weights > quantile(tempdat$ipw.weights,
1 - trunc)] <- quantile(tempdat$ipw.weights, 1 -
trunc)
}
return(tempdat)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.