R/UtilityTrial.R

#' @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)))]
    }
)
stephematician/trialcostR documentation built on May 30, 2019, 3:18 p.m.