# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#                                                                         #
#               EVALUATION OF RANDOMIZED CONTROLLED TRIALS:               #
#         A PRIMER AND TUTORIAL FOR MENTAL HEALTH RESEARCHERS             #
#                                                                         #
#                                                                         #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# 0. Dataset Emulation ----------------------------------------------------     

## 0.1 Load Dependencies --------------------------------------------------

pacman::p_load(
  protectr,   # Non-public package used to retrieve trial data 
  mice,
  miceadds,
  mitml,
  purrr,
  dplyr,
  simstudy,
  sjstats,
  ggplot2
)

set.seed(123)


## 0.2 Create Data from Warehouse RCTs -------------------------------------    

db <- get.db()
prevdep_406 <- get.data("prevdep_406", merge = "intersection", glimpse = F)
prevdep_204 <- get.data("prevdep_204", merge = "intersection", glimpse = F)

L <- list(prevdep_406, prevdep_204)
names <- Reduce(intersect, lapply(L, names))
merge(prevdep_406, prevdep_204, all = TRUE) -> dat
dat <- dat[names]
names

within(dat,{
  sex = as.factor(sex)
  degree = as.factor(degree)
  prevtraining = as.factor(prevtraining)
  prevpsychoth = as.factor(prevpsychoth)
  rel = as.factor(rel)
  child = as.factor(child)
  employ = as.factor(employ)
}) %>% 
  select(id, group, sex, age, degree,
         prevtraining, prevpsychoth, sess,
         rel, child, employ, depmed,
         ceq, atsphs.0,
         cesd.0, cesd.1, cesd.2,
         badssf.0, badssf.1, badssf.2,
         hadsa.0, hadsa.1, hadsa.2) -> dat

skimr::skim(dat)



## 0.3 Simulate Demographics ----------------------------------------------

rm(def)

### 0.3.1 Continuous Variables -------------------------------------------- 

#### 0.3.1.1 Age ----------------------------------------------------------

mean(dat$age)
var(dat$age)
def <- defData(varname = "age", dist = "normal", formula = 45,
               var = 140)


### 0.3.2 Dichotomous Variables ------------------------------------------- 

#### 0.3.2.1 Sex ----------------------------------------------------------

sjstats::prop(dat, "sex == 1", na.rm = TRUE)
def <- defData(def, varname = "sex", dist = "binary", 
               formula = "0.24", link = "logit")

#### 0.3.2.2 Previous Psychotherapy ---------------------------------------

sjstats::prop(dat, "prevpsychoth == 1", na.rm = TRUE)
def <- defData(def, varname = "prevpsychoth", dist = "binary", 
               formula = "0.42", link = "identity")

#### 0.3.2.3 Previous Health Training -------------------------------------

sjstats::prop(dat, "prevtraining == 1", na.rm = TRUE)
def <- defData(def, varname = "prevtraining", dist = "binary", 
               formula = "0.22", link = "identity")

#### 0.3.2.4 Children (yes/no) -------------------------------------------

sjstats::prop(dat, "child == 1", na.rm = TRUE)
def <- defData(def, varname = "child", dist = "binary", 
               formula = "0.57", link = "identity")

#### 0.3.2.5 Employment --------------------------------------------------

sjstats::prop(dat, "employ == 1", na.rm = TRUE)
def <- defData(def, varname = "employ", dist = "binary", 
               formula = "0.85", link = "identity")



### 0.3.3 Categorical Variables -------------------------------------------- 

#### 0.3.3.1 Relationship status -------------------------------------------

genCatFormula(0.30, 0.54, 0.15, 0.01)
sjstats::prop(dat, "rel == 0", na.rm = TRUE)

def <- defData(def, varname = "rel", dist = "categorical", 
               formula = "0.3;0.54;0.15;0.01")


#### 0.3.3.2 Education ----------------------------------------------------

genCatFormula(0, 0.03, 0.32, 0.62, 0.03)
sjstats::prop(dat, "degree == 4", na.rm = TRUE)

def <- defData(def, varname = "degree", dist = "categorical", 
               formula = "0;0.03;0.32;0.62;0.03")



## 0.4 Simulate Baseline Questionnaires ------------------------------------

### 0.4.1 Depression (CES-D) -----------------------------------------------

def <- defData(def, varname = "cesd.0", dist = "normal", formula = 27,
               variance = 58)

