Estimation

For an outcome where increase is desirable, calculate

q_ij = I(y_j^B > y_i^A) + 0.5I(y_j^B = y_i^A) for i = 1, ..., m and j = 1, ..., n. For an outcome where decrease is desirable, one would instead use

q_ij = I(y_j^B < y_i^A) + 0.5I(y_j^B = y_i^A).

The NAP effect size index is then calculated as

$$ \text{NAP} = \frac{1}{m n} \sum_{i=1}^m \sum_{j=1}^n q_{ij}. $$

Standard errors

The SingleCaseES package provides several different methods for estimating the standard error of NAP. The default method is calculated based on the exactly unbiased variance estimator described by Sen (1967; cf. Mee, 1990), which assumes that the observations are mutually independent and are identically distributed within each phase. Let

$$ \begin{aligned} Q_1 &= \frac{1}{m n^2} \sum_{i=1}^m \left[\sum_{j=1}^n \left(q_{ij} - \text{NAP}\right)\right]^2, \\ Q_2 &= \frac{1}{m^2 n} \sum_{j=1}^n \left[\sum_{i=1}^m \left(q_{ij} - \text{NAP}\right)\right]^2, \qquad \text{and} \\ Q_3 &= \frac{1}{m n} \sum_{i=1}^m \sum_{j=1}^n \left(q_{ij} - \text{NAP}\right)^2. \end{aligned} $$

The SE is then calculated as

$$ SE_{\text{unbiased}} = \sqrt{\frac{\text{NAP}(1 - \text{NAP}) + n Q_1 + m Q_2 - 2 Q_3}{(m - 1)(n - 1)}}. $$

Another method for estimating a standard error was introduced by Hanley & McNeil (1982). This standard error is calculated as

$$ SE_{\text{Hanley}} = \sqrt{\frac{1}{mn}\left(\text{NAP}(1 - \text{NAP}) + (n - 1)Q_1 + (m - 1)Q_2\right)}, $$

with Q₁ and Q₂ defined as above. This standard error is based on the same assumptions as the unbiased SE.

A limitation of SE_unbiased and SE_Hanley is that they will be equal to zero when there is complete non-overlap (i.e., when NAP is equal to zero or equal to one). In order to ensure a strictly positive standard error for NAP, the SingleCaseES package calculates SE_unbiased and SE_Hanley using a truncation of NAP. Specifically, the formulas are evaluated using

$$ \widetilde{\text{NAP}} = \text{max}\left\{\frac{1}{2 mn}, \ \text{min}\left\{\frac{2mn - 1}{2mn}, \ \text{NAP} \right\} \right\} $$ in place of NAP.

A final method for estimating a standard error is to work under the null hypothesis that there is no effect—i.e., that the data points from each phase are sampled from the same distribution. Under the null hypothesis, the sampling variance of NAP depends only on the number of observations in each phase:

$$ SE_{\text{null}} = \sqrt{\frac{m + n + 1}{12 m n}} $$ (cf. Grissom & Kim, 2001, p. 141). If null hypothesis is not true—that is, if the observations in phase B are drawn from a different distribution than the observations in phase A—then this standard error will tend to be too large.

Confidence interval

A confidence interval for θ can be calculated using a method proposed by Newcombe [Newcombe (2006); Method 5], which assumes that the observations are mutually independent and are identically distributed within each phase. Using a confidence level of 100% × (1 − α), the endpoints of the confidence interval are defined as the values of θ that satisfy the equality

$$ (\text{NAP} - \theta)^2 = \frac{z^2_{\alpha / 2} h \theta (1 - \theta)}{mn}\left[\frac{1}{h} + \frac{1 - \theta}{2 - \theta} + \frac{\theta}{1 + \theta}\right], $$

where h = (m + n)/2 − 1 and z_α/2 is 1 − α/2 critical value from a standard normal distribution. This equation is a fourth-degree polynomial in θ, solved using a numerical root-finding algorithm.

Estimation

For an outcome where increase is desirable,

$$ \text{PND} = \frac{1}{n} \sum_{j=1}^n I(y^B_j > y^A_{(m)}), $$

where y_(m)^A is the maximum value of y₁^A, ..., y_m^A. For an outcome where decrease is desirable,

