The SingleCaseES package provides R functions for calculating basic, within-case effect size indices for single-case designs, including several non-overlap measures and parametric effect size measures, and for estimating the gradual effects model (Swan & Pustejovsky, 2018). Standard errors and confidence intervals are provided for the subset of effect sizes indices with known sampling distributions.

The package also includes two graphical user interfaces for interactive use (designed using Shiny), both of which are also available as web apps hosted through shinyapps.io:

In this vignette, we introduce the package’s primary functions for carrying out effect size calculations. We demonstrate how to use the functions for calculating individual effect sizes from single data series, how to use the calc_ES() function for calculating multiple effect sizes from a single data series, and how to use batch_calc_ES() for calculating one or multiple effect sizes from multiple data series. We also highlight some further options available for the parametric effect size functions.

To start, be sure to load the package:

library(SingleCaseES)

# 1 Individual effect size functions

The SingleCaseES package includes functions for calculating the major non-overlap measures that have been proposed for use with single-case designs, as well as several parametric effect size measures. The following non-overlap measures are available (function names are listed in parentheses):

• Percentage of non-overlapping data (PND)
• Percentage of all non-overlapping data (PAND)
• Robust improvement rate difference (IRD)
• Percentage exceeding the median (PEM)
• Non-overlap of all pairs (NAP)
• Tau non-overlap (Tau)
• Tau-U, which includes baseline trend adjustment (Tau_U)

The following parametric effect sizes are available:

• Within-case standardized mean difference (SMD)
• The increasing and decreasing versions of the log response ratio (LRRi and LRRd)
• Log odds ratio (LOR)

All of the functions for calculating individual effect sizes follow the same basic syntax. For demonstration purposes, let’s take a look at the syntax for NAP(), which calculates the non-overlap of all pairs (Parker & Vannest, 2009):

args(NAP)
#> function (A_data, B_data, condition, outcome, baseline_phase = unique(condition)[1],
#>     improvement = "increase", SE = "unbiased", confidence = 0.95)
#> NULL

All of the effect sizes functions in SingleCaseES can take the data from a single SCD series in one of two formats. We will first demonstrate each of these methods, then explain the further arguments to the function.

## 1.1 Using the A_data, B_data inputs

The first input format involves providing separate vectors for the data from each phase. Here is some hypothetical data from the A phase and B phase of a single-case data series:

A <- c(20, 20, 26, 25, 22, 23)
B <- c(28, 25, 24, 27, 30, 30, 29)

We can feed these data into the NAP function as follows:

NAP(A_data = A, B_data = B)
#>    ES       Est         SE  CI_lower  CI_upper
#> 1 NAP 0.9166667 0.06900656 0.5973406 0.9860176

## 1.2 Using the condition, outcome inputs

The second input format involves providing a single vector containing all of the outcome data from the series, along with a vector that describes the phase of each observation in the data. Here we re-format the hypothetical data above to follow this structure:

phase <- c(rep("A", 6), rep("B", 7))
phase
#>  [1] "A" "A" "A" "A" "A" "A" "B" "B" "B" "B" "B" "B" "B"

outcome_dat <- c(A, B)
outcome_dat
#>  [1] 20 20 26 25 22 23 28 25 24 27 30 30 29
NAP(condition = phase, outcome = outcome_dat)
#>    ES       Est         SE  CI_lower  CI_upper
#> 1 NAP 0.9166667 0.06900656 0.5973406 0.9860176

Note that if the vector provided to condition has more than two values, the effect size function will throw an error:

phase2 <- c(rep("A", 5), rep("B", 5), rep("C",3))
NAP(condition = phase2, outcome = outcome_dat)
#> Error in calc_ES(A_data = A_data, B_data = B_data, condition = condition, : The 'condition' variable must have exactly two unique values.

By default, NAP() and all of the other effect size functions treat the first value of condition as the the baseline phase. However, in some single-case data series, the initial observation might not be in the baseline phase. For example, an SCD with four cases might use a cross-over treatment reversal design, where two of the cases follow an ABAB design and the other two cases follow a BABA design. To handle this situation, we will need to specify the baseline phase using the baseline_phase argument:

phase_rev <- c(rep("B", 7), rep("A", 6))
outcome_rev <- c(B, A)
NAP(condition = phase_rev, outcome = outcome_rev, baseline_phase = "A")
#>    ES       Est         SE  CI_lower  CI_upper
#> 1 NAP 0.9166667 0.06900656 0.5973406 0.9860176

## 1.3 Direction of improvement

All of the effect size functions in SingleCaseES are defined based on some assumption about the direction of therapeutic improvement in the outcome (e.g., improvement would correspond to increases in on-task behavior but decreases in aggressive behavior). For all of the effect size functions, it is important to specify the direction of therapeutic improvement for the data series by providing a value for the improvement argument:

NAP(A_data = A, B_data = B, improvement = "decrease")
#>    ES        Est         SE   CI_lower  CI_upper
#> 1 NAP 0.08333333 0.06900656 0.01398242 0.4026594

Note that NAP() and most of the effect size functions default to assuming that increases in the outcome correspond to improvements.

## 1.4 Further options for NAP()

The NAP function provides several possible methods for calculating the standard error. By default, the exactly unbiased standard errors are used. However, the function can also produce standard errors using the Hanley-McNeil estimator, the variance under the null hypothesis of no effect, or no standard errors at all:

NAP(A_data = A, B_data = B, SE = "unbiased")
#>    ES       Est         SE  CI_lower  CI_upper
#> 1 NAP 0.9166667 0.06900656 0.5973406 0.9860176

NAP(A_data = A, B_data = B, SE = "Hanley")
#>    ES       Est         SE  CI_lower  CI_upper
#> 1 NAP 0.9166667 0.07739185 0.5973406 0.9860176

NAP(A_data = A, B_data = B, SE = "null")
#>    ES       Est        SE  CI_lower  CI_upper
#> 1 NAP 0.9166667 0.1666667 0.5973406 0.9860176

NAP(A_data = A, B_data = B, SE = "none")
#>    ES       Est  CI_lower  CI_upper
#> 1 NAP 0.9166667 0.5973406 0.9860176

The function also produces a confidence interval for NAP. By default, a 95% CI is calculated. To calculate an interval at some other level of coverage, set the confidence argument to a value between 0 and 1:

NAP(A_data = A, B_data = B)
#>    ES       Est         SE  CI_lower  CI_upper
#> 1 NAP 0.9166667 0.06900656 0.5973406 0.9860176

NAP(A_data = A, B_data = B, confidence = .99)
#>    ES       Est         SE  CI_lower  CI_upper
#> 1 NAP 0.9166667 0.06900656 0.4875014 0.9907377

NAP(A_data = A, B_data = B, confidence = .90)
#>    ES       Est         SE  CI_lower  CI_upper
#> 1 NAP 0.9166667 0.06900656 0.6591091 0.9822249

Set confidence = NULL to omit the confidence interval calculations all together:

