Initial Value For Fitting the General Markov Model

Description

Compute the initial value for fitting the MSM illness-death general Markov model using EM type algorithm

Usage

OUT_em_weights(data,X1,X2,event1,event2,w,Trt)

Arguments

Details

See usual_illness_death_weight

Value

A list of vectors and dataframes:

See Also

usual_illness_death_weight

Examples

n <- 500
set.seed(1234)
Cens = runif(n,0.7,0.9)
set.seed(1234)
OUT1 <- sim_cox_msm_semicmrsk(beta1 = 1,beta2 = 1,beta3 = 0.5,
                              sigma_2 = 1,
                              alpha0 = 0.5, alpha1 = 0.1, alpha2 = -0.1, alpha3 = -0.2,
                              n=n, Cens = Cens)
data_test <- OUT1$data0

## Get the PS weights
vars <- c("Z1","Z2","Z3")
ps1 <- doPS(data = data_test,
            Trt = "A",
            Trt.name = 1,
            VARS. = vars,
            logistic = TRUE,w= NULL)
w <- ps1$Data$ipw_ate_stab

### Get the initial value
EM_initial <- OUT_em_weights(data = data_test,
                             X1 = "X1",
                             X2 = "X2",
                             event1 = "delta1",
                             event2 = "delta2",
                             w = w,
                             Trt = "A")

Obtaining Bayesian Bootstrap Sample for Individual Risk Difference and Risk Ratio.

Description

bayesian_boot_irrd provides the bootstrap sample for individual risk difference and risk ratio, it can be used for further inferences.

Usage

bayesian_boot_irrd(dat2,B,sigma_2_0, EM_initial, varlist, t1_star,t)

Arguments

Details

For each bootstrap sample:

1. Generate n standard exponential (mean and variance 1) random variates : u_1, u_2,..., u_n;

2. The weights for the Bayesian bootstrap are: w_{i}^{boot} = u_i / \bar{u}, where \bar{u} = n^{-1}\sum_{i=1}^{n} u_i;

3. Calculate the propensity score and IP weights w_{i}^{IPW} based on Bayesian bootstrap weighted data, and assigned the weights for fitting the MSM general Markov model as w_i = w_{i}^{boot} * w_{i}^{IPW}.

4. After obtaining \hat{\theta} and \hat{b}_i, for each individual i, calculate the IRR and IRD by plugging \hat{\theta}, \hat{b}_i and a=0, a=1 separately at time t.

The 95% prediction intervals (PI) cam be obtained by the normal approximation using bootstrap standard error.

Value

Estimating Three Cumulative Incidence Functions Using the Usual Markov Model

Description

cif_est_usual estimates the cumulative incidence function (CIF, i.e.risk) based on the MSM illness-death usual Markov model.

Usage

cif_est_usual(data,X1,X2,event1,event2,w,Trt,
              t1_star = t1_star)

Arguments

Details

After estimating the parameters in the illness-death model \lambda_{j}^a using IPW, we could estimate the corresponding CIF:

\hat{P}(T_1^a<t,\delta_1^a=1) = \int_{0}^{t} \hat{S}^a(u) d\hat{\Lambda}_{1}^a(u),

\hat{P}(T_2^a<t,\delta_1^a=0,\delta_2^a=1) = \int_{0}^{t} \hat{S}^a(u) d\hat{\Lambda}_{2}^a(u),

and

\hat{P}(T_2^a<t_2 \mid T_1^a<t_1, T_2^a>t_1) = 1- e^{- \int_{t_1}^{t_2} d \hat{\Lambda}_{12}^a(u) },

where \hat{S}^a is the estimated overall survial function for joint T_1^a, T_2^a, \hat{S}^a(u) = e^{-\hat{\Lambda}_{1}^a(u)} - \hat{\Lambda}_{2}^a(u) . We obtain three hazards by fitting the MSM illness-death model \hat\Lambda_{j}^a(u) = \hat\Lambda_{0j}(u)e^{\hat\beta_j*a} , \hat\Lambda_{12}^a(u) = \hat\Lambda_{03}(u)e^{\hat\beta_3*a} , and \hat\Lambda_{0j}(u) is a Breslow-type estimator of the baseline cumulative hazard.

