## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4,
  eval = TRUE
)

## ----load-package, eval = TRUE------------------------------------------------
library(CIMEHR)

## ----custom-column-names, eval = FALSE----------------------------------------
# fit <- CIMEHR(
#   data       = my_data,
#   id_col     = "patient_id",
#   time_col   = "visit_time",
#   y_col      = "hba1c",
#   r_col      = "measured",
#   censor_col = "followup_end",
#   covars_visit_XV               = c("age", "female"),
#   covars_outcome_fixed_XY       = c("age", "female"),
#   covars_outcome_random_link_ZY = "female",
#   covars_obs_fixed_XO           = c("age", "female"),
#   covars_obs_random_link_ZO     = "female"
# )

## ----wide-to-long-------------------------------------------------------------
library(dplyr)
library(tidyr)

wide_dat <- tibble::tibble(
  id = 1:2,
  Age = c(0.4, -0.2),
  Gender = c(1, 0),
  NSES = c(0.2, -0.4),
  C = c(60, 60),
  time_1 = c(5, 8),  time_2 = c(20, NA), time_3 = c(45, 50),
  Y_1 = c(1.2, NA), Y_2 = c(1.6, NA), Y_3 = c(2.0, 1.1)
)

long_dat <- wide_dat %>%
  pivot_longer(
    cols = matches("^(time|Y)_\\d+$"),
    names_to = c(".value", "visit"),
    names_pattern = "(time|Y)_(\\d+)"
  ) %>%
  mutate(R = as.integer(!is.na(Y))) %>%
  arrange(id, time)

head(long_dat)

## ----load-data, eval = TRUE---------------------------------------------------
data("sim_ehr_data")
nrow(sim_ehr_data)
head(sim_ehr_data)

## ----subset-data, eval = TRUE-------------------------------------------------
set.seed(1)
sub_ids <- sample(unique(sim_ehr_data$id), 500)
ehr500  <- sim_ehr_data[sim_ehr_data$id %in% sub_ids, ]
nrow(ehr500)
length(unique(ehr500$id))

## ----true-params, eval = TRUE-------------------------------------------------
true_params <- attr(sim_ehr_data, "true_params")
true_params$beta
true_params$gamma

## ----covars-vec---------------------------------------------------------------
covars_all <- c("Age", "Gender", "Marital", "Black", "Hispanic", "NSES")
out_form   <- log_HbA1c ~ Age + Gender + Marital + Black + Hispanic + NSES + time

## ----summarystat--------------------------------------------------------------
fit_ss <- Summary_stat(
  ehr500,
  outcome = "log_HbA1c",
  formula = ~ Age + Gender + Marital + Black + Hispanic + NSES
)
summary(fit_ss)

## ----lmm----------------------------------------------------------------------
fit_lmm <- Linear_mixed_model(
  ehr500,
  fixed         = out_form,
  random        = ~ 1 | id,
  obs_indicator = "R"
)
summary(fit_lmm)

## ----lmm-va-------------------------------------------------------------------
fit_va <- Linear_mixed_model(
  ehr500,
  fixed         = out_form,
  random        = ~ 1 | id,
  obs_indicator = "R",
  visit_adjust  = "VA"
)
summary(fit_va)

## ----pairwise-----------------------------------------------------------------
obs_data <- ehr500[!is.na(ehr500$R) & ehr500$R == 1, ]

fit_pl <- Pairwise_likelihood(
  data        = obs_data,
  formula     = out_form,
  y_col       = "log_HbA1c",
  pair_sample = 1000
)
summary(fit_pl)

## ----iirr---------------------------------------------------------------------
fit_iirr <- Inverse_intensity_rate_ratio_weighting(
  ehr500,
  outcome    = "log_HbA1c",
  visit_covs = covars_all
)
summary(fit_iirr)

## ----iirr-bal-----------------------------------------------------------------
fit_iirr_bal <- Inverse_intensity_rate_ratio_balancing(
  ehr500,
  outcome    = "log_HbA1c",
  visit_covs = covars_all
)
summary(fit_iirr_bal)

## ----liang--------------------------------------------------------------------
fit_liang <- Joint_modeling_visiting_and_longitudinal_Liang(
  ehr500,
  outcome     = "log_HbA1c",
  visit_covs  = covars_all,
  long_covs   = covars_all,
  random_covs = "NSES"
)
summary(fit_liang)

## ----ehrjoint-----------------------------------------------------------------
fit_ej <- EHRJoint(
  ehr500,
  outcome     = "log_HbA1c",
  visit_covs  = covars_all,
  long_covs   = covars_all,
  random_covs = "NSES"
)
summary(fit_ej)

## ----cimehr-------------------------------------------------------------------
fit_cimehr <- CIMEHR(
  data                          = ehr500,
  y_col                         = "log_HbA1c",
  covars_visit_XV               = covars_all,
  covars_outcome_fixed_XY       = covars_all,
  covars_outcome_random_link_ZY = "NSES",
  covars_obs_fixed_XO           = covars_all,
  covars_obs_random_link_ZO     = "NSES"
)
print(fit_cimehr) # Print function only prinnts the outcome process estimates
summary(fit_cimehr)

## ----stage-summary-extract----------------------------------------------------
summary_visiting(fit_cimehr)
summary_observation(fit_cimehr)
summary_outcome(fit_cimehr)

