A convenience function for calculating design comparable effect sizes

In one call, 1) clean single-case design data for treatment reversal and multiple baseline designs, 2) fit a multi-level model using restricted maximum likelihood estimation, and 3) estimate a standardized mean difference effect size.

calc_BCSMD(
  design,
  case,
  phase,
  session,
  outcome,
  cluster = NULL,
  series = NULL,
  center = 0,
  round_session = TRUE,
  treatment_name = NULL,
  FE_base = 0,
  RE_base = 0,
  RE_base_2 = NULL,
  FE_trt = 0,
  RE_trt = NULL,
  RE_trt_2 = NULL,
  corStruct = "AR1",
  varStruct = "hom",
  lmeControl = nlme::lmeControl(msMaxIter = 50, apVar = FALSE, returnObject = TRUE),
  Bayesian = FALSE,
  prior = NULL,
  chains = 4,
  iter = 2000,
  warmup = iter/2,
  thin = 10,
  cores = 1,
  seed = NA,
  A = NULL,
  B = NULL,
  D = NULL,
  cover = 95,
  bound = 35,
  symmetric = TRUE,
  summary = FALSE,
  data = NULL,
  ...
)

Arguments

design: Character string to specify whether data comes from a treatment reversal ("TR"), multiple baseline across participants ("MBP"), replicated multiple baseline across behaviors ("RMBB"), or clustered multiple baseline across participants ("CMB").
case: vector of case indicators or name of a character or factor vector within data indicating unique cases.
phase: vector of treatment indicators or name of a character or factor vector within data indicating unique treatment phases.
session: vector of measurement occasions or name of numeric vector within data of measurement times.
outcome: vector of outcome data or name of numeric vector of outcome data within data.
cluster: (Optional) vector of cluster indicators or name of a character or factor vector within data indicating clusters.
series: (Optional) vector of series indicators or name of a character or factor vector within data indicating series.
center: Numeric value for the centering value for session. Default is 0.
round_session: Logical indicating whether to round session to the nearest integer. Defaults to TRUE.
treatment_name: (Optional) character string corresponding to the name of the treatment phase.
FE_base: Vector of integers specifying which fixed effect terms to include in the baseline phase. Setting FE_base = 0 includes only a level. Setting FE_base = c(0,1) includes a level and a linear time trend.
RE_base: Vector of integers specifying which random effect terms to include in the baseline phase. Setting RE_base = 0 includes only levels (i.e., random intercepts). Setting RE_base = c(0,1) includes random levels and random linear trends.
RE_base_2: Vector of integers specifying which random effect terms to include in the baseline phase for the cluster level in clustered multiple baseline design across participants or for the case level in replicated multiple baseline across behaviors. Setting RE_base_2 = 0 includes only levels (i.e., random intercepts). Setting RE_base_2 = c(0,1) includes random levels and random linear trends.
FE_trt: Vector of integers specifying which fixed effect terms to include in the treatment phase. Setting FE_trt = 0 includes only a change in level. Setting FE_trt = c(0,1) includes a change in level and a treatment-by-linear time trend.
RE_trt: Vector of integers specifying which random effect terms to include in the treatment phase. Setting RE_trt = 0 includes only random changes in level. Setting RE_trt = c(0,1) includes random changes in level and random treatment-by-linear time trends.
RE_trt_2: Vector of integers specifying which random effect terms to include in the treatment phase for the cluster level in clustered multiple baseline design across participants or for the case level in replicated multiple baseline across behaviors. Setting RE_trt_2 = 0 includes only random changes in level. Setting RE_trt_2 = c(0,1) includes random changes in level and random treatment-by-linear time trends.
corStruct: (Optional) character string indicating the correlation structure of session-level errors. Options are "AR1" (default), "MA1", or "IID".
varStruct: (Optional) character string indicating the heteroscedasticity structure of session-level errors. Options are "hom" (default) or "het", which allows for the session-level error variances to differ by phase.
lmeControl: Specify control values for lme fit. See ?lmeControl for more detail.
Bayesian: Logical indicating whether to apply Bayesian estimation methods. If FALSE (default), apply restricted maximum likelihood estimation. If TRUE, apply Bayesian estimation methods (i.e., Markov Chain Monte Carlo) to the parameters of the specified hierarchical linear model.
prior: (Optional) One or more brmsprior objects created by "set_prior" in the "brms" R package. Default is prior = NULL when flat priors will be applied for Bayesian estimation. Ignored if Bayesian = FALSE.
chains: Number of Markov chains to use for Bayesian estimation. Default is chains = 4. Ignored if Bayesian = FALSE.
iter: Number of iterations per chain to use for Bayesian estimation. Default is iter = 2000. Ignored if Bayesian = FALSE.
warmup: Number of burn-in iterations to use for Bayesian estimation. Default is iter/2. Ignored if Bayesian = FALSE.
thin: Thinning rate to use for Bayesian estimation. Default is thin = 10. Ignored if Bayesian = FALSE.
cores: Number of cores to use for Bayesian estimation. Default is cores = 1. Ignored if Bayesian = FALSE.
seed: Random number generation seed to use for Bayesian estimation. Default is seed = NA, which will cause the seed to be set randomly. Ignored if Bayesian = FALSE.
A: The time point immediately before the start of treatment in the hypothetical between-group design.
B: The time point at which outcomes are measured in the hypothetical between-group design.
D: Numerical indicating the treatment duration across cases. Note that B = A + D and it is not allowed to specify both B and D.
cover: Confidence level.
bound: Numerical tolerance for non-centrality parameter in qt.
symmetric: If TRUE (the default), use a symmetric confidence interval. If FALSE, use a non-central t approximation to obtain an asymmetric confidence interval.
summary: Logical indicating whether to return a data frame with effect size estimates and other information. If TRUE (default), return a data.frame containing the effect size estimate, standard error, confidence interval (or credible interval if Bayesian == TRUE), and other information. If FALSE, return a list with effect size estimate, degrees of freedom, and other information.
data: (Optional) dataset to use for analysis. Must be a data.frame.
...: further arguments. Not currently used.