NAP(A_data = A, B_data = B, confidence = NULL)
#>    ES       Est         SE
#> 1 NAP 0.9166667 0.06900656

## 1.5 Other non-overlap indices

The SingleCaseES package includes functions for calculating several other non-overlap indices in addition to NAP. All of the functions accept data in either the A, B format or the condition, outcome format with optional baseline specification, and all of the functions include an argument to specify the direction of improvement. Like the function for NAP, the function for Tau (Tau) can produce unbiased standard errors, Hanley-McNeil standard errors, standard errors under the null hypothesis of no effect, or no standard errors at all. Only NAP and Tau return standard errors and confidence intervals. The remaining non-overlap measures return only a point estimate:

Tau(A_data = A, B_data = B)
#>    ES       Est        SE  CI_lower  CI_upper
#> 1 Tau 0.8333333 0.1380131 0.1946812 0.9720352

PND(A_data = A, B_data = B)
#>    ES       Est
#> 1 PND 0.7142857

PEM(A_data = A, B_data = B)
#>    ES Est
#> 1 PEM   1

PAND(A_data = A, B_data = B)
#>     ES       Est
#> 1 PAND 0.8461538

IRD(A_data = A, B_data = B)
#>    ES       Est
#> 1 IRD 0.6904762

Tau_U(A_data = A, B_data = B)
#>      ES       Est
#> 1 Tau-U 0.7380952

## 1.6 Further options for SMD()

The standardized mean difference parameter is defined as the difference between the mean level of the outcome in phase B and the mean level of the outcome in phase A, scaled by the within-case standard deviation of the outcome in phase A. As with all functions discussed so far, the SMD() function accepts data in either the A_data, B_data format or the condition, outcome format with optional baseline phase specification.

The direction of improvement can be specified with the improvement argument, with “increase” being the default. Changing the direction of the improvement does not change the magnitude of the effect size, but does change its sign:

A <- c(20, 20, 26, 25, 22, 23)
B <- c(28, 25, 24, 27, 30, 30, 29)

SMD(A_data = A, B_data = B, improvement = "increase")
#>    ES      Est        SE  CI_lower CI_upper
#> 1 SMD 1.703734 0.6370139 0.4552097 2.952258

SMD(A_data = A, B_data = B, improvement = "decrease")
#>    ES       Est        SE  CI_lower   CI_upper
#> 1 SMD -1.703734 0.6370139 -2.952258 -0.4552097

The std_dev argument controls whether the effect size estimate is based on the standard deviation of the baseline phase alone (the default, std_dev = "baseline"), or based on the standard deviation after pooling across both phases (std_dev = "pool"):

SMD(A_data = A, B_data = B, std_dev = "baseline")
#>    ES      Est        SE  CI_lower CI_upper
#> 1 SMD 1.703734 0.6370139 0.4552097 2.952258
SMD(A_data = A, B_data = B, std_dev = "pool")
#>    ES      Est        SE  CI_lower CI_upper
#> 1 SMD 1.876247 0.6374216 0.6269241 3.125571

By default the SMD() function uses the Hedges’ g bias correction for small sample sizes. The bias correction can be turned off by specifying the argument bias_correct = FALSE. The width of the confidence intervals is controlled via the confidence argument, and no confidence intervals will be produced if the argument is set to confidence = NULL.

## 1.7 Log response ratios

The response ratio parameter is the ratio of the mean level of the outcome during phase B to the mean level of the outcome during phase A. The log response ratio is the natural logarithm of the response ratio. This effect size is appropriate for outcomes measured on a ratio scale, such that zero corresponds to the true absence of the outcome.

The package includes two versions of the LRR:

• LRR-increasing (LRRi()) is defined so that positive values correspond to therapeutic improvements
• LRR-decreasing (LRRd()) is defined so that negative values correspond to therapeutic improvements.

If you are estimating an effect size for a single series, pick the version of LRR that corresponds to the therapeutic improvement expected for your dependent variable. Similarly, if you are estimating effect sizes for a set of SCD series with the same therapeutic direction, pick the version that corresponds to your intervention’s expected change.

If you are estimating effect sizes for interventions where the direction of improvement depends upon the series or study, the choice between LRRi and LRRd is slightly more involved. For example, imagine we have ten studies to meta-analyze. For eight studies, the outcome are initiations of peer interaction, so therapeutic improvements correspond to increases in behavior. For the other two studies, the outcomes were episodes of verbal aggression towards peers, so the therapeutic direction was a decrease. In this context it would be sensible to pick the LRRi() function, because most of the outcomes are positively-valenced. For the final two studies, we would specify improvement = "decrease", which would ensure that the sign and magnitude of the outcomes were consistent with the direction of therapeutic improvement ( i.e. a larger log-ratio represents a larger change in the desired direction). Conversely, if most of the outcomes had a negative valence and only a few had a positive valence, then we would use LRRd() and we would specify improvement = "increase" for the few series that had positive-valence outcomes.

### 1.7.1 Setting the scale of the outcome

LRR differs from other effect size indices for single-case designs in that calculating it requires some further information about how the outcome variable was measured. One important piece of information to know is the scale of the outcome measurements. For outcomes that are measured by frequency counting, the scale might be expressed as a raw count (scale = "count") or as a standardized rate per minute (scale = "rate"). For outcomes that are measures of state behavior, where the main dimension of interest is the proportion of time that the behavior occurs, the scale might be expressed as a percentage (ranging from 0 to 100%; scale = "percentage") or as a proportion (ranging from 0 to 1; scale = "proportion"). For outcomes that don’t fit into any of these categories, set scale = "other".

The scale of the outcome variable has two important implications for how log response ratios are estimated. First, outcomes measured as percentages or proportions need to be coded so that the direction of therapeutic improvement is consistent with the direction of the effect size. Consequently, changing the improvement direction will alter the magnitude, in addition to the sign, of the effect size (see Pustejovsky, 2018, pp. 16–18 for further details). Here is an example:

A <- c(20, 20, 26, 25, 22, 23)
B <- c(28, 25, 24, 27, 30, 30, 29)
LRRi(A_data = A, B_data = B, scale = "percentage")
#>     ES       Est         SE   CI_lower  CI_upper
#> 1 LRRi 0.1953962 0.05557723 0.08646679 0.3043255

LRRi(A_data = A, B_data = B, improvement = "decrease", scale = "percentage")
#>     ES         Est         SE   CI_lower    CI_upper
#> 1 LRRi -0.06553504 0.01810144 -0.1010132 -0.03005687

Note that if the outcome is a count (the default for both LRR functions) or rate, changing the improvement direction merely changes the sign of the effect size:

A <- c(20, 20, 26, 25, 22, 23)
B <- c(28, 25, 24, 27, 30, 30, 29)
LRRi(A_data = A, B_data = B, scale = "count")
#>     ES       Est         SE   CI_lower  CI_upper
#> 1 LRRi 0.1953962 0.05557723 0.08646679 0.3043255
LRRi(A_data = A, B_data = B, scale = "count", improvement = "decrease")
#>     ES        Est         SE   CI_lower    CI_upper
#> 1 LRRi -0.1953962 0.05557723 -0.3043255 -0.08646679

