## ----function-signature, eval=FALSE-------------------------------------------
# MRStdCRT_fit(
#   formula, data, cluster, trt, trtprob = rep(0.5, nrow(data)),
#   method, family = gaussian(link = "identity"), corstr, scale,
#   alpha = 0.05,
#   model_weights = c("none", "inverse_cluster_size", "estimand_specific"),
#   inference_scale = c("original", "log")
# )

## ----transformed-covariate-formula, eval=FALSE--------------------------------
# data$x_sq <- data$x^2
# formula <- y ~ x + x_sq

## ----treatment-probabilities--------------------------------------------------
library(MRStdCRT)
library(dplyr)
data(ppact)

prob <- rep(0.5, nrow(ppact))

## ----estimated-treatment-probability------------------------------------------
cluster_assignments <- ppact %>%
  distinct(CLUST, INTERVENTION)

estimated_prob <- rep(
  mean(cluster_assignments$INTERVENTION == 1),
  nrow(ppact)
)

## ----hypothetical-block-probabilities-----------------------------------------
clusters <- sort(unique(ppact$CLUST))
n_clusters <- length(clusters) 
block_info <- data.frame(
  CLUST = clusters
) %>%
  mutate(
    # Create a 'block' identifier for each cluster
    block = case_when(
      row_number() <= 25 ~ 1,
      row_number() <= 66 ~ 2,
      TRUE ~ 3
    )
  ) %>%
  mutate(
    # Assign the desired hypothetical treatment probability to each block
    hypothetical_prob = case_when(
      block == 1 ~ 0.3,
      block == 2 ~ 0.5,
      block == 3 ~ 0.6
    )
  )
prob_lookup <- setNames(block_info$hypothetical_prob, block_info$CLUST)
probs_vector <- prob_lookup[as.character(ppact$CLUST)]

## ----fit-ppact----------------------------------------------------------------
example <- MRStdCRT_fit(
  formula = PEGS ~ AGE + FEMALE + comorbid + Dep_OR_Anx + pain_count + PEGS_bl +
    BL_benzo_flag + BL_avg_daily + satisfied_primary + n,
  data = ppact,
  cluster = "CLUST",
  trt = "INTERVENTION",
  trtprob = prob,
  method = "GEE",
  corstr = "independence",
  scale = "RR"
)
## To view the summary, use the following command
summary(example)
## One may also extract specific statistical feature, such as table for point and interval estimates, and standard errors
example$estimate
example$ics

## ----weighted-example-data----------------------------------------------------
set.seed(68213)
sizes <- rep(c(18, 25, 33), 4)
m <- length(sizes)
trial <- data.frame(
  id = rep(seq_len(m), sizes),
  A = rep(rep(0:1, length.out = m), sizes)
)
random_intercept <- rep(rep(c(-1.4, -0.7, 0, 0.2, 0.8, 1.5), each = 2),
                        sizes)
trial$H <- rep(rnorm(m, sd = 0.5), sizes)
trial$x <- rep(rnorm(m, sd = 0.4), sizes) + rnorm(nrow(trial), sd = 0.7)
trial$continuous <- 3 + 0.5 * trial$A + 0.3 * trial$x + 0.2 * trial$H +
  random_intercept + rnorm(nrow(trial), sd = 0.8)
trial$binary <- rbinom(
  nrow(trial), 1,
  plogis(-0.2 + 0.5 * trial$A + 0.3 * trial$x + 0.2 * trial$H +
           random_intercept)
)
trial$count <- rpois(
  nrow(trial), exp(0.2 + 0.3 * trial$A + 0.15 * trial$x + 0.1 * trial$H +
                    0.7 * random_intercept)
)

## ----weighted-gee-example-----------------------------------------------------
weighted_gee <- MRStdCRT_fit(
  continuous ~ A + x + H, data = trial, cluster = "id", trt = "A",
  method = "GEE", family = gaussian(), corstr = "exchangeable",
  scale = "RD", model_weights = "inverse_cluster_size"
)
weighted_gee$estimate