$$ \text{PND} = \frac{1}{n} \sum_{j=1}^n I(y^B_j < y^A_{(1)}), $$

where y₍₁₎^A is the minimum value of y₁^A, ..., y_m^A.

The sampling distribution of PND has not been described, and so standard errors and confidence intervals are not available.

Estimation

For an outcome where increase is desirable,

$$ \text{PEM} = \frac{1}{n}\sum_{j=1}^n \left[ I(y^B_j > m_A) + 0.5 \times I(y^B_j = m_A) \right], $$

where m_A = median(y₁^A, ..., y_m^A). For an outcome where decrease is desirable,

$$ \text{PEM} = \frac{1}{n}\sum_{j=1}^n \left[ I(y^B_j < y^A_{(1)}) + 0.5 \times I(y^B_j = m_A) \right]. $$

The sampling distribution of PEM has not been described, and so standard errors and confidence intervals are not available.

Estimation

For an outcome where increase is desirable, PAND is calculated as

$$ \text{PAND} = \frac{1}{m + n} \max \left\{\left(i + j\right) I\left(y^A_{(i)} < y^B_{(n + 1 - j)}\right)\right\}, $$

where y₍₀₎^A = −∞, y_(n + 1)^B = ∞, and the maximum is taken over the values 0 ≤ i ≤ m and 0 ≤ j ≤ n. For an outcome where decrease is desirable, PAND is calculated as

$$ \text{PAND} = \frac{1}{m + n} \max \left\{\left(i + j\right) I\left(y^A_{(m + 1 - i)} > y^B_{(j)}\right)\right\}, $$

where y_(m + 1)^A = ∞, y₍₀₎^B = −∞, and the maximum is taken over the values 0 ≤ i ≤ m and 0 ≤ j ≤ n.

The sampling distribution of PAND has not been described, and so standard errors and confidence intervals are not available.

Estimation

For notational convenience, let y₍₀₎^A = y₍₀₎^B = −∞ and y_(m + 1)^A = y_(n + 1)^B = ∞. For an outcome where increase is desirable, let ĩ and j̃ denote the values that maximize the quantity

(i + j)I(y_(i)^A < y_{(n + 1 − j)}^B) for 0 ≤ i ≤ m and 0 ≤ j ≤ n. For an outcome where decrease is desirable, let ĩ and j̃ instead denote the values that maximize the quantity

(i + j)I(y_{(m + 1 − i)}^A > y_(j)^B).

Now calculate the 2 × 2 table

$$ \begin{array}{|c|c|} \hline m - \tilde{i} & \tilde{j} \\ \hline \tilde{i} & n - \tilde{j} \\ \hline \end{array} $$

Parker, Vannest, & Brown (2009) proposed the non-robust improvement rate difference, which is equivalent to the phi coefficient from this table. Parker, Vannest, & Davis (2011) proposed to instead use the robust phi coefficient, which involves modifying the table so that the row- and column-margins are equal. Robust IRD is thus equal to

$$ \text{IRD} = \frac{n - m - \tilde{i} - \tilde{j}}{2 n} - \frac{m + n - \tilde{i} - \tilde{j}}{2 m}. $$

Robust IRD is algebraically related to PAND as

$$ \text{IRD} = 1 - \frac{(m + n)^2}{2mn}\left(1 - \text{PAND}\right). $$ Just as with PAND, the sampling distribution of robust IRD has not been described, and so standard errors and confidence intervals are not available.

Estimation

For an outcome where increase is desirable, calculate

w_ij = I(y_j^B > y_i^A) − I(y_j^B < y_i^A)

For an outcome where decrease is desirable, one would instead use

w_ij = I(y_j^B < y_i^A) − I(y_j^B > y_i^A).

The Tau effect size index is then calculated as

$$ \text{Tau} = \frac{1}{m n} \sum_{i=1}^m \sum_{j=1}^n w_{ij} = 2 \times \text{NAP} - 1. $$

Standard errors

Standard errors and confidence intervals for Tau are calculated using transformations of the corresponding SEs and CIs for NAP. All of the methods assume that the observations are mutually independent and are identically distributed within each phase.

Standard errors for Tau are calculated as SE_Tau = 2SE_NAP, where SE_NAP is the standard error for NAP calculated based on one of the available methods (unbiased, Hanley, or null).

