What is already known?

• • Existing effect size measures like BC-SMD compare SCEDs with group designs but have limitations for count data due to the violation of normality and homoscedasticity assumptions, and sensitivity to the session length or time units when observing target behaviors.

• • A between-case incidence rate ratio (BC-IRR) has recently been proposed within the generalized linear mixed model (GLMM) framework. However, its computation is relatively complex, and the finite-sample statistical properties of the GLMM-based estimator have not been evaluated.

What is new?

• • This study introduces a generalized estimating equation (GEE)-based estimator for BC-IRR and systematically compares its performance with the existing GLMM-based estimator.

• • Large-scale Monte Carlo simulations demonstrate that both the GLMM- and GEE-based estimators provide unbiased estimates of BC-IRR and its standard errors across a wide range of realistic SCED conditions.

Potential impact for RSM readers

• • The satisfactory statistical properties of BC-IRR enable comparisons of SCED and group design studies, benefiting research syntheses with findings from both designs.

• • It enhances methodological rigor, providing researchers with two flexible frameworks for estimating effect sizes that are better suited for SCED count data.

1.1. Causal definition of BC-IRR for SCEDs

The marginal treatment effect aforementioned for SCEDs could be formally defined by utilizing causal inference framework. Let us consider an AB design with m participants and series length of n. For simplicity, assume that all participants enter the treatment phase at time T. This assumption is not required, as the definition also applies to MBDs in which participants enter the treatment phase at different time points. According to Luo et al., BC-IRR is defined as a causal estimand on a ratio scale:

(2) $$\begin{align}\mathrm{BC}-\mathrm{IRR}={\frac{E\left[{Y}_B(1)\right]}{E\left[{Y}_B(0)\right]}}^{29}\end{align}$$

The numerator represents the expected incidence rate at time B if a treatment strategy is implemented between times T and B. The denominator represents the expected incidence rate at time B if the treatment is never adopted. Both the numerator and the denominator are marginal means, representing the average incidence rates across the m participants under the treated and the untreated conditions. Thus, BC-IRR is a marginal effect.

Because the estimand is defined based on potential outcomes, it is necessary to link the potential outcomes to the observed data. The numerator can be directly linked to the observed incidence rate at time B because all participants receive the treatment. However, the denominator, the potential outcome under no treatment at time B, is unobserved and must be estimated using models based on observations from the baseline phase. A key identifiability condition, therefore, is that this potential outcome model is correctly specified. Additional details of the statistical models are presented in Section 2. BC-IRR has been shown to be equivalent to the IRR that would have been obtained in a hypothetical RCT, under the assumption that a never-treated group exists in which the distribution of individual effects, defined as the time-invariant effects unique to each individual, matches that of the treated units in the SCED. Due to space limitations, we do not repeat the detailed derivations here. Readers are referred to Luo et al. for the complete logical and mathematical arguments.Reference Luo, Li, Baek and Li ²⁸

2.1. Generalized linear mixed model

GLMMs are a set of conditional models that integrate the frameworks of GLM and linear mixed models (LMM). Hence, GLMMs could take the source of clustering or nested data structure in SCEDs into account by directly modeling between-case variation using random effects. With an extension to generalized framework, GLMMs could also deal with non-normal outcomes in SCEDs, such as binary, percentage/proportion, count, and rate data using various distributions (e.g., Poisson, negative binomial, binomial, beta-binomial).

We used the following Poisson GLMM to derive BC-IRR from count outcomes in an AB design:

(3) $$\begin{align}{y}_{ij}&\sim \mathrm{Poisson}\;\left({\lambda}_{ij}\right)\\\log \left({\lambda}_{ij}\right)&={\beta}_{0j}+{\beta}_{1j}{Time}_{ij}+{\beta}_{2j}{Phase}_{ij}+{\beta}_{3j}{Time}_{ij}^{\prime }{Phase}_{ij}\nonumber\end{align}$$

where ${Phase}_{ij}$ is an indicator variable equal to 0 for control phase and 1 for treatment phase; ${Time}_{ij}$ is a time variable representing the time (measured either in calendar units or sessions) elapsed from a reference time point. Researchers can choose any time point as a reference time point, which changes the meaning of the intercept ${\beta}_{0j}$ and its associated variance–covariance components. ${Time}_{ij}^{\prime }$ is another centered time variable which equals zero at the first observed treatment session, 1 at the second treatment session, and so forth, following the coding proposed by Huitema and McKean.Reference Huitema and Mckean ³¹ Readers interested in the alternative coding scheme for the treatment-by-time interaction ( ${Time}_{ij}^{\prime }$ ) are referred to Pustejovsky et al.Reference Pustejovsky, Hedges and Shadish ²⁶ Based on this coding scheme, ${e}^{\beta_{0j}}$ represents the expected outcome level at the reference time point for case j, ${e}^{\beta_{1j}}$ represents the proportional change per session during the baseline phase, ${e}^{\beta_{2j}}$ represents the immediate proportional change in the outcome level at the first observed treatment session relative to the projected baseline trend (i.e., the within-case IRR for immediate effect), and ${e}^{\beta_{3j}}$ represents the trend/slope change after the treatment (i.e., the within-case IRR for trend change). An example coding matrix for all the variables is shown in Table 1, where $C$ is the reference time point, the treatment is implemented at time $T=A$ , and the outcome $Y$ is measured at time point $B$ . Based on a multivariate normal assumption for all random coefficients (i.e., ${\beta}_{0j}-{\beta}_{3j}$ )

(4) $$\begin{align}\left[\begin{array}{c}\begin{array}{c}{\beta}_{0j}\\ {}{\beta}_{1j}\end{array}\\ {}{\beta}_{2j}\\ {}{\beta}_{3j}\end{array}\right]\sim \mathrm{MVN}\;\left[\begin{array}{c}\begin{array}{c}{\gamma}_{00}\\ {}{\gamma}_{10}\end{array}\\ {}{\gamma}_{20}\\ {}{\gamma}_{30}\end{array}\right]\left[\begin{array}{cc}\begin{array}{cc}{\sigma}_{u0}^2& {\sigma}_{u0u1}\\ {}{\sigma}_{u0u1}& {\sigma}_{u1}^2\end{array}& \begin{array}{cc}{\sigma}_{u0u2}& {\sigma}_{u0u3}\\ {}{\sigma}_{u1u2}& {\sigma}_{u1u3}\end{array}\\ {}\begin{array}{cc}{\sigma}_{u0u2}& {\sigma}_{u1u2}\\ {}{\sigma}_{u0u3}& {\sigma}_{u1u3}\end{array}& \begin{array}{cc}{\sigma}_{u2}^2& {\sigma}_{u2u3}\\ {}{\sigma}_{u2u3}& {\sigma}_{u3}^2\end{array}\end{array}\right]\end{align}$$

and the coding matrix in Table 1

(5) $$\begin{align}\mathrm{BC}-\mathrm{IRR}&=\frac{E\left[{Y}_B(1)\right]}{E\left[{Y}_B(0)\right]}\nonumber\\ &={e}^{\gamma_{20}+\frac{1}{2}{\sigma}_{u2}^2+{\sigma}_{u0u2}+\left(B-C\right){\sigma}_{u1u2}+{\gamma}_{30}\left(B-\mathrm{A}\right)+\frac{1}{2}{\left(B-A\right)}^2{\sigma}_{u3}^2+\left(B-A\right){\sigma}_{u0u3}+\left(B-A\right){\sigma}_{u2u3}+\left(B-C\right)\left(B-A\right){\sigma}_{u1u3}},\end{align}$$

where $B-C$ represents the time elapsed between the reference time point ( $C$ ) and the outcome measurement occasion $B$ , while $B-A$ represents the time elapsed between the introduction of treatment and the outcome measurement occasion. To compute BC-IRR, the estimates of the conditional effects ( $\gamma$ ) and variance/covariance components ( $\sigma$ ) from the GLMM are substituted into equation (5). The full derivation process for equation (5) can be found in Luo et al.Reference Luo, Li, Baek and Li ²⁸

Table 1.

Example coding matrix for a case with count data

Table 1. Long description

The table consists of six columns and twelve rows of data.

Column headers from left to right are:

1. Measurement occasions (I)

2. Time sub i j

3. Phase sub i j

4. Time prime sub i j Phase sub i j, defined as I minus (C plus 1)

5. B minus A

6. Outcome (Y)

Data rows:

- Occasion 1: Time -5, Phase 0, Time prime 0, B-A is blank, Outcome 1.

- Occasion 2: Time -4, Phase 0, Time prime 0, B-A is blank, Outcome 3.

- Occasion 3: Time -3, Phase 0, Time prime 0, B-A is blank, Outcome 2.

- Occasion 4: Time -2, Phase 0, Time prime 0, B-A is blank, Outcome 4.

- Occasion 5: Time -1, Phase 0, Time prime 0, B-A is blank, Outcome 2.

- Occasion 6 (C): Time 0, Phase 0, Time prime 0, B-A is blank, Outcome 3.

- Occasion 7 (A): Time 1, Phase 1, Time prime 0, B-A 0, Outcome 10.

- Occasion 8: Time 2, Phase 1, Time prime 1, B-A 1, Outcome 9.

- Occasion 9: Time 3, Phase 1, Time prime 2, B-A 2, Outcome 11.

- Occasion 10: Time 4, Phase 1, Time prime 3, B-A 3, Outcome 6.