The scale of the outcome has one further important implication. To account for the possibility of a sample mean of zero, the LRRd() and LRRi() functions use a truncated sample mean, where the truncation level is determined by the scale of the outcome and some further details of how the outcomes were measured. For rates, the truncated mean requires specifying the length of the observation session in minutes:

A <- c(0, 0, 0, 0)
B <- c(28, 25, 24, 27, 30, 30, 29)
LRRd(A_data = A, B_data = B, scale = "rate")
#>     ES Est  SE CI_lower CI_upper
#> 1 LRRd NaN NaN      NaN      NaN
LRRd(A_data = A, B_data = B, scale = "rate", observation_length = 30)
#>     ES      Est         SE CI_lower CI_upper
#> 1 LRRd 8.797947 0.03249549 8.734257 8.861637

If no additional information is provided and there is a sample mean of 0, the function returns a value of NaN.

For outcomes specified as percentages or proportions, the argument intervals must be supplied. For interval recording methods such as partial interval recording or momentary time sampling, provide the number of intervals. For continuous recording, set intervals equal to 60 times the length of the observation session in minutes:

LRRd(A_data = A, B_data = B, scale = "percentage")
#>     ES Est  SE CI_lower CI_upper
#> 1 LRRd NaN NaN      NaN      NaN
LRRd(A_data = A, B_data = B, scale = "percentage", intervals = 180)
#>     ES      Est         SE CI_lower CI_upper
#> 1 LRRd 5.984536 0.03249549 5.920846 6.048226

You can also specify your own value for the constant used to truncate the sample mean by supplying a value for D_const. If a vector, the mean will be used.

### 1.7.2 Additional arguments

Both LRR functions return a effect size that has been bias-corrected for small sample sizes by default. To omit the bias correction, set bias_correct = FALSE. Finally, as with the non-overlap measures, the confidence argument can be used to change the default 95% confidence interval, or set to NULL to omit confidence interval calculations.

## 1.8 Log-odds ratios

The odds ratio parameter is the ratio of the odds that the outcome occurs during phase B to the odds that the outcome occurs during phase A. The log-odds ratio (LOR) is the natural logarithm of the odds ratio. This effect size is appropriate for outcomes measured on a percentage or proportion scale. The LOR() function works almost identically to the LRRi() and LRRd() functions, but there are a few exceptions.

The LOR() function accepts only outcomes that are on proportion or percentage scales:

A_pct <- c(20, 20, 25, 25, 20, 25)
B_pct <- c(30, 25, 25, 25, 35, 30, 25)

LOR(A_data = A_pct, B_data = B_pct, scale = "percentage")
#>    ES       Est         SE   CI_lower  CI_upper
#> 1 LOR 0.2852854 0.09790282 0.09339935 0.4771713

LOR(A_data = A_pct/100, B_data = B_pct/100, scale = "proportion")
#>    ES       Est         SE   CI_lower  CI_upper
#> 1 LOR 0.2852854 0.09790282 0.09339935 0.4771713

LOR(A_data = A_pct, B_data = B_pct, scale = "count")
#> Warning: LOR can only be calculated for proportions or percentages. Will
#> return NAs for other outcome scales.
#>    ES Est SE CI_lower CI_upper
#> 1 LOR  NA NA       NA       NA

LOR(A_data = A_pct, B_data = B_pct, scale = "proportion")
#> Error in calc_LOR(A_data = c(20, 20, 25, 25, 20, 25), B_data = c(30, 25, : Proportions must be between 0 and 1!

As with the LRR functions, LOR() includes an argument to specify the direction of therapeutic improvement, with the default assumption being that a therapeutic improvement is an increase in the behavior. In contrast to LRRi and LRRd, changing the direction of therapeutic improvement only reverses the sign of the LOR, but does not change its absolute magnitude:

LOR(A_data = A_pct, B_data = B_pct,
scale = "percentage", improvement = "increase")
#>    ES       Est         SE   CI_lower  CI_upper
#> 1 LOR 0.2852854 0.09790282 0.09339935 0.4771713

LOR(A_data = A_pct, B_data = B_pct,
scale = "percentage", improvement = "decrease")
#>    ES        Est         SE   CI_lower    CI_upper
#> 1 LOR -0.2852854 0.09790282 -0.4771713 -0.09339935

Similar to the LRR functions, LOR() will be calculated using truncated sample means for cases where phase means are close to the extremes of the scale. To use truncated means, the number of intervals per observation session must be specified using the intervals argument:

LOR(A_data = c(0,0,0), B_data = B_pct,
scale = "percentage")
#>    ES Est  SE CI_lower CI_upper
#> 1 LOR NaN NaN      NaN      NaN
LOR(A_data = c(0,0,0), B_data = B_pct,
scale = "percentage", intervals = 20)
#>    ES      Est         SE CI_lower CI_upper
#> 1 LOR 3.828777 0.07398661 3.683766 3.973788

For data measured using continuous recording, set the number of intervals equal to 60 times the length of the observation session in minutes.

Like the LRR functions, it is possible to specify your own truncation constant using the D_const argument. By default the LOR() function uses a bias correction for small sample sizes, but this can be turned off by specifying the argument bias_correct = FALSE. The width of the confidence intervals is controlled via the confidence argument; set the argument to confidence = NULL to omit the confidence interval calculations.

# 2 calc_ES()

The calc_ES() function will calculate multiple effect sizes estimates for a single SCD series. As with the individual effect size functions, calc_ES() accepts data in either the A_data, B_data format or the condition, outcome format. To calculate multiple effect size estimates, provide a list of effect sizes to the ES argument. Here we use the A_data, B_data format:

A <- c(20, 20, 26, 25, 22, 23)
B <- c(28, 25, 24, 27, 30, 30, 29)
calc_ES(A_data = A, B_data = B, ES = c("NAP","PND","Tau-U"))
#>      ES       Est         SE  CI_lower  CI_upper
#> 1   NAP 0.9166667 0.06900656 0.5973406 0.9860176
#> 2   PND 0.7142857         NA        NA        NA
#> 3 Tau-U 0.7380952         NA        NA        NA

Here is the same calculation in the condition, outcome format:

phase <- c(rep("A", length(A)), rep("B", length(B)))
outcome <- c(A, B)
calc_ES(condition = phase, outcome = outcome, baseline_phase = "A",
ES = c("NAP","PND","Tau-U"))
#>      ES       Est         SE  CI_lower  CI_upper
#> 1   NAP 0.9166667 0.06900656 0.5973406 0.9860176
#> 2   PND 0.7142857         NA        NA        NA
#> 3 Tau-U 0.7380952         NA        NA        NA

The ES argument can include any of the following metrics: "LRRd", "LRRi", "LOR", "SMD", "NAP", "PND", "PEM", "PAND", "IRD", "Tau", or "Tau-U".