Confidence intervals

The CI for τ is calculated as

[L_τ, U_τ] = [2L_θ − 1, 2U_θ − 1],

where L_θ and U_θ are the lower and upper bounds of the CI for θ, calculated using a method proposed by Newcombe (2006).

Estimation

For an outcome where increase is desirable, calculate

w_ij^AB = I(y_j^B > y_i^A) − I(y_j^B < y_i^A)

and

w_ij^AA = I(y_j^A > y_i^A) − I(y_j^A < y_i^A)

For an outcome where decrease is desirable, one would instead use

w_ij^AB = I(y_j^B < y_i^A) − I(y_j^B > y_i^A)

and

w_ij^AA = I(y_j^A < y_i^A) − I(y_j^A > y_i^A).

The Tau-U effect size index is then calculated as

$$ \text{Tau-U} = \frac{1}{m n} \left(\sum_{i=1}^m \sum_{j=1}^n w^{AB}_{ij} - \sum_{i=1}^{m - 1} \sum_{j=i + 1}^m w^{AA}_{ij}\right). $$

The sampling distribution of Tau-U has not been described, and so standard errors and confidence intervals are not available.

Estimation

Estimation of τ_BC entails correcting the data series for the baseline slope β. If using the baseline trend pre-test, the null hypothesis of H₀ : β = 0 is first tested using Kendall’s rank correlation. If the test is not significant, then set β̂ = 0 and α̂ = 0. If the test is significant or if the pre-test for baseline time trend is not used, then the slope is estimated by Theil-Sen regression. Specifically, we calculate the slope between every pair of observations in the A phase: $$ s_{hi} = \frac{y_i^A - y_h^A}{i - h} $$ for i = 1, ..., m − 1 and h = i + 1, ..., m. The overall slope estimate is taken to be the median over all m(m − 1)/2 slope pairs:

β̂ = median{s₁₂, ..., s_{(m − 1)m}}. The intercept term is estimated by taking the median observation in the A phase after correcting for the estimated linear time trend:

α̂ = median{y₁^A − β̂ × 1, y₂^A − β̂ × 2, ..., y_m^A − β̂ × m}. However, the intercept estimate is irrelevant for purposes of estimating Tau-BC because the Tau estimator is a function of ranks and is invariant to a linear shift of the data series.

After estimating the phase A time trend, τ_BC is estimated by de-trending the full data series and calculating Kendall’s rank correlation or Tau (non-overlap) on the de-trended observations. Specifically, set ϵ̂_i^A = y_i^A − β̂(i) − α̂ for i = 1, ..., m and ϵ̂_j^B = y_j^B − β̂(m + j) − α̂. For an outcome where increase is desirable, calculate w_ij^ϵ = I(ϵ̂_j^B > ϵ̂_i^A) − I(ϵ̂_j^B < ϵ̂_i^A)

or, for an outcome where decrease is desirable, calculate w_ij^ϵ = I(ϵ̂_j^B < ϵ̂_i^A) − I(ϵ̂_j^B > ϵ̂_i^A).

Tau-BC (non-overlap) is then estimated by $$ \text{Tau}_{BC} = \frac{1}{m n} \sum_{i=1}^m \sum_{j=1}^n w^\epsilon_{ij}. $$

If calculated with Kendall’s rank correlation, Tau-BC is estimated as the rank correlation between {ϵ̂₁^A, …, ϵ̂_m^A, ϵ̂₁^B, …, ϵ̂_n^B} and a dummy coded variable {0₁, …, 0_m, 1₁, …, 1_n}, with an adjustment for ties (Kendall, 1970, p. 35). Specifically,

$$ \text{Tau}_{BC}^* = \frac{1}{D} \sum_{i=1}^m \sum_{j=1}^n w^\epsilon_{ij}, $$ where

$$ D = \sqrt{m \times n \times \left(\frac{(m+n)(m+n-1)}{2}-U\right)} $$ and U is the number of ties between all possible pairs of observations (including pairs within phase A, pairs within phase B, and pairs of one phase A and one phase B data point). U can be computed as