- Occasion 11: Time 5, Phase 1, Time prime 4, B-A 4, Outcome 9.

- Occasion 12: Time 6, Phase 1, Time prime 5, B-A 5, Outcome 8.

The Poisson GLMM in equation 3 assumes that the conditional mean and variance for the outcome are identical. However, when count data exhibit more variability than predicted by the Poisson model, overdispersion occurs. Failing to account for overdispersion leads to poor model fit and misleading statistical inferences due to underestimated standard errors of regression coefficients.Reference Li, Luo, Baek, Thompson and Lam ^{12 ,} Reference Agresti ³² Pustejovsky et al. found that overdispersion is not uncommon for count outcomes in SCEDs.Reference Pustejovsky, Swan and English ⁸ To accommodate overdispersion, the following Poisson-Gamma mixture model, also called the negative binomial model (NB), has been introduced to analyze overdispersed count data in SCEDs (Li et al., 2025; Li et al., 2024).

(6) $$\begin{align}{y}_{ij}&\sim Poisson\;\left({\lambda}_{ij}\right)\nonumber\\\log \left({\lambda}_{ij}\right)&={\beta}_{0j}+{\beta}_{1j}{Time}_{ij}+{\beta}_{2j}{Phase}_{ij}+{\beta}_{3j}{Time}_{ij}^{\prime }{Phase}_{ij}+{e}_{ij}\nonumber\\\exp \left({e}_{ij}\right)&\sim \mathrm{Gamma}\left(\theta, 1/\theta \right),\theta >0.\end{align}$$

With the log link function and the same multivariate assumption for random coefficients $\beta$ , the calculation for BC-IRR as shown in equation (5) is exactly the same for the NB model.

2.2. Marginal model via generalized estimating equation

Rather than explicitly modeling the clustering mechanism via random effects as is done with the conditional models of GLMM, marginal models treat within-subject association as a nuisance characteristic of data that should be accounted for to make correct inferences about changes in the population mean response.Reference Liang and Zeger ³⁰ A marginal model for longitudinal data has the following three-part specificationReference Fitzmaurice, Laird and Ware ³³ : (1) The marginal expectation $E\! \left({y}_{ij}|{\boldsymbol{x}}_{ij}\right)={\mu}_{ij}$ , depends only on the covariates ${\boldsymbol{x}}_{ij}$ , through a known link function $g\! \left({\mu}_{ij}\right)={\boldsymbol{x}}_{ij}^{\prime}\boldsymbol{\beta}$ ; (2) the variance of ${Y}_{ij}$ , conditional on the covariates, is given by $Var\! \left({y}_{ij}|{\boldsymbol{x}}_{ij}\right)=\phi v\left({\mu}_{ij}\right)$ , where $v\! \left({\mu}_{ij}\right)$ is a known variance function of the conditional mean ${\mu}_{ij}$ and $\phi$ is a scale parameter; (3) the covariance matrix for the within-subject responses, conditional on the covariates, is given by ${\boldsymbol{V}}_j=\phi {\sqrt{\boldsymbol{A}}}_j\mathrm{Corr}\left({\boldsymbol{Y}}_j\right)\sqrt{{\boldsymbol{A}}_j}$ , where ${\boldsymbol{A}}_j$ is a diagonal matrix with components $v\! \left({\mu}_{ij}\right)$ , and $\mathrm{Corr}\! \left({\boldsymbol{Y}}_j\right)$ is the correlation matrix as the function of a vector of parameters $\boldsymbol{\alpha}$ that estimate the correlations of observations within subjects.

When the outcomes in the marginal model are discrete (e.g., binary, count, or ordinal categorical data), there is no convenient specification of the joint distribution of ${\boldsymbol{Y}}_j$ . Hence, it requires an alternative to maximum likelihood estimation and leads to a method known as GEE. GEE approach essentially views the model as a single-level model and applies statistical corrections to produce standard error estimates that account for clustering.Reference Zeger and Liang ^{29 ,} Reference Liang and Zeger ³⁰ The advantage of GEE is that the covariance structure of the outcome does not have to be correctly specified, meaning that there are far fewer assumptions required compared to random effects models.Reference Zeger, Liang and Albert ³⁴ For the estimation, GEE requires researchers to specify a working correlation matrix $\mathrm{Corr}\left({\boldsymbol{Y}}_j\right)$ , a modeling “device” used to describe the association among repeated observations within the same cluster (e.g., measurements within a case in a SCED). The correlation is “working” because it does not need to be correctly specified to get consistent estimates of the marginal parameter estimates. Several working correlation structures are potentially suitable for SCED application, including independence (i.e., no within-case correlation), exchangeable (i.e., constant correlation among all observations within a case), autoregressive (i.e., correlation decays with time lag), and unstructured (i.e., each pair of observations has its own correlation).

After specifying a working correlation structure, an iterative estimation process is started using weighted least squares (quasi-likelihood) by (1) estimating the regression coefficients ( $\boldsymbol{\beta}$ ) as if data are independent, (2) estimating the parameters ( $\boldsymbol{\alpha}$ ) in the working correlation matrix to handle clustering, (3) using the estimated working correlation matrix to update the outcome variable covariance matrix ( ${\boldsymbol{V}}_j$ ), (4) using ${\boldsymbol{V}}_j$ to update estimates of $\boldsymbol{\beta}$ ; and (5) continuing this iteration process until convergence and then applying sandwich estimators for standard errors. For an approachable illustration of the full estimation procedure for GEE, readers are referred to McNeish and Stapleton.Reference McNeish and Stapleton ³⁵

In SCEDs, BC-IRR could be estimated by specifying a marginal model using the GEE approach with an appropriate working correlation matrix. The marginal model is specified as follows for count outcomes:

(7) $$\begin{align}E\! \left({y}_{ij}|{\boldsymbol{x}}_{ij}\right)&={\lambda}_{ij},\log \! \left({\lambda}_{ij}\right)={\beta}_0+{\beta}_1{Time}_{ij}+{\beta}_2{Phase}_{ij}+{\beta}_3{Time}_{ij}^{\prime }{Phase}_{ij}.\nonumber\\Var\! \left({y}_{ij}|{\boldsymbol{x}}_{ij}\right)&=\phi v\! \left({\lambda}_{ij}\right)\nonumber\\{\boldsymbol{V}}_{\! j}&=\phi {\boldsymbol{A}}_j^{\frac{1}{2}}{\boldsymbol{R}}_{\! j}\left(\boldsymbol{\alpha} \right){\boldsymbol{A}}_j^{\frac{1}{2}}\end{align}$$

where the Poisson variance function $v\left({\lambda}_{ij}\right)={\lambda}_{ij}$ is adopted, and the scale parameter $\phi$ could handle potential overdispersion in SCED count outcomes. The working correlation structure ${\boldsymbol{R}}_{\! j}\left(\boldsymbol{\alpha} \right)$ describes the correlation pattern of observations within-subject with $\boldsymbol{\alpha}$ as a vector of association parameters depending on the correlation structure. In our study, we adopt the “independence” correlation structure to avoid convergence issues encountered with very small sample sizes in SCEDs. We note that other correlation structures including, but not limited to, “unstructured”, “autoregressive”, “Toeplitz”, and “exchangeable” are available for the GEE approach. Based on equations (2) and (7),

(8) $$\begin{align}\mathrm{BC}-\mathrm{IRR}={e}^{\left({\beta}_2+{\beta}_3\left(B-A\right)\right)},\end{align}$$

where ${\beta}_2$ and ${\beta}_3$ represent the marginal estimates of the immediate effect and trend change, respectively. A key advantage of the GEE over GLMM is that it does not require specifying the level-2 variance–covariance structure of random effects. The GEE approach could provide robust estimates of standard errors if all relevant predictors are included, regardless of the proper specification of the working correlation matrix of the outcome measurements.Reference McNeish and Harring ³⁶

2.3. Delta methods for standard errors

As BC-IRR is essentially a combination or a transformation of parameter estimates related to the immediate effect and regarding trend change, we need to further calculate its standard error. As outlined in equations ((5)) and (8), the transformation is relatively complex for the GLMM and more intuitive for the GEE approach. To get standard errors for the transformed parameter, the delta method is typically used. The delta method is essentially calculating the variance of the Taylor series approximation of a function ( $G$ ) of random variables ( $\boldsymbol{x}$ = (x1, x2, …)) with known variance. The first two terms of the Taylor expansion are an approximation for the function $G$ ,

(9) $$\begin{align}G \! \left(\boldsymbol{x}\right)\cong G \! \left({\boldsymbol{\mu}}_x\right)+{\left(\frac{d \! \left(G \! \left({\boldsymbol{\mu}}_x\right)\right)}{d \! \left(\boldsymbol{x}\right)}\right)}^T\left(\boldsymbol{x}-{\boldsymbol{\mu}}_x\right),\end{align}$$

where ${\boldsymbol{\mu}}_x$ is the mean vector of $\boldsymbol{x}$ and $\frac{d \! \left(G \! \left({\boldsymbol{\mu}}_x\right)\right)}{d \! \left(\boldsymbol{x}\right)}$ is a vector of partial derivatives of $G \! \left(\boldsymbol{x}\right)$ at point ${\boldsymbol{\mu}}_x$ . We can then approximate the variance of $G \! \left(\boldsymbol{x}\right)$ ,