## ----estimand-specific-example------------------------------------------------
unweighted_gee <- MRStdCRT_fit(
  continuous ~ A + x + H, data = trial, cluster = "id", trt = "A",
  method = "GEE", family = gaussian(), corstr = "exchangeable",
  scale = "RD", model_weights = "none"
)
estimand_specific_gee <- MRStdCRT_fit(
  continuous ~ A + x + H, data = trial, cluster = "id", trt = "A",
  method = "GEE", family = gaussian(), corstr = "exchangeable",
  scale = "RD", model_weights = "estimand_specific"
)
stopifnot(
  isTRUE(all.equal(estimand_specific_gee$estimate["cATE", ],
                   weighted_gee$estimate["cATE", ])),
  isTRUE(all.equal(estimand_specific_gee$estimate["iATE", ],
                   unweighted_gee$estimate["iATE", ]))
)
estimand_specific_gee$estimate
estimand_specific_gee$ics

## ----weighted-lmm-example-----------------------------------------------------
weighted_lmm <- MRStdCRT_fit(
  continuous ~ A + x + H, data = trial, cluster = "id", trt = "A",
  method = "LMM", family = gaussian(), scale = "RD",
  model_weights = "inverse_cluster_size"
)
weighted_lmm$estimate

## ----weighted-glmm-example----------------------------------------------------
weighted_glmm <- MRStdCRT_fit(
  binary ~ A + x + H, data = trial, cluster = "id", trt = "A",
  method = "GLMM", family = binomial(), scale = "OR",
  model_weights = "inverse_cluster_size"
)
summary(weighted_glmm)

## ----log-inference-example----------------------------------------------------
ratio_original <- MRStdCRT_fit(
  binary ~ A + x + H, data = trial, cluster = "id", trt = "A",
  method = "GEE", family = binomial(), corstr = "independence",
  scale = "OR", model_weights = "inverse_cluster_size"
)
ratio_log <- MRStdCRT_fit(
  binary ~ A + x + H, data = trial, cluster = "id", trt = "A",
  method = "GEE", family = binomial(), corstr = "independence",
  scale = "OR", model_weights = "inverse_cluster_size",
  inference_scale = "log"
)
stopifnot(
  identical(ratio_original$estimate[, c("Estimate", "Std. Error")],
            ratio_log$estimate[, c("Estimate", "Std. Error")]),
  identical(ratio_original$ics, ratio_log$ics)
)
ratio_original$estimate
summary(ratio_log)
ratio_log$log_se

## ----cs-gee-example-----------------------------------------------------------
set.seed(20261013)
cs_sizes <- rep(c(12, 20, 32, 48), 6)
cs_m <- length(cs_sizes)
assignment <- sample(rep(0:1, each = cs_m / 2))
cs_data <- data.frame(
  id = rep(seq_len(cs_m), cs_sizes), A = rep(assignment, cs_sizes),
  size = rep(cs_sizes, cs_sizes)
)
cs_data$y <- 2 + 0.02 * cs_data$size +
  cs_data$A * (0.4 + 0.01 * cs_data$size) +
  rep(rnorm(cs_m, sd = 0.6), cs_sizes) + rnorm(nrow(cs_data))

cs_mrs <- MRStdCRT_fit(
  y ~ A * size, data = cs_data, cluster = "id", trt = "A",
  method = "GEE", family = gaussian(), corstr = "exchangeable",
  scale = "RD", model_weights = "inverse_cluster_size"
)
cs_data$weight <- 1 / cs_data$size
cs_model <- geepack::geeglm(
  y ~ A * size, data = cs_data, id = id, weights = weight,
  family = gaussian(), corstr = "exchangeable"
)
stopifnot(cs_model$geese$error == 0)
ybar <- tapply(cs_data$y, cs_data$id, mean)
cluster_treatment <- tapply(cs_data$A, cs_data$id, mean)
g_mean <- correction <- numeric(2)
for (a in 0:1) {
  newdata <- cs_data
  newdata$A <- a
  eta <- tapply(predict(cs_model, newdata, type = "response"),
                cs_data$id, mean)
  g_mean[a + 1] <- mean(eta)
  correction[a + 1] <- mean((cluster_treatment == a) / 0.5 * (ybar - eta))
}
stopifnot(max(abs(correction)) < 1e-9,
          abs(cs_mrs$estimate[1, "Estimate"] - diff(g_mean)) < 1e-9)
rbind(g_computation = g_mean, residual_correction = correction)