$$ U = \sum_{i=1}^{m - 1} \sum_{j = i+1}^m I\left(\hat\epsilon^A_i = \hat\epsilon^A_j\right) + \sum_{i=1}^{n - 1} \sum_{j = i+1}^n I\left(\hat\epsilon^B_i = \hat\epsilon^B_j\right) + \sum_{i=1}^m \sum_{j=1}^n I\left(\hat\epsilon^A_i = \hat\epsilon^B_j\right). $$

We prefer and recommend to use the Tau-AB form, which divides by m × n rather than by D, because it leads to a simpler interpretation. Furthermore, using D means that Tau_BC^* may be sensitive to variation in phase lengths. To see this sensitivity, consider a scenario where there are no tied values and so every value {ϵ̂₁^A, …, ϵ̂_m^A, ϵ̂₁^B, …, ϵ̂_n^B} is unique. In this case, U = 0 and $$ D = \sqrt{\frac{1}{2} m n (m + n)(m + n - 1)} = m n \sqrt{1 + \frac{m - 1}{2n} + \frac{n - 1}{2m}}. $$ Thus, the denominator will always be larger than mn, meaning that Tau_BC^* will always be smaller than Tau_BC. Further, the largest and smallest possible values of Tau_BC^* will be ±mn/D, or about $1 / \sqrt{2}$ when m and n are close to equal. In contrast, the largest and smallest possible values of Tau_BC are always -1 and 1, respectively.

Standard errors and confidence intervals

The exact sampling distribution of Tau_BC^* (Kendall, adjusted for ties) has not been described. Tarlow (2017) proposed to approximate its sampling variance using $$ SE_{Kendall} = \sqrt{\frac{2 (1 - \text{Tau}_{BC}^2)}{m + n}}, $$ arguing that this would generally be conservative (in the sense of over-estimating the true sampling error). When Tau-BC is calculated using Kendall’s rank correlation, the SingleCaseES package reports a standard error based on this approximation.

When calculated without adjustment for ties, the SingleCaseES package takes a different approach for estimating the standard error for Tau_BC (non-overlap), reporting approximate standard errors and confidence intervals for Tau_BC based on the methods described above for Tau (non-overlap, without baseline trend correction). An important limitation of this approach is that it does not account for the uncertainty introduced by estimating the phase A time trend (i.e., the uncertainty in β̂).

Baseline SD

Gingerich (1984) and Busk & Serlin (1992) originally suggested scaling by the SD from phase A only, due to the possibility of non-constant variance across phases. Without assuming constant SDs, an estimate of the standardized mean difference is

$$ d_A = \left(1 - \frac{3}{4m - 5}\right) \frac{\bar{y}_B - \bar{y}_A}{s_A}. $$

The term in parentheses is a small-sample bias correction term (cf. Hedges, 1981; Pustejovsky, 2019). The standard error of this estimate is calculated as

$$ SE_{d_A} = \left(1 - \frac{3}{4m - 5}\right)\sqrt{\frac{1}{m} + \frac{s_B^2}{n s_A^2} + \frac{d_A^2}{2(m - 1)}}. $$

Pooled SD

If it is reasonable to assume that the SDs are constant across phases, then one can use the pooled sample SD, defined as

$$ s_p = \sqrt{\frac{(m - 1)s_A^2 + (n - 1) s_B^2}{m + n - 2}}. $$

The SMD can then be estimated as

$$ d_p = \left(1 - \frac{3}{4(m + n) - 9}\right) \frac{\bar{y}_B - \bar{y}_A}{s_p}, $$

with approximate standard error

$$ SE_{d_A} = \left(1 - \frac{3}{4(m + n) - 9}\right)\sqrt{\frac{1}{m} + \frac{1}{n} + \frac{d_p^2}{2(m + n - 2)}}. $$

Confidence intervals

Whether the estimator is based on the baseline or pooled standard deviation, an approximate confidence interval for δ is given by

[d − z_α/2 × SE_d,  d + z_α/2 × SE_d].

LRR-decreasing and LRR-increasing

There are two variants of the LRR (Pustejovsky, 2018), corresponding to whether therapeutic improvements correspond to negative values of the index (LRR-decreasing or LRRd) or positive values of the index (LRR-increasing or LRRi). For outcomes measured as frequency counts or rates, LRRd and LRRi are identical in magnitude but have opposite sign. However, for outcomes measured as proportions (ranging from 0 to 1) or percentages (ranging from 0% to 100%), LRRd and LRRi will differ in both sign and magnitude because the outcomes are first transformed to be consistent with the selected direction of therapeutic improvement.