### 0.4.2 Help-Seeking Attitudes (ATSPHS) ----------------------------------

def <- defData(def, varname = "atsphs.0", dist = "normal", formula = 14,
               variance = 6)

### 0.4.3 Credibility & Expectancy (CEQ) -----------------------------------

def <- defData(def, varname = "ceq", dist = "normal", formula = 35,
               variance = 74)

### 0.4.4 Behavioral Activation (BADS-SF) ----------------------------------

def <- defData(def, varname = "badssf.0", dist = "normal", formula = 26,
               variance = 39)

### 0.4.5 Anxiety (HADS-A) ------------------------------------------------

def <- defData(def, varname = "hadsa.0", dist = "normal", formula = 10,
               variance = 10)
rm(dd)



## 0.5 Power Calculation --------------------------------------------------

# Power for minimal important differnce
pwr::pwr.t.test(d = 0.24, sig.level = 0.05 , power = 0.8, 
                type = c("two.sample"))

# n = 273 per group 
dd <- genData(546, def)


## 0.6 Simulate Outcomes ---------------------------------------------------

dd <- trtAssign(dd, n = 2, balanced = TRUE,
                grpName = "group")


dat.cg <- dat %>% filter(group == 0)
dat.ig <- dat %>% filter(group == 1)


### 0.6.1 Depression ------------------------------------------------------

dc <- defCondition(condition = "group == 0", formula = "23",
                   variance = 87, dist = "normal")

dc <- defCondition(dc, condition = "group == 1", formula = "17",
                   variance = 74, dist = "normal")

dd <- addCondition(dc, dd, newvar = "cesd.1")

dc <- defCondition(condition = "group == 0", formula = "22",
                   variance = 94, dist = "normal")

dc <- defCondition(dc, condition = "group == 1", formula = "17",
                   variance = 84, dist = "normal")

dd <- addCondition(dc, dd, newvar = "cesd.2")


### 0.6.2 Anxiety ---------------------------------------------------------


dc <- defCondition(condition = "group == 0", formula = "9",
                   variance = 12, dist = "normal")

dc <- defCondition(dc, condition = "group == 1", formula = "7",
                   variance = 12, dist = "normal")

dd <- addCondition(dc, dd, newvar = "hadsa.1")

dc <- defCondition(condition = "group == 0", formula = "9",
                   variance = 14, dist = "normal")

dc <- defCondition(dc, condition = "group == 1", formula = "7",
                   variance = 13, dist = "normal")

dd <- addCondition(dc, dd, newvar = "hadsa.2")


### 0.6.3 Behavioral Activation -------------------------------------------

dc <- defCondition(condition = "group == 0", formula = "27",
                   variance = 46, dist = "normal")

dc <- defCondition(dc, condition = "group == 1", formula = "30",
                   variance = 52, dist = "normal")

dd <- addCondition(dc, dd, newvar = "badssf.1")

dc <- defCondition(condition = "group == 0", formula = "27",
                   variance = 50, dist = "normal")

dc <- defCondition(dc, condition = "group == 1", formula = "28",
                   variance = 63, dist = "normal")

dd <- addCondition(dc, dd, newvar = "badssf.2")



## 0.7 Data Amputation ---------------------------------------------------

### 0.7.1 MAR Mechanism --------------------------------------------------

dd$miss.cond.1 <- NULL
dd$miss.cond.2 <- NULL

amp <- ampute(dd, mech = "MAR")
amp.data <- amp$data

prop <- amp$prop
prop <- 0.3

pattern <- amp$patterns
pattern[, names(pattern) %in% c("cesd.1", "cesd.2", "hadsa.1", "hadsa.2",
                                 "badssf.1", "badssf.2")] <- 0
pattern[, !names(pattern) %in% c("cesd.1", "cesd.2", "hadsa.1", "hadsa.2",
                                 "badssf.1", "badssf.2")] <- 1

weights <- amp$weights
rownames(weights) <- colnames(weights)
weights[c("age", "sex"), c("cesd.1", "cesd.2", "hadsa.1",
                  "hadsa.2", "badssf.1", "badssf.2")] <- 2
diag(weights) <- 0


amp2.mar <- ampute(dd, mech = "MAR", prop = prop, 
                   patterns = pattern, weights = weights)

data <- amp2.mar$amp


