Example data

In the following sections the theory behind the methods in singcar is reviewed. But at the end of each subsection the usage of the implementations will be exemplified. To aid this the package comes with the dataset size_weight_illusion, which is a dataset from a neuropsychological study investigating the size-weight illusion in a well known patient “DF” (Hassan et al. 2020). DF incurred damage to the lateral occipital complex which gave rise to visual form agnosia, potentially making her less susceptible to the size-weight illusion (Dijkerman et al. 2004). This illusion is a perceptual phenomenon making objects of small size be perceived as heavier during lifting than objects of larger size but equal in weight (Buckingham 2014).

However, due to the location of DF’s brain damage one would only suspect less susceptibility to the illusion for visual cues since she other sensory processing should be unaffected of her damage. In the study by Hassan et al. (2020) the illusion is tested given visual cues in one condition and kinaesthetic cues in another. It was predicted that she would be abnormally less susceptible to visual size-weight illusion and that there would be an abnormally large discrepancy between visual and kinaesthetic size-weight illusion. The control sample consisted of 28 participants and their age and sex was collected along with their performance in the two conditions. The size-weight illusion is given in the dataset as a scaled measure of the number of grams weight difference perceived per cubic cm of volume change. To use this data, start by installing and loading the package:

install.packages("singcar")
library("singcar")

head(size_weight_illusion) 

##   GROUP PPT    SEX YRS      V_SWI      K_SWI
## 1    SC  DF Female  65 0.02814921 0.10012712
## 2    HC E01 Female  70 0.21692634 0.21792930
## 3    HC E02   Male  64 0.17223226 0.26338899
## 4    HC E03 Female  63 0.07138049 0.09331700
## 5    HC E04 Female  65 0.10186453 0.25938045
## 6    HC E05 Female  66 0.15911439 0.08922615

V_SWI and K_SWI are the measurements for visual and kinaesthetic size-weight illusion respectively. Most functions in singcar can take summary data as input, but for demonstrative purposes the raw data of this dataset will be used. To illustrate usage of the methods more generally, the scores of the variates of interest as well as of the covariate YRS is extracted from the data and renamed:

caseA <- size_weight_illusion[1, "V_SWI"]
contA <- size_weight_illusion[-1, "V_SWI"]
caseB <- size_weight_illusion[1, "K_SWI"]
contB <- size_weight_illusion[-1, "K_SWI"]
caseAGE <- size_weight_illusion[1, "YRS"]
contAGE <- size_weight_illusion[-1, "YRS"]

Here, caseA and caseB refers to the case scores on some task A and B respectively. Similarly, contA and contB refers to the scores of the control sample on the same tasks.

Example usage in singcar

The test of deficit is called in the singcar package with the function TD(). Testing for a deficit using the example data from Section @ref(exdata) we have:

TD(case = caseA, controls = contA, alternative = "less", conf_level = 0.95)

## 
##  Test of Deficit
## 
## data:  Case = 0.03, Controls (m = 0.16, sd = 0.08, n = 28)
## t = -1.7243, df = 27, p-value = 0.04804
## alternative hypothesis: true diff. between case and controls is less than 0
## sample estimates:
##   Std. case score (Z-CC), 95% CI [-2.34, -1.15] 
##                                       -1.754857 
## Proportion below case (%), 95% CI [0.95, 12.47] 
##                                        4.804003

Where the argument alternative specifies the alternative hypothesis, i.e. whether we are performing two or one sided test and if the latter, which direction. If one instead wants to test for a dissociation the most common alternative hypothesis is two sided. Using the RSDT for this test, call:

RSDT(case_a = caseA, case_b = caseB, controls_a = contA, controls_b = contB,
      alternative = "two.sided")

## 
##  Revised Standardised Difference Test
## 
## data:  Case A: 0.03, B: 0.10, Ctrl. A (m, sd): (0.16, 0.08), B: (0.18, 0.10)
## approx. abs. t = 1.015, df = 27, p-value = 0.3191
## alternative hypothesis: true difference between tasks is not equal to 0
## sample estimates:
## Std. case score, task A (Z-CC) Std. case score, task B (Z-CC) 
##                     -1.7548574                     -0.7836956 
##  Std. task discrepancy (Z-DCC)      Proportion below case (%) 
##                     -1.0647889                     15.9560625

Notably, this function does not produce confidence intervals. If the two tasks are directly comparable without the need of standardisation (which in fact is the case for this specific dataset), call instead UDT() with the same syntax and confidence intervals will be given.

The Bayesian test of deficit

The Bayesian analogue of the test of deficit (BTD) produces asymptotically identical output as the frequentist test of deficit. It produces an estimate of p with accompanying credible intervals, i.e. the percentage of controls that would exhibit a more extreme score than the case. The method is included here for completeness but for actual analysis there is no advantage to using it over the TD, (@ref(eq:TD)).

From a sample of size n we obtain measurements of some normally distributed variate Y with unknown mean μ and unknown variance σ². Let $\overline{y}$ and s² denote the sample mean and variance and y^* the case score. The prior distribution of μ is conditioned upon the prior distribution of σ² since both parameters are unknown. The prior considered here is f(μ, σ²) ∝ (σ²)⁻¹, which is the standard non-informative prior for normal data. The marginal posterior distribution of $\frac{(n-1)s^2}{\sigma^2} \sim \chi^2_{n-1}$ and the posterior distribution $\mu|\sigma^2 \sim \mathcal{N}(\overline{y}, \sigma^2/n)$, see e.g., Gelman et al. (2013) (p. 45 and 65).

When we have a posterior distribution for the control population parameters, p can be estimated by iterating the following steps.

1. Estimate σ̂_(i)² as $\hat{\sigma}^2_{(i)} = \frac{(n-1)s^2}{\psi}$ where ψ is generated from a χ² distribution on n − 1 degrees of freedom. The subscript (i) indicates that the statistic is an estimate of σ² for this iteration.
2. The estimate for μ for this iteration is given by $\hat{\mu}_{(i)}=\overline{y}+z\sqrt{(\hat{\sigma}_{(i)}^2/n)}$, where z is generated from a standard normal distribution.
3. The estimate of p is then calculated as if σ̂_(i)² and μ̂_(i) are the true population parameters. That is, put $z^*_{(i)}=\frac{y^* - \hat{\mu}_{(i)}}{\sqrt{\hat{\sigma}_{(i)}^2}}$ as the standardised case score and obtain an iteration estimate of p as p̂_(i) = Φ(z_(i)^*) or p̂_(i) = 1 − Φ(z_(i)^*) depending on alternative hypothesis.

