--- title: "Introduction" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Introduction} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` ```{r setup} library(ClusterRandSSAdj) ``` # ClusterRandSSAdj The goal of ClusterRandSSAdj is to perform small-sample adjustment to generalized linear model based statistics when applied to one observation per cluster data from a cluster randomized trial. The small-sample adjustment method is based on the following four papers: Ford WP, Westgate PM. Improved standard error estimator for maintaining the validity of inference in cluster randomized trials with a small number of clusters. Biom J. 2017 May;59(3):478-495. doi: 10.1002/bimj.201600182. Epub 2017 Jan 27. PMID: 28128854. Westgate PM, Cheng DM, Feaster DJ, Fernández S, Shoben AB, Vandergrift N. Marginal modeling in community randomized trials with rare events: Utilization of the negative binomial regression model. Clin Trials. 2022 Apr;19(2):162-171. doi: 10.1177/17407745211063479. Epub 2022 Jan 6. PMID: 34991359; PMCID: PMC9038610. Kauermann G and Carroll RJ. A note on the efficiency of sandwich covariance matrix estimation. Journal of the American Statistical Association 2001; 96: 1387–1396. Mancl LA and DeRouen TA. A covariance estimator for GEE with improved small-sample properties. Biometrics 2001; 57: 126–134. ## Example Use The following example demonstrates how to use ClusterRandSSAdj. This example uses the simulated cluster randomized trial dataset included in ClusterRandSSAdj. The simulated trial dataset, named cluster_rct_data, contains the following variables: * clusterID * Population_size * Log_population_size * Treatment * Baserate * Log_baserate * Covar1 (0/1) * Covar2 (1/2/3) * Outcome Let's get started on demonstrating how ClusterRandSSAdj can be used for analysis by first fitting a Negative Binomial model to the example cluster randomized trial data. ```{r 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) ``` This package is designed to closely replicate the analysis results obtained from fitting a generalized linear model to a one-observation per cluster dataset using PROC GLIMMIX in SAS. As such, the following code estimates the SAS scale parameter for the Negative Binomial model and additionally calculates a Wald-type 95% confidence interval for the scale parameter: ```{r 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 ``` Note: The SAS scale parameter reported for a Negative Binomial model fit in PROC GLIMMIX is the dispersion parameter for the Negative Binomial model. In this example analysis, we apply small-sample adjustment to various statistics derived from the model fit. An initial task would be to make a small sample adjustment to the model parameters using the function ClusterRandSSAdj::ModelParmAdjEst(). The function applies a small-sample adjustment to the standard errors of the model parameters, then uses the adjusted SE to perform a t-test and construct 100(1-alpha)% confidence intervals. The following call to ModelParmAdjEst() will return 95% confidence intervals. ```{r paramAdj} ModelParmAdjEst(nb_model_test, alpha=0.05) ``` A next step in the analysis would be to apply small-sample adjustment to least square means (lsmeans) of interest. For example, consider the lsmeans from combinations of variables Treatment and Covar1. For lsmeans calculations, a right-handed formula construction is used to specify the desired lsmeans. For example, to request lsmeans for Treatment by Covar1, the right-handed formula should be ~Treatment:Covar1. As with ModelParmAdjEst(), LSMeansAdj() applies a small-sample adjustment to the standard error of each lsmean, then uses that adjusted SE to perform a t-test and construct 100(1-alpha)% confidence intervals. ```{r lsmeansAdj} LSMeansAdj(nb_model_test, formula_prt=~Treatment:Covar1, alpha=0.05) ``` Note: The lsmeans are reported on the log scale when using the Negative Binomial model. In the case of a Negative Binomial, the data scale is expressed as log of the expected count per analysis unit. For the cluster_rct_data dataset, the unit of analysis is a single person as indicated by the offset variable Log_population_size. To report these lsmeans in terms of events per population size, additional programming is needed. To calculate expected counts for each combination of Treatment by Covar 1, the log(expected count) values need to be exponentiated. The following code demonstrate a way to transform the log(expected count) to expected counts per 10,000 persons: ```{r 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)])) ``` In addition to performing small-sample adjustment on the lsmeans, ClusterRandSSAdj allows users to perform small-sample adjustment on pairwise comparisons of lsmeans using ClusterRandSSAdj::LSMeansPairwiseCompAdj(). The function's agruments are similar to those of LSMeansAdj(), with the addition of the reverse argument. For example, suppose we are interested in comparing log(expected counts) for Treatment = 1 minus Treatment = 0, when Covar1 = 0. If reverse = FALSE produces the comparison of log(expected counts) for Treatment = 0 minus Treatment = 1, when Covar1 = 0, setting reverse = TRUE will reverse the subtraction order so that the difference is oriented as desired. ```{r lsmeansAdjPWC} lsmean_comps <- LSMeansPairwiseCompAdj(nb_model_test, ~Treatment:Covar1, 0.05, reverse=TRUE) lsmean_comps ``` The lsmean_comps object contains several pairwise comparisons. In this analysis we focus on two specific comparisons: 1) log(expected counts) for Treatment = 1 minus Treatment = 0, when Covar1 = 0, and 2) log(expected counts) for Treatment = 1 minus Treatment = 0, when Covar1 = 1. In the R object, lsmeans_comps, these two comparisons are identified in the contrast variable as 'Treatment1 Covar10 - Treatment0 Covar10' and 'Treatment1 Covar11 - Treatment0 Covar11'. The following code extracts the two pairwise comparisons of interest and exponentiates the results so that they are reported on the scale of rate ratios. ```{r 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)])) ``` Now that our analysis has made small sample adjustments to model parameter estimates, lsmeans of interest, as well as pairwise comparisons of important lsmeans, let us consider how we might make small sample adjustments to user defined statistics. In SAS, procedures like PROC GLIMMIX allow the user to define estimate statements, which estimate important and allowable linear combinations of model parameters. This feature is helpful when a user wants to calculate a function of model parameters that is not equal to any estimate produced by the lsmeans statement. The function ClusterRandSSAdj::EstimateAdj() provides analgous functionality in R enabling the user to estimate a linear combination of model parameters and apply small-sample adjustment on the estimate's standard error. One illustrative use of ClusterRandSSAdj::EstimateAdj() is to replicate the lsmeans calculations returned by a call to ClusterRandSSAdj::LSMeansAdj(). For example, the following call to LSMeansAdj() will calculate the lsmeans for Treatment by Covar2: ```{r lsmAdjCov2} LSMeansAdj(nb_model_test, formula_prt=~Treatment:Covar2, alpha=0.05) ``` The following code containing a call to EstimateAdj() will replicate the lsmeans calculations for each combination of Treatment by Covar2 using a call to ClusterRandSSAdj::EstimateAdj(). ```{r 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 ``` While this may suggest that ClusterRandSSAdj::EstimateAdj() is redundant and replicates lsmeans calculations for Treatment by Covar2 that are produced by a call to ClusterRandSSAdj::LSMeansAdj(), let us continue to explore user defined linear combinations of model parameters using ClusterRandSSAdj::EstimateAdj(). For instance, consider pairwise comparisons of Treatment by Covar2 combinations, focusing on the Treatment = 1 minus Treatment = 0 comparisons within each level of Covar2. Covar2 has 3 levels, so we have 3 treatment comparisons, one in each level of Covar2. We can calculate these small-sample adjusted treatment effect for each level of Covar2 using ClusterRandSSAdj::EstimateAdj(). ```{r 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 ``` Note, similar to the individual lsmeans calculations for Tretment by Covar2, these calculation could have been calculated by a call to ClusterRandSSAdj::LSMeansPairwiseCompAdj() ```{r 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 ``` Once again, do not forget the results returned by LSMeansPairwiseCompAdj and by EstimateAdj() are on the log scale; therefore, in order to transform these statistics to effect ratios, the results in EstAdj2 (or lsmean_comps_sub2) will need to be exponentiated. ```{r estimates3b} cbind(EstAdj2[,1], exp(EstAdj2[,c(2, 7, 8)])) ``` At this point, we have reached an end of functions of linear combinations of model parameters which can be estimated using lsmeans and pairwise lsmean comparison functions such as LSMeansAdj() and LSMeansPairwiseCompAdj(). Now, let's look at differences between a specific set of pairwise comparisons, namely differences between Treatment 1 minus Treatment 0 log(expected count) estimates between different levels of Covar2. Since Covar2 has 3 levels, we have 3 treatment comparisons within each level of Covar2. 1) Estimates of Treatment 1 log(expected count) minus Treatment 0 log(expected count) when Covar2 is equal to 1 2) Estimates of Treatment 1 log(expected count) minus Treatment 0 log(expected count) when Covar2 is equal to 2 3) Estimates of Treatment 1 log(expected count) minus Treatment 0 log(expected count) when Covar2 is equal to 3 Given these three treatment effect estimates, we can compare levels of Covar2 with respect to these treatment effect estimates. For example, how different is the treatment effect when Covar2 is equal to 2, relative to when Covar2 is equal to 1? These estimates are difference in differences esimates and the estimates can be used to better understand the interaction between Treatment and Covar2. Given we have 3 treatment effects, we have 3 difference in difference estimates * Covar2 = 1 treatment effect vs. Covar2 = 2 treatment effect * Covar2 = 1 treatment effect vs. Covar2 = 3 treatment effect * Covar2 = 2 treatment effect vs. Covar2 = 3 treatment effect As such the contrast matrix is going to be a set of three differences in lsmean differences. We can estimate the parameter combination vector using the vectors in contrast_matrix2 above. logT1vT0_1 minus logT1vT0_2 will be: * c(0, 1, 0, 0, 0, 0, 0, 0, 0.5) - c(0, 1, 0, 0, 0, 0, 1, 0, 0.5) = c(0, 0, 0, 0, 0, 0, -1, 0, 0) logT1vT0_1 minus logT1vT0_3 will be: * c(0, 1, 0, 0, 0, 0, 0, 0, 0.5) - c(0, 1, 0, 0, 0, 0, 0, 1, 0.5) = c(0, 0, 0, 0, 0, 0, 0, -1, 0) logT1vT0_2 minus logT1vT0_3 will be: * c(0, 1, 0, 0, 0, 0, 1, 0, 0.5) - c(0, 1, 0, 0, 0, 0, 0, 1, 0.5) = c(0, 0, 0, 0, 0, 0, 1, -1, 0) These estimates cannot be calculated in a call to any lsmeans ot pairwise lsmeans comparison function. A linear combination of parameter estimates function is necessary, and this is where EstimateAdj() becomes very useful. ```{r 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 ``` The estimates in dodEst are differences in differences on the log(expected counts) scale. These differences in differences can be transformed into ratios of expected count ratios, providing an expression of interaction on the ratio scale. Specifically, if exp(differences in log(expected count)), is a ratio of expected counts, then exp(differences in log(expected count1) - differences in log(expected count2)) yields a ratio of expected count ratios. ```{r estimates4} new_lab <- c("T1vT0_1_2", "T1vT0_1_3", "T1vT0_2_3") cbind(new_lab, exp(dodEst[,c(2, 7, 8)])) ``` Next, let us consider another function within ClusterRandSSAdj, TypeIIIAdj(). In SAS procedures, Type III tests are produced for predictors included in a model. When PROC GLIMMIX is used, the tests are F-tests. The car::Anova() function in the car package can produce Type III tests for generalized linear models, however the tests are Chi-square tests. The TypeIIIAdj() function applies small-sample adjustment to the Type III tests as implemented in the car::Anova() function. Within a call to TypeIIIAdj(), a call to car::Anova() is made and it uses "type = 'III'" and "test.statistic = "Wald" to approximate the F-test in SAS Type III test results. For type III tests similar to those reported in SAS, the user needs to create a list of contrast structures for each categorical variable, where the contrast structure is "contr.sum". Once these are defined a call to TypeIIIAdj() will produce small-sample adjusted type III test results analogous to what SAS produces however these results are based on chi-square tests instead of F-tests. ```{r typeIII} type3_lst <- list(Treatment = contr.sum, Covar1 = contr.sum, Covar2 = contr.sum) TypeIIIAdj(nb_model_test, type3_lst) ``` These small-sample adjusted Type III tests can be used as a gate-keeping step in which statistics for combinations of multi-level covariates are reported only if the type III test for the interaction term is statistically significant. Such a test may involve at least two degrees of freedom, thus a type III test may be appropriate. In this case, small-sample adjustment can be applied to the test using TypeIIIAdj(). The chi-square test that is used in TypeIIIAdj() makes it impossible to replicate the Type III test produced by GLIMMIX. To address this, two additional functions were developed in ClusterRandSSAdj. These functions allow the users to replicate the SAS Type III test and to estimate and apply small-sample adjustment to any linear combination of the model parameters using an F-test. The function SAS_Lincom_Ftest() allows the user to complete an F-test on any linear combination of model parameters and to adjust the standard error of the estimate using HC3 (FIRORES), HC2 (Root), or no adjustment. One requirement for use of SAS_Lincomp_Ftest() is that the user must specify the desired linear combination. For example, the following code replicates the Type III test from PROC GLIMMIX when analyzing cluster_rct_data in PROC GLIMMIX using the empirical=FIRORES option. ```{r 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") ``` Finally, a small-sample adjustment can be applied on any desired linear combination of model parameters. For example, the Treatment by Covar1 interaction: ```{r typeIII_2B} Lincom_FtestAdj(nb_model_test, L_value) ```