(10) $$\begin{align}\operatorname{var} \! \left(G \! \left(\boldsymbol{x}\right)\right)\cong {\left(\frac{d \! \left(G \! \left({\boldsymbol{\mu}}_x\right)\right)}{d \! \left(\boldsymbol{x}\right)}\right)}^T\operatorname{cov} \! \left(\boldsymbol{x}\right)\left(\frac{d \! \left(G \! \left({\boldsymbol{\mu}}_x\right)\right)}{d \! \left(\boldsymbol{x}\right)}\right),\end{align}$$

where $\operatorname{cov} \! \left(\boldsymbol{x}\right)$ is the known variance–covariance matrix of $\boldsymbol{x}$ . The square root of the outcome based on equation 10 is the standard error of the transformed parameter estimate. To calculate the standard errors of BC-IRR estimated by the GEE approach as shown in equation (8), $G \! \left(\boldsymbol{x}\right)=G \! \left(\boldsymbol{\beta} \right)={e}^{\left({\beta}_2+{\beta}_3\left(B-A\right)\right)}$ and $\operatorname{cov} \! \left(\boldsymbol{x}\right)=\operatorname{cov} \! \left(\boldsymbol{\beta} \right)=\operatorname{cov} \! \left({\beta}_2,{\beta}_3\right)$ . Similarly, to obtain the standard error of BC-IRR estimated by the GLMM as shown in equation ((5)), $G \! \left(\boldsymbol{x}\right)=G \! \left(\boldsymbol{\gamma}, \boldsymbol{\sigma} \right)$ and $\operatorname{cov} \! \left(\boldsymbol{x}\right)=\operatorname{cov} \! \left(\boldsymbol{\gamma}, \boldsymbol{\sigma} \right)$ . The calculations via the delta method can be conducted using the deltamethod function in the R package msm.Reference Jackson ³⁷

2.4. Performance of GLMM with small sample size

Recent studies have systematically evaluated the performance of GLMMs in the SCED context.Reference Li, Luo, Baek, Thompson and Lam ^{12 ,} Reference Li, Luo and Baek ^{13 ,} Reference Declercq, Jamshidi and Fernández-Castilla ³⁸ They have demonstrated that GLMMs are promising methods to deal with overdispersed and zero-inflated SCED count data. These studies also provide tailored guidelines for model selection based on the properties of SCED count data. For example, Declercq et al. found that when count data are not overdispersed, Poisson GLMMs can produce accurate estimates of treatment effects evaluated by WC-IRR.Reference Declercq, Jamshidi and Fernández-Castilla ³⁸ In recent simulation studies, Li et al. also found that GLMMs can provide accurate estimates and reliable inferential statistics for treatment effects measured by WC-IRR in MBDs, when researchers select appropriate distributions (e.g., Poisson, negative binomial, zero-inflated Poisson, and zero-inflated negative binomial).Reference Li, Luo and Baek ¹³ In addition, there are few studies directly comparing the performance of GLMM and GEE under small sample sizes in more general contexts.Reference McNeish and Stapleton ^{35 ,} Reference McNeish and Harring ³⁶ They found that GLMM outperforms GEE in terms of standard error bias, coverage rate, and statistical power for both continuous and binary outcomes. More specifically, for both continuous and binary outcomes with a small number of clustering units (≤ 50), GLMM with degrees of freedom adjusted method (i.e., Kenward-Roger adjustment) showed minimal bias in the standard error estimates of fixed effects and had the coverage rate close to the nominal level of 95%.

2.5. Performance of GEE with small sample size

Unlike solid evidence provided for the statistical properties of GLMMs, there have been no direct studies of GEE applied to the SCED context. However, the performance of GEE for correlated binary and continuous data has been well documented. In the large sample framework, it is known that estimates of regression parameters of GEE are generally unbiased regardless of the working correlation matrix, but misspecifications of the working correlation matrix leads to biased estimates of the standard errors of the regression estimates. Thus, it is typical practice to obtain standard errors based on the robust or sandwich estimator proposed by Liang and ZegerReference Liang and Zeger ³⁰ for the variance–covariance matrix of regression coefficients, which is consistent even under misspecification of the working correlation matrix.Reference Zorn ³⁹

With small sample size, previous research found that GEE does not perform well.Reference McNeish and Stapleton ^{35 ,} Reference McNeish and Harring ^{36 ,} Reference Thompson, Hemming, Forbes, Fielding and Hayes ⁴⁰ Specifically, McNeish and StapletonReference McNeish and Stapleton ³⁵ found that GEEs had relatively large bias for fixed effects estimates regarding level-2 predictors and cross-level interactions when dealing with a very small number of clustering units (e.g., n = 4). The 95% confidence interval based on the classical GEE with the sandwich estimator for standard errors often produced significant under-coverage across different levels of number of units (n = 4 to 14).

To enhance both parameter estimates and standard errors of the estimates under small sample sizes, several bias correction methods have been proposed.Reference Thompson, Hemming, Forbes, Fielding and Hayes ^{40 ,} Reference Paul and Zhang ⁴¹ For example, Paul and Zhang demonstrated that by using bias-adjusted approaches, unbiased parameter estimates can be obtained for small to moderate sample sizes (n = 15 to 50).Reference Paul and Zhang ⁴¹ Some other researchers have proposed various corrections for the sandwich estimator to obtain accurate statistical inferences for the GEE regression parameter estimators under small sample sizes, including FG (Fay and Graubard),Reference Fay and Graubard ⁴² GST (Gosho et al),Reference Gosho, Sato and Takeuchi ⁴³ KC (Kauermann and Carroll),Reference Kauermann and Carroll ⁴⁴ MK (MacKinnon),Reference MacKinnon and White ⁴⁵ MD (Mancl and DeRouen),Reference Mancl and DeRouen ⁴⁶ MBN (Morel et al.),Reference Morel and Bokossa ⁴⁷ PAN (Pan),Reference Pan ⁴⁸ and WL (Wang and Long).Reference Wang and Long ⁴⁹ Recently, Thompson et al. compared the performance of several methods to improve the standard error estimates of GEE with binary outcomes in small samples (n = 6 to 54).Reference Thompson, Hemming, Forbes, Fielding and Hayes ⁴⁰ They found that some bias-corrected standard error methods (i.e., KC and FG standard errors) outperformed others to mitigate this bias when used with an independent working correlation matrix and were robust to the misspecification of the working correlation matrix. However, Thompson et al. noted significant overcorrections for the FG standard error with a very small number of clusters (n = 6),Reference Thompson, Hemming, Forbes, Fielding and Hayes ⁴⁰ which is also consistent with other findings.Reference Scott, Decamp, Juraska, Fay and Gilbert ^{50 ,} Reference Leyrat, Morgan, Leurent and Kahan ⁵¹

2.6. Small-sample correction methods

For GLMMs, researchers could use Kenward-RogerReference Kenward and Roger ^{52 ,} Reference Kenward and Roger ⁵³ or SatterthwaiteReference Satterthwaite ⁵⁴ approaches to adjust for standard errors and degrees of freedom of fixed-effects estimates with small sample sizes. Numerous studies have provided illustrations and/or evaluations of these small-sample corrections for multilevel models in psychological and educational research.Reference McNeish and Stapleton ^{35 ,} Reference McNeish and Harring ^{36 ,} Reference Stroup ^{55 –} Reference Ferron, Bell, Hess, Rendina-Gobioff and Hibbard ⁵⁷ However, researchers in the social sciences are not familiar with small-sample corrections for the GEE approach, which is an important issue to consider in the calculation of the standard error of BC-IRR.

As mentioned in the review of performance, the sandwich estimator for the standard error in GEE is downwardly biased when the sample size is small. In our study, we chose four corrections that could potentially improve small-sample performance and are able to account for unbalanced designs (e.g., different numbers of measurements across cases) based on previous studies. To understand these corrections, the sandwich estimator of the standard error of $\boldsymbol{\beta}$ based on equation (7) is provided:

(11) $$\begin{align}{\boldsymbol{\Phi}}_{LZ}={\left(\sum \limits_{j=1}^J{\boldsymbol{D}}_j^{\prime }{\boldsymbol{V}}_j^{-1}{\boldsymbol{D}}_j\right)}^{-1}\boldsymbol{M}\;{\left(\sum \limits_{j=1}^J{\boldsymbol{D}}_j^{\prime }{\boldsymbol{V}}_j^{-1}{\boldsymbol{D}}_j\right)}^{-1}\end{align}$$

where $\boldsymbol{M}$ is the “meat” of the “sandwich”: $\boldsymbol{M}=\sum_{j=1}^J{\boldsymbol{D}}_j^{\prime }{\boldsymbol{V}}_j^{-1}{\widehat{\boldsymbol{r}}}_j{\widehat{\boldsymbol{r}}}_j^{\prime }{\boldsymbol{V}}_j^{-1}{\boldsymbol{D}}_j$ , with ${\boldsymbol{D}}_j={\boldsymbol{X}}_j^{\prime }{\boldsymbol{A}}_j$ , ${\boldsymbol{V}}_j=\phi {\boldsymbol{A}}_j^{\frac{1}{2}}{\boldsymbol{R}}_{\! j}\left(\boldsymbol{\alpha} \right){\boldsymbol{A}}_j^{\frac{1}{2}}$ , ${\widehat{\boldsymbol{r}}}_j{\widehat{\boldsymbol{r}}}_j^{\prime }=\left({\boldsymbol{Y}}_j-{\widehat{\boldsymbol{\mu}}}_j\right){\left({\boldsymbol{Y}}_j-{\widehat{\boldsymbol{\mu}}}_j\right)}^{\prime }$ . When sample size is small, ${\widehat{\boldsymbol{r}}}_j{\widehat{\boldsymbol{r}}}_j^{\prime }$ is downwardly biased and leads to underestimated standard error of $\boldsymbol{\beta}$ .