Value

Returns a table containing the estimated CIF for the event of interest for control and treated group.

References

Meira-Machado, Luis and Sestelo, Marta (2019). “Estimation in the progressive illness-death model: A nonexhaustive review,” Biometrical Journal 61(2), 245–263.

Estimating Three Conditional Cumulative Incidence Functions Using the General Markov Model Conditional on Random Effect

Description

conditional_cif_b estimates the cumulative incidence function based on the MSM illness-death general Markov model conditional on the fixed random effect b.

Usage

conditional_cif_b(res1,
                  t1_star,
                  b)

Arguments

Details

Similar as cif_est_usual, after estimating the parameters in the illness-death model \lambda_{j}^a using IPW, we could estimate the corresponding conditional CIF under fixed b:

\hat{P}(T_1^a<t,\delta_1^a=1 \mid b) = \int_{0}^{t} \hat{S}^a(u \mid b) d\hat{\Lambda}_{1}^a(u \mid b ),

\hat{P}(T_2^a<t,\delta_1^a=0,\delta_2^a=1 \mid b) = \int_{0}^{t} \hat{S}^a(u \mid b) d\hat{\Lambda}_{2}^a(u \mid b),

and

\hat{P}(T_2^a<t_2 \mid T_1^a<t_1, T_2^a>t_1 \mid b) = 1- e^{- \int_{t_1}^{t_2} d \hat{\Lambda}_{12}^a(u \mid b) },

where \hat{S}^a is the estimated overall survial function for joint T_1^a, T_2^a, \hat{S}^a(u) = e^{-\hat{\Lambda}_{1}^a(u)} - \hat{\Lambda}_{2}^a(u) . We obtain three hazards by fitting the MSM illness-death model \hat\Lambda_{j}^a(u) = \hat\Lambda_{0j}(u)e^{\hat\beta_j*a} , \hat\Lambda_{12}^a(u) = \hat\Lambda_{03}(u)e^{\hat\beta_3*a} , and \hat\Lambda_{0j}(u) is a Breslow-type estimator of the baseline cumulative hazard.

where S(t \mid b;a) = \exp[- \int_0^{t} \{ \lambda_{01} (u)e^{\beta_1a + b} + \lambda_{02} (u )e^{\beta_2a + b} \} d u ] = \exp \{- e^{\beta_1a + b} \Lambda_{01}(t) - e^{\beta_2a + b} \Lambda_{02} (t ) \}

Value

See Also

cif_est_usual

Generate the Inverse Probability Treatment Weights

Description

doPS calculates the unstabilized and stabilized inverse probability treatment weights (IPW) for average treatment effect using propensity score. The propensity score is calculated by twang package using the boosted logistic regression.

Usage

doPS(data,Trt,Trt.name,VARS.,logistic = FALSE,w=NULL)

Arguments

Details

The treatment variable should only contain 2 levels of treatment, and one should be viewed as treated group and another is control group.

For stabilized weights:

For the treated individuals, we assign the IPW: w = Pr(T=1)/Pr(T=1|X=x), for control individuals, the stabilized weight is: w = (1-Pr(T=1))/(1-Pr(T=1|X=x)).

Value

doPS returns an object of class "PS". An object of class "PS" is a list containing the following components:

See Also

ps

Examples

n <- 500
set.seed(1234)
Cens = runif(n,0.7,0.9)
set.seed(1234)
OUT1 <- sim_cox_msm_semicmrsk(beta1 = 1,beta2 = 1,beta3 = 0.5,
                              sigma_2 = 1,
                              alpha0 = 0.5, alpha1 = 0.1, alpha2 = -0.1, alpha3 = -0.2,
                              n=n, Cens = Cens)