Setting ES = "all" will return all available effect sizes:

calc_ES(A_data = A, B_data = B, ES = "all")
#>       ES        Est         SE    CI_lower    CI_upper
#> 1   LRRd -0.1953962 0.05557723 -0.30432554 -0.08646679
#> 2   LRRi  0.1953962 0.05557723  0.08646679  0.30432554
#> 3    LOR  0.2609312 0.07356710  0.11674235  0.40512007
#> 4    SMD  1.7037340 0.63701390  0.45520971  2.95225831
#> 5    NAP  0.9166667 0.06900656  0.59734061  0.98601758
#> 6    IRD  0.6904762         NA          NA          NA
#> 7   PAND  0.8461538         NA          NA          NA
#> 8    PND  0.7142857         NA          NA          NA
#> 9    PEM  1.0000000         NA          NA          NA
#> 10   Tau  0.8333333 0.13801311  0.19468122  0.97203517
#> 11 Tau-U  0.7380952         NA          NA          NA

ES = "NOM" will return all of the non-overlap measures:

calc_ES(A_data = A, B_data = B, ES = "NOM")
#>      ES       Est         SE  CI_lower  CI_upper
#> 1   NAP 0.9166667 0.06900656 0.5973406 0.9860176
#> 2   IRD 0.6904762         NA        NA        NA
#> 3  PAND 0.8461538         NA        NA        NA
#> 4   PND 0.7142857         NA        NA        NA
#> 5   PEM 1.0000000         NA        NA        NA
#> 6   Tau 0.8333333 0.13801311 0.1946812 0.9720352
#> 7 Tau-U 0.7380952         NA        NA        NA

and ES = "parametric" will return all of the parametric effect sizes:

calc_ES(A_data = A, B_data = B, ES = "parametric")
#>     ES        Est         SE    CI_lower    CI_upper
#> 1 LRRd -0.1953962 0.05557723 -0.30432554 -0.08646679
#> 2 LRRi  0.1953962 0.05557723  0.08646679  0.30432554
#> 3  LOR  0.2609312 0.07356710  0.11674235  0.40512007
#> 4  SMD  1.7037340 0.63701390  0.45520971  2.95225831

If the ES argument is omitted, calc_ES() will return LRRd, LRRi, SMD, and Tau by default.

## 2.1 Further arguments

All of the individual effect size functions have the further argument improvement, and several of them also have further optional arguments. Include these arguments in calc_ES() in order to pass them on to the individual effect size calculation functions. For example, we can set the direction of improvement to decrease:

calc_ES(A_data = A, B_data = B, ES = "NOM", improvement = "decrease")
#>      ES         Est         SE    CI_lower   CI_upper
#> 1   NAP  0.08333333 0.06900656  0.01398242  0.4026594
#> 2   IRD  0.07142857         NA          NA         NA
#> 3  PAND  0.53846154         NA          NA         NA
#> 4   PND  0.00000000         NA          NA         NA
#> 5   PEM  0.00000000         NA          NA         NA
#> 6   Tau -0.83333333 0.13801311 -0.97203517 -0.1946812
#> 7 Tau-U -0.73809524         NA          NA         NA

To omit the confidence interval calculations for NAP and Tau, we can include the argument confidence = NULL:

calc_ES(A_data = A, B_data = B, ES = "NOM", improvement = "decrease", confidence = NULL)
#>      ES         Est         SE
#> 1   NAP  0.08333333 0.06900656
#> 2   IRD  0.07142857         NA
#> 3  PAND  0.53846154         NA
#> 4   PND  0.00000000         NA
#> 5   PEM  0.00000000         NA
#> 6   Tau -0.83333333 0.13801311
#> 7 Tau-U -0.73809524         NA

Details such as the measurement scale can also be passed on to functions that will make use of them:

calc_ES(A_data = A, B_data = B, ES = "parametric", scale = "count")
#> Warning: LOR can only be calculated for proportions or percentages. Will
#> return NAs for other outcome scales.
#>     ES        Est         SE    CI_lower    CI_upper
#> 1 LRRd -0.1953962 0.05557723 -0.30432554 -0.08646679
#> 2 LRRi  0.1953962 0.05557723  0.08646679  0.30432554
#> 3  LOR         NA         NA          NA          NA
#> 4  SMD  1.7037340 0.63701390  0.45520971  2.95225831

Any additional arguments will be used in the calculation of effect sizes for which they are relevant and will be ignored if they are not relevant.

## 2.2 Long vs. wide format

Finally, calc_ES() includes an option to change the format of the output. The function defaults to format = "long"; setting format = "wide" will return all of the results as a single line, rather than one line per effect size:

calc_ES(A_data = A, B_data = B, ES = c("NAP","PND","SMD"))
#>    ES       Est         SE  CI_lower  CI_upper
#> 1 NAP 0.9166667 0.06900656 0.5973406 0.9860176
#> 2 PND 0.7142857         NA        NA        NA
#> 3 SMD 1.7037340 0.63701390 0.4552097 2.9522583

calc_ES(A_data = A, B_data = B, ES = c("NAP","PND","SMD"), format = "wide")
#>     NAP_Est     NAP_SE NAP_CI_lower NAP_CI_upper   PND_Est  SMD_Est
#> 1 0.9166667 0.06900656    0.5973406    0.9860176 0.7142857 1.703734
#>      SMD_SE SMD_CI_lower SMD_CI_upper
#> 1 0.6370139    0.4552097     2.952258

# 3 batch_calc_ES()

Most single-case studies include multiple cases, and many also include multiple dependent variables measured on each case. Thus, it will often be of interest to calculate effect size estimates for multiple data series from a study, or even from multiple studies. The batch_calc_ES() function does exactly this—calculating any of the previously detailed effect sizes for each of several data series. Its syntax is a bit more involved than the previous functions, and so we provide several examples here. In what follows, we will assume that you are already comfortable using the es_calc() function as well as the other individual effect size functions in the package.

## 3.1 Data organization

Unlike with the other functions in the package, the input data for batch_calc_ES() must be organized in a data frame, with one line corresponding to each observation within a series, and columns corresponding to different variables (e.g. outcome, phase, session number). One or more variables must be included that uniquely identify every data series. Let’s look at two examples.

### 3.1.1 McKissick

The McKissick dataset is data drawn from McKissick, Hawkins, Lentz, Hailley, & McGuire (2010), a single-case design study of a group contingency intervention. The study used a multiple baseline design across three classrooms. The outcome data are event counts of disruptive behaviors observed at the classroom level.

data(McKissick)

Here are the first few rows of the data:

Case_pseudonym Session_number Condition Outcome Session_length Procedure
Period 1 1 A 13.62 20 count
Period 1 2 A 12.57 20 count
Period 1 3 A 15.76 20 count
Period 1 4 B 5.97 20 count
Period 1 5 B 4.63 20 count
Period 1 6 B 5.82 20 count
Period 1 7 B 3.72 20 count
Period 1 8 B 8.07 20 count
Period 1 9 B 2.95 20 count
Period 1 10 B 11.86 20 count