Value

If summary == TRUE, a data frame containing the design-comparable effect size estimate, standard error, confidence interval (or credible interval if Bayesian == TRUE), and other information. If summary == FALSE, a list containing all elements of a `g_mlm()` object, plus the fitted `lme()` model.

Examples

data(Laski)

# Change-in-levels model with fixed treatment effect
calc_BCSMD(design = "MBP",
           case = case, phase = treatment,
           session = time, outcome = outcome,
           FE_base = 0, RE_base = 0, FE_trt = 0,
           data = Laski)
#> $g_AB
#> [1] 1.404887
#> 
#> $SE_g_AB
#> [1] 0.2863249
#> 
#> $CI_L
#> [1] 0.8046221
#> 
#> $CI_U
#> [1] 2.005152
#> 
#> $nu
#> [1] 18.5524
#> 
#> $phi
#> [1] 0.2527692
#> 
#> $var_param
#> [1] NA
#> 
#> $rho
#> case.var((Intercept)) 
#>             0.5606065 
#> 
#> $model
#> Linear mixed-effects model fit by REML
#>   Data: dat 
#>   Log-restricted-likelihood: -519.1424
#>   Fixed: outcome ~ 1 + trt 
#> (Intercept)         trt 
#>    39.07612    30.68366 
#> 
#> Random effects:
#>  Formula: ~1 | case
#>         (Intercept) Residual
#> StdDev:    15.68279 13.88421
#> 
#> Correlation Structure: AR(1)
#>  Formula: ~session | case 
#>  Parameter estimate(s):
#>       Phi 
#> 0.2527692 
#> Number of Observations: 128
#> Number of Groups: 8 
#> 
#> $converged
#> [1] TRUE
#> 
#> $CI_cover
#> [1] 95
#> 
#> $A
#> [1] 4
#> 
#> $B
#> [1] 13
#> 
#> attr(,"class")
#> [1] "bcsmd"

if (FALSE) { # \dontrun{
# Bayesian estimation
calc_BCSMD(design = "MBP",
           case = case, phase = treatment,
           session = time, outcome = outcome,
           FE_base = 0, RE_base = 0, FE_trt = 0,
           Bayesian = TRUE,
           data = Laski)
} # }

# Model with linear time trends in baseline and treatment phases,
# random baseline slopes, fixed treatment effects
calc_BCSMD(design = "MBP",
           case = case, phase = treatment,
           session = time, outcome = outcome, center = 4,
           FE_base = c(0,1), RE_base = c(0,1), FE_trt = c(0,1),
           data = Laski)