data_test <- OUT1$data0

## Get the PS weights
vars <- c("Z1","Z2","Z3")
ps1 <- doPS(data = data_test,
            Trt = "A",
            Trt.name = 1,
            VARS. = vars,
            logistic = TRUE,w=NULL)
w <- ps1$Data$ipw_ate_stab

Using EM Type Algorithm for MSM Illness-death General Markov Model

Description

Under the general Markov illness-death model, with normal frailty term which is a latent variable. We use the EM type algorithm to estimate the coefficient in the MSM illness-death general Markov model.

Usage

em_illness_death_phmm_weight(data,X1,X2,event1,event2,w,Trt,
                            EM_initial,sigma_2_0)

Arguments

Details

Similar as the usual Markov model. We postulate the semi-parametric Cox models with a frailty term for three transition rates in marginal structural illness-death model:

\lambda_{1}(t_1 ; a) = \lambda_{01}(t)e^{\beta_1 a + b}, t_1 > 0 ;

\lambda_{2}(t_2 ; a) = \lambda_{02}(t)e^{\beta_2 a + b}, t_2 > 0 ;

and

\lambda_{12}(t_2 \mid t_1 ; a) = \lambda_{03}(t_2)e^{\beta_3 a + b}, 0 < t_1 < t_2 ,

where b \sim N(0,1). Since b is not observed in the data, we use the IP weighted EM type algorithm to estimate all the parameters in the MSM illness-death general Markov model.

Value

Examples

n <- 500
set.seed(1234)
Cens = runif(n,0.7,0.9)
set.seed(1234)
OUT1 <- sim_cox_msm_semicmrsk(beta1 = 1,beta2 = 1,beta3 = 0.5,
                              sigma_2 = 1,
                              alpha0 = 0.5, alpha1 = 0.1, alpha2 = -0.1, alpha3 = -0.2,
                              n=n, Cens = Cens)
data_test <- OUT1$data0

## Get the PS weights
vars <- c("Z1","Z2","Z3")
ps1 <- doPS(data = data_test,
            Trt = "A",
            Trt.name = 1,
            VARS. = vars,
            logistic = TRUE,w=NULL)
w <- ps1$Data$ipw_ate_stab

### Fit the General Markov model
EM_initial <- OUT_em_weights(data = data_test,
                             X1 = "X1",
                             X2 = "X2",
                             event1 = "delta1",
                             event2 = "delta2",
                             w = w,
                             Trt = "A")

res1 <- em_illness_death_phmm_weight(data = data_test,
                                     X1 = "X1",
                                     X2 = "X2",
                                     event1 = "delta1",
                                     event2 = "delta2",
                                     w = w,
                                     Trt = "A",
                                     EM_initial = EM_initial,
                                     sigma_2_0 = 2)

print(paste("The estimated value for beta1 is:", round(res1$beta1[res1$em.n],5) ) )

Compute the (Cumulative) Baseline Hazard from Cox Model

Description

Compute the Breslow type baseline hazard and cumulative baseline hazard at each event time from a Cox model.

Usage

get_hazard(fit)

Arguments

Details

See also basehaz, we only extract the estimated baseline hazard and baseline cumulative hazard from the results of a coxph fit.

Value

A list contains two dataframes.

See Also

basehaz

Compute the (Cumulative) Baseline Hazard from Cox Model with Offsets

Description

Compute the Breslow type baseline hazard and cumulative baseline hazard at each event time from a weighted Cox model with offsets.

Usage

get_hazard_offset_weights(fit,data,time1= NULL,time2,w)

Arguments

Details

See also get_hazard, handles the offset term in coxph for predicting the baseline hazard.

Value

A list contains two dataframes.

See Also

get_hazard, basehaz

Estimating Three Individual Risk Difference and Risk Ratio Using the General Markov Model Conditional on Predicted Random Effect