Iterating these steps a large number of times yields a distribution of p̂, the mean of which is the point estimate of p, which can be used for significance testing. The 2.5th and 97.5th percentile of the distribution of p̂ as well as of z^* gives the boundaries for the 95% credible interval of p̂ and Ẑ_CC respectively. The point estimate of the latter is that of (@ref(eq:zcc)).

The Bayesian standardised difference test

In contrast to the BTD, the Bayesian analogue of the difference tests does not produce asymptotically identical results to the RSDT.
For the Bayesian standardised difference test (BSDT) we now obtain measurements on Y₁ and Y₂ following a bivariate normal distribution representing two tasks of interest, from a control sample of size n. The goal is to estimate the percentage p of the control population that would be expected to show a greater difference Y₁ − Y₂ than the case.

Let $\overline{y}_1$ and $\overline{y}_2$ be the sample means and
the sample variance-covariance matrix. Then S = A(n − 1) is the sums of squares and cross products (SSCP) matrix. The case scores are denoted y₁^* and y₂^*.

The population parameters are both unknown. To estimate p we must first obtain a posterior of these parameters. Since both are unknown the prior distribution of μ is again conditioned upon the prior distribution of Σ.

The prior considered for f(μ, Σ⁻¹) in Crawford and Garthwaite (2007) was f(μ, Σ⁻¹) ∝ |Σ|. This is a common non-informative prior, however, it was chosen in favour of the perhaps even more commonly used f(μ, Σ⁻¹) ∝ |Σ|^{(k + 1)/2} for multivariate normal data, where k= the number of variates. This was because the former prior was shown to give asymptotically identical interval estimates as the UDT for unstandardised case scores. Even though the latter perhaps is more commonly applied, Crawford and Garthwaite (2007) refer to the former as the “standard theory” prior.

For this prior we have that Σ⁻¹ ∼ 𝒲_k(S⁻¹, n). That is, the posterior marginal distribution for Σ⁻¹ is a Wishart distribution parametrised with scale matrix S⁻¹ and n degrees of freedom. The Wishart distribution is a multivariate generalisation of the χ² distribution. It is possible to view S ∼ 𝒲_p(Σ, n) as a distribution of SSCP-matrices associated with Σ. Inversely, one can view the inverse Wishart distribution parametrised with some SSCP matrix and n degrees of freedom as a distribution of variance-covariance matrices related to that SSCP matrix.

We then have that $\pmb\mu|\Sigma \sim \mathcal{N}_k(\pmb{\overline{y}}, \Sigma/n)$ where $\mathcal{N}_k(\pmb{\overline{y}}, \Sigma/n)$ denotes a multivariate normal distribution with mean $\pmb{\overline{y}}$ and variance Σ/n, see e.g., Gelman et al. (2013) (p. 72-73).

The posterior marginal distribution of Σ is then constructed by iterating draws of random observations from 𝒲⁻¹(S, n), each observation, Σ̂_(i), being the estimation of Σ for the ith iteration.

As previously mentioned, frequentist properties are desired if the estimations are to be used for significance testing. Berger and Sun (2008) investigated convergence to frequentist estimates for bivariate normal distributions using various priors. They showed that the “standard theory” prior produced converging estimates of σ₁ and σ₂ but that the convergence was not as good for ρ and differences in means. As a “general purpose” prior they instead recommend using $f(\pmb\mu, \Sigma^{-1}) \propto \frac{1}{\sigma_1\sigma_2(1-\rho^2)}$. Posteriors with this prior is constructed by using rejection sampling. Random draws from an inverse Wishart on n − 1 degrees of freedom are generated, but only a subset of the the observations are retained.

Crawford, Garthwaite, and Ryan (2011) considered this prior but noticed that it gave rise to too narrow credible intervals, overestimating the certainty of the estimations. However, they showed that if the sample was treated as being of size n − 1 when estimating Σ, i.e. so that the generated intervals were somewhat more conservative, frequentist properties were retained. This “calibrated” prior is recommended by Crawford, Garthwaite, and Ryan (2011), and they refer to it as such.

To construct the posterior using the calibrated prior the following steps are iterated:

1. Since we are reducing the sample size we must put S^* = (n − 2)S/(n − 1) to retain the unbiased property of S/(n − 1) as an estimate of Σ.

2. Draw a random observation from 𝒲⁻¹(S^*, n − 2). This draw is denoted: $$ \hat{\Sigma} = \begin{bmatrix} \hat{\sigma}^2_{1} & \hat{\sigma}_{12} \\ \hat{\sigma}_{12} & \hat{\sigma}^2_{2} \end{bmatrix} $$ and $$ \hat{\rho}= \frac{\hat{\sigma}_{12}}{\sqrt{\hat{\sigma}^2_{1}\hat{\sigma}^2_{2}}}. $$

3. If u is a random draw from a uniform distribution u ∼ U(0, 1), Σ̂ is only accepted as the estimation of Σ if u² ≤ 1 − ρ̂², if not we iterate the procedure until we have an accepted draw of Σ̂. Once accepted we have an estimation of Σ for this iteration and denote the draw:

Even though this latter calibrated prior is recommended, some might argue that frequentist properties are of no concern and therefore choose the “standard theory” prior. But regardless of the prior chosen, when we have an estimate of Σ, the following steps are iterated to obtain an estimate of p, the percentage of controls expected to exhibit a more extreme difference between the variates than the case.

1. Generate two draws from a standard normal distribution and denote them z_r1 and z_r2. Find the lower triangular matrix T such that TT′ = Σ̂_(i) using Choleksy decomposition. The estimate of μ for this iteration is then: See for example Gelman et al. (2013) (p. 580).