#> Warning: nlminb problem, convergence error code = 1
#>   message = iteration limit reached without convergence (10)
#> $g_AB
#> [1] 1.885295
#> 
#> $SE_g_AB
#> [1] 0.5846749
#> 
#> $CI_L
#> [1] 0.675189
#> 
#> $CI_U
#> [1] 3.0954
#> 
#> $nu
#> [1] 22.79118
#> 
#> $phi
#> [1] 0.1313209
#> 
#> $var_param
#> [1] NA
#> 
#> $rho
#> case.var((Intercept)) 
#>             0.7188766 
#> 
#> $model
#> Linear mixed-effects model fit by REML
#>   Data: dat 
#>   Log-restricted-likelihood: -513.0308
#>   Fixed: outcome ~ 1 + I(session^1) + trt + I(session_trt^1) 
#>      (Intercept)     I(session^1)              trt I(session_trt^1) 
#>       38.8505064       -0.1372739       30.9354157        0.3454063 
#> 
#> Random effects:
#>  Formula: ~1 + I(session^1) | case
#>  Structure: General positive-definite, Log-Cholesky parametrization
#>              StdDev     Corr  
#> (Intercept)  20.6729955 (Intr)
#> I(session^1)  0.9944003 -1    
#> Residual     12.9278045       
#> 
#> Correlation Structure: AR(1)
#>  Formula: ~session | case 
#>  Parameter estimate(s):
#>       Phi 
#> 0.1313209 
#> Number of Observations: 128
#> Number of Groups: 8 
#> 
#> $converged
#> <simpleWarning in lme.formula(fixed = fixed, random = random, correlation = cor_struct,     weights = var_struct, data = dat, control = lmeControl): nlminb problem, convergence error code = 1
#>   message = iteration limit reached without convergence (10)>
#> 
#> $CI_cover
#> [1] 95
#> 
#> $A
#> [1] 4
#> 
#> $B
#> [1] 13
#> 
#> attr(,"class")
#> [1] "bcsmd"

if (FALSE) { # \dontrun{
# Bayesian estimation
calc_BCSMD(design = "MBP",
           case = case, phase = treatment,
           session = time, outcome = outcome, center = 4,
           FE_base = c(0,1), RE_base = c(0,1), FE_trt = c(0,1),
           Bayesian = TRUE,
           data = Laski)
} # }

data(Anglesea)
calc_BCSMD(design = "TR",
           case = case, phase = condition,
           session = session, outcome = outcome,
           treatment_name = "treatment",
           FE_base = 0, RE_base = 0, FE_trt = 0,
           data = Anglesea)
#> $g_AB
#> [1] 1.490965
#> 
#> $SE_g_AB
#> [1] 0.9876571
#> 
#> $CI_L
#> [1] -1.470357
#> 
#> $CI_U
#> [1] 4.452286
#> 
#> $nu
#> [1] 3.359203
#> 
#> $phi
#> [1] 0.4978687
#> 
#> $var_param
#> [1] NA
#> 
#> $rho
#> case.var((Intercept)) 
#>             0.7260842 
#> 
#> $model
#> Linear mixed-effects model fit by REML
#>   Data: dat 
#>   Log-restricted-likelihood: -241.5317
#>   Fixed: outcome ~ 1 + trt 
#> (Intercept)         trt 
#>    60.76171    89.59506 
#> 
#> Random effects:
#>  Formula: ~1 | case
#>         (Intercept) Residual
#> StdDev:    38.85317  23.8639
#> 
#> Correlation Structure: AR(1)
#>  Formula: ~session | case 
#>  Parameter estimate(s):
#>       Phi 
#> 0.4978687 
#> Number of Observations: 55
#> Number of Groups: 3 
#> 
#> $converged
#> [1] TRUE
#> 
#> $CI_cover
#> [1] 95
#> 
#> $B
#> [1] NA
#> 
#> $A
#> [1] NA
#> 
#> attr(,"class")
#> [1] "bcsmd"