### 0.7.2 MNAR Mechanism -------------------------------------------------

dd$miss.cond.1 <- sample(0:1, 546, replace = T)
dd$miss.cond.2 <- sample(0:4, 546, replace = T)

amp2.mnar <- ampute(dd, mech = "MAR", prop = prop)

weights <- amp2.mnar$weights
rownames(weights) <- colnames(weights)

pattern <- amp2.mnar$patterns
pattern[, names(pattern) %in% c("cesd.1", "cesd.2", "hadsa.1", "hadsa.2",
                                "badssf.1", "badssf.2")] <- 0
pattern[, !names(pattern) %in% c("cesd.1", "cesd.2", "hadsa.1", "hadsa.2",
                                 "badssf.1", "badssf.2")] <- 1
weights[c(!names(pattern) %in% 
            c("cesd.1", "cesd.2", "hadsa.1", "hadsa.2",
              "badssf.1", "badssf.2")),
        c("cesd.1", "cesd.2", "hadsa.1",
          "hadsa.2", "badssf.1", "badssf.2")] <- 0

weights["miss.cond.1", c("cesd.1", "cesd.2", "hadsa.1",
                         "hadsa.2", "badssf.1", "badssf.2")] <- 10

weights["miss.cond.2", c("cesd.1", "cesd.2", "hadsa.1",
                         "hadsa.2", "badssf.1", "badssf.2")] <- 10
diag(weights) <- 0

amp2.mnar <- ampute(dd, mech = "MAR", prop = prop, patterns = pattern)
data.mnar <- amp2.mnar$amp
amp2.mnar$weights


### 0.7.3 Correct implausible values and set classes ------------------------

data[data$cesd.1 < 0 & !is.na(data$cesd.1), ]$cesd.1 <- 0
data[data$cesd.2 < 0 & !is.na(data$cesd.2), ]$cesd.2 <- 0

data[data$hadsa.0 < 0 & !is.na(data$hadsa.0), ]$hadsa.0 <- 0
data[data$hadsa.1 < 0 & !is.na(data$hadsa.1), ]$hadsa.1 <- 0
data[data$hadsa.2 < 0 & !is.na(data$hadsa.2), ]$hadsa.2 <- 0

data[data$badssf.1 < 0 & !is.na(data$badssf.1), ]$badssf.1 <- 0
data[data$badssf.2 < 0 & !is.na(data$badssf.2), ]$badssf.2 <- 0

data[data$age < 18, ]$age <- 18

data.mnar[data.mnar$cesd.1 < 0 & !is.na(data.mnar$cesd.1), ]$cesd.1 <- 0
data.mnar[data.mnar$cesd.2 < 0 & !is.na(data.mnar$cesd.2), ]$cesd.2 <- 0

data.mnar[data.mnar$hadsa.0 < 0 & !is.na(data.mnar$hadsa.0), ]$hadsa.0 <- 0
data.mnar[data.mnar$hadsa.1 < 0 & !is.na(data.mnar$hadsa.1), ]$hadsa.1 <- 0
data.mnar[data.mnar$hadsa.2 < 0 & !is.na(data.mnar$hadsa.2), ]$hadsa.2 <- 0

data.mnar[data.mnar$badssf.1 < 0 & !is.na(data.mnar$badssf.1), ]$badssf.1 <- 0
data.mnar[data.mnar$badssf.2 < 0 & !is.na(data.mnar$badssf.2), ]$badssf.2 <- 0

data.mnar[data.mnar$age < 18, ]$age <- 18

within(data,{
  sex = as.factor(sex)
  degree = as.factor(degree)
  prevtraining = as.factor(prevtraining)
  prevpsychoth = as.factor(prevpsychoth)
  rel = as.factor(rel)
  child = as.factor(child)
  employ = as.factor(employ)
}) -> data

within(data.mnar,{
  sex = as.factor(sex)
  degree = as.factor(degree)
  prevtraining = as.factor(prevtraining)
  prevpsychoth = as.factor(prevpsychoth)
  rel = as.factor(rel)
  child = as.factor(child)
  employ = as.factor(employ)
}) -> data.mnar



## 0.8 Save data ---------------------------------------------------------------

save(data, file = "data/data.mar.rda")
save(data.mnar, file = "data/data.mnar.rda") # Used in the tutorial