## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup--------------------------------------------------------------------
library(ClusterRandSSAdj)

## ----model_run----------------------------------------------------------------
library(ClusterRandSSAdj)
library(MASS)
nb_model_test <- glm.nb(Outcome ~ Treatment + Log_baserate + Covar1 + Covar2 +
                      Treatment:Covar2 + Treatment:Covar1 + offset(Log_population_size), data = cluster_rct_data)

# Summarize the model
summary(nb_model_test)

## ----scale--------------------------------------------------------------------
scale_est <- 1/nb_model_test$theta
scale_95CI <- 1/(c(nb_model_test$theta - 1.96*nb_model_test$SE.theta, nb_model_test$theta + 1.96*nb_model_test$SE.theta))
scale_results <- data.frame(SAS_scale = scale_est, Lower95 = scale_95CI[2], Upper95 = scale_95CI[1])
scale_results

## ----paramAdj-----------------------------------------------------------------
ModelParmAdjEst(nb_model_test, alpha=0.05)

## ----lsmeansAdj---------------------------------------------------------------
LSMeansAdj(nb_model_test, formula_prt=~Treatment:Covar1, alpha=0.05)

## ----lsmeansAdjB--------------------------------------------------------------
rates <- LSMeansAdj(nb_model_test, ~Treatment:Covar1, 0.05)[, c(1,2,4,8,9)]
cbind(rates[,c(1,2)], 10000*exp(rates[,c(3,4,5)]))

## ----lsmeansAdjPWC------------------------------------------------------------
lsmean_comps <- LSMeansPairwiseCompAdj(nb_model_test, ~Treatment:Covar1, 0.05, reverse=TRUE)
lsmean_comps

## ----lsmAdjPWCs---------------------------------------------------------------
lsmean_comps_sub <- lsmean_comps |> dplyr::filter(contrast %in% c("Treatment1 Covar10 - Treatment0 Covar10",
                                                           "Treatment1 Covar11 - Treatment0 Covar11"))
cbind(lsmean_comps_sub[,1], exp(lsmean_comps_sub[,c(3, 7, 8)]))

## ----lsmAdjCov2---------------------------------------------------------------
LSMeansAdj(nb_model_test, formula_prt=~Treatment:Covar2, alpha=0.05)

## ----estimates----------------------------------------------------------------
contrast_matrix <- rbind(
  "logT0_1" = c(1, 0,1.937217, 0.5, 0, 0, 0, 0, 0),
  "logT1_1" = c(1, 1,1.937217, 0.5, 0, 0, 0, 0, 0.5),
  "logT0_2" = c(1, 0,1.937217, 0.5, 1, 0, 0, 0, 0),
  "logT1_2" = c(1, 1,1.937217, 0.5, 1, 0, 1, 0, 0.5),
  "logT0_3" = c(1, 0,1.937217, 0.5, 0, 1, 0, 0, 0),
  "logT1_3" = c(1, 1,1.937217, 0.5, 0, 1, 0, 1, 0.5)
)
lsmean_replicate <- EstimateAdj(nb_model_test, contrast_matrix, 0.05)
lsmean_replicate

## ----estimates2---------------------------------------------------------------
contrast_matrix2 <- rbind(
  "logT1vT0_1" = c(0, 1, 0, 0, 0, 0, 0, 0, 0.5),
  "logT1vT0_2" = c(0, 1, 0, 0, 0, 0, 1, 0, 0.5),
  "logT1vT0_3" = c(0, 1, 0, 0, 0, 0, 0, 1, 0.5)
)
EstAdj2 <- EstimateAdj(nb_model_test, contrast_matrix2, 0.10)
EstAdj2

## ----estimates3---------------------------------------------------------------
lsmean_comps2 <- LSMeansPairwiseCompAdj(nb_model_test, ~Treatment:Covar2, 0.05, reverse=TRUE)
lsmean_comps_sub2 <- lsmean_comps2 |> dplyr::filter(contrast %in% c("Treatment1 Covar21 - Treatment0 Covar21",
                                                                    "Treatment1 Covar22 - Treatment0 Covar22",
																	"Treatment1 Covar23 - Treatment0 Covar23"))
lsmean_comps_sub2

## ----estimates3b--------------------------------------------------------------
cbind(EstAdj2[,1], exp(EstAdj2[,c(2, 7, 8)]))

## ----estimates2B--------------------------------------------------------------
contrast_matrix3 <- rbind(
"logT1vT0_1mlogT1vT0_2" = c(0, 0, 0, 0, 0, 0, -1,  0, 0),
"logT1vT0_1mlogT1vT0_3" = c(0, 0, 0, 0, 0, 0,  0, -1, 0),
"logT1vT0_2mlogT1vT0_3" = c(0, 0, 0, 0, 0, 0,  1, -1, 0)
)
dodEst <- EstimateAdj(nb_model_test, contrast_matrix3, 0.10)
dodEst

## ----estimates4---------------------------------------------------------------
new_lab <- c("T1vT0_1_2", "T1vT0_1_3", "T1vT0_2_3")
cbind(new_lab, exp(dodEst[,c(2, 7, 8)]))

## ----typeIII------------------------------------------------------------------
type3_lst <- list(Treatment = contr.sum, Covar1 = contr.sum, Covar2 = contr.sum)
TypeIIIAdj(nb_model_test, type3_lst) 

## ----typeIII_2----------------------------------------------------------------

# Type III F-test for treatment parameter
L_value <- matrix(c(0, 1, 0, 0, 0, 0, 0.333333333333, 0.333333333333, 0.5), nrow=1, ncol=9)
SAS_Lincom_Ftest(nb_model_test, L_value, "HC3")

#The SAS type III test for Log_baserate?
L_value <- matrix(c(0, 0, 1, 0, 0, 0, 0, 0, 0), nrow=1, ncol=9)
SAS_Lincom_Ftest(nb_model_test, L_value, "HC3")

# Type III F-test for covar1 parameter
L_value <- matrix(c(0, 0, 0, 1, 0, 0, 0, 0, 0.5), nrow=1, ncol=9)
SAS_Lincom_Ftest(nb_model_test, L_value, "HC3")

# Type III F-test for covar2 parameter
L_value <- rbind(
  c(0, 0, 0, 0, 1, 0, 0.5,   0, 0),
  c(0, 0, 0, 0, 0, 1, 0,   0.5, 0)
)
SAS_Lincom_Ftest(nb_model_test, L_value, "HC3")

#The SAS type III test for Treatment:Covar2?
L_value <- rbind(
  c(0, 0, 0, 0, 0, 0, 1, 0, 0),
  c(0, 0, 0, 0, 0, 0, 0, 1, 0)
)
SAS_Lincom_Ftest(nb_model_test, L_value, "HC3")

#The SAS type III test for Treatment:Covar1?
L_value <- matrix(c(0, 0, 0, 0, 0, 0, 0, 0, 1), nrow=1, ncol=9)
SAS_Lincom_Ftest(nb_model_test, L_value, "HC3")


## ----typeIII_2B---------------------------------------------------------------
  Lincom_FtestAdj(nb_model_test, L_value)