ExampleWorkflow.Rmd

---
title: "Recommended Workflow for using devMSMs with Longitudinal Data"
author: "Isabella C. Stallworthy", "Meriah L. DeJoseph", "Emily R. Padrutt", "Noah Greifer", "Daniel Berry"
date: "`r Sys.Date()`"
output: html_document
---
This example workflow provides more detail on the application of *devMSMs* to the empirical example in the manuscript, xxxx, which can be found here: xx.

Please review the accompanying manuscript for a full conceptual and practical introduction to MSMs in the context of developmental data. Please also see the vignettes on the *devMSMs* website for step-by-step guidance on the use of this code: https://istallworthy.github.io/devMSMs/index.html. 

We highly recommend first implementing the *Data Requirements* and *Specify Core Inputs* Vignettes prior to using this workflow. 

Headings denote accompanying website vignettes and steps that correspond to the manuscript Section 7. We suggest using the interactive outline tool (located above the Console) for ease of navigation.

The code in each code chunk is set up to showcase all possible inputs to each function (both required and optional) to aid the user's use of the full range of package functionality. Example possible values for the optional input are shown for each function.

Please see the website vignettes and/or type `?functionName` into the console for more guidance on the arguments for each function. 

Functions from devMSMs: https://github.com/istallworthy/devMSMs


# *Installation*
https://istallworthy.github.io/devMSMs/index.html 