### 3.1.2 Schmidt (2007)

The Schmidt2007 dataset are data drawn from Schmidt (2007). This data set is somewhat more complicated. It has two outcomes for each participant, and the outcomes differ in directions of therapeutic improvement and measurement scale. The study used an ABAB design, replicated across three participants. Each series therefore has four phases: a baseline phase, a treatment phase, a return to baseline phase, and a second treatment phase.

data(Schmidt2007)

Here are the first few rows of the data

Case_pseudonym Behavior_type Session_number Outcome Condition Phase_num Metric Session_length direction n_Intervals
Albert Disruptive Behavior 1 22.944464 A 1 count 10 decrease NA
Albert Disruptive Behavior 2 22.431292 A 1 count 10 decrease NA
Albert Disruptive Behavior 3 27.785380 A 1 count 10 decrease NA
Albert Disruptive Behavior 4 16.928954 A 1 count 10 decrease NA
Albert Disruptive Behavior 5 21.838294 A 1 count 10 decrease NA
Albert Disruptive Behavior 6 3.780363 A 1 count 10 decrease NA
Albert Disruptive Behavior 7 18.137758 A 1 count 10 decrease NA
Albert Disruptive Behavior 8 11.774433 A 1 count 10 decrease NA
Albert Disruptive Behavior 9 22.083476 A 1 count 10 decrease NA
Albert Disruptive Behavior 10 4.986945 B 1 count 10 decrease NA

The Schmidt (2007) dataset contains many variables, but for now let’s focus on the following:

• Case_Pseudonym uniquely identifies each of the three participants
• Behavior_type specifies whether the outcome is disruptive behavior or on-task behavior
• Session_number specifies the order of the sessions within each data series
• Outcome contains the dependent variable measurements
• Condition specifies whether the outcome is in a baseline (“A”) condition or a treatment (“B”) condition
• Phase_num specifies whether the session is in the first or second pair of phases in the design
• Metric specifies whether the dependent variable is percentage or count data
• Session_length specifies the length of the observation session
• direction specifies the direction of therapeutic improvement
• n_Intervals specifies the number of intervals per session for the dependent variable measured using partial interval recording.

## 3.2 Main arguments of batch_calc_ES()

Here are the arguments for the batch calculator function:

args(batch_calc_ES)
#> function (dat, grouping, condition, outcome, session_number = NULL,
#>     baseline_phase = NULL, ES = c("LRRd", "LRRi", "SMD", "Tau"),
#>     improvement = "increase", scale = "other", intervals = NA,
#>     observation_length = NA, confidence = 0.95, format = "long",
#>     warn = TRUE, ...)
#> NULL

This function has a bunch of arguments, but many of them are optional and only used for certain effect size metrics. For the moment, let’s focus on the first few arguments, which are all we need to get going:

• The argument dat should be a dataframe containing all of the observations for all of the data series of interest.

• The grouping argument should specify the set of variables that uniquely identify each series. For a single study consisting of several series, like the McKissick dataset, this might simply be a variable name that identifies the participant pseudonym. Specify using bare variable names (i.e., without quotes).

• The condition argument should be the variable that identifies the treatment condition for each observation in the series. Specify using a bare variable name. The values for the baseline and treatment phases should be uniform across all of the series within a dataset. That is, if some series are coded as “0” for baseline and “1” for treatment, whereas other series had “A” as baseline and “B” as treatment, you’ll first need to clean your data and standardize the coding.

• The outcome argument should be the variable that contains the outcomes of interest. Specify using a bare variable name.

• The ES argument allows you to specify which effect sizes are desired. By default, the batch calculator generates estimates of LRRd, LRRi, SMD, and Tau. However, you’re probably going to want to specify your own effect sizes. Just like calc_ES, you request your desired effect sizes as a character vector, with the individual options of "LRRd", "LRRi", "LOR", "SMD", "NAP", "PND", "PEM", "PAND", "IRD", "Tau", or "Tau-U", in addition to "all" for all effect sizes, "NOM" for all non-overlap measures, and "parametric" for all parametric effect sizes.

All of the other arguments are truly optional, and we’ll introduce them as we go along.

## 3.3 Using grouping variables

Let’s try applying the function to the McKissick data. Remember that these data contains an identifier for each case (Case_pseudonym), a variable (Condition) identifying the baseline (“A”) and treatment (“B”) phases, and an outcome variable containing the values of the outcomes. The outcomes are disruptive behaviors, so a decrease in the behavior corresponds to therapeutic improvement. Just as with the calc_ES() function, we’ll need to specify that using the improvement argument. In the example below, we will calculate estimates of NAP and PND, to keep things simple.

mckissick_ES <- batch_calc_ES(dat = McKissick,
grouping = Case_pseudonym,
condition = Condition,
outcome = Outcome,
improvement = "decrease",
ES = c("NAP", "PND"))

Note that all of the inputs related to variable names are bare (i.e., no quotes). Let’s take a look at a table of the output.

Case_pseudonym ES Est SE CI_lower CI_upper
Period 1 NAP 1.0000000 0.0000000 1.0000000 1.0000000
Period 1 PND 1.0000000 NA NA NA
Period 2 NAP 0.7714286 0.1538619 0.4305321 0.9322444
Period 2 PND 0.4285714 NA NA NA
Period 3 NAP 0.9166667 0.0833333 0.5676324 0.9874545
Period 3 PND 0.7500000 NA NA NA

The output will always start with one or more columns corresponding to each unique combination of values from the grouping argument, followed by a column for each effect size. If any of the requested effect sizes have standard errors and confidence intervals, there will also be columns corresponding to the standard error and the upper and lower limit. Here, PND has NA for each of those, because it does not have a known standard error or confidence interval.

Now let’s look at an example using the Schmidt data. Remember that these data contain a pseudonym that uniquely identifies each of the three participants (Case_Pseudonym) as well as a variable that specifies whether the outcome is disruptive behavior or on-task behavior (Behavior_type). Furthermore, these data come from a treatment reversal design with two pairs of AB phases for each combination of case and behavior type (Phase_num). We’re going to want an effect size for each combination of pseudonym, behavior, and phase pair. The data also have an outcome variable (Outcome) and a variable identifying whether it was in the baseline (“A”) or treatment (“B”) phase (Condition). Finally, the the two different behavior types have different direction therapeutic improvement, so there is a variable called direction that specifies "increase" for on-task behavior or "decrease" for disruptive behavior.

Here’s an example of how to calculate NAP and PND for these data:

schmidt_ES <- batch_calc_ES(dat = Schmidt2007,
grouping = c(Case_pseudonym, Behavior_type, Phase_num),
condition = Condition,
outcome = Outcome,
improvement = direction,
ES = c("NAP", "PND"))

