#' @import methods
#' @import dplyr
#' @import tidyr
#' @include ExtendTrial.R
NULL
# Unable to find a better work-around for the expand(nesting_(select_(.
utils::globalVariables('.')
#' Trial data with utility-weighted beyond-trial extension
#'
#' An S4 class of the data from a trial with hazards and utility data supplied
#' for each outcome in the trial
#'
#' Hazard data will be used to extend the life-year measurements beyond the
#' follow-up times of the trial. Utility data will be used to weight the
#' life-years.
#'
#' @slot df_Ti a tbl_df of survival times intervals for health state regimes
#' @slot df_tj a tbl_df of health state utilities and associated time interval
#' data
#' @slot regime a formula describing the dependency of survival time intervals
#' for each health state regime
#' @slot state a formula describing the dependency of time intervals each
#' individual health state
#' @slot w a numeric vector for the utility in health state data
#'
#' @name UtilityTrial-class
#' @rdname UtilityTrial-class
#' @exportClass UtilityTrial
setClass('UtilityTrial',
contains='ExtendTrial',
slots = list(df_Ti='tbl_df',
df_tj='tbl_df',
regime='formula',
state='formula',
censor='logical',
w='numeric'),
validity = function(.Object) {
msg <- character(0)
if(!valid_response_formula(.Object@regime))
msg <- c(msg,
paste('Invalid data-specification expression for health',
'regimes.'))
if(!valid_interval_spec_formula(.Object@state))
msg <- c(msg,
paste('Invalid data-specification expression for',
'health-state time intervals.'))
if(!df_has_all_terms(.Object@df_Ti,
.Object@regime))
msg <- c(msg,
paste('Health-state regime data does not contain',
'required variables given the data specification'))
if(!df_has_all_terms(.Object@df_tj,
.Object@state))
msg <- c(msg,
paste('Health-state interval data does not contain',
'required variables given the data specification'))
if(!df_has_linear_terms(.Object@df,
.Object@regime))
msg <- c(msg,
paste('Trial data does not contain required variables',
'given the data specification for the health',
'state regime.'))
if(!df_has_linear_terms(.Object@df,
.Object@state))
msg <- c(msg,
paste('Trial data does not contain required variables',
'given the data specification for the health',
'health-state.'))
if(length(.Object@w) != nrow(.Object@df_tj))
msg <- c(msg,
paste('Lengths of utility and health-state interval data',
'must be equal.'))
# TODO:
# further checks on Ti and tj - e.g
# - formulas make sense
# - data is complete - i.e. all levels of an object are accounted
# for?
# - levels of all data make sense and are comparable
if (length(msg) == 0 &&
any(!(all.vars(.Object@regime)[-1] %in%
all.vars(.Object@state)[-(1:3)])))
msg <- c(msg,
paste('Data specification for the health state regimes',
'cannot depend upon values that are not also',
'included in health state interval data',
'specification.'))
if (length(msg) == 0 &&
any(select_(.Object@df_Ti,
.dots=all.vars(.Object@regime)[1]) %>%
unlist() < 0))
msg <- c(msg,
paste('All health state interval time values must be',
'non-negative'))
if (length(msg) == 0 &&
any(select_(.Object@df_tj,
.dots=all.vars(.Object@state)[1]) %>%
unlist() < 0))
msg <- c(msg,
paste('All health state interval time values must be',
'non-negative'))
if (length(msg) == 0)
# Comprehensive checks
msg <- c(msg,
check_Ti_tj_data(.Object))
if (length(msg) == 0)
TRUE
else
msg
}
)
#' Create a UtilityTrial
#'
#' TODO documentation
#'
#' @param df a data.frame of the trial data (long format)
#' @param df_mu a data.frame of beyond-trial hazard data
#' @param df_Ti a data.frame of health state regime survival time data
#' @param df_tj a data.frame of health state interval data
#' @param participant a formula describing the variables related to
#' participants
#' @param outcome a formula describing the variables related to outcomes
#' @param hazard a formula describing the variables of the hazard
#' @param regime a formula describing the variables of the health state regime
#' @param state a formula describing the variables of health state intervals
#' @param followup an expression for follow-up in trial data (NSE)
#' @param censor an expression for censoring status in trial data (NSE)
#' @param discount numeric (scalar) value of lifeyear discount rate
#' @param s an expression for the time variable in the hazard data (NSE)
#' @param s_max expression for extended follow-up time variable in participant
#' data (NSE)
#' @param w expression for the utility variable in the health state data (NSE)
#'
#' @name UtilityTrial
#' @rdname UtilityTrial-class
#'
#' @export
UtilityTrial <- function(df,
df_mu,
df_Ti,
df_tj,
participant,
outcome,
hazard,
regime,
state,
followup,
censor,
s,
s_max,
w,
discount) {
new('UtilityTrial',
df=select_(tbl_df(df),
.dots=as.list(
unlist(
lapply(c(terms(participant),
terms(outcome)),
labels)
)
)),
df_mu=select_(tbl_df(df_mu),
.dots=as.list(all.vars(hazard))),
df_Ti=select_(tbl_df(df_Ti),
.dots=as.list(all.vars(regime))),
df_tj=select_(tbl_df(df_tj),
.dots=as.list(all.vars(state))),
participant=participant,
outcome=outcome,
hazard=hazard,
regime=regime,
state=state,
followup=eval(substitute(followup), df, parent.frame()),
censor=(eval(substitute(censor), df, parent.frame())),
s=eval(substitute(s), df_mu, parent.frame()),
s_max=eval(substitute(s_max), df, parent.frame()),
w=eval(substitute(w), df_tj, parent.frame()),
discount=discount
)
}
#' Quality adjusted life-years from trial data
#' @param trial a trial object
#' @name dqaly
#' @rdname dqaly-methods
#' @exportMethod dqaly
setGeneric('dqaly',
function(trial) {
standardGeneric('dqaly')
}
)
prepare_tj_df <- function(df_tj, tj_vars) {
# Determine time intervals for health states
rbind(
df_tj %>%
filter_(~!T_relative) %>%
arrange_(~t_j) %>%
rename_(t.b=~t_j) %>%
group_by_(.dots=tj_vars) %>%
mutate_(t.a=~c(0, t.b[-n()])),
df_tj %>%
filter_(~T_relative) %>%
arrange_(~desc(t_j)) %>%
rename_(t.b=~t_j) %>%
group_by_(.dots=tj_vars) %>%
mutate_(t.a=~c(t.b[-1], 0))
) %>%
arrange_(~t.a) %>%
mutate_(r_j=~row_number(),
t.n=~max(ifelse(T_relative, 0, t.b))) %>%
select_(.dots=c(names(df_tj)[!(names(df_tj) %in% 't_j')],
'r_j', 't.a', 't.b', 't.n')) %>%
ungroup()
}
prepare_Ti_df <- function(df_Ti, Ti_vars) {
# Determine time intervals for regimes
df_Ti %>%
arrange_(~T_i) %>%
rename_(T.b=~T_i) %>%
group_by_(.dots=Ti_vars) %>%
mutate_(.dots=c(R_i=~row_number(),
T.a=~c(0, T.b[-n()]))) %>%
select_(.dots=c(names(df_Ti)[!(names(df_Ti) %in% 'T_i')],
'R_i', 'T.a', 'T.b')) %>%
ungroup()
}
check_Ti_tj_data <- function(trial) {
msg <- character(0)
name_lookup <- c(setNames(all.vars(trial@state)[-(1:3)],
all.vars(trial@state)[-(1:3)]),
setNames('T_i',
all.vars(trial@regime)[1]),
setNames('t_j',
all.vars(trial@state)[1]),
setNames('R_i',
all.vars(trial@state)[2]),
setNames('T_relative',
all.vars(trial@state)[3]))
Ti_vars <- all.vars(trial@regime)[-1]
tj_vars <- c('R_i', all.vars(trial@state)[-(1:3)])
df_Ti <- setNames(trial@df_Ti[,all.vars(trial@regime)],
name_lookup[all.vars(trial@regime)])
df_tj <- setNames(trial@df_tj[,all.vars(trial@state)],
name_lookup[all.vars(trial@state)])
if (length(Ti_vars) == 0) {
Ti_vars <- 'k'
tj_vars <- c(tj_vars, 'k')
df_Ti <- df_Ti %>% mutate(k=1)
df_tj <- df_tj %>% mutate(k=1)
}
# Check for duplicated interval data
if(any(df_tj %>%
group_by_(.dots=c(tj_vars, 'T_relative')) %>%
summarise_(test_duplicate=~anyDuplicated(t_j)) %>%
ungroup() %>%
select_(~test_duplicate) %>%
unlist(use.names=F)))
msg <- c(msg,
'All health state intervals must have unique break-points')
# Check for duplicated regime data
if(any(df_Ti %>%
group_by_(.dots=Ti_vars) %>%
summarise_(test_duplicate=~anyDuplicated(T_i)) %>%
ungroup() %>%
select_(~test_duplicate) %>%
unlist(use.names=F)))
msg <- c(msg,
'All health state regimes must have unique break-points')
df_Ti <- prepare_Ti_df(df_Ti, Ti_vars)
df_tj <- prepare_tj_df(df_tj, tj_vars)
# Check that there are no competing infinite spans
if (!all(
df_Ti %>%
inner_join(df_tj,
by=c(Ti_vars, 'R_i')) %>%
group_by_(.dots=tj_vars) %>%
summarise_(test_infinite=~sum(is.infinite(t.b)) <= 1) %>%
ungroup() %>%
select_(~test_infinite) %>%
unlist(use.names=F)
))
msg <- c(msg,
paste('Cannot have competing infinite span health state',
'intervals within a regime'))
# Check that all states are reachable
if (!all(
df_Ti %>%
inner_join(df_tj,
by=c(Ti_vars, 'R_i')) %>%
group_by_(.dots=tj_vars) %>%
mutate_(t.n=~max(ifelse(T_relative, 0, t.b)),
test_reach=~ifelse(T_relative,
ifelse(is.infinite(t.n),
F,
(T.b - t.a - t.n) > 0),
t.a < T.b)) %>%
ungroup() %>%
select_(~test_reach) %>%
unlist(use.names=F)
))
msg <- c(msg,
paste('All health state intervals must be reachable within',
'the survival time interval specified by the regime'))
# Check that entire regime interval is covered
if (!all(
df_Ti %>%
inner_join(df_tj,
by=c(Ti_vars, 'R_i')) %>%
group_by_(.dots=tj_vars) %>%
summarise_(t.n=~max(ifelse(T_relative, 0, t.b)),
t.m=~min(ifelse(T_relative,
ifelse(is.infinite(t.b),
-Inf,
T.b - t.b),
T.b))) %>%
mutate_(test_complete=~t.n>=t.m) %>%
ungroup() %>%
select_(~test_complete) %>%
unlist(use.names=F)
))
msg <- c(msg,
paste('Health state intervals must cover entire the entire',
'survival time interval specified by the regime'))
msg
}
# Calculate the weighted integral given piece-wise weight information
weighted_uncensored_integral <- function(w,
followup,
t.a,
t.b,
T_relative,
t.n,
d) {
if (d > 1e-14) {
ifelse(
T_relative,
(w / d) * (exp(-d * pmax(followup - t.b, t.n)) -
exp(-d * pmax(followup - t.a, 0))),
(w / d) * (exp(-d * t.a) - exp(-d * pmin(t.b, followup)))
)
} else {
ifelse(
T_relative,
pmax(followup - t.a, 0) - pmax(followup - t.b, t.n),
pmin(t.b, followup) - t.a
)
}
}
# Discounted QALY for uncensored data e.g. no beyond-trial extension.
uncensored_dqaly <- function(df, df_Ti, df_tj, d) {
# Check if df is trivially empty
if (nrow(df) == 0) {
return(c())
}
# Some helpful terms
Ti_vars <- names(df_Ti)[!(names(df_Ti) %in% c('T_i'))]
tj_vars <- names(df_tj)[!(names(df_tj) %in% c('w',
't_j',
'T_relative'))]
# Dummy join variable when regime is not dependent upon any trial variable
if (length(Ti_vars) == 0) {
Ti_vars <- 'k'
df_Ti <- df_Ti %>% mutate(k=1)
}
prepare_Ti_df(df_Ti, Ti_vars) %>%
# join onto trial data
inner_join(df %>%
mutate(id=row_number(),
k=1),
by=Ti_vars,
suffix=c('.a', '.b')) %>%
select_(~-k) %>%
filter_(.dots=c(~followup > T.a,
~followup <= T.b)) %>%
# Join onto health states within health-state regime
inner_join(prepare_tj_df(df_tj, tj_vars),
by=c(tj_vars)) %>%
# Does this free up some space? I dno.
select_(~id, ~followup, ~w, ~T_relative, ~t.a, ~t.b, ~t.n) %>%
# Discard final states that overlap initial states and any initial
# state that can't occur within follow up time
group_by_(~id) %>%
filter_(~ifelse(T_relative,
(followup - t.a) > t.n,
followup > t.a)) %>%
# calculate exact QALY
summarise_(qale=~sum(weighted_uncensored_integral(w,
followup,
t.a,
t.b,
T_relative,
t.n,
d))) %>%
# Complete any empty data
complete(id=1:nrow(df), fill=list(qale=0)) %>%
select_(~qale) %>%
unlist(use.names=F)
}
integral_survival_dt <- function(mu, ds, S) {
sum(ifelse(abs(mu) < 1E-14,
ds,
(c(1, S[-length(S)]) - S) / mu))
}
integral_decay_only <- function(d, t.a, t.b) {
if (abs(d) < 1e-14) {
pmax(t.b - t.a, 0)
} else {
(1 / d) * (exp(-d * t.a) - exp(-d * pmax(t.a, t.b)))
}
}
# Calculate width of the interval (s.a, s.b) restricted to interval (t.a, t.b)
ds_r_interval <- function(t.a, t.b, s.a, s.b) {
pmax(t.a,
pmin(t.b, s.b)) -
pmax(t.a,
pmin(t.b, s.a))
}
survival_final_only <- function(relative, mu, ds) {
if(first(relative)) {
exp(-cumsum(ifelse(ds > 0,
mu * ds,
0)))
} else {
0
}
}
# Discounted QALY for censored data, e.g. where beyond-trial extension needed.
censored_dqaly <- function(df, df_mu, df_Ti, df_tj, d) {
# Check if df is trivially empty
if (nrow(df) == 0) {
return(c())
}
# Joining variables
mu_vars <- names(df_mu)[!(names(df_mu) %in% c('s', 'mu'))]
Ti_vars <- names(df_Ti)[!(names(df_Ti) %in% c('T_i'))]
tj_vars <- names(df_tj)[!(names(df_tj) %in% c('w',
't_j',
'T_relative'))]
# Dummy join variable when regime is not dependent upon any trial variable
if (length(Ti_vars) == 0) {
Ti_vars <- 'k'
df_Ti <- df_Ti %>% mutate(k=1)
}
# New format - calculate the two integrals that each type of health state
# requires, then sum over all health states and all health regimes. Should
# be quite quick!
prepare_Ti_df(df_Ti, Ti_vars) %>%
# Join onto trial data
inner_join(df %>%
mutate(id=row_number(),
k=1),
by=Ti_vars) %>%
select_(~-k) %>%
filter_(.dots=c(~followup <= T.b)) %>%
# Minimum possible survival time given regime+followup
mutate_(t.alpha=~pmax(T.a, followup)) %>%
# Join onto health states within health-state regime
inner_join(prepare_tj_df(df_tj, tj_vars),
by=c(tj_vars)) %>%
mutate_(.dots=c(# Earliest possible entry to state (kinda; if
# it is reachable?)
t.hat.a=~ifelse(T_relative,
pmin(s_max,
pmax(t.n, # we require t.n >= 0
t.alpha - t.b)),
pmin(s_max, t.a)),
# Latest possible exit from state (kinda; if
# it is reachable?)
t.hat.b=~ifelse(T_relative,
pmin(s_max,
pmax(t.n, # we require t.n >= 0
T.b - t.a)),
pmin(s_max, t.b)))) %>%
# Filter out any initial states which are not reachable due to
# time-horizon, or any final states which are unreachable due to
# follow-up time, time-horizon, etc.
filter_(.dots=c(~t.hat.a < t.hat.b)) %>%
mutate_(
.dots=c(
# smallest possible survival time given the regime and followup
# restricted to the interval of the health state
s0.a=~pmax(t.hat.a,
pmin(t.hat.b,
t.alpha - ifelse(T_relative, t.a, 0))),
# Need to be careful here
s1.b=~ifelse(T_relative,
pmax(t.hat.a,
pmin(t.hat.b,
ifelse(is.infinite(T.b),
Inf,
T.b - t.b))),
t.hat.a)
)
) %>%
# Join onto hazard data to calculate non-constant integrand domains
inner_join(df_mu %>%
arrange_(~s) %>%
rename_(s.a=~s) %>%
group_by_(.dots=mu_vars) %>%
mutate_(s.b=~c(s.a[-1], Inf)) %>%
ungroup(),
by=mu_vars) %>%
# Filter out hazard data that is not within integral range
filter_(.dots=c(~s.b > followup,
~s.a < pmax(ifelse(T_relative,
s_max + t.b,
t.hat.b),
T.b))) %>%
mutate_(
.dots=c(
# delta-s for integral from followup to the value of s.hat.a
ds0.a=~ds_r_interval(followup,
s0.a + ifelse(T_relative, t.a, 0),
s.a,
s.b),
# delta-s for integral from followup to the end of the health
ds0.b=~ds_r_interval(followup,
t.hat.b + ifelse(T_relative, t.a, 0),
s.a,
s.b),
ds1.a=~ds_r_interval(followup, t.hat.a + t.b, s.a, s.b),
ds1.b=~ds_r_interval(followup, s1.b + t.b, s.a, s.b),
# delta-s for the survival (conditioned on followup) to the
# minimum possible survival time given the regime and followup
ds0.t.alpha=~ds_r_interval(followup, t.alpha, s.a, s.b),
# delta-s for the survival (conditioned on followup) to the
# maximum possible survival time given the regime.
ds1.T.b=~ds_r_interval(followup, T.b, s.a, s.b)
)
) %>%
group_by_(.dots=c(~id, ~R_i, ~r_j)) %>%
mutate_(
.dots=c(
# survival values for integrand
S0.t.a=~exp(-cumsum(ifelse(ds0.a > 0,
(mu + d) * ds0.a,
0))),
S0.t.b=~exp(-cumsum(ifelse(ds0.b > 0,
(mu + d) * ds0.b,
0))),
# survival values for integrand
S1.t.a=~survival_final_only(T_relative,
mu + d,
ds1.a),
S1.t.b=~survival_final_only(T_relative,
mu + d,
ds1.b)
)
) %>%
summarise_(
.dots=c(
# Calculate 'constant' survival terms that can be brought out
# the front of a simple discount (exp) integral
S0.t.alpha=~exp(-sum(ifelse(ds0.t.alpha > 0,
mu * ds0.t.alpha,
0))),
S1.T.b=~exp(-sum(ifelse(ds1.T.b > 0,
mu * ds1.T.b,
0))),
# Calculate discounted integral of non-constant survival
# integrands
IS0.t.hat.a=~first(w) * exp(-d*(first(followup) -
ifelse(first(T_relative),
first(t.a),
0))) *
integral_survival_dt(mu + d,
ds0.a,
S0.t.a),
IS0.t.hat.b=~first(w) * exp(-d*(first(followup) -
ifelse(first(T_relative),
first(t.a),
0))) *
integral_survival_dt(mu + d,
ds0.b,
S0.t.b),
# Need to deal with infinite t.b somehow? (integral should be
# zero)
IS1.t.hat.a=~first(w) * exp(-d*(first(followup) -
ifelse(first(T_relative),
first(t.b),
0))) *
integral_survival_dt(mu + d,
ds1.a,
S1.t.a),
IS1.t.hat.b=~first(w) * exp(-d*(first(followup) -
ifelse(first(T_relative),
first(t.b),
0))) *
integral_survival_dt(mu + d,
ds1.b,
S1.t.b),
s0.a=~first(s0.a),
s1.b=~first(s1.b),
t.hat.a=~first(t.hat.a),
t.hat.b=~first(t.hat.b),
w=~first(w)
)
) %>%
ungroup() %>%
mutate_(
.dots=c(
# Integral of the survival (conditional on followup) at entry
# to health state
# 1. Time-varying domain
IS0=~IS0.t.hat.b - IS0.t.hat.a,
# 2. Discount/decay-only domain with constant survival factor.
IS0.t.alpha=~w * S0.t.alpha *
integral_decay_only(d,
t.hat.a,
s0.a),
# Integral of the survival (conditional on follow up) at exit of
# the health state.
# 1.
IS1=~ifelse(is.infinite(IS1.t.hat.b) &
is.infinite(IS1.t.hat.a),
0,
IS1.t.hat.b - IS1.t.hat.a),
# 2. Discount/decay-only domain with constant survival factor.
IS1.T.b=~w * S1.T.b *
integral_decay_only(d,
s1.b,
t.hat.b)
)
) %>%
group_by_(~id) %>%
# Calculate life-expectancy from sum of integrals
summarise_(qale=~sum((IS0 + IS0.t.alpha) - (IS1 + IS1.T.b))) %>%
complete(id=1:nrow(df), fill=list(qale=0)) %>%
select_(~qale) %>%
unlist(use.names=F)
}
#' @rdname dqaly-methods
#' @aliases dqaly,UtilityTrial-method
setMethod('dqaly',
c(trial = 'UtilityTrial'),
function(trial) {
name_lookup <- c(standard_covariate_names(trial),
setNames('mu', all.vars(trial@hazard)[1]),
setNames('T_i', all.vars(trial@regime)[1]),
setNames('t_j', all.vars(trial@state)[1]),
setNames('R_i', all.vars(trial@state)[2]),
setNames('T_relative', all.vars(trial@state)[3]))
# Calculate no-extension and extension separately then reorder
c(uncensored_dqaly(
setNames(trial@df[!trial@censor,],
name_lookup[c(all.vars(trial@participant),
all.vars(trial@outcome))]) %>%
mutate(followup=trial@followup[!trial@censor]),
setNames(trial@df_Ti[,all.vars(trial@regime)],
name_lookup[all.vars(trial@regime)]),
setNames(trial@df_tj[,all.vars(trial@state)],
name_lookup[all.vars(trial@state)]) %>%
mutate(w=trial@w),
trial@discount
),
censored_dqaly(
setNames(trial@df[trial@censor,],
name_lookup[c(all.vars(trial@participant),
all.vars(trial@outcome))]) %>%
mutate(followup=trial@followup[trial@censor],
s_max=trial@s_max[trial@censor]),
setNames(trial@df_mu[,all.vars(trial@hazard)],
name_lookup[all.vars(trial@hazard)]) %>%
mutate(s=trial@s),
setNames(trial@df_Ti[,all.vars(trial@regime)],
name_lookup[all.vars(trial@regime)]),
setNames(trial@df_tj[,all.vars(trial@state)],
name_lookup[all.vars(trial@state)]) %>%
mutate(w=trial@w),
trial@discount
))[order(c(which(!trial@censor),
which(trial@censor)))]
}
)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.