Estimation

To account for the possibility that the sample means may be equal to zero, even if the mean levels are strictly greater than zero, the LRR is calculated using truncated sample means, given by $$ \tilde{y}_A = \text{max} \left\{ \bar{y}_A, \ \frac{1}{2 D m}\right\} \qquad \text{and} \qquad \tilde{y}_B = \text{max} \left\{ \bar{y}_B, \ \frac{1}{2 D n}\right\}, $$ where D is a constant that depends on the scale and recording procedure used to measure the outcomes (Pustejovsky, 2018). To ensure that the standard error of LRR is strictly positive, it is calculated using truncated sample variances, given by $$ \tilde{s}_A^2 = \text{max}\left\{s_A^2, \ \frac{1}{D^2 m^3}\right\} \qquad \text{and} \qquad \tilde{s}_B^2 = \text{max} \left\{ s_B^2, \ \frac{1}{D^2 n^3}\right\}. $$

A basic estimator of the LRR is then given by

R₁ = ln (ỹ_B) − ln (ỹ_A).

However, R₁ will be biased when one or both phases include only a small number of observations. A bias-corrected estimator is given by

$$ R_2 = \ln\left(\tilde{y}_B\right) + \frac{\tilde{s}_B^2}{2 n \tilde{y}_B^2} - \ln\left(\tilde{y}_A\right) - \frac{\tilde{s}_A^2}{2 m \tilde{y}_A^2}. $$ The bias-corrected estimator is the default option in SingleCaseES.

Standard errors

Under the assumption that the outcomes in each phase are mutually independent, an approximate standard error for R₁ or R₂ is given by

$$ SE_R = \sqrt{\frac{\tilde{s}_A^2}{m \tilde{y}_A^2} + \frac{\tilde{s}_B^2}{n \tilde{y}_B^2}}. $$

Confidence intervals

Under the same assumptions, an approximate confidence interval for ψ is

[R − z_α/2 × SE_R,  R + z_α/2 × SE_R].

Estimation

To account for the possibility that the sample means may be equal to zero or one, even if the mean levels are strictly between zero and one, the LOR is calculated using truncated sample means, given by

$$ \tilde{y}_A = \text{max} \left\{ \text{min}\left[\bar{y}_A, 1 - \frac{1}{2 D m}\right], \frac{1}{2 D m}\right\} $$ and

$$ \tilde{y}_B = \text{max} \left\{ \text{min}\left[\bar{y}_B, 1 - \frac{1}{2 D n}\right], \frac{1}{2 D n}\right\}, $$

where D is a constant that depends on the scale and recording procedure used to measure the outcomes (Pustejovsky, 2018).To ensure that the corresponding standard error is strictly positive, it is calculated using truncated sample variances, given by

$$ \tilde{s}_A^2 = \text{max}\left\{s_A^2, \ \frac{1}{D^2 m^3}\right\} \qquad \text{and} \qquad \tilde{s}_B^2 = \text{max} \left\{ s_B^2, \ \frac{1}{D^2 n^3}\right\}. $$

A basic estimator of the LOR is given by

LOR₁ = ln (ỹ_B) − ln (1 − ỹ_B) − ln (ỹ_A) + ln (1 − ỹ_A).

However, like the LRR, this estimator will be biased when the one or both phases include only a small number of observations. A bias-corrected estimator of the LOR is given by

$$ LOR_2 = \ln\left(\tilde{y}_B\right) - \ln\left(1-\tilde{y}_B\right) - \frac{\tilde{s}_B^2(2 \tilde{y}_B - 1)}{2 n_B (\tilde{y}_B)^2(1-\tilde{y}_B)^2} - \ln\left(\tilde{y}_A\right) + \ln\left(1-\tilde{y}_A\right) + \frac{\tilde{s}_A^2(2 \tilde{y}_A - 1)}{2 n_A (\tilde{y}_A)^2(1-\tilde{y}_A)^2}. $$ This estimator uses a small-sample correction to reduce bias when one or both phases include only a small number of observations.