The syntax is similar to the example with the McKissick dataset, except for two things. Here, we’ve provided a vector of variable names for grouping that identify each series for which we want an effect size. Instead of providing a uniform direction of improvement to the improvement variable, we’ve provided a variable name, direction, which will account for the fact that the two behavior types have different directions of therapeutic improvement. Here is a table of the output:

Case_pseudonym Behavior_type Phase_num ES Est SE CI_lower CI_upper
Albert Disruptive Behavior 1 NAP 0.9583333 0.0416667 0.7041919 0.9947601
Albert Disruptive Behavior 1 PND 0.6250000 NA NA NA
Albert Disruptive Behavior 2 NAP 1.0000000 0.0000000 1.0000000 1.0000000
Albert Disruptive Behavior 2 PND 1.0000000 NA NA NA
Albert On Task Behavior 1 NAP 0.7708333 0.1269686 0.4880323 0.9162355
Albert On Task Behavior 1 PND 0.0000000 NA NA NA
Albert On Task Behavior 2 NAP 0.9333333 0.0666667 0.4950786 0.9942060
Albert On Task Behavior 2 PND 0.8000000 NA NA NA
Faith Disruptive Behavior 1 NAP 0.7767857 0.1468576 0.5206466 0.9116290
Faith Disruptive Behavior 1 PND 0.0714286 NA NA NA
Faith Disruptive Behavior 2 NAP 1.0000000 0.0000000 1.0000000 1.0000000
Faith Disruptive Behavior 2 PND 1.0000000 NA NA NA
Faith On Task Behavior 1 NAP 0.5803571 0.1353858 0.3399107 0.7843053
Faith On Task Behavior 1 PND 0.0000000 NA NA NA
Faith On Task Behavior 2 NAP 0.8666667 0.1333333 0.4329439 0.9799281
Faith On Task Behavior 2 PND 0.6000000 NA NA NA
Lilly Disruptive Behavior 1 NAP 1.0000000 0.0000000 1.0000000 1.0000000
Lilly Disruptive Behavior 1 PND 1.0000000 NA NA NA
Lilly Disruptive Behavior 2 NAP 0.8611111 0.1438020 0.4426084 0.9770460
Lilly Disruptive Behavior 2 PND 0.5000000 NA NA NA
Lilly On Task Behavior 1 NAP 0.7350427 0.1451105 0.4843292 0.8849638
Lilly On Task Behavior 1 PND 0.0000000 NA NA NA
Lilly On Task Behavior 2 NAP 0.4444444 0.2939724 0.1530430 0.7830060
Lilly On Task Behavior 2 PND 0.1666667 NA NA NA

The first three columns are the unique values from the variables supplied to grouping, followed by the effect size information.

## 3.4 Dealing with outcome scales.

By default, the batch calculator assumes the outcome scale is "other", which means that if a phase mean is equal to zero, the logs odd ratio or the log response ratio will not be calculated. Just like calc_ES(), when calculating parametric effect sizes, you may need to specify the outcome scales as well as things like the length of the observation session or the number of intervals in each observation session. If these values are the same for all observations in the dataset, it can be quite simple, like this example using the McKissick dataset:

mckissick_ES <- batch_calc_ES(dat = McKissick,
grouping = Case_pseudonym,
condition = Condition,
outcome = Outcome,
improvement = "decrease",
scale = "count",
observation_length = 20,
ES = "parametric")
#> Warning: LOR can only be calculated for proportions or percentages. Will
#> return NAs for other outcome scales.

#> Warning: LOR can only be calculated for proportions or percentages. Will
#> return NAs for other outcome scales.

#> Warning: LOR can only be calculated for proportions or percentages. Will
#> return NAs for other outcome scales.

Here, we’ve specified the outcome scale is a count, and that each session lasted 20 minutes. Note that we get a warning about the log odds ratio. Let’s take a look at the output:

Case_pseudonym ES Est SE CI_lower CI_upper
Period 1 LRRd -0.8070640 0.1977514 -1.1946496 -0.4194785
Period 1 LRRi 0.8070640 0.1977514 0.4194785 1.1946496
Period 1 LOR NA NA NA NA
Period 1 SMD 3.5061471 1.2334987 1.0885340 5.9237602
Period 2 LRRd -0.6096109 0.3485842 -1.2928234 0.0736015
Period 2 LRRi 0.6096109 0.3485842 -0.0736015 1.2928234
Period 2 LOR NA NA NA NA
Period 2 SMD 1.2745438 0.6734201 -0.0453353 2.5944229
Period 3 LRRd -0.7478387 0.3534855 -1.4406575 -0.0550199
Period 3 LRRi 0.7478387 0.3534855 0.0550199 1.4406575
Period 3 LOR NA NA NA NA
Period 3 SMD 2.9206280 1.0781761 0.8074417 5.0338143

Once again, we have a column specifying the case to which the effect sizes correspond, as well as a column specifying the effect size metric. The log odds ratio returns all NAs, because the log odds ratio can’t be estimate for count outcomes.

Let’s assume that you are interested in estimating effect sizes for data where the measurement scale— as well as perhaps measurement details like the observation length or the number of intervals —varies depending on the data series. The Schmidt data is one example of this. Remember that the Schmidt data has a variable specifying the measurement scale of the outcome (Metric) which is "percentage" for desirable behavior and "count" for disruptive behaviors. It also has a variable that specifies the length of the observation session (Session_length), and a variable that specifies the number of intervals per session for the dependent variable measured using partial interval recording (n_Intervals). The value of Session_length is NA for the percentage outcomes and the value of n_Intervals is NA for the count outcomes because those details are not relevant for those outcome measurement scales. Let’s try it out:

schmidt_ES <- batch_calc_ES(dat = Schmidt2007,
grouping = c(Case_pseudonym, Behavior_type, Phase_num),
condition = Condition,
outcome = Outcome,
improvement = direction,
scale = Metric,
observation_length = Session_length,
intervals = n_Intervals,
ES = c("parametric"))
#> Warning: LOR can only be calculated for proportions or percentages. Will
#> return NAs for other outcome scales.

#> Warning: LOR can only be calculated for proportions or percentages. Will
#> return NAs for other outcome scales.

#> Warning: LOR can only be calculated for proportions or percentages. Will
#> return NAs for other outcome scales.

#> Warning: LOR can only be calculated for proportions or percentages. Will
#> return NAs for other outcome scales.

#> Warning: LOR can only be calculated for proportions or percentages. Will
#> return NAs for other outcome scales.

#> Warning: LOR can only be calculated for proportions or percentages. Will
#> return NAs for other outcome scales.

Unlike the previous example, where we specified a uniform value for the scale and observation_length, we now have to specify variable names for scale, observation_length, and the number of intervals. Note that we get some warnings again about the LOR effect size. Let’s take a look at the output:

