# Load packages
library(nlme)
library(scdhlm)

# Load data
data(Bryant2018)
dat <- Bryant2018


# Clean data

dat <- dat[,c("school", "case", "treatment", "session", "outcome")]
names(dat) <- c("cluster","case","phase","session","outcome")

dat <- preprocess_SCD(design = "CMB", 
                      case = case, 
                      phase = phase, 
                      session = session, 
                      outcome = outcome, 
                      cluster = cluster,
                      center = 0,
                      data = dat)

# Fit the model
fit_RML <- lme(fixed = outcome ~ 1 + trt, 
               random = list(cluster = ~ 1, case = ~ 1), 
               correlation = corAR1(0.01, ~ session | cluster / case), 
               data = dat,
               control = lmeControl(msMaxIter = 50, apVar = FALSE, returnObject = TRUE))
summary(fit_RML)

# Calculate effect size with g_mlm()
A <- 0
B <- 1
p_const <- c(rep(0L, length(c(0))), (B - A)^as.integer(c(0))) 
p_const <- c(0,1)
r_const <- c(c(1), c(), 
             c(1), c(), 
             c(0), c(), 1L) # specify whether using random effects, cor struct, var struct, and level-1 errors
r_const <- c(1,1,0,1)

ES_RML <- g_mlm(fit_RML, p_const = p_const, r_const = r_const, infotype = "expected")