2. We now have estimates of both μ and Σ and calculate p as if these were the true population parameters. The method to do so differ slightly depending on whether a standardised or unstandardised test is required. For an unstandardised test, the size of the difference score is calculated as:
   More commonly we would want to calculate the difference score from standardised variates. If this is the case, put: and

3. The estimate of p, p̂_(i), for this iteration is then the tail area of a standard normal distribution more extreme than z_(i)^*.

Iterate the procedure a large number of times and save the estimations for each iteration. This will yield a distribution of p̂ and ẑ^*, the quantiles of which again are used to calculate the boundaries of the credible intervals for p and Z_DCC. The point estimate of the former is the mean of the estimated distribution, while the point estimate of the latter is that of (@ref(eq:PJ)).

Example usage in singcar

Testing for a deficit using the Bayesian test of deficit, call BTD():

BTD(case = caseA, controls = contA, alternative = "less",
     iter = 10000, int_level = 0.95)

## 
##  Bayesian Test of Deficit
## 
## data:  Case = 0.03, Controls (m = 0.16, sd = 0.08, n = 28)
## df = 27, p-value = 0.04787
## alternative hypothesis: true diff. between case and controls is less than 0
## sample estimates:
##   Std. case score (Z-CC), 95% CI [-2.35, -1.15] 
##                                       -1.754857 
## Proportion below case (%), 95% CI [0.94, 12.57] 
##                                        4.787144

As can be seen, compared to TD() this function has the additional argument iter, which indicates the number of iterations to be used for the posterior approximation. This number should be based on desired accuracy of the analysis.

Testing for a dissociation using the BSDT (recommended over RSDT), call:

BSDT(case_a = caseA, case_b = caseB, controls_a = contA, controls_b = contB,
      alternative = "two.sided", iter = 10000, unstandardised = FALSE,
      calibrated = TRUE)

## 
##  Bayesian Standardised Difference Test
## 
## data:  Case A: 0.03, B: 0.10, Ctrl. A (m, sd): (0.16,0.08), B: (0.18,0.10)
## df = 26, p-value = 0.3212
## alternative hypothesis: true difference between tasks is not equal to 0
## sample estimates:
##                  Std. case score, task A (Z-CC) 
##                                      -1.7548574 
##                  Std. case score, task B (Z-CC) 
##                                      -0.7836956 
## Std. discrepancy (Z-DCC), 95% CI [-1.69, -0.42] 
##                                      -1.0647889 
## Proportion below case (%), 95% CI [4.60, 33.67] 
##                                      16.0600481

This function has several additional arguments compared to RSDT(). Instead of having a separate function for an unstandardised Bayesian difference test one can specify whether to standardise the variates or not within the function. If unstandardised is set to TRUE, the results from BSDT() will converge to that of UDT(). To use the “standard theory” prior instead of the calibrated, set calibrated to FALSE.

Tests using Bayesian regression to allow for covariates

It is seldom possible to design perfect experiments. This is especially true in psychological science where causal relationships can be confounded by a great variety of factors. When comparing a case to a control sample it is therefore important to reduce noise by matching the control group to the case as well as possible on variates beyond the ones of interest, such as age or education level, making the collection of control samples cumbersome.

However, this matching of covariates can instead be achieved statistically. Crawford, Garthwaite, and Ryan (2011) describe methods where they extend the tests presented in @ref(sec22) to allow for the inclusion of covariates using Bayesian regression. These methods thus allow you to compare the case’s score on the task(s) of interest against the control population conditioned on them having the same score as the case on the covariate. This both reduces noise and alleviate the collection of control data. But because one degree of freedom is lost for each included covariate it is advisable to only include covariates with a correlation to the variate(s) of interest larger than 0.3, as a rule of thumb (Crawford, Garthwaite, and Ryan 2011).

We assume that the variates of interest are normally distributed, but no assumption of the distribution of the covariates is made. The procedure to get an estimate of p in the presence of covariates is presented below for the general case, that is either when we are interested in a deficit or a dissociation. The main difference between the methods is the recommended prior.

From a sample of size n we obtain measurements on k = 1, 2 variates of interest an m number of covariates, denoted X = (X₁, …, X_m).

We wish to estimate the regression coefficients for each covariate on the variate(s) of interest: $$ \boldsymbol{B} = \begin{bmatrix} \boldsymbol{\beta}_1 & \cdots & \boldsymbol{\beta}_k \end{bmatrix}. $$ Here B is an m + 1 × 1, 2 matrix where each column contains the regression coefficients for the ith variate of interest, the first row being the intercept(s). Further, we wish to estimate the covariance matrix or variance for the variate(s) of interest conditioned upon the covariates, we assume these statistics not to vary with the covariates. $$ \Sigma \ | \ \boldsymbol{X} = \begin{bmatrix} \sigma^2_1 \ | \ \boldsymbol{X} & \rho\sigma_1\sigma_2 \ | \ \boldsymbol{X} \\ \rho\sigma_1\sigma_2 \ | \ \boldsymbol{X} & \sigma^2_2 \ | \ \boldsymbol{X} \end{bmatrix} \ \text{or} \ \left[\sigma^2 | \boldsymbol{X} \right]. $$