kable(schmidt_ES) %>%
kable_styling() %>%
scroll_box(width = "100%", height = "800px")
Case_pseudonym Behavior_type Phase_num ES Est SE CI_lower CI_upper
Albert Disruptive Behavior 1 LRRd -1.6062336 0.3238081 -2.2408858 -0.9715813
Albert Disruptive Behavior 1 LRRi 1.6062336 0.3238081 0.9715813 2.2408858
Albert Disruptive Behavior 1 LOR NA NA NA NA
Albert Disruptive Behavior 1 SMD 1.9197796 0.5318026 0.8774657 2.9620935
Albert Disruptive Behavior 2 LRRd -1.6511675 0.3758734 -2.3878658 -0.9144692
Albert Disruptive Behavior 2 LRRi 1.6511675 0.3758734 0.9144692 2.3878658
Albert Disruptive Behavior 2 LOR NA NA NA NA
Albert Disruptive Behavior 2 SMD 1.2574608 0.5673053 0.1455628 2.3693587
Albert On Task Behavior 1 LRRd -1.7159973 0.3933614 -2.4869715 -0.9450231
Albert On Task Behavior 1 LRRi 0.3226186 0.1293019 0.0691916 0.5760456
Albert On Task Behavior 1 LOR 2.0386159 0.5003548 1.0579385 3.0192934
Albert On Task Behavior 1 SMD 0.9324019 0.3696794 0.2078436 1.6569602
Albert On Task Behavior 2 LRRd -1.1652267 0.9169379 -2.9623919 0.6319385
Albert On Task Behavior 2 LRRi 0.2412181 0.2007565 -0.1522573 0.6346935
Albert On Task Behavior 2 LOR 1.4064448 1.0777694 -0.7059443 3.5188340
Albert On Task Behavior 2 SMD 0.6457036 0.4938273 -0.3221800 1.6135873
Faith Disruptive Behavior 1 LRRd -1.1680541 0.2268083 -1.6125902 -0.7235180
Faith Disruptive Behavior 1 LRRi 1.1680541 0.2268083 0.7235180 1.6125902
Faith Disruptive Behavior 1 LOR NA NA NA NA
Faith Disruptive Behavior 1 SMD 1.1042564 0.4090346 0.3025634 1.9059494
Faith Disruptive Behavior 2 LRRd -1.4268305 0.3100491 -2.0345154 -0.8191455
Faith Disruptive Behavior 2 LRRi 1.4268305 0.3100491 0.8191455 2.0345154
Faith Disruptive Behavior 2 LOR NA NA NA NA
Faith Disruptive Behavior 2 SMD 1.3269609 0.5818037 0.1866466 2.4672753
Faith On Task Behavior 1 LRRd -0.0090416 0.3639177 -0.7223071 0.7042239
Faith On Task Behavior 1 LRRi 0.0147192 0.1143270 -0.2093577 0.2387961
Faith On Task Behavior 1 LOR 0.0237609 0.4782220 -0.9135371 0.9610588
Faith On Task Behavior 1 SMD 0.0620767 0.5873491 -1.0891064 1.2132597
Faith On Task Behavior 2 LRRd -2.6984424 0.5341596 -3.7453761 -1.6515088
Faith On Task Behavior 2 LRRi 0.6044398 0.6716160 -0.7119033 1.9207829
Faith On Task Behavior 2 LOR 3.3028823 1.1778944 0.9942516 5.6115130
Faith On Task Behavior 2 SMD 0.8084812 0.4838922 -0.1399301 1.7568925
Lilly Disruptive Behavior 1 LRRd -1.7491983 0.2097252 -2.1602520 -1.3381445
Lilly Disruptive Behavior 1 LRRi 1.7491983 0.2097252 1.3381445 2.1602520
Lilly Disruptive Behavior 1 LOR NA NA NA NA
Lilly Disruptive Behavior 1 SMD 1.9391889 0.5215124 0.9170434 2.9613344
Lilly Disruptive Behavior 2 LRRd -0.9465255 0.5379386 -2.0008657 0.1078148
Lilly Disruptive Behavior 2 LRRi 0.9465255 0.5379386 -0.1078148 2.0008657
Lilly Disruptive Behavior 2 LOR NA NA NA NA
Lilly Disruptive Behavior 2 SMD 0.8023088 0.5242376 -0.2251780 1.8297956
Lilly On Task Behavior 1 LRRd -1.5600664 0.5111279 -2.5618586 -0.5582741
Lilly On Task Behavior 1 LRRi 0.4210015 0.1621196 0.1032530 0.7387500
Lilly On Task Behavior 1 LOR 1.9810679 0.6295731 0.7471273 3.2150084
Lilly On Task Behavior 1 SMD 1.0477764 0.3963880 0.2708701 1.8246826
Lilly On Task Behavior 2 LRRd -0.8708337 0.9403393 -2.7138649 0.9721975
Lilly On Task Behavior 2 LRRi 0.0524781 0.1169389 -0.1767180 0.2816742
Lilly On Task Behavior 2 LOR 0.9233118 1.0551020 -1.1446500 2.9912736
Lilly On Task Behavior 2 SMD 0.2230707 0.4396028 -0.6385349 1.0846764

In this case, LOR is all NA for the outcomes that are disruptive behaviors because those are counts and therefore the LOR isn’t an appropriate effect size. However, for the percentage of on task behavior, the LOR was estimated.

## 3.5 Further arguments

### 3.5.1 Output format

We can also request the effect sizes in a wide format:

mckissick_wide_ES <- batch_calc_ES(dat = McKissick,
grouping = Case_pseudonym,
condition = Condition,
outcome = Outcome,
improvement = "decrease",
ES = c("NAP", "PND"),
format = "wide")

The default argument for the batch calculator is format = "long", but if you want each case to be on a single line, specifying format = "wide" will provide the output that way, just like calc_ES(). Here’s the output:

Case_pseudonym NAP_Est NAP_SE NAP_CI_lower NAP_CI_upper PND_Est
Period 1 1.0000000 0.0000000 1.0000000 1.0000000 1.0000000
Period 2 0.7714286 0.1538619 0.4305321 0.9322444 0.4285714
Period 3 0.9166667 0.0833333 0.5676324 0.9874545 0.7500000

In this case there is a column for NAP, NAP’s standard error, and the upper and lower bound of the confidence interval. PND only has a column for the estimate, but remember that the values for SE and upper and lower CI were all NA in the long format. Columns that would have all NA values are removed when specifying format = "wide".

### 3.5.2 Suppressing warnings

Remember how, when we asked for the LOR for counts, the calculator gave us a bunch of warning messages? If you’re asking for the LOR, and some of your outcomes are in a scale other than percentage or proportion, you can specify the argument warn = "FALSE" (by default it is set to TRUE) if you want to suppress the warning messages. You will still get NA for any series with an inappropriate outcome scale.