Description

individual_RR_RD estimates the individual risk difference and risk ratio based on the MSM illness-death general Markov model conditional on predicted random effect for each data point at a fixed time point.

Usage

individual_RR_RD(dat1,res1,t1_star ,t)

Arguments

Details

Similar as cif_est_usual, after estimating the parameters in the illness-death model \lambda_{j}^a using IPW, we could estimate the corresponding conditional CIF under the predicted b:

\hat{P}(T_1^a<t,\delta_1^a=1 \mid b) = \int_{0}^{t} \hat{S}^a(u \mid b) d\hat{\Lambda}_{1}^a(u \mid b ),

\hat{P}(T_2^a<t,\delta_1^a=0,\delta_2^a=1 \mid b) = \int_{0}^{t} \hat{S}^a(u \mid b) d\hat{\Lambda}_{2}^a(u \mid b),

and

\hat{P}(T_2^a<t_2 \mid T_1^a<t_1, T_2^a>t_1 \mid b) = 1- e^{- \int_{t_1}^{t_2} d \hat{\Lambda}_{12}^a(u \mid b) },

The frailty term, or equivalently, the random effect b represents the unobserved heterogeneity among the individuals. As such, the above conditional risk represents individual risk, and the risk contrasts the individual risk contrasts. We therefore have the individual risk difference (IRD) and the individual risk ratio (IRR).

Under the random effects model, for i = 1,2,...,n, the predicted random effect is \hat{b}_i = E(b_i \mid O_i, \hat{\theta}). We then obtain the predicted IRD and the predicted IRR.

Value

Returns a data frame that includes the individual risk difference / ratio for three type of events.

Fit the MSM Cox Model with IP Weights

Description

Fit the MSM cox model with IPW as the initial value for EM algorithm to fit the illness-death general Markov model

Usage

initial_fit_em_weights(data,X1,X2,event1,event2,w,Trt)

Arguments

Details

As initial values we use for \beta_j, j=1, 2, 3, the estimates from IP weighted Cox regression without the offsets, i.e. from the usual Markov model.

Value

A list of objects from survival package:

See Also

Surv, coxph

Compute the Initial (Cumulative) Baseline Hazard From the MSM Illness-death Model

Description

Compute the Breslow type baseline hazard and cumulative baseline hazard at each event time from the MSM illness-death model.

Usage

initial_lambda_em (OUT)

Arguments

Details

See also get_hazard

Value

A list contains six dataframes: including baseline hazard and cumulative baseline hazard for non-terminal event, terminal event without non-terminal event, and terminal event following non-terminal event.

See Also

get_hazard

Plotting Histogram of Propensity Score and Balancing Plot for Covariates in the Propensity Score Model

Description

Displays a the histogram plots for the propensity score, stratified by treated and control group and a graph of standardized mean difference of potential confounders before and after weigthing.

Usage

## S3 method for class 'PS'
plot(x,...)

Arguments

Details

Only available when logistic = FALSE in doPS. The standardized mean difference (SMD), defined as the (weighted) treatment group mean minus the (weighted) control group mean divided by the (weighted) pooled sample (treatment and control) standard deviation. SMD between -0.1 and 0.1 typically indicates good balance.

Value

Histogram of propensity score and balancing plot for covariates in the propensity score model corresponding to the output from doPS.

See Also

bal.table

Simulating Semi-competing Risks with Right-censored Survival Data under Marginal Structural Illness-death Cox Model

Description

The function to simulate semi-competing risk with right-censored survival data under marginal structural illness-death Cox model.

Usage

sim_cox_msm_semicmrsk(beta1,beta2,beta3,sigma_2,
        alpha0,alpha1,alpha2,alpha3,
        n,Cens)

Arguments

Details