Two classes of corrected sandwich estimators have been proposed to deal with small samples: residual-based correction and design-based corrections.Reference McNeish and Stapleton ³⁵ Residual-based corrections account for the small-sample bias by adding matrixes to the innermost part of $\boldsymbol{M}$ in the sandwich estimators, which could inflate the standard errors. For example, ${\boldsymbol{\Phi}}_{KC}$ is a bias-corrected sandwich estimator proposed by Kauermann and Carroll,Reference Kauermann and Carroll ⁴⁴ which rewrites the estimator from equation 11 such that

(12) $$\begin{align}{\boldsymbol{\Phi}}_{KC}&={\left(\sum \limits_{j=1}^J{\boldsymbol{D}}_j^{\prime }{\boldsymbol{V}}_j^{-1}{\boldsymbol{D}}_j\right)}^{-1}{\boldsymbol{M}}_{KC}{\left(\sum \limits_{j=1}^J{\boldsymbol{D}}_j^{\prime }{\boldsymbol{V}}_j^{-1}{\boldsymbol{D}}_j\right)}^{-1}\nonumber\\{\boldsymbol{M}}_{KC}&=\sum \limits_{j=1}^J{\boldsymbol{D}}_j^{\prime }{\boldsymbol{V}}_j^{-1}{\left(\boldsymbol{I}-{\boldsymbol{H}}_j\right)}^{-\frac{1}{2}}{\widehat{\boldsymbol{r}}}_j{\widehat{\boldsymbol{r}}}_j^{\prime }{\left(\boldsymbol{I}-{\boldsymbol{H}}_j^{\prime}\right)}^{-\frac{1}{2}}{\boldsymbol{V}}_j^{-1}{\boldsymbol{D}}_j,\end{align}$$

where ${\boldsymbol{H}}_j={\boldsymbol{D}}_j{\left(\sum_{j=1}^J{\boldsymbol{D}}_j^{\prime }{\boldsymbol{V}}_j^{-1}{\boldsymbol{D}}_j\right)}^{-1}{\boldsymbol{D}}_j^{\prime }{\boldsymbol{V}}_j^{-1}$ .

Other commonly used residual-based corrections include ${\boldsymbol{\Phi}}_{MD}$ proposed by Mancl and DerouedReference Mancl and DeRouen ⁴⁶ :

(13) $$\begin{align}{\boldsymbol{\Phi}}_{MD}&={\left(\sum \limits_{j=1}^J{\boldsymbol{D}}_j^{\prime }{\boldsymbol{V}}_j^{-1}{\boldsymbol{D}}_j\right)}^{-1}{\boldsymbol{M}}_{\! MD}{\left(\sum \limits_{j=1}^J{\boldsymbol{D}}_j^{\prime }{\boldsymbol{V}}_j^{-1}{\boldsymbol{D}}_j\right)}^{-1}\nonumber\\{\boldsymbol{M}}_{\! MD}&=\sum \limits_{j=1}^J{\boldsymbol{D}}_j^{\prime }{\boldsymbol{V}}_j^{-1}{\left(\boldsymbol{I}-{\boldsymbol{H}}_j\right)}^{-1}{\widehat{\boldsymbol{r}}}_j{\widehat{\boldsymbol{r}}}_j^{\prime }{\left(\boldsymbol{I}-{\boldsymbol{H}}_j^{\prime}\right)}^{-1}{\boldsymbol{V}}_j^{-1}{\boldsymbol{D}}_j,\end{align}$$

and ${\boldsymbol{\Phi}}_{FG}$ proposed by Fay and GraubardReference Fay and Graubard ⁴² :

(14) $$\begin{align}{\boldsymbol{\Phi}}_{FG}&={\left(\sum \limits_{j=1}^J{\boldsymbol{D}}_j^{\prime }{\boldsymbol{V}}_j^{-1}{\boldsymbol{D}}_j\right)}^{-1}{\boldsymbol{M}}_{\! FG}{\left(\sum \limits_{j=1}^J{\boldsymbol{D}}_j^{\prime }{\boldsymbol{V}}_j^{-1}{\boldsymbol{D}}_j\right)}^{-1}\nonumber\\{\boldsymbol{M}}_{\! FG}&=\sum \limits_{j=1}^J Diag\left\{{\left(1-\mathit{\min}\left\{b,{\boldsymbol{Q}}_{jj}\right\}\right)}^{-1/2}\right\}{\boldsymbol{D}}_j^{\prime }{\boldsymbol{V}}_j^{-1}{\left(\boldsymbol{I}-{\boldsymbol{H}}_j\right)}^{-\mathbf{1}}{\widehat{\boldsymbol{r}}}_j{\widehat{\boldsymbol{r}}}_j^{\prime }{\left(\boldsymbol{I}-{\boldsymbol{H}}_j^{\prime}\right)}^{-1}\nonumber\\&\quad\times {\boldsymbol{V}}_j^{-1}{\boldsymbol{D}}_j Diag\left\{{\left(1-\mathit{\min}\left\{b,{\boldsymbol{Q}}_{jj}\right\}\right)}^{-1/2}\right\},\end{align}$$

where $\boldsymbol{Q}={\boldsymbol{D}}_j^{\prime }{\boldsymbol{V}}_j^{-1}{\boldsymbol{D}}_j{\left(\sum_{j=1}^J{\boldsymbol{D}}_j^{\prime }{\boldsymbol{V}}_j^{-1}{\boldsymbol{D}}_j\right)}^{-1}$ and $b$ is an upper bound for the correction; by default, it has a value of 0.75.

On the other hand, design-based corrections take a different form and include additional additive terms to the sandwich estimator rather than modifying the “meat” of the “sandwich.” The Morel-Bokossa-Neerchal correction ( ${\boldsymbol{\Phi}}_{MBN}$ )Reference Morel and Bokossa ⁴⁷ is a commonly adopted design-based correction in applied studies, which comprises the sandwich estimator and an additive term:

(15) $$\begin{align}\begin{aligned}{\boldsymbol{\Phi}}_{MBN}&={\boldsymbol{\Phi}}_{LZ}+\delta \phi {\left(\sum \limits_{j=1}^J{\boldsymbol{D}}_j^{\prime }{\boldsymbol{V}}_j^{-1}{\boldsymbol{D}}_j\right)}^{-1}\\\delta &=\left\{\begin{array}{c}\!\! P/ \! \left(J-P\right)\;\mathrm{if}\;J>\left(d+1\right)P\\ {}1/d\;\mathrm{otherwise}\end{array}\right.,\end{aligned}\end{align}$$

where $P$ is the number of predictors, $J$ equals the number of clusters, and $\phi =\max \left(r, trace\left({\left(\sum_{j=1}^J{\boldsymbol{D}}_j^{\prime }{\boldsymbol{V}}_j^{-1}{\boldsymbol{D}}_j\right)}^{-1}\boldsymbol{M}\right)/P\right)$ . The default values for $d$ and $r$ are 2 and 1, respectively. A summary of the sandwich estimator and discussed corrections is provided in Table 2.

Table 2.

Small-sample corrections to the GEE sandwich variance estimator

Table 2. Long description

The table contains four columns: Category, Estimator, Expression, and Key modification.

* Sandwich Category:

- Estimator: L Z (Liang and Zeger).

- Expression: Phi sub L Z = B super minus 1 M B super minus 1.

- Key modification: M = summation from j = 1 to J of D sub j prime V sub j super minus 1 r-hat sub j r-hat sub j prime V sub j super minus 1 D sub j.

* Residual-based Category:

- Estimator: K C (Kauermann and Carroll).

- Expression: Phi sub K C = B super minus 1 M sub K C B super minus 1.

- Key modification: M sub K C = summation from j = 1 to J of D sub j prime V sub j super minus 1 open parenthesis I minus H sub j close parenthesis super minus 1/2 r-hat sub j r-hat sub j prime open parenthesis I minus H sub j prime close parenthesis super minus 1/2 V sub j super minus 1 D sub j.

- Estimator: M D (Mancl and DeRouen).

- Expression: Phi sub M D = B super minus 1 M sub M D B super minus 1.

- Key modification: M sub M D = summation from j = 1 to J of D sub j prime V sub j super minus 1 open parenthesis I minus H sub j close parenthesis super minus 1 r-hat sub j r-hat sub j prime open parenthesis I minus H sub j prime close parenthesis super minus 1 V sub j super minus 1 D sub j.

- Estimator: F G (Fay and Graubard).

- Expression: Phi sub F G = B super minus 1 M sub F G B super minus 1.

- Key modification: M sub F G = summation from j = 1 to J of Diag open brace open parenthesis 1 minus min open brace b, Q sub j j close brace close parenthesis super minus 1/2 close brace D sub j prime V sub j super minus 1 open parenthesis I minus H sub j close parenthesis super minus 1 r-hat sub j r-hat sub j prime open parenthesis I minus H sub j prime close parenthesis super minus 1 V sub j super minus 1 D sub j Diag open brace open parenthesis 1 minus min open brace b, Q sub j j close brace close parenthesis super minus 1/2 close brace.

* Design-based Category:

- Estimator: M B N (Morel et al).

- Expression: Phi sub M B N = Phi sub L Z plus delta phi open parenthesis summation from j = 1 to J of D sub j prime V sub j super minus 1 D sub j close parenthesis super minus 1.

- Key modification: delta = P all over open parenthesis J minus P close parenthesis if J is greater than open parenthesis d plus 1 close parenthesis P; otherwise 1 all over d.

3.1. Data generation

Data were simulated using a NB GLMM as shown in equation 6. We varied the following design factors, including series length ( $I$ ), number of cases ( $J$ ), immediate effect ( ${\gamma}_{20}$ ), trend change ( ${\gamma}_{30}$ ), variance components ( ${\sigma}_{u0}^2$ , ${\sigma}_{u1}^2$ , ${\sigma}_{u2}^2$ , ${\sigma}_{u3}^2$ ), dispersion ( $\theta$ ), and the time elapsed between intervention and measurement sessions ( $B-A$ ). The chosen conditions and associated rationale are similar to the previous simulation studies with SCED count data.Reference Li, Luo and Baek ^{13 ,} Reference Declercq, Jamshidi and Fernández-Castilla ³⁸ Specifically, the number of cases was 4 or 8, and the series length was 10 or 20. This represented a common setting based on reviews of single-case studies.Reference Pustejovsky, Swan and English ^{8 ,} Reference Shadish and Sullivan ⁵⁸ Following the MBD, the start of the treatment was staggered across cases: $I$ = 10 (starting points of the intervention: 3, 4, 6, 7) and $I$ = 20 (starting points of the intervention: 6, 8, 12, 14). These intervention starting points were selected to reflect the range of baseline lengths commonly observed in published SCEDs,Reference Shadish and Sullivan ⁵⁸ including shorter baseline phases that may not fully satisfy current phase-length recommendations according to the What Works Clearinghouse (WWC). The trend in the baseline phase ${\gamma}_{10}$ = log (1.00), and thus there is no trend in the baseline phase. The immediate effect ${\gamma}_{20}$ was set to log (1.00), log (2.00), log (6.00), or log (12.00), representing the zero, small, medium, and large treatment conditions, respectively. Although the effect size benchmark has not been established for immediate effect, the chosen conditions from zero to large effects are consistent with the previous simulation study with SCED count outcomes.Reference Li, Luo, Baek, Thompson and Lam ^{12 ,} Reference Declercq, Jamshidi and Fernández-Castilla ^{38 ,} Reference Li and Luo ⁵⁹ For trend change ${\gamma}_{30}$ , the zero, medium, and large effects were set to log(1.00), log(1.05), and log(1.10), which leads to the frequency count increasing by 0%, 5%, and 10% for each measurement in the treatment phase, respectively. The values chosen for trend change are based on reanalysis of meta-analyses of SCEDs using GLMMs.Reference Verschuur, Didden, Lang, Sigafoos and Huskens ^{60 –} Reference Dart, Collins, Klingbeil and McKinley ⁶²

The between-case variance components for the intercept ( ${\sigma}_{u0}^2$ ), baseline trend ( ${\sigma}_{u1}^2$ ), immediate effect ( ${\sigma}_{u2}^2$ ), and trend change ( ${\sigma}_{u3}^2$ ) were paired as [0.1, 0, 0.1,0], [0.3, 0.0003, 0.3, 0.0003], or [0.5, 0.0005, 0.5, 0.0005]. Following the same practice in Moeyaert et al.Reference Moeyaert, Ugille, Ferron, Beretvas and Van den Noortgate ⁶³ the correlations between all random effects were set as −0.5. The 1:1000 ratio of ${\sigma}_{u1}^2$ and ${\sigma}_{u3}^2$ versus ${\sigma}_{u0}^2$ and ${\sigma}_{u2}^2$ was also determined based on the reanalysis of SCED count data in previous meta-analyses.Reference Verschuur, Didden, Lang, Sigafoos and Huskens ^{60 –} Reference Dart, Collins, Klingbeil and McKinley ⁶² Table 3 summarizes all design conditions. The full factorial design resulted in 864 conditions, each with 2,000 independent replications.

Table 3.

Summary of factors and conditions for simulated count data

Table 3. Long description

The table consists of two columns: Parameter and Value.

* Series length (I): 10 (starting points of the intervention: 3, 4, 6, 7) or 20 (starting points of the intervention: 6, 8, 12, 14).

* Number of cases (J): 4 or 8.

* Intercept (gamma sub 00): log 2.

* Baseline trend (gamma sub 10): log 1.00.

* Immediate effect (gamma sub 20): log 1.00, log 2.00, log 6.00, or log 12.00.

* Trend change (gamma sub 30): log 1.00, log 1.05, or log 1.10.

* Between-case variance for intercept (sigma sub u0 super 2), baseline trend (sigma sub u1 super 2), immediate effect (sigma sub u2 super 2), and trend change (sigma sub u3 super 2): [0.1, 0, 0.1, 0], [0.3, 0.0003, 0.3, 0.0003], or [0.5, 0.0005, 0.5, 0.0005].

* Correlations between all random effects: minus 0.5.

* Dispersion parameter (theta): 5 or 10.

* Time elapsed (B minus A): 0, 1, or 5.

To analyze the simulated data, we fitted the NB GLMMs by using the R package glmmTMB Reference Brooks, Kristensen and van ⁶⁴ and the GEE models with independent error structures via the R package gee.Reference Carey, Lumley, Moler and Ripley ⁶⁵ In GEE models, we applied four small-sample corrections in addition to the sandwich estimator, GEE_LZ ( ${\boldsymbol{\Phi}}_{LZ}$ ), to get standard errors of coefficients, including GEE_KC ( ${\boldsymbol{\Phi}}_{KC}$ ), GEE_MD ( ${\boldsymbol{\Phi}}_{MD}$ ), GEE_FG ( ${\boldsymbol{\Phi}}_{FG}$ ), and GEE_MBN ( ${\boldsymbol{\Phi}}_{MBN}$ ). All corrections were implemented using the R package geesmv.Reference Wang ⁶⁶ For GLMMs, we allowed all variance components and only one covariance between random intercept and immediate effect to be freely estimated to avoid convergence issues. Hence, the fitted GLMMs were under-specified regarding the level-2 variance–covariance structure, which allows us to examine the robustness of GLMMs to misspecification. For both GLMM and GEE, each model was estimated using two reference time points (C), namely, the first and last sessions of the baseline phase, whereas the data were consistently generated with C at the last session of the baseline phase. This allows us to examine whether BC-IRR is robust to different coding strategies adopted for ${Time}_{ij}$ . Standard errors of BC-IRRs for both models were computed via the delta method using the R package msm.Reference Jackson ³⁷ The simulation code is available at https://osf.io/eqa4h/?view_only=98358eeef01f4d20b043d902e8d16c0a.

3.2. Performance measure

For each simulated condition, estimator performance was evaluated with four metrics. Bias of the BC-IRR was computed as the mean difference between the estimated and true BC-IRR across the 2,000 replications, as well as the mean-squared error (MSE), which is the sum of the estimator’s variance and squared bias. We calculated the bias of the model-based standard error (SE) from the GEE approach as the difference between the mean reported SE and the empirical standard deviation of BC-IRR estimates. Finally, we evaluated the 95% confidence-interval coverage for the small-sample corrections found effective in Section 2, defined as the proportion of replications whose Wald interval for BC-IRR contains the true value. Sampling variability in coverage was calculated with exact binomial error bounds: $.95\pm 1.96\sqrt{\frac{0.95\left(1-0.95\right)}{2000}}=\left(.94,.96\right)$ , so deviations outside this range showed under- or over-coverage beyond Monte Carlo error. To quantify simulation uncertainty, we calculated Monte Carlo standard errors (MCSEs) for each performance measure following the framework outlined by Morris et al.Reference Morris, White and Crowther ⁶⁷

3.3. Results

To facilitate comparison of the two estimation models, we summarized the differences between the GLMM and GEE across the simulation conditions. Overall, both methods produced negligible bias. Across all design conditions, the median bias was −0.0005 for the GLMM and 0.0004 for the GEE, whereas the mean bias was 0.001 for both methods. Bias ranged from −0.093 to 0.167 for the GLMM and from −0.137 to 0.137 for the GEE. Furthermore, 83.6% of the GLMM conditions and 92.2% of the GEE conditions produced absolute bias smaller than 0.05, indicating that differences in point estimate bias between the two approaches were negligible. The smallest differences in bias occurred when the time elapsed between treatment introduction and outcome measurement was short (particularly when $B-A=0$ ) and the between-case intercept variance was small ( ${\sigma}_{u0}^2=0.10$ ). Larger differences emerged as the follow-up duration increased and the between-case intercept variance became moderate or large ( ${\sigma}_{u0}^2=0.30$ or $0.50$ ). Even under these conditions, however, the practical differences between the two estimators remained modest, with bias generally ranging from approximately −0.03 to 0.04 across the most influential design factors.

For estimation accuracy, the GLMM was consistently more efficient than the GEE. The median MSEs were 0.19 for the GLMM and 0.26 for the GEE, with corresponding interquartile ranges of 0.21 and 0.27, respectively. Differences in MSE were smallest under the conditions with longer series ( $I=20$ ), larger numbers of cases ( $J=8$ ), and shorter follow-up durations after treatment introduction. In contrast, the largest differences occurred under the conditions when the series length was short ( $I=10$ ), the number of cases was small ( $J=4$ ), and the follow-up duration was longest ( $B-A=5$ ). Under these conditions, both estimators exhibited increased MSE, although the GEE consistently showed lower efficiency than the GLMM.

In Figure 1, we compare the bias of the SE of the BC-IRR using six different methods. The GLMM, along with the KC and MD corrections for sandwich estimators in GEE, yielded the lowest bias. The SE bias ranged from −0.728 to 0.013 for the GLMM (Median = −0.035, Mean = −0.047, The interquartile range (IQR) = 0.03, MCSE = 0.008), from −0.444 to 0.412 for the KC correction (Median = −0.048, Mean = −0.035, IQR = 0.07, MCSE = 0.01), and from −0.044 to 0.467 for the MD correction (Median = 0.021, Mean = 0.045, IQR = 0.07, MCSE = 0.011). In contrast, the GEE_FG correction showed the most extreme bias, ranging from −0.054 to 12.100 (Median = 0.50, Mean = 0.86, MCSE = 0.023), with the largest variability across conditions (IQR = 0.87), reflecting systematic over-correction of the SE. The GEE_LZ (Median = −0.130, Mean = −0.161, MCSE = 0.01) and GEE_MBN (Median = 0.058, Mean = 0.079, MCSE = 0.01) methods showed moderate bias but were not stable across conditions, exhibiting wide ranges and IQRs (0.13 and 0.15, respectively).

Figure 1.

Boxplot of the standard error bias of BC-IRR.

Figure 1. Long description

The Y-axis is labeled Bias of the standard error with numerical increments from -0.5 to 1.0. The X-axis is labeled S E methods and lists six categories. From left to right, the data shows.

* G L M M. A very narrow box centered slightly below 0.0 with a long tail of outliers extending down to -0.5.

* G E E underscore L Z. A box positioned between -0.1 and -0.2 with a dense cluster of outliers extending below -0.5.

* G E E underscore F G. The largest box in the set, spanning from approximately 0.2 to 1.0 with a median line at 0.5.

* G E E underscore K C. A box centered just below 0.0 with symmetrical outlier whiskers extending both up to 0.4 and down to -0.4.

* G E E underscore M D. A small box positioned just above 0.0 with a vertical line of outliers reaching up to 0.5.

* G E E underscore M B N. A box centered near 0.1 with significant outlier distributions extending both up toward 0.9 and down toward -0.5.

As GLMM, GEE_KC, and GEE_MD yielded more accurate estimates and higher stability across all conditions for SEs than other estimators, we only calculated the coverage rate of BC-IRR for the three estimators as shown in Figure 2. The median coverage rates for BC-IRR under the GLMM and the MD-corrected SEs are nearly identical (.92 versus .91; MCSE = 0.006 and 0.005, respectively), with tight IQRs of .03 and .02, respectively. In comparison, the KC adjustment yields a slightly lower median coverage rate (.87; MCSE = 0.006) and wider IQR (.06).

Figure 2.

Boxplot of the coverage rate of BC-IRR.

Figure 2. Long description

The y-axis is labeled Coverage rate with values from 0.8 to 1.0. The x-axis is labeled S E methods with three categories.

* G L M M: The box is positioned between 0.9 and 0.94 with a median around 0.92. It has a long lower whisker extending to 0.85 and several outliers reaching down to 0.83.

* G E E underscore M D: The box is narrower and positioned similarly to G L M M with a median around 0.91. It has a shorter lower whisker and a few outliers between 0.85 and 0.88.

* G E E underscore K C: This box is much larger, spanning from 0.85 to 0.90 with a median around 0.87. It has a very long upper whisker reaching near 1.0 and a dense cluster of outliers extending down to 0.73.

Next, to evaluate the impact of design factors on performance measures, we conducted mixed two-way ANOVA in which the between-condition factors were series length $I$ , number of cases  $J$ , immediate-effect magnitude  ${\gamma}_{20}$ , trend change magnitude  ${\gamma}_{30}$ , between-case variance of the intercept  ${\sigma}_{00}^2$ , and dispersion  $\theta$ , while the within-condition factor included the estimation strategy (only in ANOVAs with bias and MSE), standard error estimators (only in ANOVAs with SE bias and coverage rate), time elapsed ( $B-A$ ) and reference time point ( $C$ ). Main effects and all two-way interactions were tested, with ${\eta}^2$ statistics reported. All effects that were reported in the Results section had statistical significance at the α = .05 level. To concentrate on the most consequential sources of variation, we used scree-type plots of the ordered ${\eta}^2$ values to find main effects and interaction terms that fall to the left of the “elbow” point, which is analogous to using the scree plots to determine the number of factors in exploratory factor analyses (Li, 2024).

The bias of BC-IRR was mainly impacted by the time elapsed ( $F\left(2,554\right)=2518.41$ , ${\eta}^2=.19$ ) and the interaction of the intercept variance and the time elapsed ( $F\left(4,554\right)=1344.96$ , ${\eta}^2=.20$ ). As depicted in Figure 3, when the between-case variance of the intercept was set at its smallest level ( ${\sigma}_{00}^2=0.1)$ , the bias remained at the lowest level across different conditions of time elapsed. For the medium and large intercept-variance levels ( ${\sigma}_{00}^2=0.3\;\mathrm{or}\;0.5$ ), the mean bias increases as time elapsed increases. Although the main and interaction effects were statistically significant, there were no meaningful differences, with the bias ranging from −0.03 to 0.04, indicating that the BC-IRR was practically unbiased throughout the design space.

Figure 3.

Bias of BC-IRR as a function of time elapsed and variance of the intercept.

Figure 3. Long description

The figure consists of three vertical panels arranged from left to right. All panels share a common y-axis labeled Mean bias ranging from -0.02 to 0.04 and a common x-axis labeled Time elapsed (B minus A) with discrete points at 0, 1, and 5.

Panel 1 (Left): Titled Variance of Intercept (sigma sub u 0 squared) = 0.1. The data shows a slight downward trend. At time 0, mean bias is approximately 0.005. It decreases slightly at time 1 and drops to approximately -0.005 at time 5.

Panel 2 (Center): Titled Variance of Intercept (sigma sub u 0 squared) = 0.3. The data shows an upward trend. At time 0, mean bias is approximately -0.012. It rises to approximately -0.005 at time 1 and increases sharply to approximately 0.025 at time 5.

Panel 3 (Right): Titled Variance of Intercept (sigma sub u 0 squared) = 0.5. The data shows a steep upward trend. At time 0, mean bias is at its lowest point of approximately -0.028. It rises to approximately -0.012 at time 1 and increases very sharply to its peak of approximately 0.042 at time 5.

The MSE of BC-IRR was mainly impacted by time elapsed ( $F\left(2,554\right)=3059.73$ , ${\eta}^2=.24$ ), series length ( $F\left(1,237\right)=4458.61$ , ${\eta}^2=.18$ ), and the number of cases ( $F\left(1,237\right)=3567.22$ , ${\eta}^2=.14$ ). As shown in Figure 4, MSE increases when the time elapsed increases from zero to five. Regarding series length, MSE tends to decrease as the series length increased from 10 to 20. Similarly, increasing the number of cases consistently lowers the MSE across all design conditions.

Figure 4.

MSE of BC-IRR as a function of time elapsed, number of cases, and series length.

Figure 4. Long description

The grid is organized into two columns and three rows. The y-axis for all panels is Mean squared error with increments of 0.3, 0.6, and 0.9. The x-axis for all panels is Number of cases (J) with values 4 and 8.

Column 1: Series length (I) = 10.

* Top panel (Time-elapsed (B minus A) = 0): A line starts at approximately 0.35 for 4 cases and decreases to 0.2 for 8 cases.

* Middle panel (Time-elapsed (B minus A) = 1): A line starts at approximately 0.4 for 4 cases and decreases to 0.2 for 8 cases.

* Bottom panel (Time-elapsed (B minus A) = 5): A steep line starts at approximately 1.1 for 4 cases and decreases to 0.5 for 8 cases.

Column 2: Series length (I) = 20.

* Top panel (Time-elapsed (B minus A) = 0): A line starts at approximately 0.25 for 4 cases and decreases to 0.15 for 8 cases.

* Middle panel (Time-elapsed (B minus A) = 1): A line starts at approximately 0.25 for 4 cases and decreases to 0.15 for 8 cases.

* Bottom panel (Time-elapsed (B minus A) = 5): A line starts at approximately 0.35 for 4 cases and decreases to 0.2 for 8 cases.

Overall, M S E decreases as the number of cases increases, and the error is significantly higher when the series length is shorter and the time-elapsed interval is larger.

As shown in Figure 5, the bias of the standard error for BC-IRR was impacted mostly by standard error methods ( $F\left(5,1385\right)=1127.53$ , ${\eta}^2=.30$ ) and its interaction with immediate effect ( $F\left(15,1385\right)=171.78$ , ${\eta}^2=.14$ ). The GEE_FG correction exhibited the largest positive bias, which deteriorated exponentially as the immediate effect increased from 0 to 2.48. GEE_MBN correction showed overestimation, while GEE_LZ showed underestimation across different settings. The GLMM estimator, as well as the KC and MD corrections for GEE’s SE, showed minimal bias and were relatively insensitive to the magnitude of the immediate effect.

Figure 5.

Bias of SE as a function of SE methods and immediate effect.

Figure 5. Long description

The figure consists of six panels arranged in a two-by-three grid. The y-axis for all panels is Bias of the standard error, ranging from 0.0 to 1.5. The x-axis is Immediate effect gamma sub 2 0, with discrete values at 0, 0.69, 1.79, and 2.48.

Top row from left to right:

* G L M M: Shows a flat horizontal line near 0.0 across all x-axis values.

* G E E underscore L Z: Shows a flat horizontal line slightly below 0.0.

* G E E underscore F G: Shows a sharp exponential increase, starting near 0.2 at x equals 0 and rising to approximately 1.8 at x equals 2.48.

Bottom row from left to right:

* G E E underscore K C: Shows a mostly flat line near 0.0 with a slight upward tick at the final data point 2.48.

* G E E underscore M D: Shows a flat horizontal line near 0.0.

* G E E underscore M B N: Shows a very slight linear increase, starting near 0.0 and ending just above 0.1 at x equals 2.48.

From the preceding analysis of SE bias, we only include the standard errors based on GLMM, as well as KC and MD corrections for GEE, all of which showed minimal bias, for calculating the coverage rate. The SE methods ( $F\left(2,554\right)=6905.16$ , ${\eta}^2=.34$ ) and the number of cases ( $F\left(1,237\right)=3750.37$ , ${\eta}^2=.19$ ) were largely associated with coverage rate. In Figure 6, with four cases, the GLMM’s SE and MD-based intervals had coverage rates close to the nominal .95 band (mean = .91 and .90, respectively), whereas KC-based interval showed under-coverage (mean = .85). Increasing the number of cases to eight improved coverage rates for each method, with GLMM to .93, MD to .92, and KC to .90. Overall, all methods showed under-coverage across conditions.

Figure 6.

Coverage rate of BC-IRR as a function of se methods and number of cases.

Figure 6. Long description

A multi-panel figure with three vertical panels. The y-axis is shared across all panels, labeled Coverage rate, with percentage increments from 86% to 92%. The x-axis for each panel is labeled Number of cases (J) with two discrete points: 4 and 8.

* The first panel, G L M M, shows a linear increase from approximately 91% at J equals 4 to nearly 93% at J equals 8.

* The second panel, G E E underscore K C, shows the steepest linear increase, starting at approximately 84.5% at J equals 4 and rising to just below 90% at J equals 8.

* The third panel, G E E underscore M D, shows a linear increase from approximately 90% at J equals 4 to approximately 92% at J equals 8.

All three methods demonstrate a positive correlation between the number of cases and the coverage rate.

4.1. Data

Ota and DuPaulReference Ota and DuPaul ⁶⁸ evaluated a computer-assisted mathematics program with a game format for three fourth- to sixth-grade students who have attention-deficit hyperactivity disorder (ADHD). An MBD across participants was used: each student completed a series of baseline typical mathematics instruction sessions followed by an intervention phase in which students used 20-minute computer sessions three to four times per week. For every session, the investigators recorded two curriculum-based measures, including Digits Correct per Minute (DCPM) and Problems Correct per Minute (NPCPM), For illustrative purposes, we used the DCPM outcome only, and raw data of the three participants were extracted from the graphs presented in the original paper, as shown in Figure 7.

Figure 7.

Performance in DCPM for three students with ADHD during baseline and a computer-assisted intervention, adapted from Ota and DePaul. Reference Ota and DuPaul ⁶⁸ .

Figure 7. Long description

The graph consists of three vertically stacked panels. The y-axis for all panels is labeled Digits Correct Per Minute D C P M with increments of 10 from 20 to 40. The x-axis is labeled Session with increments from 1 to 16. Vertical dashed lines indicate the transition from baseline to intervention.

* Top panel Matthew. Baseline sessions 1 to 5 show scores fluctuating between 15 and 23. Intervention begins at session 6. Scores immediately jump to 33 and remain relatively stable between 26 and 33 until session 15, where the score rises to 38.

* Middle panel David. Baseline sessions 1 to 9 show a variable trend between 20 and 31. Intervention begins at session 10. Scores start at 36, dip slightly to 32 at session 12, then peak at 45 at session 14 before dropping to 32 at session 15.

* Bottom panel Andrew. Baseline sessions 1 to 12 show a stable but low trend fluctuating between 18 and 24. Intervention begins at session 13. Scores show a modest increase, peaking at 26 at session 15 and ending at 23 at session 16.

4.2. Considerations for effect size selection

For the choice of design-comparable effect size in research synthesis, the BC-SMD is not ideal for this study. Within the primary study, math performance was assessed using two different session lengths (DCPM based on 2-minute sessions and NPCPM based on 5-minute sessions). Prior work has shown that SMDs are sensitive to observation session length, and scaling the scores to a per-minute metric does not fix the problem, where the standard deviation remains a function of session length (Pustejovsky, 2019). As a result, SMDs may not be strictly comparable even across outcomes within the same study.

Beyond the primary study, in a subsequent meta-analysis in which this study was included (DuPaul et al., 2012), outcomes across SCED studies encompassed a wide range of measurement types, including curriculum-based measures, direct observations, percentage of work completed, and counts of words written. Because these measures differ substantially in scale and variance, BC-SMD is not directly comparable across studies and may misrepresent relative intervention impacts. This measurement heterogeneity motivates the adoption of the alternative effect size, BC-IRR, which is less sensitive to between-outcome and between-study differences in outcome variability and observation procedures.

4.3. Model specification

The equations (5) and 8 define BC-IRR within fully specified GLMM and GEE models that account for the baseline trend, immediate effect, and trend change. However, analysts often fit parsimonious models in practice. A reduced specification may reflect substantive knowledge about real data patterns or be imposed by non-convergence issues for complicated models. For demonstration purposes, we illustrated the computation of BC-IRR using GLMM and GEE under three realistic scenarios with increasing complexities. For GEE models, we used two different SE methods, including GEE_KC and GEE_MD, which showed more favorable statistical properties in the simulation study. In Scenario 1, we assume that there is no baseline trend or trend change, while the immediate effect varies across individuals. The GLMM is specified as follows (Model 1):

(16) $$\begin{align}{Y}_{ij}&\sim Poisson\;\left({\lambda}_{ij},\theta \right)\; or\;{Y}_{ij}\sim Negative\kern0.17em Binomial\;\left({\lambda}_{ij},\theta \right)\nonumber\\\mathrm{Level}\;1&:\log \left({\lambda}_{ij}\right)={\beta}_{0j}+{\beta}_{1j}{Phase}_{ij}\nonumber\\\mathrm{Level}\;2&:\left\{\begin{array}{c}{\beta}_{0j}={\gamma}_{00}+{u}_{0j}\\ {}{\beta}_{1j}={\gamma}_{10}+{u}_{1j}\end{array}\right.\;\left[\begin{array}{c}{u}_{0j}\\ {}{u}_{1j}\end{array}\right]\sim MVN\left(\left[\begin{array}{c}0\\ {}0\end{array}\right],\left[\begin{array}{cc}{\sigma}_{u0}^2& {\sigma}_{u0u1}\\ {}{\sigma}_{u1u0}& {\sigma}_{u1}^2\end{array}\right]\right).\end{align}$$

Based on this GLMM, the formula for BC-IRR is

(17) $$\begin{align}\mathrm{BC}-\mathrm{IRR}={e}^{\gamma_1+\frac{1}{2}{\sigma}_{u1}^2+{\sigma}_{u0u1}}.\end{align}$$

The GEE is specified as follows, with ${\boldsymbol{R}}_{\! j}\left(\boldsymbol{\alpha} \right)$ being the independence working correlation structure (Model 2):

(18) $$\begin{align}E\! \left({y}_{ij}|{\boldsymbol{x}}_{\boldsymbol{ij}}\right)&={\lambda}_{ij},\log \left({\lambda}_{ij}\right)={\beta}_0+{\beta}_1{Phase}_{ij}.\nonumber\\Var\! \left({y}_{ij}|{\boldsymbol{x}}_{\boldsymbol{ij}}\right)&=\phi v\left({\lambda}_{ij}\right)\nonumber\\{\boldsymbol{V}}_{\boldsymbol{j}}&=\phi {\boldsymbol{A}}_j^{\frac{1}{2}}{\boldsymbol{R}}_{\! j} \! \left(\boldsymbol{\alpha} \right){\boldsymbol{A}}_j^{\frac{1}{2}}\end{align}$$

Based on the specified GEE, the formula for BC-IRR is

(19) $$\begin{align}\mathrm{BC}-\mathrm{IRR}={e}^{\beta_1},\end{align}$$

In Scenario 2, we assume that there is a fixed baseline trend, a fixed trend change, and varying immediate effects. The GLMM model is specified as follows (Model 3):

${Y}_{ij}\sim Poisson\;\left({\lambda}_{ij},\theta \right)\; or\;{Y}_{ij}\sim Negative\kern0.17em Binomial\;\left({\lambda}_{ij},\theta \right)$

$\mathrm{Level}\;1:\log \left({\lambda}_{ij}\right)={\beta}_{0j}+{\beta}_{1j}{Time}_{ij}+{\beta}_{2j}{Phase}_{ij}+{\beta}_{3j}{Time}_{ij}^{\prime }{Phase}_{ij}$

(20) $$\begin{align}\mathrm{Level}\;2:\left\{\begin{array}{l}{\beta}_{0j}={\gamma}_{00}+{u}_{0j}\\ {}{\beta}_{1j}={\gamma}_{10}\kern1.2em \\ {}{\beta}_{2j}={\gamma}_{20}+{u}_{2j}\\ {}{\beta}_{3j}={\gamma}_{30\kern1.56em }\kern0.24em \end{array}\right.\left[\begin{array}{c}{u}_{0j}\\ {}{u}_{2j}\end{array}\right]\sim MVN\left(\left[\begin{array}{c}0\\ {}0\end{array}\right],\left[\begin{array}{cc}{\sigma}_{u0}^2& {\sigma}_{u0u2}\\ {}{\sigma}_{u2u0}& {\sigma}_{u2}^2\end{array}\right]\right).\end{align}$$

Based on this GLMM, the formula for BC-IRR is simplified to

(21) $$\begin{align}\mathrm{BC}-\mathrm{IRR}={e}^{\gamma_2+\frac{1}{2}{\sigma}_{u2}^2+{\sigma}_{u0u2}+{\gamma}_3\left(B-\mathrm{A}\right)}.\end{align}$$

The GEE is specified as follows, with ${\boldsymbol{R}}_{\boldsymbol{j}}\left(\boldsymbol{\alpha} \right)$ being the independence working correlation structure (Model 4):

(22) $$\begin{align}E\! \left({y}_{ij}|{\boldsymbol{x}}_{\boldsymbol{ij}}\right)&={\lambda}_{ij},\log \left({\lambda}_{ij}\right)={\beta}_0+{\beta}_1{Time}_{ij}+{\beta}_2{Phase}_{ij}+{\beta}_3{Time}_{ij}^{\prime }{Phase}_{ij}.\nonumber\\Var\! \left({y}_{ij}|{\boldsymbol{x}}_{\boldsymbol{ij}}\right)&=\phi v\left({\lambda}_{ij}\right)\nonumber\\{\boldsymbol{V}}_{\boldsymbol{j}}&=\phi {\boldsymbol{A}}_j^{\frac{1}{2}}{\boldsymbol{R}}_{\! j}\left(\boldsymbol{\alpha} \right){\boldsymbol{A}}_j^{\frac{1}{2}}\end{align}$$

Based on the specified GEE, the formula for BC-IRR is

(23) $$\begin{align}\mathrm{BC}-\mathrm{IRR}={e}^{\left({\beta}_2+{\beta}_3\left(B-A\right)\right)},\end{align}$$

In Scenario 3, we assume that there are varying baseline trends, immediate effects, and trend changes, and only the covariance between the random intercept and the random immediate effect is modeled. The model is specified as follows (Model 5):

${Y}_{ij}\sim Poisson\;\left({\lambda}_{ij},\theta \right)\;\mathrm{or} \ {Y}_{ij}\sim Negative\kern0.17em Binomial\;\left({\lambda}_{ij},\theta \right)$

$\mathrm{Level}\;1:\log \left({\lambda}_{ij}\right)={\beta}_{0j}+{\beta}_{1j}{Time}_{ij}+{\beta}_{2j}{Phase}_{ij}+{\beta}_{3j}{Time}_{ij}^{\prime }{Phase}_{ij}$

(24) $$\begin{align}\mathrm{Level}\;2:\left\{\begin{array}{l}\!\! {\beta}_{0j}={\gamma}_{00}+{u}_{0j}\\ {}\!\! {\beta}_{1j}={\gamma}_{10}+{u}_{1j}\\ {}\begin{array}{l}\!\! {\beta}_{2j}={\gamma}_{20}+{u}_{2j}\\ {}\!\! {\beta}_{3j}={\gamma}_{30}+{u}_{3j}\end{array}\end{array}\;\left[\begin{array}{c}\begin{array}{c}{u}_{0j}\\ {}{u}_{1j}\end{array}\\ {}\begin{array}{c}{u}_{2j}\\ {}{u}_{3j}\end{array}\end{array}\right]\sim MVN\left(\left[\begin{array}{c}\begin{array}{c}0\\ {}0\end{array}\\ {}\begin{array}{c}0\\ {}0\end{array}\end{array}\right],\left[\begin{array}{cc}\begin{array}{cc}{\sigma}_{u0}^2& 0\\ {}0& {\sigma}_{u1}^2\end{array}& \begin{array}{cc}{\sigma}_{u0u2}& 0\\ {}0& 0\end{array}\\ {}\begin{array}{cc}{\sigma}_{u2u0}& 0\\ {}0& 0\end{array}& \begin{array}{cc}{\sigma}_{u2}^2& 0\\ {}0& {\sigma}_{u3}^2\end{array}\end{array}\right]\right)\right. .\end{align}$$

Based on this model, the formula for BC-IRR is simplified to

(25) $$\begin{align}\mathrm{BC}-\mathrm{IRR}={e}^{\gamma_2+\frac{1}{2}{\sigma}_{u2}^2+{\sigma}_{u0u2}+{\gamma}_3\left(B-\mathrm{A}\right)+\frac{1}{2}{\left(B-A\right)}^2{\sigma}_{u3}^2}.\end{align}$$

The GEE (Model 6) for Scenario 3 is the same as Scenario 2 (see equation 22) and the resulting BC-IRR formula is the same as well (see equation 23).

The coding of the variables used in these models is the same as that shown in Table 1. For the second and third scenarios where $B-A$ (i.e., time elapsed between the introduction of treatment and the outcome measurement occasion) should be identified for BC-IRR calculation, multiple durations were applied including $B-A=3$ (shortest treatment duration of this study) and $B-A=9$ (longest treatment duration of this study). Models were fitted and analyzed using glmmTMB package version 1.1.11,Reference Brooks, Kristensen and van ⁶⁴ gee package version 4.13.29,Reference Carey, Lumley, Moler and Ripley ⁶⁵ geesmv version 1.3,Reference Wang ⁶⁶ and saws version 0.9.7.0Reference Fay and Graubard ⁴² in R version 4.5.0. Standard errors of BC-IRRs were computed based on the delta method, using the msm package in R .Reference Jackson ³⁷

4.4. Convergence and test of overdispersion

To determine whether a NB model should be used in GLMM-based BC-IRR calculation, overdispersion was tested after fitting the Poisson model in each scenario for the GLMM. Then the Pearson’s Chi-squared test was conducted to test overdispersion. For Scenarios 1 and 2, the test results indicated that there was no significant overdispersion in the count outcome based on the Poisson model (Model 1: ${\chi}^2=19.98, df=41,p=.99$ , Model 3: ${\chi}^2=19.71, df=39,p=.99)$ . Thus, Poisson models were selected for Scenarios 1 and 2. For Scenario 3, the Poisson GLMM (Model 6) failed to converge due to a non-positive-definite Hessian matrix when estimating fixed-effect SEs. All GEE models with Poisson mean specification converged successfully across the three scenarios.

4.5. Results

Table 4 summarizes the parameter estimates from both the conditional GLMMs and the marginal GEEs, together with their corresponding BC-IRRs. For brevity, we do not interpret coefficients of the GLMM and GEE models. Our focus is on the calculation and interpretation of the final BC-IRR values. Detailed explanations of the GLMM parameters can be found in Li et al. (2025). The GEE coefficients are interpreted and used to obtain BC-IRR.

Table 4.

Parameter estimates and BC-IRR for models 1 to 6

Table 4. Long description

The table is organized into five columns: Model, Parameter, Estimate, S E, and 95% C I.

Scenario 1 includes:

* Model 1-G L M M: gamma sub 00 (3.11), gamma sub 10 (0.30), sigma squared sub u0 (0.006), sigma squared sub u1 (0.004), sigma sub u1 u2 (0.005), and B C - I R R (1.36).

* Model 2-G E E (K C-adjusted): beta sub 0 (3.11), beta sub 1 (0.33), and B C - I R R (1.39).

* Model 2-G E E (M D-adjusted): beta sub 0 (3.11), beta sub 1 (0.33), and B C - I R R (1.39).

Scenario 2 includes:

* Model 3-G L M M: gamma sub 00 (3.12), gamma sub 10 (0.002), gamma sub 20 (0.27), gamma sub 30 (0.008), sigma squared sub u0 (0.006), sigma squared sub u2 (0.003), sigma sub u0 u2 (0.005), B C - I R R (B minus A equals 3) (1.34), and B C - I R R (B minus A equals 9) (1.41).

* Model 4-G E E (K C-adjusted): beta sub 0 (3.14), beta sub 1 (0.007), beta sub 2 (0.27), beta sub 3 (0.002), B C - I R R (B minus A equals 3) (1.31), and B C - I R R (B minus A equals 9) (1.32).

* Model 4-G E E (M D-adjusted): beta sub 0 (3.14), beta sub 1 (0.007), beta sub 2 (0.27), beta sub 3 (0.002), B C - I R R (B minus A equals 3) (1.31), and B C - I R R (B minus A equals 9) (1.32).

Scenario 3 includes:

* Model 5-G L M M: Listed as Non-convergence.

* Model 6-G E E: Listed as having the same results as Model 4.