extract_coefficient(fit_cimehr, stage = "outcome",     parameter = "NSES")
extract_coefficient(fit_cimehr, stage = "visiting",    parameter = "NSES")
extract_coefficient(fit_cimehr, stage = "observation", parameter = "NSES")

## ----cimehr-integral,eval = FALSE---------------------------------------------
# fit_integral <- CIMEHR_timevarying_integral(
#   data                          = ehr500,
#   y_col                         = "log_HbA1c",
#   covars_visit_XV               = covars_all,
#   covars_outcome_fixed_XY       = covars_all,
#   covars_outcome_random_link_ZY = "NSES",
#   covars_obs_fixed_XO           = covars_all,
#   covars_obs_random_link_ZO     = "NSES",
#   gh_points_U = 7,
#   gh_points_q = 5
# )

## ----cimehr-ou,eval = FALSE---------------------------------------------------
# fit_ou <- CIMEHR_timevarying_ou(
#   data                          = ehr500,
#   y_col                         = "log_HbA1c",
#   covars_visit_XV               = covars_all,
#   covars_outcome_fixed_XY       = covars_all,
#   covars_outcome_random_link_ZY = "NSES",
#   covars_obs_fixed_XO           = covars_all,
#   covars_obs_random_link_ZO     = "NSES",
#   pair_type = "adjacent"
# )

## ----method-comparisons-available---------------------------------------------
available_comparison_methods()

## ----method-comparisons-------------------------------------------------------
cmp_subset <- method_comparisons(
  data    = ehr500,
  methods = c("Summary_stat", "Linear_mixed_model",
              "Inverse_intensity_rate_ratio_weighting",
              "Joint_modeling_visiting_and_longitudinal_Liang",
              "CIMEHR"),
  method_args = list(
    Summary_stat = list(
      outcome = "log_HbA1c",
      formula = ~ Age + Gender + Marital + Black + Hispanic + NSES
    ),
    Linear_mixed_model = list(
      fixed = out_form, random = ~ 1 | id, obs_indicator = "R"
    ),
    Inverse_intensity_rate_ratio_weighting = list(
      outcome = "log_HbA1c", visit_covs = covars_all
    ),
    Joint_modeling_visiting_and_longitudinal_Liang = list(
      outcome = "log_HbA1c",
      visit_covs = covars_all, long_covs = covars_all, random_covs = "NSES"
    ),
    CIMEHR = list(
      y_col                         = "log_HbA1c",
      covars_visit_XV               = covars_all,
      covars_outcome_fixed_XY       = covars_all,
      covars_outcome_random_link_ZY = "NSES",
      covars_obs_fixed_XO           = covars_all,
      covars_obs_random_link_ZO     = "NSES"
    )
  ),
  printSE = TRUE,
  print95CI = TRUE,
  true_value_verbose = TRUE
)

print(cmp_subset)
cmp_subset$tables$outcome

## ----bootstrap, eval = FALSE--------------------------------------------------
# bs <- bootstrap(
#   fit_iirr,
#   data = ehr500,
#   B    = 200,
#   seed = 1
# )
# print(bs)
# summary(bs)

## ----bootstrap-cimehr, eval = FALSE-------------------------------------------
# bs_cimehr <- bootstrap(
#   fit_cimehr,
#   data = ehr500,
#   B    = 200,
#   seed = 1,
#   covars = list(
#     covars_visit_XV               = covars_all,
#     covars_outcome_fixed_XY       = covars_all,
#     covars_outcome_random_link_ZY = "NSES",
#     covars_obs_fixed_XO           = covars_all,
#     covars_obs_random_link_ZO     = "NSES"
#   )
# )

## ----custom-sim---------------------------------------------------------------
ehr_covs <- list(
  Age      = list(type = "continuous", dist = "uniform",
                  min = 18, max = 99, standardize = TRUE),
  Gender   = list(type = "binary", prob = 0.5),
  Marital  = list(type = "binary", prob = 0.4),
  Black    = list(type = "binary", prob = 0.30),
  Hispanic = list(type = "binary", prob = 0.20),
  NSES     = list(type = "continuous", dist = "normal",
                  mean = 0, sd = 1, standardize = FALSE)
)

dat <- sim_data_gen(
  n          = 500,
  time.start = 0,
  time.end   = 120,
  seed       = 42,
  covariates = ehr_covs,
  gamma = c(intercept = -2.8, Age = 0.5, Gender = -0.10, Marital = -0.05,
            Black = 0.15, Hispanic = -0.05, NSES = 0.3),
  sigma_zeta = 0.5,
  alpha = c(intercept = 0.2, Age = -0.3, Gender = -0.10, Marital = 0.00,
            Black = 0.05, Hispanic = -0.05, NSES = -0.6),
  obs_random_covs = "NSES",
  delta   = c(intercept = -0.3, NSES = -0.3),
  sigma_q = c(intercept =  0.5, NSES =  0.2),
  phi     = 0,
  beta = c(intercept = -2.0, Age = 0.5, Gender = 0.10,
           Marital = -0.03, Black = 0.15, Hispanic = 0.05,
           NSES = -0.5),
  beta_t              = 0.01,
  outcome_random_covs = "NSES",
  theta   = c(intercept = 0.4, NSES = 0.0),
  Sigma_b = matrix(c(0.5, 0, 0, 0.8), 2, 2),
  sigma_y = 0.7,
  verbose = FALSE
)

head(dat$long_data)