Until *devMSMs* is available on CRAN, you will need to install it directly from Github (https://github.com/istallworthy/devMSMs), as shown below.  

```{r}
install.packages("devtools", quiet = TRUE)
library(devtools)

install_github("istallworthy/devMSMs", quiet = TRUE)
library(devMSMs)

install_github("istallworthy/devMSMsHelpers", quiet = TRUE)
library(devMSMsHelpers)
```
<br>

# *Specify Core Inputs* and *Data Requirements & Inspection* Vignettes
## Phase 0: Preliminary Conceptual & Data Preparation
Please see the accompanying manuscript for steps P1 (creating hypotheses) and P2 (creating a DAG).

### STEP P3. Specify Core Inputs
The first step is to specify core inputs for use with the package. 

In addition to the core inputs supplied directly to `initMSM()` in Step P4.4, we specify the following (See the Specify Core Inputs vignette for more information including guidance on specifying other core inputs (P3.1-P3.8)) : 

#### P3.4. Required for Continuous Exposures: Identify High and Low Cutoff Values 
See the Specify Core Inputs vignette for more information. 

Below, we specify the 60th and 30th percentiles to demarcate high and low levels of economic strain exposure, respectively.  
```{r}
hi_lo_cut <- c(0.6, 0.3)
```

#### P3.6. Recommended: Specify Hypotheses-Relevant Exposure Histories 
See the Specify Core Inputs vignette for more information. 

Below, we specify all reference histories and all comparison histories associated with our two hypotheses. We will specify more specific reference and comparison histories corresponding to each hypothesis in Step 5.   
```{r}
reference <- c("h-l-l", "l-h-h", "l-l-l", "l-h-l", "l-l-h", "h-h-h")

comparison <- c("l-l-h", "l-h-l", "l-l-l", "h-l-h", "h-h-l", "h-l-l", "l-h-h")
```

#### P3.7 Required: Specify Outcome
```{r}
outcome <-  "StrDif_Tot.58"
```
<br>

### STEP P4. Data Preparation & Inspection
We highly recommend first implementing the *Data Preparation Vignette* https://istallworthy.github.io/devMSMs/articles/Preliminary_Steps.html for preparing data starting in different formats (P4.1 - P4.3a).

#### P4.3b. Required: Read in Wide Data  
Please see the *Data Requirements* vignette for guidance on preparing complete data in wide format 
  
*devMSMs* requires data in one of the following wide data formats (see Figure 1) for use in the package:  

* a single data frame of data in wide format with no missing data  

* a mids object (output from mice::mice()) of data imputed in wide format  

* a list of data imputed in wide format as data frames  
<br>

We first read in data that has been formatted and multiply imputed 5 times using the *mice* package using the random forest imputation technique. 
```{r}
data <- readRDS("/Users/isabella/Library/CloudStorage/Box-Box/BSL General/MSMs/testing/isa/imputations/ESETA1-StrDif_Tot.58_all_imp.rds") 
```
<br>

#### P4.4 Required: Create MSM Object
Below, we create an MSM object by supplying our imputed data, exposure, time invariant and time-varying confounders, exposure epochs, the separator preceding the time point for time-varying variables, and an optional home directory for saving output. 
```{r}
set.seed(1234)

obj <- initMSM(
  data,
  exposure = c("ESETA1.6", "ESETA1.15", "ESETA1.24", "ESETA1.35", "ESETA1.58"), 
  ti_conf =  c("state", "BioDadInHH2", "PmAge2", "PmBlac2", "TcBlac2", "PmMrSt2", 
               "PmEd2", "KFASTScr", "RMomAgeU", "RHealth", "HomeOwnd", "SWghtLB", 
               "SurpPreg", "SmokTotl", "DrnkFreq",
               "peri_health", "caregiv_health", "gov_assist"),
  tv_conf = c("SAAmylase.6","SAAmylase.15", "SAAmylase.24", 
              "MDI.6", "MDI.15",                                            
              "RHasSO.6", "RHasSO.15", "RHasSO.24","RHasSO.35",                                         
              "WndNbrhood.6","WndNbrhood.24", "WndNbrhood.35",                                      
              "IBRAttn.6", "IBRAttn.15", "IBRAttn.24",                                   
              "B18Raw.6", "B18Raw.15", "B18Raw.24",                                           
              "HOMEETA1.6", "HOMEETA1.15", "HOMEETA1.24", "HOMEETA1.35",                              
              "InRatioCor.6", "InRatioCor.15", "InRatioCor.24", "InRatioCor.35",                         
              "CORTB.6", "CORTB.15", "CORTB.24",                                                                  
              "EARS_TJo.24", "EARS_TJo.35",                                        
              "LESMnPos.24", "LESMnPos.35",                                  
              "LESMnNeg.24", "LESMnNeg.35",       
              "StrDif_Tot.35", 
              "fscore.35"
              ),
  epoch = c("Infancy", "Infancy", "Toddlerhood", "Toddlerhood", "Childhood"),
  sep = "\\.",
  home_dir = '/Users/isabella/Library/CloudStorage/Box-Box/BSL General/MSMs/testing/isa' 
)
```

View the msm object.
```{r}
print(obj)
```


#### P4.5 Recommended: Inspect Inspect Exposure Histories and Data
The helper `inspectData()` function outputs the following files into the home directory: a correlation plot of all variables in the dataset, tables of exposure and outcome descriptive statistics, and two summary tables of the confounders considered at each time point.

```{r, eval = FALSE}
inspectData(data = data, 
            obj = obj, 
            outcome = outcome, 
            hi_lo_cut = hi_lo_cut,
            reference = reference, 
            comparison = comparison, 
            verbose = TRUE,
            save.out = TRUE)
```
Here, we see summaries of the data types as well as reasonable cell counts in each of our specified histories, for each imputed dataset. 
<br>
<br>

# *Workflow: Continuous Exposure* Vignette
https://istallworthy.github.io/devMSMs/articles/Workflow_Continuous_Exposure.html

## PHASE 1: Confounder Adjustment
The first phase of the MSM process is focused on eliminating confounding of the relation between exposure and outcome.  

### STEP 1: Create Full Balancing Formulas & Conduct Pre-Balance Checking
The first step is to create full balancing formulas that reflect all measured confounders at each exposure time point.  

#### 1a. Create Full Balancing Formulas at each Exposure Time Point
Users have the option to specify concurrent confounders to retain and we recommend doing so consistently throughout this workflow. 

Please see the Customize Balancing Formulas Vignette at the link below for more detail on custom formulas. 
https://istallworthy.github.io/devMSMs/articles/Customize_Balancing_Formulas.html  

```{r}
type <- "full"

full_formulas <- createFormulas(obj = obj, 
                                type = type, 
                                save.out = TRUE)
```

Inspect
```{r}
print(full_formulas)
```


#### 1b. Conduct Exploratory Pre-Balance Assessment
The next step is to examine initial imbalance between confounders and exposure prior to IPTW weighting. Users have the option to specify balance threshold(s), which we recommend doing consistently throughout this workflow.  

```{r}
balance_thresh <- c(0.05, 0.1) 

imp_conf <- c("InRatioCor.6", "InRatioCor.15", "InRatioCor.24", "InRatioCor.35",
              "PmEd2") 

prebalance_stats <- assessBalance(obj = obj, 
                                  data = data,
                                  balance_thresh = balance_thresh, 
                                  imp_conf = imp_conf, 
                                  save.out = TRUE)
```

Inspect and save out 
```{r}
# summarize prebalance stats averaging across imputed datasets 
summary(prebalance_stats, 
        save.out = TRUE)

# view and save out balance stats averaged across imputed datasets
print(prebalance_stats,
      save.out = TRUE)

# save out plots for all time points and imputed datasets 
for (i in 1:5){
  for (t in 1:5){
    plot(prebalance_stats, 
         i = i, 
         t = t, 
         save.out = TRUE)
  }
}

# save out balance stats for all imputed datasets
for (i in 1:5){
    print(prebalance_stats, 
          i = i, 
          save.out = TRUE) 
}
```


### STEP 2: Create Simplified Balancing Formulas & Determine Optimal Weighting Method
The next step is to specify shorter, simplified balancing formula for the purposes of determining the weighting method optimal for the data.   

#### 2a. Create Simplified Balancing Formulas
First, create shorter, simplified balancing formulas at each exposure time point.  

```{r}
type <- "short" 

short_formulas <- createFormulas(obj = obj, 
                                 type = type, 
                                 save.out = TRUE) 
```

Inspect
```{r}
print(short_formulas)
```


#### 2b. Create IPTW Balancing Weights Using Multiple Weighting Methods
We recommend users use the short formulas to create IPTW weights using all the available weighting methods: "glm", "gbm", "bart", "super", and "cbps." 

```{r}
formulas <- short_formulas

method <- "cbps" 

weights.cbps <- createWeights(obj = obj, 
                              data = data, 
                              formulas = formulas, 
                              method = method, 
                              # maxit = 1, # for testing purposes only; makes it run faster 
                              save.out = TRUE)
# 
# weights.cbps <- readRDS("/Users/isabella/Library/CloudStorage/Box-Box/BSL General/MSMs/testing/isa/weights/type_short-exposure_ESETA1-method_cbps.rds")
```

Inspect by imputed dataset
```{r}
print(weights.cbps, 
      i = 1)

plot(weights.cbps, 
     i = 1, 
     save.out = TRUE)

summary(weights.cbps[[1]]) 
```

We then fit weights for all other available weighting methods. 
```{r}
method <- "glm"
weights.glm <- createWeights(obj = obj, 
                             data = data, 
                             formulas = formulas, 
                             method = method,
                             save.out = TRUE)
print(weights.glm, 
      i = 2)
plot(weights.glm, 
     i = 1,
     save.out = TRUE)


method <- "gbm"
weights.gbm <- createWeights(obj = obj, 
                             data = data, 
                             formulas = formulas, 
                             method = method,
                             save.out = TRUE)
print(weights.gbm, 
      i = 2)
plot(weights.gbm, 
     i = 1, 
     save.out = TRUE)


method <- "bart"
weights.bart <- createWeights(obj = obj, 
                              data = data, 
                              formulas = formulas, 
                              method = method,
                              save.out = TRUE)
print(weights.bart, 
      i = 2)
plot(weights.bart, 
     i = 1, 
     save.out = TRUE)



method <- "super"
weights.super <- createWeights(obj = obj, 
                               data = data, 
                               formulas = formulas, 
                               method = method,
                               save.out = TRUE)
print(weights.super, 
      i = 2)
plot(weights.super, 
     i = 1, 
     save.out = TRUE)
```


#### 2c. Assess All Weighting Methods to Determine Optimal Method
The next step is to assess balance for each weighting method and for the user to determine the optimal weighting method.  

```{r}
balance_thresh <- c(0.05, 0.1) 

imp_conf <- c("InRatioCor.6", "InRatioCor.15", "InRatioCor.24", 
              "InRatioCor.35", "PmEd2") 

weights <- weights.cbps 

balance_stats.cbps <- assessBalance(data = data, 
                                    obj = obj, 
                                    weights = weights,
                                    imp_conf = imp_conf, 
                                    balance_thresh = balance_thresh,
                                    save.out = TRUE)

# to load in from file (output from saving above!)
# balance_stats.cbps <- readRDS("/Users/isabella/Library/CloudStorage/Box-Box/BSL General/MSMs/testing/isa/balance/weighted-exposure_A-method_cbps.rds")

weights <- weights.glm
balance_stats.glm <- assessBalance(data = data, 
                                   obj = obj, 
                                   weights = weights,
                                   imp_conf = imp_conf, 
                                   balance_thresh = balance_thresh,
                                   save.out = TRUE)

weights <- weights.gbm
balance_stats.gbm <- assessBalance(data = data, 
                                   obj = obj, 
                                   weights = weights,
                                   imp_conf = imp_conf, 
                                   balance_thresh = balance_thresh,
                                   save.out = TRUE)


weights <- weights.bart
balance_stats.bart <- assessBalance(data = data, 
                                    obj = obj, 
                                    weights = weights,
                                    imp_conf = imp_conf, 
                                    balance_thresh = balance_thresh,
                                    save.out = TRUE)


weights <- weights.super 
balance_stats.super <- assessBalance(data = data, 
                                     obj = obj, 
                                     weights = weights,
                                     imp_conf = imp_conf, 
                                     balance_thresh = balance_thresh,
                                     save.out = TRUE)
```

We then compare the balance of the different weighting methods. 

```{r}
summary(balance_stats.cbps, 
        save.out = TRUE)
print(balance_stats.cbps, 
      save.out = TRUE)

summary(balance_stats.glm, 
        save.out = TRUE)
print(balance_stats.glm, 
      save.out = TRUE)

summary(balance_stats.gbm, 
        save.out = TRUE)
print(balance_stats.gbm, 
      save.out = TRUE)

summary(balance_stats.bart, 
        save.out = TRUE)
print(balance_stats.bart, 
      save.out = TRUE)

summary(balance_stats.super, 
        save.out = TRUE)
print(balance_stats.super, 
      save.out = TRUE)
```



We identify CBPS as the weighting method that most successfully reduces imbalance. 



### STEP 3: Create Updated Formulas & Re-Specify Weights Using Optimal Weighting Method
The next step is to update the short balancing formulas with any imbalanced confounders, and re-specify the IPTW weights.  

#### 3a. Inspect Balance of Best-Performing Weighting Method
First, assess how well each of the IPTW achieve balance for all measured confounders. 

```{r}
# save all time points averaged across imps
for (t in 1:5){
  plot(balance_stats.glm, 
       t = t, 
       save.out = TRUE)
}

# save for all imps 
for (t in 1:5){
  for (i in 1:5){
  plot(balance_stats.glm, 
       t = t, 
       i = i, 
       save.out = TRUE)
  }
}


# save all time points averaged across imps
for (t in 1:5){
  plot(balance_stats.gbm, 
       t = t, 
       save.out = TRUE)
}

# save for all imps 
for (t in 1:5){
  for (i in 1:5){
  plot(balance_stats.gbm, 
       t = t, 
       i = i, 
       save.out = TRUE)
  }
}


# save all time points averaged across imps
for (t in 1:5){
  plot(balance_stats.super, 
       t = t, 
       save.out = TRUE)
}

# save for all imps 
for (t in 1:5){
  for (i in 1:5){
  plot(balance_stats.super, 
       t = t, 
       i = i, 
       save.out = TRUE)
  }
}

# save all time points averaged across imps
for (t in 1:5){
  plot(balance_stats.bart, 
       t = t, 
       save.out = TRUE)
}

# save for all imps 
for (t in 1:5){
  for (i in 1:5){
  plot(balance_stats.bart, 
       t = t, 
       i = i, 
       save.out = TRUE)
  }
}




summary(balance_stats.cbps)

plot(balance_stats.cbps, 
     t = 1, 
     save.out = FALSE)

# save all time points averaged across imps
for (t in 1:5){
  plot(balance_stats.cbps, 
       t = t, 
       save.out = TRUE)
}

# save for all imps 
for (t in 1:5){
  for (i in 1:5){
  plot(balance_stats.cbps, 
       t = t, 
       i = i, 
       save.out = TRUE)
  }
}

print(balance_stats.cbps) 
```


#### 3b. Update Simplified Formulas 
Next, update the short formulas with any imbalanced confounders.  

```{r}
type <- "update"

bal_stats <- balance_stats.cbps

updated_formulas <- createFormulas(obj = obj, 
                                   type = type,
                                   bal_stats = bal_stats,
                                   save.out = TRUE)
```

Inspect
```{r}
print(updated_formulas)
```


#### 3c. Create Final Balancing Weights
Next, create final balancing weights using the optimal weighting method and updated balancing formulas.  

```{r}
formulas <- updated_formulas

method <- "cbps" 

final_weights <- createWeights(data = data, 
                               obj = obj, 
                               method = method, 
                               formulas = formulas,
                               # maxit = 1, # just for testing purposes 
                               save.out = TRUE)
# final_weights <- readRDS("/Users/isabella/Library/CloudStorage/Box-Box/BSL General/MSMs/testing/isa/weights/type_update-exposure_ESETA1-method_cbps.rds")
```

Inspect
```{r}
print(final_weights, 
      i = 1)

summary(final_weights[[3]])

plot(final_weights, 
     i = 1,
     save.out = TRUE)
```


#### 3d. Trim Final Balancing Weights
Next, trim the final balancing weights to reduce the heavy right tail.  

##### Main
First, trim the main weights.  

```{r}
quantile <- 0.95 

weights <- final_weights
trim_weights <- trimWeights(weights = weights, 
                            at = quantile, 
                            save.out = TRUE)
```

Inspect
```{r}
print(trim_weights, 
      i = 1)

plot(trim_weights, 
     i = 1, 
     save.out = TRUE)
```


##### Sensitvity Analyses
Next, conduct sensitivity analyses using two other quantile values.  

```{r}
quantile <- 0.92 
trim_weights.s1 <- trimWeights(weights = weights, 
                               at = quantile,
                               save.out = TRUE)
plot(trim_weights.s1, 
     i = 1, 
     save.out = TRUE)


quantile <- 0.98
trim_weights.s2 <- trimWeights(weights = weights, 
                               at = quantile,
                               save.out = TRUE)
plot(trim_weights.s2, 
     i = 1, 
     save.out = TRUE)
```



### STEP 4: Conduct Final Balance Assessment
Next, conduct a final balance assessment using all measured confounders (i.e., the full balancing formulas).  

#### Main
First, conduct the main balance assessment.  

```{r}
balance_thresh <- c(0.05, 0.1) 

imp_conf <- c("InRatioCor.6", "InRatioCor.15", "InRatioCor.24", 
              "InRatioCor.35", "PmEd2") 

weights <- trim_weights
final_balance_stats <- assessBalance(data = data, 
                                     obj = obj, 
                                     balance_thresh = balance_thresh,
                                     imp_conf = imp_conf,
                                     weights = weights,
                                     save.out = TRUE)
```

Inspect
```{r}
summary(final_balance_stats, 
        save.out = TRUE)

plot(final_balance_stats, 
     t = 2, 
     save.out = TRUE)

print(final_balance_stats, 
      save.out = TRUE)
```


#### Sensitvity Analyses
Next, conduct the recommended specifying sensitivity analyses to match the main analyses above.  

```{r}
weights <- trim_weights.s1
final_balance_stats.s1 <- assessBalance(data = data, 
                                        obj = obj,
                                        weights = weights,
                                        save.out = "balance-weighted-exposure_ESETA1-method_cbps-form_update-s1.rds")
summary(final_balance_stats.s1, 
        save.out = TRUE)


weights <- trim_weights.s2
final_balance_stats.s2 <- assessBalance(data = data, 
                                        obj = obj, 
                                        weights = weights,
                                        save.out = "balance-weighted-exposure_ESETA1-method_cbp-form_update-s2.rds")
summary(final_balance_stats.s2, 
        save.out = TRUE)
```



## PHASE 2: Assess Substantive Associations between Exposure and Outcome
Lastly, having attenuated confounder associations, we model substantive associations.  

### STEP 5: Fit Weighted Model & Summarize & Visualize Results

#### 5a. Select and Fit a Weighted Outcome Model
First, select and fit a marginal outcome model.  

##### Main
First, fit the main model.  

```{r}
family <- gaussian

link <- "identity" 

weights <- trim_weights

m <- "m0"

outcome <- "StrDif_Tot.58"

models <- fitModel(data = data, 
                   obj = obj, 
                   weights = weights, 
                   outcome = outcome, 
                   model = m,
                   save.out = TRUE) 

models_unweighted <- fitModel(data = data,
                   obj = obj,
                   weights = NULL,
                   outcome = outcome,
                   model = m,
                   save.out = TRUE)
```

Inspect
```{r}
print(models, 
      save.out = TRUE)

print(models, 
      i = 3,
      save.out = TRUE)


# # unweighted
# print(models_unweighted, 
#       save.out = "fit_model_summary-outcome_StrDif_Tot_58-exposure_ESETA1-model_m0-all_avg_UNWEIGHTED.docx")
```



##### Sensitvity Analyses
Next, fit the recommended  specifying sensitivity analyses to match the main analyses above.   

```{r}
weights <- trim_weights.s1
models.s1 <- fitModel(data = data, 
                      obj = obj, 
                      weights = weights, 
                      outcome = outcome, 
                      model = m,
                      save.out = sprintf("outcome_StrDif_Tot-exposure_ESETA1-model_%s_s1.rds", m))
print(models.s1, 
      save.out = sprintf("fit_model_summary-outcome_StrDif_Tot_58-exposure_ESETA1-model_%s-all_avg_s1.docx", m))


weights <- trim_weights.s2
models.s2 <- fitModel(data = data, 
                      obj = obj, 
                      weights = weights, 
                      outcome = outcome, 
                      model = m,
                      save.out = sprintf("outcome_StrDif_Tot-exposure_ESETA1-model_%s_s2.rds", m))
print(models.s2,
      save.out = sprintf("fit_model_summary-outcome_StrDif_Tot_58-exposure_ESETA1-model_%s-all_avg_s2.docx", m))
```


#### 5b. Estimate, Compare, & Visualize Model-Predicted Outcome as a Function of History
Lastly, estimate and compare user-specified exposure histories.   

```{r}
hi_lo_cut <- c(0.6, 0.3) 

mc_comp_method <- "BH" 

dose_level <- "h"

exp_lab <- "Economic Strain" 

out_lab <- "Behavior Problems"

colors <- c("Dark2") 

model <- models 
```


##### Timing hypotheses
Assess timing hypotheses about the role of exposure during infancy.  

##### 1 dose timing
```{r}
reference_t1 <- "h-l-l" 
comparison_t1 <- c("l-l-h", "l-h-l", "l-l-l")

results_t1 <- compareHistories(fit = model, 
                            comparison = comparison_t1,
                            reference = reference_t1,
                            hi_lo_cut = hi_lo_cut,
                            mc_comp_method = mc_comp_method, 
                            dose = dose_level,
                            save.out = TRUE)
summary(results_t1)

print(results_t1, 
      save.out = "comparisons_table-outcome_StrDif_Tot_58-exposure_ESETA1-t1.docx")
plot(results_t1,
     save.out = "comparisons_plot-outcome_StrDif_Tot_58-exposure_ESETA1-t1.jpeg")


# 
# unweighted
results_t1_uw <- compareHistories(fit = models_unweighted,
                            comparison = comparison_t1,
                            reference = reference_t1,
                            hi_lo_cut = hi_lo_cut,
                            mc_comp_method = mc_comp_method,
                            dose = dose_level,
                            save.out = TRUE)
summary(results_t1_uw)

print(results_t1,
      save.out = "comparisons_table-outcome_StrDif_Tot_58-exposure_ESETA1-t1-uw.docx")
plot(results_t1,
     save.out = "comparisons_plot-outcome_StrDif_Tot_58-exposure_ESETA1-t1-uw.jpeg")
```

##### 2 dose timing
```{r}
reference_t2 <- "l-h-h" 
comparison_t2 <- c("h-l-h", "h-h-l") 

results_t2 <- compareHistories(fit = model, 
                            comparison = comparison_t2,
                            reference = reference_t2,
                            hi_lo_cut = hi_lo_cut,
                            mc_comp_method = mc_comp_method, 
                            dose = dose_level,
                            save.out = TRUE)
summary(results_t2)

print(results_t2, 
      save.out = "comparisons_table-outcome_StrDif_Tot_58-exposure_ESETA1-t2.docx")
plot(results_t2,
     save.out = "comparisons_plot-outcome_StrDif_Tot_58-exposure_ESETA1-t2.jpeg")



# unweighted
results_t2_uw <- compareHistories(fit = models_unweighted,
                            comparison = comparison_t2,
                            reference = reference_t2,
                            hi_lo_cut = hi_lo_cut,
                            mc_comp_method = mc_comp_method,
                            dose = dose_level,
                            save.out = FALSE)
summary(results_t2_uw)

print(results_t2_uw,
      save.out = "comparisons_table-outcome_StrDif_Tot_58-exposure_ESETA1-t2-uw.docx")
plot(results_t2_uw,
     save.out = "comparisons_plot-outcome_StrDif_Tot_58-exposure_ESETA1-t2-uw.jpeg")
```


##### Dose hypotheses
Assess dose hypotheses.

##### 0 vs 1
```{r}
reference_d0 <- "l-l-l"
comparison_d0 <- c("h-l-l", "l-h-l", "l-l-h")

results_d0 <- compareHistories(fit = model, 
                            comparison = comparison_d0,
                            reference = reference_d0,
                            hi_lo_cut = hi_lo_cut,
                            mc_comp_method = mc_comp_method, 
                            dose = dose_level,
                            save.out = TRUE)
summary(results_d0)

print(results_d0, 
      save.out = "comparisons_table-outcome_StrDif_Tot_58-exposure_ESETA1-d0.docx")
plot(results_d0,
     save.out = "comparisons_plot-outcome_StrDif_Tot_58-exposure_ESETA1-d0.jpeg")



# unweighted
results_d0_uw <- compareHistories(fit = models_unweighted,
                            comparison = comparison_d0,
                            reference = reference_d0,
                            hi_lo_cut = hi_lo_cut,
                            mc_comp_method = mc_comp_method,
                            dose = dose_level,
                            save.out = TRUE)
summary(results_d0_uw)

print(results_d0_uw,
      save.out = "comparisons_table-outcome_StrDif_Tot_58-exposure_ESETA1-d0-uw.docx")
plot(results_d0_uw,
     save.out = "comparisons_plot-outcome_StrDif_Tot_58-exposure_ESETA1-d0-uw.jpeg")
```

##### 1 vs 2
```{r}
reference_d1 <- c("h-l-l", "l-h-l", "l-l-h")
comparison_d1 <- c("h-h-l", "l-h-h", "h-l-h")

results_d1 <- compareHistories(fit = model, 
                            comparison = comparison_d1,
                            reference = reference_d1,
                            hi_lo_cut = hi_lo_cut,
                            mc_comp_method = mc_comp_method, 
                            dose = dose_level,
                            save.out = TRUE)
summary(results_d1)

print(results_d1, 
      save.out = "comparisons_table-outcome_StrDif_Tot_58-exposure_ESETA1-d1.docx")
plot(results_d1,
     save.out = "comparisons_plot-outcome_StrDif_Tot_58-exposure_ESETA1-d1.jpeg")



# unweighted
results_d1_uw <- compareHistories(fit = models_unweighted,
                            comparison = comparison_d1,
                            reference = reference_d1,
                            hi_lo_cut = hi_lo_cut,
                            mc_comp_method = mc_comp_method,
                            dose = dose_level,
                            save.out = TRUE)
summary(results_d1_uw)

print(results_d1_uw,
      save.out = "comparisons_table-outcome_StrDif_Tot_58-exposure_ESETA1-d1-uw.docx")
plot(results_d1_uw,
     save.out = "comparisons_plot-outcome_StrDif_Tot_58-exposure_ESETA1-d1-uw.jpeg")
```

##### 2 vs 3
```{r}
reference_d2 <- "h-h-h"
comparison_d2 <- c("h-h-l", "l-h-h", "h-l-h")

results_d2 <- compareHistories(fit = model, 
                            comparison = comparison_d2,
                            reference = reference_d2,
                            hi_lo_cut = hi_lo_cut,
                            mc_comp_method = mc_comp_method, 
                            dose = dose_level,
                            save.out = TRUE)
summary(results_d2)

print(results_d2, 
      save.out = "comparisons_table-outcome_StrDif_Tot_58-exposure_ESETA1-d2.docx")
plot(results_d2,
     save.out = "comparisons_plot-outcome_StrDif_Tot_58-exposure_ESETA1-d2.jpeg")



# unweighted
results_d2_uw <- compareHistories(fit = models_unweighted,
                            comparison = comparison_d2,
                            reference = reference_d2,
                            hi_lo_cut = hi_lo_cut,
                            mc_comp_method = mc_comp_method,
                            dose = dose_level,
                            save.out = TRUE)
summary(results_d2_uw)

print(results_d2_uw,
      save.out = "comparisons_table-outcome_StrDif_Tot_58-exposure_ESETA1-d2-uw.docx")
plot(results_d2_uw,
     save.out = "comparisons_plot-outcome_StrDif_Tot_58-exposure_ESETA1-d2-uw.jpeg")
```

##### 3 vs 1
```{r}
reference_d3 <- "h-h-h"
comparison_d3 <- c("h-l-l", "l-h-l", "l-l-h")

results_d3 <- compareHistories(fit = model, 
                            comparison = comparison_d3,
                            reference = reference_d3,
                            hi_lo_cut = hi_lo_cut,
                            mc_comp_method = mc_comp_method, 
                            dose = dose_level,
                            save.out = TRUE)

summary(results_d3)

print(results_d3, 
      save.out = "comparisons_table-outcome_StrDif_Tot_58-exposure_ESETA1-d3.docx")
plot(results_d3,
     colors = c("#7570B3", "#D95F02", "#D95F02", "#D95F02"),
     save.out = "comparisons_plot-outcome_StrDif_Tot_58-exposure_ESETA1-d3.jpeg")


# sensiivity analyses
results_d3.s1 <- compareHistories(fit = models.s1, 
                            comparison = comparison_d3,
                            reference = reference_d3,
                            hi_lo_cut = hi_lo_cut,
                            mc_comp_method = mc_comp_method, 
                            dose = dose_level,
                            save.out = TRUE)
summary(results_d3.s1)

print(results_d3.s1, 
      save.out = "comparisons_table-outcome_StrDif_Tot_58-exposure_ESETA1-d3.s1.docx")
plot(results_d3.s1,
     colors = c("#7570B3", "#D95F02", "#D95F02", "#D95F02"),
     save.out = "comparisons_plot-outcome_StrDif_Tot_58-exposure_ESETA1-d3.s1.jpeg")


results_d3.s2 <- compareHistories(fit = models.s2, 
                            comparison = comparison_d3,
                            reference = reference_d3,
                            hi_lo_cut = hi_lo_cut,
                            mc_comp_method = mc_comp_method, 
                            dose = dose_level,
                            save.out = TRUE)
summary(results_d3.s2)

print(results_d3.s2, 
      save.out = "comparisons_table-outcome_StrDif_Tot_58-exposure_ESETA1-d3.s2.docx")
plot(results_d3.s2,
     colors = c("#7570B3", "#D95F02", "#D95F02", "#D95F02"),
     save.out = "comparisons_plot-outcome_StrDif_Tot_58-exposure_ESETA1-d3.s2.jpeg")



# unweighted
results_d3_uw <- compareHistories(fit = models_unweighted,
                            comparison = comparison_d3,
                            reference = reference_d3,
                            hi_lo_cut = hi_lo_cut,
                            mc_comp_method = mc_comp_method,
                            dose = dose_level,
                            save.out = FALSE)

summary(results_d3_uw)

print(results_d3_uw,
      save.out = "comparisons_table-outcome_StrDif_Tot_58-exposure_ESETA1-d3-uw.docx")
plot(results_d3_uw,
     colors = c("#7570B3", "#D95F02", "#D95F02", "#D95F02"),
     save.out = "comparisons_plot-outcome_StrDif_Tot_58-exposure_ESETA1-d3-uw.jpeg")
```

##### 3 vs 0
```{r}
reference_d4 <- "h-h-h"
comparison_d4 <- "l-l-l"

results_d4 <- compareHistories(fit = model, 
                            comparison = comparison_d4,
                            reference = reference_d4,
                            hi_lo_cut = hi_lo_cut,
                            mc_comp_method = mc_comp_method, 
                            dose = dose_level,
                            save.out = TRUE)
summary(results_d4)

print(results_d4, 
      save.out = "comparisons_table-outcome_StrDif_Tot_58-exposure_ESETA1-d4.docx")
plot(results_d4,
     colors = c("#1B9E77", "#E7298A","#66A61E", "#E6AB02") ,#"#A6761D" "#666666"
     save.out = "comparisons_plot-outcome_StrDif_Tot_58-exposure_ESETA1-d4.jpeg")


# sensitivity analyses
results_d4.s1 <- compareHistories(fit = models.s1, 
                            comparison = comparison_d4,
                            reference = reference_d4,
                            hi_lo_cut = hi_lo_cut,
                            mc_comp_method = mc_comp_method, 
                            dose = dose_level,
                            save.out = TRUE)
summary(results_d4.s1)

print(results_d4.s1, 
      save.out = "comparisons_table-outcome_StrDif_Tot_58-exposure_ESETA1-d4.s1.docx")
plot(results_d4.s1,
     save.out = "comparisons_plot-outcome_StrDif_Tot_58-exposure_ESETA1-d4.s1.jpeg")


results_d4.s2 <- compareHistories(fit = models.s2, 
                            comparison = comparison_d4,
                            reference = reference_d4,
                            hi_lo_cut = hi_lo_cut,
                            mc_comp_method = mc_comp_method, 
                            dose = dose_level,
                            save.out = TRUE)
summary(results_d4.s2)

print(results_d4.s2, 
      save.out = "comparisons_table-outcome_StrDif_Tot_58-exposure_ESETA1-d4.s2.docx")
plot(results_d4.s2,
     save.out = "comparisons_plot-outcome_StrDif_Tot_58-exposure_ESETA1-d4.s2.jpeg")



# unweighted
results_d4_uw <- compareHistories(fit = models_unweighted,
                            comparison = comparison_d4,
                            reference = reference_d4,
                            hi_lo_cut = hi_lo_cut,
                            mc_comp_method = mc_comp_method,
                            dose = dose_level,
                            save.out = TRUE)
summary(results_d4_uw)

print(results_d4_uw,
      save.out = "comparisons_table-outcome_StrDif_Tot_58-exposure_ESETA1-d4-uw.docx")
plot(results_d4,
     save.out = "comparisons_plot-outcome_StrDif_Tot_58-exposure_ESETA1-d4-uw.jpeg")
```


##### 2 vs 0
```{r}
reference_d5 <- "l-l-l"
comparison_d5 <- c("h-h-l", "l-h-h", "h-l-h")

results_d5 <- compareHistories(fit = model, 
                            comparison = comparison_d5,
                            reference = reference_d5,
                            hi_lo_cut = hi_lo_cut,
                            mc_comp_method = mc_comp_method, 
                            dose = dose_level,
                            save.out = TRUE)
summary(results_d5)

print(results_d5, 
      save.out = "comparisons_table-outcome_StrDif_Tot_58-exposure_ESETA1-d5.docx")
plot(results_d5,
     save.out = "comparisons_plot-outcome_StrDif_Tot_58-exposure_ESETA1-d5.jpeg")



# unweighted
results_d5_uw <- compareHistories(fit = models_unweighted,
                            comparison = comparison_d5,
                            reference = reference_d5,
                            hi_lo_cut = hi_lo_cut,
                            mc_comp_method = mc_comp_method,
                            dose = dose_level,
                            save.out = TRUE)
summary(results_d5_uw)

print(results_d5_uw,
      save.out = "comparisons_table-outcome_StrDif_Tot_58-exposure_ESETA1-d5-uw.docx")
plot(results_d5_uw,
     save.out = "comparisons_plot-outcome_StrDif_Tot_58-exposure_ESETA1-d5-uw.jpeg")
```






##### More stringent p-value correction considering all comparisons as one family. 
```{r}
all_results <- rbind(results_t1$comps,
                     results_t2$comps,
                     results_d0$comps,
                     results_d1$comps,
                     results_d2$comps,
                     results_d3$comps,
                     results_d4$comps,
                     results_d5$comps)
all_results$p_corr_all <- p.adjust(all_results$p.value, method = mc_comp_method)
all_results[] <- lapply(all_results, function(x) if(is.numeric(x)) round(x, 3) else x)

tinytable::tt(all_results) |> 
  tinytable::save_tt(file.path(obj[["home_dir"]], "histories/all_comps.html"), 
                               overwrite = TRUE)


# unweighted
all_results_uw <- rbind(results_t1_uw$comps,
                     results_t2_uw$comps,
                     results_d0_uw$comps,
                     results_d1_uw$comps,
                     results_d2_uw$comps,
                     results_d3_uw$comps,
                     results_d4_uw$comps,
                     results_d5_uw$comps)
all_results_uw$p_corr_all <- p.adjust(all_results_uw$p.value, method = mc_comp_method)
all_results_uw[] <- lapply(all_results_uw, function(x) if(is.numeric(x)) round(x, 3) else x)

tinytable::tt(all_results_uw) |>
  tinytable::save_tt(file.path(obj[["home_dir"]], "/histories/all_comps_uw.html"),
                               overwrite = TRUE)
```


# Package Citations
We are grateful to the authors of many existing packages that *devMSMs* draws from!

```{r}
grateful::cite_packages(out.dir = obj[["home_dir"]], 
                        omit = c("devMSMs", "devMSMsHelpers"), 
                        out.format = "docx")
```