If X is the n × (m + 1) design matrix and Y, the n × k response matrix $$ \boldsymbol{X} = \begin{bmatrix} 1 & x_{11} & \cdots & x_{1m} \\ \vdots & \vdots & \ddots & \vdots \\ 1 & x_{n1} & \cdots & x_{nm} \end{bmatrix},\ \ \boldsymbol{Y} = \begin{bmatrix} y_{11} & \cdots & y_{1k} \\ \vdots & \ddots & \vdots \\ y_{n1} & \cdots & y_{nk} \end{bmatrix} $$ the data estimates of B and Σ are $$ \boldsymbol{B}^* = (\boldsymbol{X}'\boldsymbol{X})^{-1}\boldsymbol{X}'\boldsymbol{Y} \ \text{and} \ \Sigma^* = \frac{1}{n-m-1}(\boldsymbol{Y}-\boldsymbol{X}\boldsymbol{B}^*)'(\boldsymbol{Y}-\boldsymbol{X}\boldsymbol{B}^*). $$

When k = 1 or the “standard theory” prior is used for k = 2 the posterior of Σ is an inverse Wishart with scale matrix (n − m − 1)Σ^* on n − m − 1 or n − m degrees of freedom for k = 1 and k = 2 respectively (Crawford, Garthwaite, and Ryan 2011; Tiao and Zellner 1964). Hence, to get an estimate Σ̂_(i) of Σ we generate a draw from these distributions.

For the “calibrated” prior we instead follow the procedure outlined in Section @ref(BSDT) and use S^* = (n − m − 2)S/(n − m − 1) as scale matrix to sample from an inverse Wishart on n − m − 2 degrees of freedom, S = (n − m − 1)Σ^*. The estimates are then:

   k=2.
(#eq:sigma) \end{equation}

We turn the (m + 1) × k matrix B^* into a k(m + 1) × 1 vector by concatenating the columns in B^*. B_vec^* = (β₁^*′, ..., β_k^*′)′ is the estimate of B_vec. The posterior of B_vec|Σ̂_(i) is a multivariate normal distribution with mean B_vec^* and covariance matrix Λ_(i) = Σ̂_(i) ⊗ (X′X)⁻¹, where ⊗ denotes the Kronecker product. To get an estimate of $\hat{\boldsymbol{B}}_{(i)}$ we generate a random observation from this distribution where the k(m + 1) × 1 vector of obtained values $\hat{\boldsymbol{B}}^*_{\text{vec(i)}} =(\hat{\boldsymbol{\beta}}'_{1(i)}, ...,\hat{\boldsymbol{\beta}}'_{k(i)})'$ is turned back into the m + 1 × k matrix $\hat{\boldsymbol{B}}_{(i)} = (\hat{\boldsymbol{\beta}}_{1(i)}, ...,\hat{\boldsymbol{\beta}}_{k(i)})$. If x^* is a vector of the case values on the covariates we obtain estimates of the conditional case scores on the task of interest as so: $$ \hat{\boldsymbol{\mu}}_{(i)} = \hat{\boldsymbol{B}}_{(i)}\boldsymbol{x}^*. $$

With these values we can now estimate the conditional effect size of the case. Crawford, Garthwaite, and Ryan (2011) denote these effect measures Z_CCC and Z_DCCC depending on whether k = 1 or k = 2. As the names indicate they are related to Z_CC (@ref(eq:zcc)) and Z_DCC (@ref(eq:PJ)), but calculated from statistics conditioned upon the covariates, hence the extra C in the subscript.

If y₁^* is the case score on a variate Y₁, we would calculate the conditional deficit (Z_DCCC) of the case on variate Y₁ like so: If y₂^* is the case score on variate Y₂ and we are interested in a dissociation between Y₁ and Y₂, we denote the two conditional means on these variates μ̂_1(i) and μ̂_2(i). Z_DCCC is then calculated like so: The conditional standard deviations ŝ_1(i) and ŝ_2(i) are obtained from Σ̂_(i) in (@ref(eq:sigma)), the conditional correlation between the variates (ρ̂_12(i)) is given by $$\hat{\rho}_{12(i)} = \frac{\hat{s}_{12(i)}}{\sqrt{\hat{s}^2_{1(i)}\hat{s}^2_{2(i)}}}.$$

To get estimates of p we find the tail area under the standard normal distribution with more extreme values on Ẑ_CCC(i) or Ẑ_DCCC(i), depending on the problem and alternative hypothesis specified. If testing for a defict we would for example have: p̂_(i) = Φ(Ẑ_CCC(i)).

Iterate the procedure a large number of times and save the estimations for each iteration. The mean of the resultant distribution of p̂ is taken at the point estimate of p. Point estimates of Z_CCC and Z_DCCC are obtained by using (@ref(eq:zccc)) and (@ref(eq:zdccc)) with the conditional means, standard deviations and correlation calculated directly from the control sample. The 1 − α credible interval boundaries for each estimate is obtained from the quantiles of the estimated distributions as described for the previous methods.

Example usage in singcar

To specify a Bayesian test of deficit with one or more covariates, call:

BTD_cov(case_task = caseA, case_covar = caseAGE, control_task = contA,
         control_covar = contAGE, alternative = "less", int_level = 0.95,
         iter = 10000, use_sumstats = FALSE)

## 
##  Bayesian Test of Deficit with Covariates
## 
## data:  Case = 0.03, Controls (m = 0.16, sd = 0.08, n = 28)
## df = 26, p-value = 0.05202
## alternative hypothesis: true diff. between case and controls is less than 0
## sample estimates:
##  Std. case score (Z-CCC), 95% CI [-2.31, -1.11] 
##                                       -1.749556 
## Proportion below case (%), 95% CI [1.05, 13.40] 
##                                        5.202005

The covariates should for the case be specified as a single value or vector of values containing the case scores for each covariate included. For the controls the covariates should be specified as a vector or matrix where each column represents one of the covariates.

To specify a Bayesian test for a dissociation with one or more covariates present, call:

BSDT_cov(case_tasks = c(caseA, caseB), case_covar = caseAGE, 
          control_tasks = cbind(contA, contB), control_covar = contAGE, 
          alternative = "two.sided", int_level = 0.95, iter = 10000,
          calibrated = TRUE, use_sumstats = FALSE)

## 
##  Bayesian Standardised Difference Test with Covariates
## 
## data:  Case A = 0.03, B = 0.10, Ctrl. (m, sd) A: (0.16,0.08), B: (0.18,0.10)
## df = 25, p-value = 0.33
## alternative hypothesis: true difference between tasks is not equal to 0
## sample estimates:
##                   Std. case score, task A (Z-CC) 
##                                        -1.754857 
##                   Std. case score, task B (Z-CC) 
##                                        -0.783696 
## Std. discrepancy (Z-DCCC), 95% CI [-1.67, -0.41] 
##                                        -1.064152 
##  Proportion below case (%), 95% CI [4.79, 34.18] 
##                                        16.500000

Here, the case scores as well as the control scores from variates of interest must be concatenated. The covariates are specified in the same manner as for BTD_cov(). The prior used is again specified with the argument calibrated.

If using summary statistics for these functions, the argument use_sumstats must be set to TRUE. Information on how to specify the data in such a case consult the documentation.

Example usage in R

Begin by installing and loading required packages:

install.packages(c("lme4", "lmerTest", "MASS"))
library("lme4")
library("lmerTest")
library("MASS") # For multivariate simulation

Because the example data is not suited for repeated measures analysis a simulated dataset will be used. Say that we have 30 measurements from four normally distributed variates Y_i, i = 1, 2, 3, 4. Call these variates items. They are thought to converge on some cognitive ability of interest and have the distribution: $$ \begin{pmatrix} Y_1 \\ Y_2 \\ Y_3 \\ Y_4 \end{pmatrix} \sim \mathcal{N}\left( \begin{bmatrix} 100 \\ 80 \\ 50 \\ 20 \end{bmatrix}, \begin{bmatrix} 15^2 & 108 & 78 & 30 \\ 108 & 10^2 & 12 & 15 \\ 78 & 12 & 10^2 & 25 \\ 30 & 15 & 25 & 5^2 \end{bmatrix} \right). $$ The case has incurred a deficit of two standard deviations on all four variates. These data could then be simulated as follows:

simdata <-  expand.grid(case = c(1, rep(0, 29)), item = factor(1:4))
simdata$ID <- factor(rep(1:30, 4)) 

set.seed(123) # For replicability
mu <- c(100, 80, 50, 20)  
sigma <- matrix(c(15^2, 108, 78, 30,
                  108, 10^2, 12, 15,
                  78,  12, 10^2, 25,
                  30,  15,  25,  5^2), nrow = 4, byrow = T)
dv <- mvrnorm(30, mu = mu, Sigma = sigma) # Dependent variable

dv[1, ] <- dv[1, ] + c(-30, -20, -20,  -10) # Case scores
simdata$dv <- c(dv)

The linear mixed model to test the deficit of the case would then be

model1 <- lmer(dv ~ case + (1|ID) + (1|item), data = simdata)
summary(model1)[["coefficients"]]

##              Estimate Std. Error        df   t value    Pr(>|t|)
## (Intercept)  62.15587  17.512303  3.030078  3.549269 0.037506019
## case        -25.07172   7.793911 27.999869 -3.216834 0.003263006

Where (1|ID) and (1|item) specifies the random effects structure, in this case the random intercepts for participant and item. This is approximately the same t statistic that would be produced by the test of deficit if averaging the scores for each participant:

agg_data <- rowMeans(dv)
TD(agg_data[1], agg_data[-1], alternative = "two.sided")[["statistic"]]

##        t 
## -3.21684

Now, say that we want to model a dissociation using LMM. We have collected 30 measurements from two items during two different conditions. It is thought that these two different conditions are tapping two independent cognitive abilities but that the items collected in each condition are converging on the same. Hence we have four variates Y_ij, i = 1, 2, j = 1, 2, distributed according to the multivariate distribution: where e.g. Y₁₂ is the random variable corresponding to the first item in the second condition. Given that the two conditions in fact taps different cognitive structures, the incurred impairment of the case would potentially affect performance in one of the conditions more than the other. Hence, we impose a deficit of two standard deviations on the case scores. The data can then be simulated as follows:

simdata <-  expand.grid(case = c(1, rep(0, 29)),
                        item = factor(1:2), condition = factor(1:2))
simdata$ID <- factor(rep(1:30, 4)) 
contrasts(simdata$condition) <- contr.sum(2) # Effect coding

set.seed(123) # For replicability
mu <- c(100, 80, 50, 20)  
sigma <- matrix(c(15^2, 120,    45,  7,
                  120, 10^2,    11, 15,
                  45,    11,  10^2, 40,
                  7,     15,    40, 5^2), nrow = 4, byrow = T)
dv <- mvrnorm(30, mu = mu, Sigma = sigma) # Dependent variable
dv[1, ] <- dv[1, ] + c(-30, -20, -0,  -0) # Case scores
simdata$dv <- c(dv)

The fully crossed linear mixed model used to test a dissociation for these data would then be specified like so:

model2 <- lmer(dv ~ case*condition + (condition|ID) + (1| item), data = simdata)
round(summary(model2)[["coefficients"]], 5)

##                  Estimate Std. Error       df  t value Pr(>|t|)
## (Intercept)      61.87447   12.18065  1.02161  5.07973  0.11992
## case            -17.24035    7.66804 28.00044 -2.24834  0.03261
## condition1       28.08855    0.98494 28.00011 28.51813  0.00000
## case:condition1 -13.82757    5.39472 28.00011 -2.56317  0.01603

Comparing the p value of the interaction effect to the p value obtained from RSDT() we see that that they are far from equivalent.

agg_data1 <- rowMeans(dv[ , 1:2])
agg_data2 <- rowMeans(dv[ , 3:4])
RSDT(agg_data1[1], agg_data2[1], agg_data1[-1], agg_data2[-1])[["p.value"]]

## [1] 0.06737687

However, a small simulation study conducted with the data described above revealed that the RSDT and the LMM yields comparable results regarding power and Type I error. It takes the specified impairments of the case for each variate and returns the ratio of significant interaction effects from (@ref(eq:lmm4)) divided by the total number of simulations run, as well as the same ratio obtained from either the RSDT or the BSDT. The distribution of the data is the same in this small simulation study as (@ref(eq:lmmdist)).

R> powerLMM(nsim = 1000, case_impairment = c(-30, -20, -0, -0),
+ compare_against = "RSDT", alpha = 0.05)
RSDT LMM
1 0.353 0.462

It seems like the method based on the linear mixed model has quite a good power advantage over the RSDT for data simulated from this particular distribution. To compare the error rate of two tests when can use the same ratio as for power, when there is no true effect present. In this case that would be either when the case has no incurred impairment at all or when the impairments are equal in both conditions.

R> powerLMM(nsim = 1000, case_impairment = c(-0, -0, -0, -0),
+ compare_against = "RSDT", alpha = 0.05)
RSDT LMM
1 0.032 0.039

Even though a power advantage of the mixed model was observed it seems to produce an error rate under the nominal level of 0.05 as well. If this was a consistent behaviour it would eliminate the usefulness of the RSDT. However, if the case exhibits extreme impairments in both conditions we would have another null scenario. This was discussed previously in Section @ref(sec22) where it was noted that one of the reasons Crawford and Garthwaite (2007) developed a Bayesian approach to test for dissociations was the highly increasing error rate of the RSDT, for increasing deficits on both tasks with no discrepancy between them. Say that the case has impairments of six standard deviations on both items in both conditions and let us compare the LMM to the BSDT instead of the RSDT:

R> powerLMM(nsim = 1000, case_impairment = c(-90, -60, -60, -30),
+ compare_against = "BSDT", alpha = 0.05)
BSDT LMM
1 0.053 0.61

The BSDT seems to keep the error rate almost at its nominal level while it is enormous for the LMM. Even though a more extensive simulation study with various types of data structures should be performed to solidify these findings, the preliminary results shown here indicate some limitations of using LMM for dissociations and reiterates the recommendation of Crawford and Garthwaite (2007) to use the BSDT to test for dissociations, even when using aggregated repeated measures data. That is, because the effects of brain damage can be severe I would recommend against using LMM to model dissociations.

Estimation based on confidence intervals

Reiser (2001) suggested a method to construct confidence intervals on δ². Using this method Elfadaly, Garthwaite, and Crawford (2016) propose a way to construct confidence intervals on p as well. We have: $$ F_0 = \frac{n\nu_2}{(n-1)\nu_1}\lambda_0 \sim F_{\nu_1, \nu_2}(n\lambda), $$ that is, the sample statistic F₀ has a non-central F distribution on ν₁ and ν₂ degrees of freedom and non-centrality parameter nλ. We call the value of nλ, where F₀ is the α-quantile of this distribution, L_α. We can construct the lower and upper boundaries of a confidence interval for p using the CDF for the χ²-distribution on ν₁ degrees of freedom, denoted G(.), as such:

Similarly, they suggest a point estimator of p based on the median of F_ν1, ν2(nλ). So that if we let L_0.5 be the value of nλ where F₀ is the median of F_ν1, ν2(nλ) the estimator is:

Similarly we would construct a confidence interval around λ by choosing λ_u and λ_l such that

However, this method has a drawback. It exactly matches the nominal level, i.e. all of the say 95% intervals (assuming α = 0.05) will in fact contain the true value of δ 95% of the time, proved by Wunderlich et al. (2015). But when F₀ is small enough so that $\mathbb{P}[F_{\nu_1,n- \nu_1}(0)\leq F_0] \leq 1- \frac{\alpha}{2}$, we cannot find a λ_u so that $\mathbb{P}\left[F_{\nu_1,n- \nu_1}(n\lambda_u)\leq F_0\right] =1- \frac{\alpha}{2}$ and it is instead set to 0. Similarly we set λ_l to 0 when $\mathbb{P}[F_{\nu_1,n- \nu_1}(0)\leq F_0] \leq \frac{\alpha}{2}$. But we know with full certainty that an interval of [0, 0] will not contain λ, since y^* is continuous with ℙ[y^* = μ] = 0.

This problem follows in the construction of confidence intervals for p, (@ref(eq:DCI)) as well as for the estimator p̂_D, (@ref(eq:pmed)). The median of the distribution F_ν1, ν2(nλ) decreases of course as the non-centrality parameter nλ decreases. We have that λ ≥ 0, therefore the lower limit of the median for F_ν1, ν2(nλ) is the median of of this distribution when the non-centrality parameter nλ = 0, i.e. the median of a central F distribution on ν₁ and ν₂ degrees of freedom. So, using the approach of Reiser (2001), if F₀ is smaller than this limit one set nλ to zero.

We see, however, a bias in this estimator when F₀ is small even if it is not smaller than the lower limit. Elfadaly, Garthwaite, and Crawford (2016) gives as an example of this scenario where we have ν₁ = 4, ν₂ = 20 and λ₀ = 0.4, hence n = 24 and we have $F_0 = \frac{24* 20}{23 * 4}0.4 = 2.09$. The value of the non-centrality parameter nλ where the median of F_ν1, ν2(nλ) equals 2.09 is ~5.015 and p̂_D = 1 − G₄(5.015/24) = 0.9949 or 99.49%. So when using the estimated Mahalanobis index of 0.4 we would infer that only 0.51% of the population would exhibit a lower Mahalanobis index. If instead calculating the true percentage of the population that would exhibit a smaller Mahalanobis index when 0.4 is the true case value, that is p = 1 − G₄(0.4) = 0.9825 or 98.25%, which is quite a considerable discrepancy.

Garthwaite, Elfadaly, and Crawford (2017) proposed a solution to this in the construction of confidence intervals of δ by modifying the method by Reiser (2001). They suggested to form a Bayesian credible interval instead of the frequentist confidence interval, whenever F₀ is small. Credible intervals have the more intuitive interpretation that they contain the true value of interest with a probability of 1 − α. As such this interval will not give rise to unreasonable values such as [0, 0]. But to form a credible interval, a prior must be specified.

Let y be a vector of values from an observation of the control population and not necessarily the case. We denote the prior for δ when δ is the Mahalanobis distance of the case as π(δ) and when δ is the distance between the mean of the controls and a randomly selected observation as ψ(δ). A priori δ should be larger when it is calculated from the observation of the case than a randomly selected observation of the controls. Garthwaite, Elfadaly, and Crawford (2017) assume that for any δ̂² and any specified quantile π_q(δ|δ̂) ≥ ψ_q(δ|δ̂), where π_q(.) and ψ_q(.) are the q:th quantile of the prior distributions.

Note that they do not aim to specify π(δ), instead they want to put limits on the quantiles of the resultant posterior through ψ(δ). This is because we know neither π(δ) nor π_q(δ|δ̂) but it is possible to establish ψ(δ|δ̂) and the α/2 and 1 − α/2 quantiles of this distribution are then the upper and lower boundaries of a 1 − α credible interval of δ.

The logic behind this modification is intuitive: treat the case like a randomly drawn observation from the control distribution whenever the case observation is more similar to the average of the controls than the majority of the control observations. That is, we take the boundaries that are lowest of the confidence and credible intervals.

Again, we set λ = δ² and λ̂ = δ̂² and let ψ^*(λ) be the prior corresponding to ψ(δ). This is a χ² distribution on ν₁ degrees of freedom. If we set m = n/(n − 1) the likelihood of λ, L(λ; λ̂), is given by Combining the prior with the likelihood, ψ^*(λ)L(λ; λ̂) yields the posterior ψ^*(λ|λ̂):

Here B(., .) and Γ(.) are the beta and gamma function respectively, c is the norming constant which is obtained through numeric integration of @ref(eq:post) over 0 < λ < ∞.

Sums over infinite ranges are problematic for practical computation. One way do deal with this is to set reasonable endpoints instead of infinity. For the implementation in singcar the endpoint of the infinite sum was set by noting that the denominator in (@ref(eq:post)) goes to infinity faster than the numerator. In fact, in R and many other languages r! = ∞ for r > 170, so a reasonable boundary for the infinite sum would then be 170, because anything divided by infinity is 0. However, for some parameter values infinity will be approached in the numerator before it is approached in the denominator. Infinity divided by any real value is still infinity but infinity divided by infinity is undefined. Therefore, to be able to sum, we remove any infinite or undefined value if present. As such, the integral over the sum in (@ref(eq:post)) will always be convergent and the normalising constant c can be calculated.

singcar uses the function integrate in the base R package stats to do this. This numerical integration procedure, based on routines from the QUADPACK package (Piessens et al. 2012), often yields correct integrals. However, it uses subdivision of the interval to be integrated and when integrating to infinity the area of interest will be undersampled. Hence, in singcar we set: $$ c = \int_0^{0.5\sqrt{\hat\lambda}}\psi^*(\lambda | \hat{\lambda}) d\lambda + \int_{0.5\sqrt{\hat\lambda}}^{\sqrt{\hat\lambda}}\psi^*(\lambda | \hat{\lambda}) d\lambda + \int_{\sqrt{\hat\lambda}}^{1.5\sqrt{\hat\lambda}}\psi^*(\lambda | \hat{\lambda}) d\lambda + \int_{1.5\sqrt{\hat\lambda}}^{2\sqrt{\hat\lambda}}\psi^*(\lambda | \hat{\lambda}) d\lambda + \int_{2\sqrt{\hat\lambda}}^{\infty}\psi^*(\lambda | \hat{\lambda}) \ d\lambda. $$ This ensures heavier sampling around $\sqrt{\hat\lambda}$, which of course is the area we are most interested in. Further, for large values of λ̂ some discontinuities will be introduced in the posterior which can be problematic for some numerical integration techniques as well but the above subdivision will alleviate this problem somewhat.

The q:th quantile of the posterior ψ^*(λ|λ̂), λ_q is then found by an optimisation search procedure (nlminb in stats) that aims to minimise |∫_λ0^λ_qψ^*(λ|λ̂)dλ − q|.

To get an estimate of p we use the same procedure to find the median of the posterior, λ_0.5. Elfadaly, Garthwaite, and Crawford (2016)’s modified estimator of p, p̂_MD, is then: p̂_MD = min[p̂_D, 1 − G(λ_0.5)]. Note however that due to the fact that c will sometimes be very small and hence 1/c very large, computations of the normalised posterior can be very slow, but this mostly happens when δ̂ is so large that p̂_D would be greater than 1 − G(λ_0.5) anyway. In singcar it is therefore recommended to use p̂_D but to calculate p̂_MD if the Mahalanobis distance of the case seems suspiciously small.

The difference between the estimators p̂_D and p̂_MD is small and Elfadaly, Garthwaite, and Crawford (2016) show that their bias and means square error are almost identical, but they recommend p̂_MD. This is because, as δ̂ → 0, p̂_D → 1 much faster than it, by common sense, should.

It should be noted that the implementation of p̂_MD in singcar is unstable when compared to the implementation by Garthwaite, Elfadaly, and Crawford (2017) written in C. I have not been able to find the cause of this discrepancy, but it only happens for relatively large δ̂. It is likely that the numerical integration procedures differ and that the procedure used in the implementation by Garthwaite, Elfadaly, and Crawford (2017) better handles discontinuous functions, which we see in the posterior (@ref(eq:post)) for high values of δ̂. Alternatively, the discontinuities cause problems for the optimisation routine. It is therefore recommended that both estimators are calculated, and the smallest of the two are chosen, but if that is p̂_MD then compare the result to the implementation by Garthwaite, Elfadaly, and Crawford (2017).

Example usage in singcar

Say that we want to know the case abnormality in the example data on both tasks simultaneously. To estimate this abnormality on a multivariate space, call:

MTD(case = c(caseA, caseB), controls = cbind(contA, contB),
    conf_level = 0.95, method = c("pd", "pchi", "pf", "pmd"),
    mahalanobis_dist = NULL, k = NULL, n = NULL)

As can be seen, p̂_D is the default method to estimate the abnormality. This is because the difference in results between p̂_D and p̂_MD is very small, but the difference in efficiency is huge. In addition, "pmd" in singcar is actually 1 − G(λ_0.5) rather than min[p̂_D, 1 − G(λ_0.5)], hence taking the minimum is left to the user. The reason for this was to minimise confusion. However, I realise that this could yield the opposite effect and might be changed in future updates.

The p value in the output of MTD() is the p value associated with the Hotelling’s T² statistic and the same as the abnormality estimate if the method pf is chosen. The three last arguments are only required if summary statistics are used.

Power calculators

A further capacity of singcaris that it can be used to calculate statistical power of the tests. The notion of power when comparing cases to control samples have been somewhat overlooked for this class of statistical tests. In recent work (McIntosh and Rittmo 2020), we argued that, even though power is inherently limited in this paradigm, a priori calculations are still useful for study design and interpretation in neuropsychological and other applications. Calculating power for the test of deficit is similar to calculating power for any t test and can be done analytically. Where T_n − 1(.|θ) is the cumulative distribution function for the non-central t distribution with n − 1 degrees of freedom and non-centrality parameter $\frac{y^* - \overline{y}}{\sigma \sqrt{\frac{n+1}{n}}}$ (i.e., TD, Equation @ref(eq:TD) and t_{α, n − 1} is the α quantile of the t distribution on n − 1 degrees of freedom (note that this is for a one-sided test). For the unstandardised difference test power is calculated in an analogous way by putting Equation @ref(eq:UDT) as the non-centrality parameter. Deriving power for the other functions in an analytic manner is however not possible (the RSDT is only approximately t distributed) and a Monte Carlo approach has been used for these tests. To call any power calculator in the package one simply uses the function names with _power added as a suffix.

So, for example, to calculate power for the test of deficit we call TD_power(). The expected case score and either sample size or desired power must be supplied. The mean and standard deviation of the control sample can also be specified with the arguments mean and sd. If not, they take the default values of 0 and 1 respectively so that the case score is interpreted as distance from the mean in standard deviations. A conventional α-level of 0.05 is assumed if nothing else is supplied. The alternative hypothesis can also be specified by the argument alternative: specify "less" (default) or "greater" for a one-tailed test, specify "two.sided" for a two-tailed test.

TD_power(case = 70,
         mean = 100,
         sd = 15,
         sample_size = 16,
         power = NULL,
         alternative = "less",
         alpha = 0.05,
         spec = 0.005)

## [1] 0.5819579

TD_power() can also be used to calculate required sample size for a desired level of power. For example, if we specify a desired power level of 0.6, leave sample_size to its default and let the rest of the arguments be as in the previous example, we see from the output of the function that power will not increase more than 0.5% for any additional participant after a sample size of 15. That is, the algorithm stops searching when this level of specificity has been reached and we are nearing the asymptotic maximum power for this effect size. We can increase the specificity by lowering the spec argument.

TD_power(case = 70,
         mean = 100,
         sd = 15,
         sample_size = NULL,
         power = 0.6,
         alternative = "less",
         alpha = 0.05,
         spec = 0.005)

## Power (0.578) will not increase more than 0.5% per participant for n > 15

##    n     power
## 1 15 0.5780555

Power calculators for the Bayesian tests of deficit cannot calculate required sample size. This is because they rely on simulation methods to estimate approximate power and deploying a search algorithm to find the required sample size for a given level of power would be computationally too intense. The syntax is otherwise relatively similar to that of TD_power(). For BTD_power() we have the two extra arguments nsim and iter, indicating the number of simulations used in the power function and by BTD(), respectively.

BTD_power(case = 70,
         mean = 100,
         sd = 15,
         sample_size = 15,
         alternative = "less",
         alpha = 0.05,
         nsim = 1000,
         iter = 1000)

## [1] 0.574

The only difference in syntax of BTD_cov_power() is due to the inclusion of covariates. The variate of interest must be specified as a vector where the first element gives the mean and the second the standard deviation in the argument control_task. The covariates can be specified similarly or as an m × 2 matrix where the first column gives the means of each covariate and the second column gives the standard deviations. The correlation matrix of the variates must be given as well. In the example below, power is evaluated for a test taking two covariates into account, both with a mean of 0 and a standard deviation of 1. The correlation is specified as a 3 × 3 matrix with pairwise correlations of 0.3. The default settings include only one covariate having a 0.3 correlation with the variate of interest. This function is computationally intense and hence, the number of simulations has, for the example below, been decreased.

covars <- matrix(c(0, 1,
                   0, 1), ncol = 2, byrow = TRUE)
BTD_cov_power(case = -2,
              case_cov = c(0.2, -0.6),
              control_task = c(0, 1),
              control_covar = covars,
              cor_mat = diag(3) + 0.3 - diag(c(0.3, 0.3, 0.3)),
              sample_size = 15,
              alternative = "less",
              alpha = 0.05,
              nsim = 100,
              iter = 100)

## [1] 0.48

For the difference tests one must supply the expected case scores on both variates as well as sample size. The means and standard deviations of the control sample can also be specified. If unspecified, they take on the default values of 0 and 1 respectively, so that the expected case scores are interpreted as distances from the means in standard deviations. RSDT_power(), BSDT_power() and UDT_power() additionally require an estimate of the sample correlation between the variates of interest, r_ab. If this is not specified a correlation of 0.5 is assumed by default. For BSDT_cov_power() the correlation matrix between the variates of interest and the covariates must instead be supplied (i.e., at least a 3 × 3 matrix where the first correlation is the correlation between the variates of interest).

The alternative hypothesis is by default assumed to be "two.sided" since the direction of the effect is dependent on the order of the inputs, but can be specified to be "less" or "greater" as well. The syntax is similar for all three functions but with small differences. For UDT_power() one can request required sample size for a desired power, as for TD_power(). Calculators for the Bayesian tests have the extra argument calibrated as to be able to specify the prior. BSDT_cov_power() requires input in the same format as BTD_cov_power() for both control_tasks and control_covar. The two examples below demonstrate usage for RSDT_power() and BSDT_cov_power().

RSDT_power(case_a = 70,
           case_b = 55,
           mean_a = 100,
           mean_b = 50,
           sd_a = 15,
           sd_b = 10,
           r_ab = 0.5,
           sample_size = 15,
           alternative = "two.sided",
           alpha = 0.05,
           nsim = 1000)

## [1] 0.604

cor_mat <- matrix(c(1,   0.5, 0.6,
                    0.5,   1, 0.3,
                    0.6, 0.3,   1), ncol = 3, byrow = TRUE)
BSDT_cov_power(case_tasks = c(70, 55),
               case_cov = 65,
               control_tasks = matrix(c(100, 15,
                                        50, 10), ncol = 2, byrow = TRUE),
               control_covar = c(50, 25),
               cor_mat = cor_mat,
               sample_size = 15,
               alternative = "two.sided",
               alpha = 0.05,
               nsim = 100,
               iter = 100,
               calibrated = TRUE)

## [1] 0.75