batch_calc_ES(dat = McKissick,
grouping = Case_pseudonym,
condition = Condition,
outcome = Outcome,
improvement = "decrease",
scale = "count",
observation_length = 20,
ES = c("LRRi","LOR"),
warn = FALSE)
#> # A tibble: 6 x 6
#>   Case_pseudonym ES       Est     SE CI_lower CI_upper
#>   <chr>          <chr>  <dbl>  <dbl>    <dbl>    <dbl>
#> 1 Period 1       LRRi   0.807  0.198   0.419      1.19
#> 2 Period 1       LOR   NA     NA      NA         NA
#> 3 Period 2       LRRi   0.610  0.349  -0.0736     1.29
#> 4 Period 2       LOR   NA     NA      NA         NA
#> 5 Period 3       LRRi   0.748  0.353   0.0550     1.44
#> 6 Period 3       LOR   NA     NA      NA         NA

### 3.5.3 Arguments to individual effect size functions

The ... argument allows you to specify arguments particular to an individual function such as std_dev for the SMD() function. For instance, compare the results of calculating a pooled SMD versus the default, baseline phase only SMD:

batch_calc_ES(dat = McKissick,
grouping = Case_pseudonym,
condition = Condition,
outcome = Outcome,
ES = "SMD",
improvement = "decrease")
#> # A tibble: 3 x 6
#>   Case_pseudonym ES      Est    SE CI_lower CI_upper
#>   <chr>          <chr> <dbl> <dbl>    <dbl>    <dbl>
#> 1 Period 1       SMD    3.51 1.23    1.09       5.92
#> 2 Period 2       SMD    1.27 0.673  -0.0453     2.59
#> 3 Period 3       SMD    2.92 1.08    0.807      5.03

batch_calc_ES(dat = McKissick,
grouping = Case_pseudonym,
condition = Condition,
outcome = Outcome,
ES = "SMD",
improvement = "decrease",
std_dev = "pool")
#> # A tibble: 3 x 6
#>   Case_pseudonym ES      Est    SE CI_lower CI_upper
#>   <chr>          <chr> <dbl> <dbl>    <dbl>    <dbl>
#> 1 Period 1       SMD    2.58 0.853   0.909      4.25
#> 2 Period 2       SMD    1.12 0.588  -0.0345     2.27
#> 3 Period 3       SMD    2.34 0.727   0.920      3.77

Arguments common to several functions will be used when calculating any of the effect sizes for which they are relevant. For example, the bias_correct argument applies to all of the parametric effect sizes:

batch_calc_ES(dat = McKissick,
grouping = Case_pseudonym,
condition = Condition,
outcome = Outcome,
ES = "parametric",
improvement = "decrease",
scale = Procedure,
observation_length = Session_length,
bias_correct = FALSE,
warn = FALSE)
#> # A tibble: 12 x 6
#>    Case_pseudonym ES        Est     SE CI_lower CI_upper
#>    <chr>          <chr>   <dbl>  <dbl>    <dbl>    <dbl>
#>  1 Period 1       LRRd   -0.822  0.198  -1.21    -0.435
#>  2 Period 1       LRRi    0.822  0.198   0.435    1.21
#>  3 Period 1       LOR    NA     NA      NA       NA
#>  4 Period 1       SMD     4.82   2.17    0.571    9.07
#>  5 Period 2       LRRd   -0.650  0.349  -1.33     0.0328
#>  6 Period 2       LRRi    0.650  0.349  -0.0328   1.33
#>  7 Period 2       LOR    NA     NA      NA       NA
#>  8 Period 2       SMD     1.51   0.840  -0.133    3.16
#>  9 Period 3       LRRd   -0.807  0.353  -1.50    -0.114
#> 10 Period 3       LRRi    0.807  0.353   0.114    1.50
#> 11 Period 3       LOR    NA     NA      NA       NA
#> 12 Period 3       SMD     3.19   1.22    0.807    5.58

The bias_correct argument cannot be specified differently for different effect size functions. If you wanted bias corrected values for the LRRd effect size and not for the SMD effect size, you would need to call batch_calc_ES() separately for the two different effect sizes.

### 3.5.4 Order of observations

The session_number argument orders the data within each series by the specified variable. This argument is only important if Tau-U is being calculated, because it adjusts for trend in the baseline and so the ordering of the baseline phase is important. This argument is irrelevant for all of the other effect sizes.

### 3.5.5 Specifying a baseline phase

The baseline_phase argument works in a similar way the calc_ES() argument. If nothing is specified, the first phase in each series will be treated as the baseline phase. However, if the baseline phase is not always the first phase in each series, such as an SCD with four cases that use a cross-over treatment reversal design, where two of the cases follow an ABAB design and the other two cases follow a BABA design, you will have to specify the baseline_phase in the same way as the calc_ES() function.

### 3.5.6 Confidence levels

The confidence argument controls the confidence intervals in the same way as all the other functions. To skip calculating confidence intervals, specify confidence = NULL:

batch_calc_ES(dat = McKissick,
grouping = Case_pseudonym,
condition = Condition,
outcome = Outcome,
ES = "parametric",
improvement = "decrease",
scale = Procedure,
observation_length = Session_length,
confidence = NULL,
warn = FALSE)
#> # A tibble: 12 x 4
#>    Case_pseudonym ES        Est     SE
#>    <chr>          <chr>   <dbl>  <dbl>
#>  1 Period 1       LRRd   -0.807  0.198
#>  2 Period 1       LRRi    0.807  0.198
#>  3 Period 1       LOR    NA     NA
#>  4 Period 1       SMD     3.51   1.23
#>  5 Period 2       LRRd   -0.610  0.349
#>  6 Period 2       LRRi    0.610  0.349
#>  7 Period 2       LOR    NA     NA
#>  8 Period 2       SMD     1.27   0.673
#>  9 Period 3       LRRd   -0.748  0.353
#> 10 Period 3       LRRi    0.748  0.353
#> 11 Period 3       LOR    NA     NA
#> 12 Period 3       SMD     2.92   1.08

# References

McKissick, C., Hawkins, R. O., Lentz, F. E., Hailley, J., & McGuire, S. (2010). Randomizing multiple contingency components to decrease disruptive behaviors and increase student engagement in an urban second-grade classroom. Psychology in the Schools, 47(9), 944–959. https://doi.org/10.1002/pits.20516

Parker, R. I., & Vannest, K. (2009). An improved effect size for single-case research: Nonoverlap of all pairs. Behavior Therapy, 40(4), 357–367.

Pustejovsky, J. E. (2015). Measurement-comparable effect sizes for single-case studies of free-operant behavior. Psychological Methods, 20(3), 342.

Pustejovsky, J. E. (2018). Using response ratios for meta-analyzing single-case designs with behavioral outcomes. Journal of School Psychology, 68, 99–112.

Schmidt, A. C. (2007). The effects of a group contingency on group and individual behavior in an urban first-grade classroom (Doctoral dissertation). University of Kansas.

Scruggs, T. E., Mastropieri, M. A., & Casto, G. (1987). The quantitative synthesis of single-subject research: Methodology and validation. Remedial and Special Education, 8(2), 24–33.

Swan, D. M., & Pustejovsky, J. E. (2018). A Gradual Effects Model for Single-Case Designs. Multivariate Behavioral Research. https://doi.org/10.1080/00273171.2018.1466681