data(Thiemann2001)
calc_BCSMD(design = "RMBB",
           case = case, series = series, phase = treatment,
           session = time, outcome = outcome,
           FE_base = 0, RE_base = 0, RE_base_2 = 0, FE_trt = 0,
           data = Thiemann2001)
#> $g_AB
#> [1] 1.598701
#> 
#> $SE_g_AB
#> [1] 0.1803463
#> 
#> $CI_L
#> [1] 1.239027
#> 
#> $CI_U
#> [1] 1.958375
#> 
#> $nu
#> [1] 70.16912
#> 
#> $phi
#> [1] 0.19949
#> 
#> $var_param
#> [1] NA
#> 
#> $rho
#> [1] "Level2:0.322  Level3:0.109"
#> 
#> $model
#> Linear mixed-effects model fit by REML
#>   Data: dat 
#>   Log-restricted-likelihood: -1138.817
#>   Fixed: outcome ~ 1 + trt 
#> (Intercept)         trt 
#>    5.929056   10.795439 
#> 
#> Random effects:
#>  Formula: ~1 | case
#>         (Intercept)
#> StdDev:    2.203097
#> 
#>  Formula: ~1 | series %in% case
#>         (Intercept) Residual
#> StdDev:    3.088218 5.498575
#> 
#> Correlation Structure: AR(1)
#>  Formula: ~session | case/series 
#>  Parameter estimate(s):
#>     Phi 
#> 0.19949 
#> Number of Observations: 363
#> Number of Groups: 
#>             case series %in% case 
#>                5               14 
#> 
#> $converged
#> [1] TRUE
#> 
#> $CI_cover
#> [1] 95
#> 
#> $A
#> [1] 9
#> 
#> $B
#> [1] 15
#> 
#> attr(,"class")
#> [1] "bcsmd"


data(Bryant2018)
calc_BCSMD(design = "CMB",
           cluster = group, case = case, phase = treatment,
           session = session, outcome = outcome, center = 49,
           treatment_name = "treatment",
           FE_base = c(0,1), RE_base = 0, RE_base_2 = 0,
           FE_trt = c(0,1), RE_trt = NULL, RE_trt_2 = NULL,
           data = Bryant2018)
#> $g_AB
#> [1] -2.245586
#> 
#> $SE_g_AB
#> [1] 0.5876192
#> 
#> $CI_L
#> [1] -3.425588
#> 
#> $CI_U
#> [1] -1.065584
#> 
#> $nu
#> [1] 50.45912
#> 
#> $phi
#> [1] 0.6305068
#> 
#> $var_param
#> [1] NA
#> 
#> $rho
#> [1] "Level2:0.771  Level3:0"
#> 
#> $model
#> Linear mixed-effects model fit by REML
#>   Data: dat 
#>   Log-restricted-likelihood: -1964.54
#>   Fixed: outcome ~ 1 + I(session^1) + trt + I(session_trt^1) 
#>      (Intercept)     I(session^1)              trt I(session_trt^1) 
#>      133.3943765        1.4149293        0.2407699       -0.8734667 
#> 
#> Random effects:
#>  Formula: ~1 | cluster
#>         (Intercept)
#> StdDev: 0.009365135
#> 
#>  Formula: ~1 | case %in% cluster
#>         (Intercept) Residual
#> StdDev:    16.39448 8.930915
#> 
#> Correlation Structure: ARMA(1,0)
#>  Formula: ~session | cluster/case 
#>  Parameter estimate(s):
#>      Phi1 
#> 0.6305068 
#> Number of Observations: 536
#> Number of Groups: 
#>           cluster case %in% cluster 
#>                12                33 
#> 
#> $converged
#> [1] TRUE
#> 
#> $CI_cover
#> [1] 95
#> 
#> $A
#> [1] 5
#> 
#> $B
#> [1] 54
#> 
#> attr(,"class")
#> [1] "bcsmd"
 
if (FALSE) { # \dontrun{          
# Bayesian estimation
calc_BCSMD(design = "CMB",
           cluster = group, case = case, phase = treatment,
           session = session, outcome = outcome, center = 49,
           treatment_name = "treatment",
           FE_base = c(0,1), RE_base = c(0,1), RE_base_2 = 0,
           FE_trt = c(0,1), RE_trt = 1, RE_trt_2 = NULL,
           Bayesian = TRUE,
           data = Bryant2018)
} # }