Standard errors

Under the assumption that the outcomes in each phase are mutually independent, an approximate standard error for LOR is given by

$$ SE_{LOR} = \sqrt{\frac{\tilde{s}^2_A}{n_A \tilde{y}_A^2 (1 - \tilde{y}_A)^2} + \frac{\tilde{s}^2_B}{n_B \tilde{y}_B^2 (1 - \tilde{y}_B)^2}}. $$

Confidence intervals

Under the same assumption, an approximate confidence interval for ψ is

[LOR − z_α/2 × SE_LOR,  LOR + z_α/2 × SE_LOR],

Estimation

A natural estimator of the λ is given by

LRM = ln (m_B) − ln (m_A), where m_B and m_A are the sample medians for phase B and phase A, respectively. Note that the sample median might be zero for either phase B and phase A in some single-case design data, resulting in infinite LRM.

Standard errors

Standard errors and confidence intervals for LRM can be obtained under the assumption that the outcome data within each phase are mutually independent and follow a common distribution. Using the fact that the logarithm of the median is the same or close to the median of the log-transformed outcomes, the standard error for LRM can be calculated using the order statistics within each phase (Bonett & Price, 2020a). Let $$ \begin{aligned} l_A &= \text{max}\left\{1, \ \frac{m}{2} - \sqrt{m}\right\}, \quad &u_A &= m - l_A + 1, \\ l_B &= \text{max}\left\{1, \ \frac{n}{2} - \sqrt{n}\right\}, \quad &u_B &= n - l_B + 1, \end{aligned} $$ and find $$ q_A = \Phi^{-1}\left(\frac{1}{2^m}\sum_{i=0}^{l_A - 1} \frac{m!}{i!(m - i)!}\right) \quad \text{and} \quad q_B = \Phi^{-1}\left(\frac{1}{2^n}\sum_{j=0}^{l_B - 1} \frac{n!}{j!(n - j)!}\right). $$ The standard error of LRM is then $$ SE_{LRM} = \sqrt{\left(\frac{\ln\left(y^B_{(u_B)}\right)-\ln\left(y^B_{(l_B)}\right)}{2\ q_B}\right)^2 + \left(\frac{\ln\left(y^A_{(u_A)}\right)-\ln\left(y^A_{(l_A)}\right)}{2\ q_A}\right)^2}, $$ (Bonett & Price, 2020a) where y_(lA)^A, y_(uA)^A are the l_A and u_A order statistics of the phase A outcomes and y_(lB)^B, y_(uB)^B are the l_B and u_B order statistics of the phase B outcomes.

Confidence intervals

An approximate confidence interval for λ is [LRM − z_α/2 × SE_LRM,  LRM + z_α/2 × SE_LRM], where z_α/2 is 1 − α/2 critical value from a standard normal distribution.

Estimation

Approaches for estimation of PoGO depend on one’s assumption about the stability of the observations in phases A and B. Under the assumption that the observations are temporally stable, a natural estimator of PoGO is $$ PoGO = \frac{\bar{y}_B - \bar{y}_A}{\gamma - \bar{y}_A} \times 100\%. $$

Standard errors

Patrona, Ferron, Olszewski, Kelley, & Goldstein (2022) proposed a method for calculating a standard error for the PoGO estimator under the assumptions that the observations within each phase are mutually independent. The standard error uses an approximation for the standard error of two independent, normally distributed random variables due to Dunlap & Silver (1986). It is calculated as $$ SE_{PoGO} = \frac{1}{\gamma - \bar{y}_A} \sqrt{\frac{s_A^2}{n_A} + \frac{s_B^2}{n_B} + \left(\frac{\bar{y}_B - \bar{y}_A}{\gamma - \bar{y}_A}\right)^2 \frac{s_A^2}{n_A}}. $$ Patrona et al. (2022) also provided a more general approximation, which can be applied when PoGO is estimated using regressions that control for time trends or auto-correlation. However, these methods are not implemented in SingleCaseES.

Confidence intervals

An approximate confidence interval for PoGO is given by [PoGO − z_α/2 × SE_PoGO,  PoGO + z_α/2 × SE_PoGO], where z_α/2 is the 1 − α/2 critical value from a standard normal distribution (Patrona et al., 2022).