We simulate data followed by Xu(2010) to generate semi-competing risk data under illness-death model, where we have baseline hazard \lambda_{01}(t) = \lambda_{02}(t) = 2exp(-t)I(0 \le t \le 3) + 2exp(-3)I(t \ge 3), and \lambda_{03}(t) = 2\lambda_{01}(t).

We also have the propensity score model to generate treatment assignment P_A = logit^{-1}(\alpha_0 + \alpha_1 Z_1 + \alpha_2 Z_2 + \alpha_3 Z3).

Value

Returns a data frame that contains time to non-terminal event, T1, terminal event, T2 and censoring time C with their event indicator, delta1 and delta2. Three covariates Z1, Z2, Z3, and treatment assignment A are also included.

Examples

n <- 500
set.seed(1234)
Cens = runif(n,0.7,0.9)
set.seed(1234)
OUT1 <- sim_cox_msm_semicmrsk(beta1 = 1,beta2 = 1,beta3 = 0.5,
                              sigma_2 = 1,
                              alpha0 = 0.5, alpha1 = 0.1, alpha2 = -0.1, alpha3 = -0.2,
                              n=n, Cens = Cens)
data_test <- OUT1$data0

Fit MSM Illness-death Usual Markov Model For Semi-competing Risks Data

Description

Fit the marginal structural three-state illness-death model with Cox representation and IP weights for semi-competing risks data. Inference under this model can be carried out using estimating equations with IP weights.

Usage

usual_illness_death_weight(data,X1,X2,event1,event2,w,Trt)

Arguments

Details

Let T_1, T_2 be the time to non-terminal event and terminal event, A be the treatment assignment. We postulate the semi-parametric Cox models for three transition rates in marginal structural illness-death model:

\lambda_{1}(t_1 ; a) = \lambda_{01}(t)e^{\beta_1 a}, t_1 > 0 ;

\lambda_{2}(t_2 ; a) = \lambda_{02}(t)e^{\beta_2 a}, t_2 > 0 ;

and

\lambda_{12}(t_2 \mid t_1 ; a) = \lambda_{03}(t_2)e^{\beta_3 a}, 0 < t_1 < t_2 .

The coefficients as well as Breslow type baseline hazards can be estimated by fitting the IP weights Cox proportional hazards models. Meanwhile, if we assume the estimated weights as known, then the robust sandwich variance estimator can be used to obtain the estimated variance.

The usual Markov model is also the same as the initial value for the general Markov model.

Value

A list of values and dataframes:

Examples

n <- 500
set.seed(1234)
Cens = runif(n,0.7,0.9)
set.seed(1234)
OUT1 <- sim_cox_msm_semicmrsk(beta1 = 1,beta2 = 1,beta3 = 0.5,
                              sigma_2 = 1,
                              alpha0 = 0.5, alpha1 = 0.1, alpha2 = -0.1, alpha3 = -0.2,
                              n=n, Cens = Cens)
data_test <- OUT1$data0

## Get the PS weights
vars <- c("Z1","Z2","Z3")
ps1 <- doPS(data = data_test,
            Trt = "A",
            Trt.name = 1,
            VARS. = vars,
            logistic = TRUE,w=NULL)
w <- ps1$Data$ipw_ate_stab

### Fit the Usual Markov model
res1 <- usual_illness_death_weight(data = data_test,
                                   X1 = "X1",
                                   X2 = "X2",
                                   event1 = "delta1",
                                   event2 = "delta2",
                                   w = w,
                                   Trt = "A")
print(paste("The estimated value for beta1 is:", round(res1$beta1,5) ) )

Variance of parameters in MSM Illness-death General Markov Model

Description

Use bootstrap to obtain the variance estimator for parameters in MSM illness-death general markov model.

Usage

var_em_illness_death_phmm(data,sigma_2_0,VARS.)

Arguments

Details

See em_illness_death_phmm_weight. In each bootstrap, the propensity score model needs to be re-fitted, and fit the MSM illness-death general markov model with new IP weights.

Value

List of bootstrap SE for all the parameters in the general Markov model