1vanderbilt university, department of biostatistics

2vanderbilt university medical center, department of psychiatry and behavioral
sciences

Please address correspondence to:
Kaidi Kang
2525 West End Ave., #1136
Department of Biostatistics
Vanderbilt University
Nashville, TN 37203
<EMAIL_ADDRESS>
Abstract

We recently proposed a robust effect size index (RESI) that is related to the
non-centrality parameter of a test statistic. RESI is advantageous over common
indices because (1) it is widely applicable to many types of data; (2) it can
rely on a robust covariance estimate; (3) it can accommodate the existence of
nuisance parameters. We provided a consistent estimator for the RESI, however,
there is no established confidence interval (CI) estimation procedure for the
RESI. Here, we use statistical theory and simulations to evaluate several CI
estimation procedures for three estimators of the RESI. Our findings show (1)
in contrast to common effect sizes, the robust estimator is consistent for the
true effect size; (2) common CI procedures for effect sizes that are non-
centrality parameters fail to cover the true effect size at the nominal level.
Using the robust estimator along with the proposed bootstrap CI is generally
accurate and applicable to conduct consistent estimation and valid inference
for the RESI, especially when model assumptions may be violated. Based on the
RESI, we propose a general framework for the analysis of effect size (ANOES),
such that effect sizes and confidence intervals can be easily reported in an
analysis of variance (ANOVA) table format for a wide range of models.

  * Key words: effect size; robust effect size index; confidence intervals; analysis of effect sizes (ANOES); non-centrality parameter; bootstrap; non-central Chi-squared distribution; non-central F distribution.

##  1 Introduction

Effect size indices are measures quantifying the strength of association
between a covariate and an outcome of interest that are unaffected by the
sample size. They play an important role in power analyses, sample size
planning and meta-analyses (Cohen,, 1988; Chinn,, 2000; Morris and DeShon,,
2002). Given the recent criticism of the misuse and misinterpretation of null
hypothesis significance testing by the American Statistical Association (ASA)
(Wasserstein and Lazar,, 2016; Wasserstein et al.,, 2019), there is an
increasing interest in seeking alternatives for communicating study findings.
Reporting effect size estimates with confidence intervals (CIs) can be an
excellent way to simultaneously communicate the strength of evidence (sample
size independent) as well as the confidence one should have in the evidence
(sample size dependent). The American Psychological Association (APA) has also
been calling for reporting of effect sizes and their CIs for almost two
decades (American Psychological Association,, 2001, 2010), but they are still
not routinely reported by psychological studies. Fritz et al., (2012) reviewed
articles published in 2009 and 2010 in the _Journal of Experimental
Psychology: General_ , and noted that less than half of the articles they
reviewed reported effect sizes and no article reported a confidence interval
for an effect size. The barriers stopping researchers from easily reporting
effect sizes along with their CIs not only lies in their unfamiliarity of
different effect size indices but also in the lack of guidance of how to
correctly estimate the CI for a specific effect size index.

We recently proposed a robust effect size index (RESI) (Vandekar et al.,,
2020), which has several advantages over previously proposed indices (Cohen,,
1988; Hedges and Olkin,, 1985; Long and Freese,, 2006; Zhang and Schoeps,,
1997; Rosenthal,, 1994) because (1) it is widely applicable to many types of
data since it is constructed from M-estimators, which are generally defined;
(2) it can rely on a robust covariance estimate; (3) it can accommodate the
existence of nuisance covariates/parameters. We also proposed a simple
consistent estimator for the RESI that is a function of the Chi-squared test
statistic (Vandekar et al.,, 2020). The RESI is defined on the basis of the
Wald test statistic – it is related to the non-centrality parameter (NCP) of
the test statistic under the alternative hypothesis, therefore, it has the
generality that it can be estimated in different types of data. Researchers
can use the RESI to report their observed effect sizes regardless of the
model. Furthermore, the studies on the same scientific topic, but using
different types of data can easily communicate their observed effect sizes
without translating between different effect size indices. While the RESI can
currently be used to report the strength of a finding through the RESI
estimate, we did not establish a CI estimation procedure for the RESI, which
quantifies the amount of certainty of the estimate.

The goal of this paper is to establish an accurate CI estimation procedure for
the RESI and establish a framework for the analysis of effect sizes (ANOES)
based on the RESI. Because the RESI is related to the NCP of the Chi-squared
statistic, an intuitive approach is to use existing methods of constructing
intervals for the NCP of a noncentral Chi-squared or F distribution (Kent and
Hainsworth,, 1995; Harlow et al.,, 2013). Here, we use statistical theory and
simulations to show that the Chi-squared CI provides low coverage for the NCP
when the variance must be estimated. In fact, the coverage gets lower with
increasing sample and effect size. Similarly, the F CI has decreasing coverage
with increasing sample size when the robust covariance estimate is used. We
use theory to show that this occurs because the variance of the estimators is
not asymptotically equivalent. As a solution, we propose several bootstrap CI
construction procedures and evaluate their coverage performance through
simulations in various scenarios. We propose a bootstrap CI for the robust
estimator of the RESI that is generally accurate and applicable, even when
model assumptions are violated. Based on the RESI, we propose a general
framework for ANOES, such that effect sizes and CIs can be easily reported in
an analysis of variance (ANOVA) table format. We use this framework to study
the effect of early psychosis and schizophrenia on relational memory for
illustrative purpose. Our early stage RESI R package is available to install
through github (https://github.com/statimagcoll/RESI).

##  2 Statistical theory

###  2.1 Estimators for the Robust Effect Size Index (RESI)

In this section, we define the RESI and describe three estimators for the
parameter. Let $W=(W_{1},\ldots,W_{n})$ denote the full dataset, where $W_{i}$
is independent of $W_{j}$ for all $i\neq j$. Assume
$\theta=(\alpha,\beta)\in\Theta\subset\mathbb{R}^{m}$ are model parameters,
where $\alpha\in\mathbb{R}^{m_{0}}$ is a vector of nuisance parameters,
$\beta\in\mathbb{R}^{m_{1}}$ is a vector of target parameter, and
$m=m_{0}+m_{1}$. We assume
$\Psi(\theta;W)=n^{-1}\sum_{i=1}^{n}\psi(\theta;W)\in R$ is an estimating
equation, where $\psi$ is a known function and $\Psi$ can be maximized to
obtain the M-estimator $\hat{\theta}$

$\hat{\theta}=\operatorname*{arg\,max}_{\theta^{*}\in\Theta}\Psi(\theta^{*};W).$

If $\Psi$ is a likelihood function then $\hat{\theta}$ corresponds to the
maximum likelihood estimator. Denote the target value of the parameter
$\theta=\arg\max_{\theta^{*}}\mathbb{E}\Psi(\theta^{*};W)$.

We define the components of the asymptotic covariance of
$\sqrt{n}(\hat{\theta}-\theta)$:

$\mathbf{J}_{jk}(\theta)=-\left.\lim_{n\to\infty}\mathbb{E}\frac{\partial^{2}\Psi(\theta^{*},W)}{\partial\theta^{*}_{j}\partial\theta^{*}_{k}}\right\rvert_{\theta}$
$\mathbf{K}_{jk}(\theta)=\left.\lim_{n\to\infty}\mathbb{E}\frac{\partial\Psi(\theta^{*},W)}{\partial\theta^{*}_{j}}\frac{\partial\Psi(\theta^{*},W)}{\partial\theta^{*}_{k}}\right\rvert_{\theta}$

Under mild conditions (Van der Vaart,, 2000; Vandekar et al.,, 2020), the
variance of $\sqrt{n}(\hat{\theta}-\theta)$ is

$\Sigma_{\theta}=\mathbf{J}(\theta)^{-1}\mathbf{K}(\theta)\mathbf{J}(\theta)^{-1}$
(1)

If $\Psi$ is a correctly specified likelihood function, then
$\mathbf{J}(\theta)=\mathbf{K}(\theta)$, and the asymptotic covariance matrix
of $\sqrt{n}(\hat{\theta}-\theta)$ is $\mathbf{J}(\theta)^{-1}$.

We defined the RESI from the test statistic for $H_{0}:\beta=\beta_{0}$, where
$\beta_{0}$ is a vector-valued reference point. We, previously, suggested that
the typical Wald-style statistic for the test of the null hypothesis follows a
Chi-squared distribution with $m_{1}$ degrees of freedom and non-centrality
parameter $n(\beta-\beta_{0})^{T}\Sigma_{\beta}(\beta-\beta_{0})$,

$\displaystyle
T^{2}=n(\hat{\beta}-\beta_{0})^{T}\Sigma_{\beta}^{-1}(\hat{\theta})(\hat{\beta}-\beta_{0})\sim\chi^{2}_{m_{1}}\\{n(\beta-\beta_{0})^{T}\Sigma_{\beta}^{-1}(\theta)(\beta-\beta_{0})\\},$
(2)

where the $\hat{\beta}$ is the estimated value of $\beta$ and
$\Sigma_{\beta}(\theta)$ is the asymptotic covariance matrix of
$\sqrt{n}(\hat{\beta}-\beta)$.

The RESI was defined as the square root of the component of the NCP that is
due to the deviation of $\beta$ from the reference value,

$\displaystyle
S_{\beta}=\sqrt{(\beta-\beta_{0})^{T}\Sigma_{\beta}^{-1}(\theta)(\beta-\beta_{0})}.$

The estimator for the RESI is defined as (Vandekar et al.,, 2020)

$\hat{S}_{\beta}=\left\\{\max\left[0,\frac{T^{2}-m_{1}}{n}\right]\right\\}^{\frac{1}{2}}$
(3)

where $T^{2}$ is as defined in (2). This estimator was derived by setting the
observed statistic $T^{2}$ equal to the expected value of the non-central Chi-
squared distribution and solving for $S_{\beta}$

$\displaystyle\hat{S}_{\beta}^{2}$ $\displaystyle=\frac{T^{2}-m_{1}}{n}.$

Because $S_{\beta}$ must be nonnegative, the estimator (3) has lower mean
square error (Vandekar et al.,, 2020; Neff and Strawderman,, 1976; Kubokawa et
al.,, 1993).

There are several ways of constructing the statistic $T^{2}$ and
$\hat{S}_{\beta}$ through the estimation of $\Sigma_{\beta}(\hat{\theta})$
(White,, 1980; MacKinnon and White,, 1985; Long and Ervin,, 2000). First, the
matrices $\mathbf{J}$ and $\mathbf{K}$ are estimated by

$\displaystyle\hat{\mathbf{J}}_{jk}(\hat{\theta})$
$\displaystyle=-\left.n^{-1}\sum_{i=1}^{n}\frac{\partial^{2}\psi(\theta^{*},W)}{\partial\theta^{*}_{j}\partial\theta^{*}_{k}}\right\rvert_{\hat{\theta}}$
(4) $\displaystyle\hat{\mathbf{K}}(\hat{\theta})$
$\displaystyle=\left.n^{-1}\sum_{i=1}^{n}\frac{\partial\psi(\theta^{*},W)}{\partial\theta^{*}_{j}}\frac{\partial\psi(\theta^{*},W)}{\partial\theta^{*}_{k}}\right\rvert_{\hat{\theta}}$
$\displaystyle\hat{\Sigma}(\hat{\theta})$
$\displaystyle=\hat{\mathbf{J}}(\hat{\theta})^{-1}\hat{\mathbf{K}}(\hat{\theta})\hat{\mathbf{J}}(\hat{\theta})^{-1}.$
(5)

Then we can estimate the RESI using three versions of test statistics:

  1. 1.

Oracle test statistics: when the true covariance of
$\sqrt{n}(\hat{\beta}-\beta)$ is known,

$T_{(o)}^{2}=n(\hat{\beta}-\beta_{0})^{T}\Sigma_{\beta}^{-1}(\theta)(\hat{\beta}-\beta_{0}),$

where $\Sigma_{\beta}(\theta)$ is obtained from the block diagonal of (1)
corresponding to $\beta$. This is called an “oracle” statistic because it
depends on the true covariance matrix, which is not known in practice.

  2. 2.

Parametric test statistics: when we believe the working model is correctly
specified,

$T_{(p)}^{2}=n(\hat{\beta}-\beta_{0})^{T}\hat{\mathbf{J}}_{\beta}(\hat{\theta})(\hat{\beta}-\beta_{0}),$

where $\hat{\mathbf{J}}_{\beta}(\hat{\theta})$ is obtained from the block
diagonal component of (4) corresponding to $\beta$.

  3. 3.

Robust test statistics: when the working model is not assumed to be correct,

$T_{(r)}^{2}=n(\hat{\beta}-\beta_{0})^{T}\hat{\Sigma}_{\beta}^{-1}(\hat{\theta})(\hat{\beta}-\beta_{0}),$

where $\hat{\Sigma}_{\beta}(\hat{\theta})$ is obtained from the block diagonal
component of (5) corresponding to $\beta$.

With the oracle, parametric, or robust test statistics, we can derive the
corresponding oracle, parametric or robust estimator for the RESI using (3).
These different versions of the estimator for $S_{\beta}$ have different
sampling distributions depending on which test statistic is used.

###  2.2 Chi-squared CIs underestimate the variance

When the covariance of $\sqrt{n}(\hat{\beta}-\beta)$ is known, the test
statistics follows the non-central Chi-squared distribution (2), where its
variance depends on the true effect size. However, in practice, the covariance
needs to be estimated. In this section, we will show that the large sample
distribution of the test statistic deviates from the theoretical non-central
Chi-squared distribution when an estimator for the covariance is used. As a
result, the CIs constructed using Chi-squared distribution (Kent and
Hainsworth,, 1995) will fail to provide nominal coverage. This is in contrast
to the null case, where the asymptotic distribution of the Chi-square
statistic is valid whether or not $\Sigma_{\beta}(\theta)$ is known or
estimated. To illustrate this problem here, we compare the asymptotic
equivalence of the variance of the oracle, parametric, and robust estimators
in a linear model. Two functions $f(n)$, $g(n)$ are said to be asymptotically
equivalent if $\lim_{n\to\infty}\frac{f(n)}{g(n)}=1$. We need asymptotic
equivalence here, because a standard asymptotic approach cannot be used; under
the alternative, the mean and variance of the test statistics depend on the
sample size.

Throughout this section, we assume $Y=X\beta+\epsilon$, where $X$ is full rank
and $\epsilon_{i}$ are independent with mean 0 and variance $\sigma^{2}_{i}$.
Note that $\sigma^{2}_{i}$ may not be equal. The ordinary least squares (OLS)
estimator of $\beta$ is $\hat{\beta}_{OLS}=(X^{T}X)^{-1}X^{T}Y$.

##### Oracle Estimator:

Assuming known covariance, the test statistic (2) is approximately central
Chi-squared by the central limit theorem under the null. Under the
alternative, a similar argument suggests its approximation to a non-central
Chi-squared distribution.

$T_{(o)}^{2}\sim\chi^{2}_{m_{1}}(nS^{2}_{\beta}).$

The expectation and variance of the oracle test statistics are:

$\displaystyle\mathbb{E}T_{(o)}^{2}$ $\displaystyle=m_{1}+nS^{2}_{\beta}$
$\displaystyle Var(T_{(o)}^{2})$ $\displaystyle=2(m_{1}+2nS^{2}_{\beta}).$ (6)

Thus, the estimator (3) is consistent and the variance of $T_{(o)}$ is linear
in the sample size.

##### Parametric Estimator:

Assuming homoskedasticity (i.e., $\sigma^{2}_{i}=\sigma^{2}$, $\forall i$) and
normality of $Y$ makes finding the distribution of the test statistic
tractable. The covariance matrix of $\hat{\beta}_{OLS}$ can be estimated as
$(X^{T}X)^{-1}\hat{\sigma}^{2}$, where
$\hat{\sigma}^{2}=(n-m)^{-1}Y^{T}(I-H)Y$ and $H$ is the hat matrix
$H=X(X^{T}X)^{-1}X^{T}$. In this situation, the parametric version of test
statistic for $H_{0}:\beta=\beta_{0}$ is

$\displaystyle T_{(p)}^{2}$
$\displaystyle=n(\hat{\beta}_{OLS}-\beta_{0})^{T}\hat{\Sigma}_{\beta_{OLS}}^{-1}(\hat{\beta}_{OLS}-\beta_{0})$
$\displaystyle=\hat{\sigma}^{-2}(\hat{\beta}_{OLS}-\beta_{0})^{T}X^{T}X(\hat{\beta}_{OLS}-\beta_{0})$
(7)

Additional variability is introduced into the test statistic when plugging in
the estimated covariance of $\hat{\beta}_{OLS}$ and, consequently, the
parametric test statistic does not follow the non-central Chi-squared
distribution (2). For linear regression models, it can be shown that the test
statistics divided by its own degrees of freedom $m_{1}$ follows the non-
central F-distributions (see Appendix)

$T_{(p)}^{2}/m_{1}\sim F(m_{1},n-m;nS^{2}_{\beta}).$ (8)

By plugging the parametric test statistic (2.2) into (3), we can derive a
parametric version of the estimator for the RESI. The expectation and variance
of the parametric test statistic are derived from the moments of the
F-distribution,

$\displaystyle\mathbb{E}T_{(p)}^{2}$
$\displaystyle=\frac{(n-m)(m_{1}+nS^{2}_{\beta})}{n-m-2}$ $\displaystyle
Var\left(T_{(p)}^{2}\right)$
$\displaystyle=2\frac{(m_{1}+nS^{2}_{\beta})^{2}+(m_{1}+2nS^{2}_{\beta})(n-m-2)}{(n-m-2)^{2}(n-m-4)}\left(n-m\right)^{2}$

$\mathbb{E}T_{(p)}^{2}$ and $\mathrm{Var}\left(T_{(p)}^{2}\right)$ are
asymptotically equivalent to

$\displaystyle\mathbb{E}T_{(p)}^{2}$ $\displaystyle\simeq
m_{1}+nS^{2}_{\beta}$ $\displaystyle\mathrm{Var}\left(T_{(p)}^{2}\right)$
$\displaystyle\simeq
2(m_{1}+2nS^{2}_{\beta})+2\left(2m_{1}S^{2}_{\beta}+nS^{4}_{\beta}\right)$ (9)
$\displaystyle=\mathrm{Var}\left(T_{(o)}^{2}\right)+2\left(2m_{1}S^{2}_{\beta}+nS^{4}_{\beta}\right).$

This parametric test statistic, which is simplest unknown variance case, has a
variance that is not asymptotically equal to the variance of the oracle test
statistic, which follows the non-central Chi-squared distribution. The
consequence of this fact is that a CI using the Chi-squared distribution will
have lower than nominal coverage for the NCP and the effect sizes that are a
function of the NCP. In fact, the CI will become increasingly inaccurate as
the effect size or sample size get larger. We demonstrate this with
simulations in Section 3.2.

##### Robust Estimator:

When there is a suspected unknown heteroskedasticity, a robust version of
covariance estimator of $\hat{\beta}_{OLS}$ can be applied instead,

$\displaystyle\hat{\mathrm{V}}\mathrm{ar}(\hat{\beta}_{OLS})=(X^{T}X)^{-1}X^{T}\hat{G}X(X^{T}X)^{-1}$

where $\hat{G}_{ij}=\frac{e_{i}^{2}}{(1-h_{i})^{2}}$ if $i=j$ and
$\hat{G}_{ij}=0$ if $i\neq j$, $e_{i}=Y_{i}-X_{i}\hat{\beta}$ is the $i$th
residual, and $h_{i}$ is the $i$-th elements on the diagonal of the hat matrix
$H$, i.e., $h_{i}=x_{i}(X^{T}X)^{-1}x_{i}^{T}$ (Long and Ervin,, 2000). Many
versions of this robust “sandwich” covariance estimator have been proposed,
this version is discussed by Long and Ervin, (2000) as a jackknife
approximation.

The robust version of the test statistic in this case is

$\displaystyle T_{(r)}^{2}$
$\displaystyle=n(\hat{\beta}_{OLS}-\beta_{0})^{T}\hat{\Sigma}_{\beta}^{-1}(\hat{\theta})(\hat{\beta}_{OLS}-\beta_{0})$
$\displaystyle=(\hat{\beta}_{OLS}-\beta_{0})^{T}(X^{T}X)(X^{T}\hat{G}X)^{-1}(X^{T}X)(\hat{\beta}_{OLS}-\beta_{0}).$
(10)

Then the robust version of the estimator for RESI can be constructed by
plugging (2.2) into (3).

Similar to the situation of parametric test statistic, the covariance of
parameter estimator $\hat{\beta}_{OLS}$ is estimated. The estimator is still
consistent

$\displaystyle\mathbb{E}T_{(r)}^{2}$
$\displaystyle=\mathrm{tr}\left(\left(n^{-1}\hat{\Sigma}_{\hat{\beta}_{OLS}}\right)^{-1}n^{-1}\Sigma_{\hat{\beta}_{OLS}}\right)+n(\mathbb{E}\hat{\beta}_{OLS}-\beta_{0})^{T}\hat{\Sigma}_{\hat{\beta}_{OLS}}^{-1}(\mathbb{E}\hat{\beta}_{OLS}-\beta_{0})$
$\displaystyle\simeq\mathrm{tr}\left(\Sigma_{\hat{\beta}_{OLS}}^{-1}\Sigma_{\hat{\beta}_{OLS}}\right)+n(\beta-\beta_{0})^{T}\Sigma_{\hat{\beta}_{OLS}}^{-1}(\beta-\beta_{0})$
$\displaystyle=\mathrm{tr}\left(\mathbf{I}_{m_{1}}\right)+nS^{2}$
$\displaystyle=m_{1}+nS^{2}=\mathbb{E}T^{2}_{(o)},$

but using the robust covariance estimator increases variance of the robust
test statistic above the non-central Chi-squared and the F distributions. We
found an analytical form for the variance of $T_{(r)}^{2}$ to be intractable,
so estimated it using simulations.

Figure 1 shows the simulated variance of each test statistic as a linear
function of $n$. When the effect size is large (e.g., $S=1$) and with a fixed
covariate, it’s obvious that the variance of oracle and parametric statistics
equals the variance of non-central Chi-squared and F distributions,
respectively. This is expected since these two distributions are the
asymptotic distributions of these two statistics, respectively (Section 2.2).
Because the test statistic using the robust covariance has larger variance
than the oracle and parametric statistics, the robust test statistic is not
non-central Chi-squared or non-central F distributed. This implies that the
CIs constructed using either of these two distributions will not have accurate
coverage for the robust RESI estimator and that the coverage will get worse
with increasing effect size or sample size. When the covariate is random, the
variance of all three statistics further grows, and both of the oracle and
parametric statistics deviate from the non-central Chi-squared and F
distributions. This implies that when the covariate(s) is random by designs
(such as observational studies), neither of these two distributions will
provide CIs with accurate coverage for either of these three estimators. While
we only studied linear models in this section, we expect that the robust test
statistic will not follow the Chi-squared or F distribution under the
alternative in general.

Figure 1: A comparison of simulated variance of the three test statistics from
a simple linear regression model with a binary covariate and assumptions of
homoskedasticity and symmetric errors. The grey lines are the variance of
theoretical non-central F or Chi-squared distribution using formulas (2.2) and
(9). The black lines are the simulated variance of test statistics with random
or fixed covariate. With fixed covariate, the variance of oracle and
parametric statistics equals theoretical variance of the non-central Chi-
squared or F distribution, respectively. With random covariates, the true
variance of the test statistic deviates from the theoretical value. $S$
denotes the RESI, Chisq and F denote the non-central Chi-squared and F
distributions, respectively.

##  3 Methods

###  3.1 Confidence Intervals for RESI

Here, we discuss several potential procedures to construct CIs using non-
central Chi-squared and F distributions and bootstrapping.

####  3.1.1 Theoretical methods: using non-central distributions

As discussed in the previous section, the test statistic may follow a non-
central Chi-squared or F distribution with a non-centrality parameter (NCP)
not equal to 0. Since the squared RESI was defined as the NCP divided by
sample size $n$, there is a one-to-one relationship between the NCP and RESI.
If a CI for the NCP can be constructed, the CI for RESI can be derived from
the CI for NCP. The confidence interval construction for NCP has been
discussed for non-central Chi-squared and F distributions (Kent and
Hainsworth,, 1995; Harlow et al.,, 2013).

Suppose the test statistics from the sample $T_{obs}^{2}$ is observed and its
degrees of freedom is $m_{1}$. Let $F(T_{obs}^{2};m_{1},\lambda)$ denote the
cumulative distribution function (CDF) of the non-central Chi-squared
distribution the test statistic $T^{2}_{(o)}$ follows. Then
$F(T_{obs}^{2};m_{1},\lambda$) is a monotonic and strictly decreasing function
of the NCP, $\lambda$.

The lower bound ($\ell$) and upper bound ($u$) of central $(1-\alpha)\times
100\%$ CI for $\lambda$ can be chosen as the values that satisfy the
equalities (Kent and Hainsworth,, 1995)

$\displaystyle F(T_{obs}^{2};m_{1},\lambda=\ell)=1-\alpha/2$ $\displaystyle
F(T_{obs}^{2};m_{1},\lambda=u)=\alpha/2$

For the non-central F distribution, the procedure to construct a CI for its
NCP is similar, except there are 2 degree of freedom parameters. The lower
($\ell$) and upper ($u$) bounds for the NCP of the F distribution can be
determined from the observed statistic, $T_{obs}^{2}$, with degrees of freedom
being $m_{1}$ and $(n-m)$ are chosen to satisfy

$\displaystyle F(T_{obs}^{2};m_{1},n-m,\lambda=\ell)=1-\alpha/2$
$\displaystyle F(T_{obs}^{2};m_{1},n-m,\lambda=u)=\alpha/2$

We provide code to compute these intervals in the RESI R package
(https://github.com/statimagcoll/RESI). Our code is adapted from Kent and
Hainsworth, (1995) and Harlow et al., (2013).

####  3.1.2 Resampling methods: bootstrap CIs

As we discussed in the previous section, the distribution of the test
statistic may deviate from the theoretical non-central Chi-squared or F
distribution if the covariance must be estimated and/or the covariate(s) is
random instead of fixed by designs. In this situation, the CIs built through
non-central distributions may have bad coverage of the true effect size.
Bootstrapping is a good alternative that can be used to approximate the actual
distribution of the test statistic and can be used to construct the CIs for
test statistic and the estimated effect size (Efron,, 1979; Hall,, 1992).

In the body of this paper, we consider the standard nonparametric bootstrap by
sampling the data with replacement $R$ times and estimating the corresponding
RESI for each resampled data. The lower and upper $\alpha/2\times 100\%$
percentiles of the $R$ estimated RESIs are the bootstrapped lower and upper
bounds of the $(1-\alpha)\times 100\%$ bootstrap CI for the estimated RESI.

In addition to the non-parametric bootstrap, we also consider several variants
of the wild bootstrap, which was originally proposed by Wu, (1986). The wild
bootstrap is a kind of residual bootstrap which leaves the covariates at the
sample value but resamples the outcome values based on the residual values.
That is, in each replicate, the resampled outcome value of the $i$-th
observation is

$y_{i}^{*}=\hat{y}_{i}+\hat{\epsilon}_{i}\times w_{i},$

where $\hat{y}_{i}$ and $\hat{\epsilon}_{i}$ are the expected value and
estimated residual of the $i$-th observation from the analysis model,
respectively. $w_{i}$ is a random number, called multiplier, drawn from a
distribution $W$ such that $\mathbb{E}(W)=0$ and $Var(W)=1$. There are several
different distributions from which $w_{i}$ can be drawn, such as the
Rademacher and standard Normal distributions.

The original wild bootstrap doesn’t resample the data. We propose several
variants of the original wild bootstrap: (1) resampling covariates along with
residuals with replacements; (2) fixing covariates and only resampling
residuals with replacements; (3) independently resampling covariates and
residuals with replacements. We also considered 3 types of multipliers: (1) no
multipliers, i.e., $w_{i}=1$, $\forall i$; (2) multipliers drawn from the
Rademacher distribution; (3) multipliers drawn from standard Normal
distribution. Table 1 summarizes the original and variant wild bootstraps
considered in this paper. The original wild bootstrap without using
multipliers is excluded because no resampling happens within this combination.
Therefore, there are 11 combinations of the wild bootstrapping evaluated in
the simulations. Detailed results for all the bootstrap methods are provided
in the Supplementary Material.

Wild bootstrap type | Multiplier | Bootstrap model |  Assumption(s)
---|---|---|---
Original | None | $Y^{(b)}=X\hat{\beta}+\hat{\epsilon}\times 1=Y$ |  Excluded b/c no sampling happens here
| Rademacher |  $Y^{(b)}=X\hat{\beta}+\hat{\epsilon}\times W$, $w_{i}\sim Rad$ |  Symmetric errors
| $N(0,1)$ |  $Y^{(b)}=X\hat{\beta}+\hat{\epsilon}\times W$, $w_{i}\sim N(0,1)$ |  Symmetric errors
(1) Resampling covariates along with residuals | None | $Y^{(b)}=R^{(b)}X\hat{\beta}+R^{(b)}\hat{\epsilon}\times 1$ |  Non-parametric bootstrap. not assuming homo-/hetero-skedasticity or symmetric errors
| Rademacher |  $Y^{(b)}=R^{(b)}X\hat{\beta}+R^{(b)}\hat{\epsilon}\times W$, $w_{i}\sim Rad$ |  Symmetric errors
| $N(0,1)$ |  $Y^{(b)}=R^{(b)}X\hat{\beta}+R^{(b)}\hat{\epsilon}\times W$, $w_{i}\sim N(0,1)$ |  Symmetric errors
(2) Fixing covariates and only resampling residuals | None | $Y^{(b)}=X\hat{\beta}+R^{(b)}\hat{\epsilon}\times 1$ |  Homoskedasticity
| Rademacher |  $Y^{(b)}=X\hat{\beta}+R^{(b)}\hat{\epsilon}\times W$, $w_{i}\sim Rad$ |  Homoskedasticity and symmetric errors
| $N(0,1)$ |  $Y^{(b)}=X\hat{\beta}+R^{(b)}\hat{\epsilon}\times W$, $w_{i}\sim N(0,1)$ |  Homoskedasticity and symmetric errors
(3) independently resampling covariates and residuals | None | $Y^{(b)}=R^{(b)}_{1}X\hat{\beta}+R^{(b)}_{2}\hat{\epsilon}\times 1$ |  Homoskedasticity
| Rademacher |  $Y^{(b)}=R^{(b)}_{1}X\hat{\beta}+R^{(b)}_{2}\hat{\epsilon}\times W$, $w_{i}\sim Rad$ |  Homoskedasticity and symmetric errors
| $N(0,1)$ |  $Y^{(b)}=R^{(b)}_{1}X\hat{\beta}+R^{(b)}_{2}\hat{\epsilon}\times W$, $w_{i}\sim N(0,1)$ |  Homoskedasticity and symmetric errors
Table 1: Eleven bootstrap procedures considered in this paper. The original
wild bootstrap with constant multiplier is excluded since it performs no
resampling. $R^{(b)}$ denotes the bootstrap matrix for the $b$-th replicate;
$Y^{(b)}$ is the bootstrapped outcome values; $X$ is the covariate(s);
$\hat{\beta}$ is the estimated parameters; $\hat{\epsilon}$ is the estimated
residuals; $W$ is the randomly drawn values for the multiplier.

###  3.2 Simulation setup

In the previous sections, we proposed three estimators and several ways of
constructing CIs for RESI. In this section, we use 1,000 simulations to
evaluate the performance of the proposed CI construction procedures with
respect to each estimator under 128 different scenarios. We evaluate the
influence of heteroskedasticity, data skewness and fixed/random covariate(s)
on the performance of estimators and CIs. All CIs are constructed at the
significance level of 0.05.

We simulate a simple linear regression model
$Y=\beta_{0}+\beta_{1}x+\epsilon$, where $x\in\\{0,1\\}$. We vary the sample
size $n\in\\{50,100,250,500\\}$ across 4 different true effect sizes,
$S\in\\{0,0.33,0.66,1\\}$. In the scenario of symmetric errors, the errors are
independently sampled from $N(0,\sigma_{0}^{2}=1^{2})$ under homoskedasticity;
Under heteroskedasticity, the errors are independently sample from
$N(0,\sigma_{0}^{2}=0.5^{2})$ if $x_{i}=0$ and from
$N(0,\sigma_{1}^{2}=1.5^{2})$ if $x_{i}=1$. In the scenario of skewed errors,
$\epsilon+\sqrt{0.1}$ $\sim$ Gamma($\alpha=0.1$, $\beta=\sqrt{0.1}$) under
homoskedasticity and $\epsilon+\sqrt{0.1}\cdot(x+0.5)$ $\sim$
Gamma($\alpha=0.1$, $\beta=\sqrt{0.1}/(x+0.5)$) under heteroskedasticity.
These errors are very heavily right-skewed, with Pearson’s moment coefficient
of skewness equal to $2/\sqrt{0.1}\approx 6.32$. To illustrate the difference
in the performance of CI between a randomized controlled and observational
design where the covariate(s) $X$ is treated as fixed or random, the values of
the covariate were generated in two ways: (1) when the covariate is fixed,
$\pi=0.3$ and ceiling($n\pi$) of the $n$ individuals have their covariate with
value 1 and the remaining have value 0; (2) when the covariate is random,
$X\sim\mathrm{Bern}(0.3)$ is sampled from a Bernoulli distribution with
parameter $\pi=0.3$. For each bootstrap CI, 1,000 bootstraps are used.

##  4 Simulation Results

We ran simulations to assess the bias of the estimators and the coverage of
the confidence intervals. We considered eight possible cases where there are
homo- or hetero-skedaticity errors, symmetric or skewed errors, and fixed or
random covariate in each simulation.

In small samples ($n=50$) the estimators are positively biased for $S=0$, but
negatively biased for all other values of $S$ (Figure 2). For all other sample
sizes $\hat{S}$ has small bias. As expected, under heteroskedasticity, the
parametric estimator is heavily biased (Vandekar et al.,, 2020). When the
error distribution is heavily skewed and $S>0$, the robust estimator is
biased, but the bias goes to zero in large samples (Figure 3). The randomness
of covariate doesn’t have an effect on the consistency of all these 3
estimators.

Figure 2: A comparison of bias for three estimators of the RESI under the
simulation settings described in Section 3.2. Under homoskedasticity the
estimators have small bias (top). Under heteroskedasticity, the parametric
estimator for effect size is heavily biased. $S$ denotes the RESI, Oracle
estimator assumes known variance, parametric assumes homoskedasticity, Robust
is the sandwich covariance estimator (Section 2.2). Figure 3: A comparison of
bias for three estimators of the RESI under the simulation settings described
in Section 3.2. Under homoskedasticity and heavily skewed residuals, the
estimators have bias (top), but the bias goes to 0 in large samples. Under
heteroskedasticity and heavily skewed residuals, the parametric estimator for
effect size is heavily biased, the bias of the robust estimator approaches 0
in large samples. $S$ denotes the RESI, Oracle estimator assumes known
variance, parametric assumes homoskedasticity, Robust is the sandwich
covariance estimator (Section 2.2).

The pattern of differences between the oracle, parametric, and robust
estimators is most obvious when the effect size is larger (e.g. $S=1$) with
fixed covariate (Figure 4). In this case, the Chi-squared CI should have
nominal coverage for the oracle estimator, the F CI should have nominal
coverage for the parametric estimator and the bootstrap CI should should have
nominal coverage for all (asymptotically; Figure 4). As expected, larger
effect sizes have worse coverage for the Chi-squared and F CIs when the
variance is estimated (parametric and robust statistics; Figure 4) because the
variance depends on the true effect size (see Section 2.2).

When there is random covariate, the Chi-squared and F CIs both fail to provide
nominal coverage for the oracle and parametric estimators, respectively. This
is because the extra variance introduced by random covariate(s) into the test
statistics (also shown in Figure 1).

Figure 4: A comparison of the coverage of different CIs for the three
estimators of RESI under homoskedasticity, with symmetric residuals and
random/fixed covariates. With fixed covariates, the Chi-squared CI has nominal
coverage for the oracle estimator and F CI has nominal coverage for the
parametric estimator, whereas they both fail when the covariates are random.
The bootstrapped CI has nominal coverage for all estimators, with random or
fixed covariates. $S$ denotes the RESI, Oracle estimator assumes known
variance, parametric assumes homoskedasticity, Robust is the sandwich
covariance estimator (Section 2.2)

Under heteroskedasticity, all CIs fail to provide nominal coverage for the
parametric estimator as expected because this estimator is biased as shown in
Figure 2. Both Chi-squared and F CIs fail to provide nominal coverage for the
robust estimator. The bootstrap CI has nominal coverage for both of the oracle
and robust estimators (Figure 5).

Figure 5: A comparison of the coverage of different CIs for the three
estimators of RESI under heteroskedasticity, with symmetric residuals and
random covariates. All three CIs fail to provide nomial coverage for the
parametric estimator. The Chi-squred and F CIs fail to provide nominal
coverage for the robust estimator. The bootstrapped CI has nominal coverage
for both the oracle and robust estimators. $S$ denotes the RESI, Oracle
estimator assumes known variance, parametric assumes homoskedasticity, Robust
is the sandwich covariance estimator (Section 2.2)

The skewness of the residuals has a big impact on the CIs’ performance. Under
homoskedasticity, the Chi-squared and F CIs provide nominal coverage for the
oracle estimator, but they both fail to provide nominal coverage for the
parametric and robust estimators (Figure 6). In large samples, the coverage of
the bootstrap CI for all three estimators approaches to the nominal level.
Under heteroskedasticity, all CIs fail to provide nominal coverage for all
estimators except F CI for the oracle estimator. The coverage of bootstrap CI
for the oracle and robust estimators approaches to the nominal level in large
samples.

Figure 6: A comparison of the coverage of different CIs for the three
estimators of RESI with skewed residuals and random covariates. In
homoskedasticity, the coverage of the bootstrapped CI approaches to the
nominal level in large sample for all three estimators. In heteroskedasticity,
the boostrapped CI approaches to the nominal level in large samples for the
oracle and robust estimators. $S$ denotes the RESI, Oracle estimator assumes
known variance, parametric assumes homoskedasticity, Robust is the sandwich
covariance estimator (Section 2.2)

##  5 Application of analysis of effect sizes

In this section, we use two datasets from studies of relational memory among
schizophrenia and early psychosis patients (Armstrong et al.,, 2012; Avery et
al.,, 2021) to illustrate how to conduct the analysis of effect sizes (ANOES)
with the function anoes in the RESI R package. Relational memory is the
ability to bind information into complex memories and is impaired in chronic
schizophrenia and in the early stages of psychosis (Armstrong et al.,, 2012,
2018; Avery et al.,, 2021). In both studies a relational memory paradigm was
used to compare the ability of psychosis patients and healthy individuals to
identify novel stimulus pairings (Armstrong et al.,, 2012, 2018; Avery et
al.,, 2021). The first study compared relational memory accuracy in 60
patients with schizophrenia or schizoaffective disorder to 38 healthy control
subjects Armstrong et al., (2012). The second study assessed participant’s
relational memory accuracy in 66 early psychosis patients and 64 healthy
control subjects (Avery et al.,, 2021).

In our analyses, in order to demonstrate the communicability of the RESI
across different models, we use both logistic regression models and multiple
linear models to quantify the differences in relational memory performance in
schizophrenia and early psychosis after controlling for age and gender,
respectively (See Appendix for R code). Then, ANOESs are conducted based on
the model results as the RESI and corresponding CIs are estimated to indicate
the effect size of each factor after controlling others. The results are
summarized in an ANOVA table format (see Table 2 and 3). We perform the ANOES
on the model fit object using the anoes function in R with the robust
covariance estimator, and the nonparametric bootstrap with 1,000 bootstrap
samples. We include estimates of effect size and CIs for each parameter
estimate as well as for the overall model fit. We also include the results
using the hypothesis testing framework output by the analysis of variance.

Table 2 shows that the estimated RESI of schizophrenia (after controlling for
age and gender) is 0.44 with 95% CI (0.24, 0.63) from the logistic model and
0.46 with 95% CI (0.24, 0.69) from the linear model. Table 3 shows that the
estimated RESI of early psychosis (after controlling for age and gender) is
0.49 with 95% CI (0.31, 0.72) from the logistic model and 0.50 with 95% CI
(0.32, 0.71) from linear model.

The estimated RESIs and their CIs from logistic and linear models are very
close to each other. If one model had much greater sensitivity to a given
effect, it would be represented by a larger effect size comparing across the
models. Using the RESI makes it easier to compare across logistic and linear
models.

The RESI also makes it easy to compare findings across these two studies that
have different sample sizes. While the test statistics are sample size
dependent, the RESI can be compared across the four models that come from two
different samples and two model types.

Logistic model: schizophrenia (SZ) vs. healthy controls
---
factor | estimate | robust s.e. | Chi-squared | d.f. |  $p$-value | RESI | 95% CI
group (SZ) | -0.77 | 0.18 | 18.92 | 1 |  $<$0.001 | 0.44 | (0.25, 0.63)
age | -0.03 | 0.01 | 14.18 | 1 |  $<$0.001 | 0.37 | (0.17, 0.59)
gender (female) | -0.15 | 0.14 | 1.19 | 1 | 0.27 | 0.05 | (0, 0.30)
Overall |  |  | 50.55 | 3 |  $<$0.001 | 0.71 | (0.53, 1.02)
Residual |  |  |  | 94 |  |  |
Linear model: schizophrenia (SZ) vs. healthy controls
factor | estimate | robust s.e. | Chi-squared | d.f. |  $p$-value | RESI | 95% CI
group (SZ) | -0.17 | 0.04 | 20.99 | 1 |  $<$0.001 | 0.46 | (0.24, 0.68)
age | -0.01 | 0.00 | 15.35 | 1 |  $<$0.001 | 0.39 | (0.18, 0.65)
gender (female) | -0.03 | 0.03 | 1.00 | 1 | 0.32 | 0 | (0, 0.29)
Overall |  |  | 72.76 | 3 |  $<$0.001 | 0.86 | (0.60, 1.38)
Residual |  |  |  | 94 |  |  |
Table 2: A comparison of ANOESs of the schizophrenia (SZ) from logistic and
linear models. The observed RESIs and their CIs for SZ after controlling for
age and gender are very close from the two different models, which
demonstrates the robustness of RESI across different models. The $p$-values
are calculated from the Chi-squared distribution. Logistic model: early
psychosis (EP) vs. healthy controls
---
factor | estimate | robust s.e. | Chi-squared | d.f. |  $p$-value | RESI | 95% CI
group (EP) | -1.26 | 0.22 | 31.43 | 1 |  $<$0.001 | 0.49 | (0.31, 0.72)
age | 0.04 | 0.04 | 0.79 | 1 | 0.37 | 0 | (0, 0.28)
gender (female) | 0.33 | 0.27 | 1.50 | 1 | 0.22 | 0.06 | (0, 0.26)
Overall |  |  | 45.17 | 3 |  $<$0.001 | 0.58 | (0.43, 0.82)
Residual |  |  |  | 126 |  |  |
Linear model: early psychosis (EP) vs. healthy controls
factor | estimate | robust s.e. | Chi-squared | d.f. |  $p$-value | RESI | 95% CI
group (EP) | -0.18 | 0.03 | 32.54 | 1 |  $<$0.001 | 0.50 | (0.32, 0.71)
age | 0.01 | 0.01 | 0.93 | 1 | 0.33 | 0 | (0, 0.29)
gender (female) | 0.04 | 0.03 | 1.73 | 1 | 0.19 | 0.08 | (0, 0.29)
Overall |  |  | 45.59 | 3 |  $<$0.001 | 0.58 | (0.44, 0.82)
Residual |  |  |  | 126 |  |  |
Table 3: A comparison of ANOESs of the early psychosis (EP) from logistic and
linear models. The observed RESIs and their CIs for EP after controlling for
age and gender are very close from the two different models, which
demonstrates the robustness of RESI across different models. The $p$-values
are calculated from the Chi-squared distribution.

##  6 Discussion

In this paper, we derived confidence intervals (CIs) for the robust effect
size index (RESI) and used them to describe an analysis of effect sizes
(ANOES) approach. We proposed 3 different estimators for the RESI and a
variety of ways to construct CIs. The oracle estimator assumes the covariance
matrix is known, so is not possible to compute in applications. Through
simulations, we showed that all 3 estimators are consistent under
homoskedasticity and the robust estimator is consistent under
heteroskedasticity. We used statistical theory and simulations to demonstrate
that the non-central Chi-squared CI has low coverage when the covariance in
the Wald test statistic needs to be estimated. In addition, the randomness of
covariate(s) also reduces the coverage of Chi-squared and F CIs, which is an
important implication for the observational studies where the covariates are
random instead of fixed/controlled by the experimenters. According to the
simulation results, using the robust estimator along with the (non-parametric)
bootstrap CI is generally most accurate and applicable to conduct consistent
estimation and valid inference using the RESI.

The RESI estimator and CI reduce the barriers to effect size analysis by
introducing a single framework that is widely applicable across different
models. They provide a basis for a framework of the ANOES, where effect sizes
with confidence intervals can easily be reported in summary tables alongside
$p$-values. This approach may help to address the limitations of null
hypothesis significance testing (Wasserstein and Lazar,, 2016; Wasserstein et
al.,, 2019) and may provide guidance on whether a study finding is under-
powered. We hope the ease of using RESI will broaden the use of effect sizes
and their CIs in study reporting.

We also provided functions to perform ANOES based on the RESI in our RESI R
package (https://github.com/statimagcoll/RESI). It outputs the estimated RESI
along with CIs in a summary table, with which researchers can conveniently
report their observed effect sizes along with trustworthy CIs. Coupled with
the generality of the RESI, researchers studying the same scientific questions
but using different types of data can easily communicate their observed effect
sizes and CIs instead of having to translate between different effect size
indices.

Our research focused on the performance of the CIs in a linear regression
model setting because these are widely used models. Although it adequately
illustrated the problems we wanted to discuss, theory and evaluation for other
models (e.g., mixed effects models and survival analysis) requires further
research. We formed the basis of applying the RESI in the ANOES on cross-
sectional data. It can be expected that there would be more technical
questions to be solved before making solid ANOES in longitudinal studies. In
the future, it would be interesting to work on the estimation and inference of
the RESI with longitudinal data.

## References

  * American Psychological Association, (2001) American Psychological Association (2001).  Publication manual of the American Psychological Association. American Psychological Association, Washington, DC, 5th ed. edition.
  * American Psychological Association, (2010) American Psychological Association (2010).  Publication manual of the American Psychological Association. American Psychological Association, Washington, DC, 6th ed. edition.
  * Armstrong et al., (2018) Armstrong, K., Avery, S., Blackford, J. U., Woodward, N., and Heckers, S. (2018).  Impaired associative inference in the early stage of psychosis.  Schizophrenia Research, 202:86–90.
  * Armstrong et al., (2012) Armstrong, K., Kose, S., Williams, L., Woolard, A., and Heckers, S. (2012).  Impaired associative inference in patients with schizophrenia.  Schizophrenia Bulletin, 38(3):622–629.
  * Avery et al., (2021) Avery, S. N., Armstrong, K., McHugo, M., Vandekar, S., Blackford, J. U., Woodward, N. D., and Heckers, S. (2021).  Relational memory in the early stage of psychosis: a 2-year follow-up study.  Schizophrenia bulletin, 47(1):75–86.  Publisher: Oxford University Press US.
  * Chinn, (2000) Chinn, S. (2000).  A simple method for converting an odds ratio to effect size for use in meta-analysis.  Statistics in Medicine, 19(22):3127–3131.
  * Cohen, (1988) Cohen, J. (1988).  Statistical power analysis for the behavioral sciences.  Erlbaum Associates, Hillsdale, NJ.
  * Efron, (1979) Efron, B. (1979).  Bootstrap Methods: Another Look at the Jackknife.  The Annals of Statistics, 7(1):1–26.
  * Fritz et al., (2012) Fritz, C. O., Morris, P. E., and Richler, J. J. (2012).  Effect size estimates: Current use, calculations, and interpretation.  Journal of Experimental Psychology: General, 141(1):2–18.
  * Hall, (1992) Hall, P. (1992).  The Bootstrap and Edgeworth Expansion.  Springer Series in Statistics. Springer-Verlag, New York.
  * Harlow et al., (2013) Harlow, L. L., Mulaik, S. A., and Steiger, J. H., editors (2013).  What If There Were No Significance Tests? Psychology Press, New York.
  * Hedges and Olkin, (1985) Hedges, L. V. and Olkin, I. (1985).  Statistical Methods for Meta-Analysis.  Elsevier, London, UK.
  * Kent and Hainsworth, (1995) Kent, J. T. and Hainsworth, T. J. (1995).  Confidence intervals for the noncentral chi-squared distribution.  Journal of Statistical Planning and Inference, 46(2):147–159.
  * Kubokawa et al., (1993) Kubokawa, T., Robert, C. P., and Saleh, A. K. M. E. (1993).  Estimation of Noncentrality Parameters.  The Canadian Journal of Statistics / La Revue Canadienne de Statistique, 21(1):45–57.
  * Long and Ervin, (2000) Long, J. S. and Ervin, L. H. (2000).  Using heteroscedasticity consistent standard errors in the linear regression model.  The American Statistician, 54(3):217–224.
  * Long and Freese, (2006) Long, J. S. and Freese, J. (2006).  Regression models for categorical dependent variables using Stata.  Stata press, 2 edition.
  * MacKinnon and White, (1985) MacKinnon, J. G. and White, H. (1985).  Some heteroskedasticity-consistent covariance matrix estimators with improved finite sample properties.  Journal of econometrics, 29(3):305–325.
  * Morris and DeShon, (2002) Morris, S. B. and DeShon, R. P. (2002).  Combining effect size estimates in meta-analysis with repeated measures and independent-groups designs.  Psychological methods, 7(1):105–125.
  * Neff and Strawderman, (1976) Neff, N. and Strawderman, W. E. (1976).  Further remarks on estimating the parameter of a noncentral chi-square distribution.  Communications in Statistics - Theory and Methods, 5(1):65–76.
  * Rosenthal, (1994) Rosenthal, R. (1994).  Parametric measures of effect size.  The handbook of research synthesis, 621:231–244.
  * Van der Vaart, (2000) Van der Vaart, A. W. (2000).  Asymptotic statistics, volume 3.  Cambridge university press.
  * Vandekar et al., (2020) Vandekar, S., Tao, R., and Blume, J. (2020).  A Robust Effect Size Index.  Psychometrika, 85(1):232.  Publisher: Springer.
  * Wasserstein and Lazar, (2016) Wasserstein, R. L. and Lazar, N. A. (2016).  The ASA’s statement on p-values: context, process, and purpose.  The American Statistician, 70(2):129–133.
  * Wasserstein et al., (2019) Wasserstein, R. L., Schirm, A. L., and Lazar, N. A. (2019).  Moving to a World Beyond “p < 0.05”.  The American Statistician, 73(sup1):1–19.
  * White, (1980) White, H. (1980).  A heteroskedasticity-consistent covariance matrix estimator and a direct test for heteroskedasticity.  Econometrica: Journal of the Econometric Society, pages 817–838.
  * Wu, (1986) Wu, C. F. J. (1986).  Jackknife, Bootstrap and Other Resampling Methods in Regression Analysis.  The Annals of Statistics, 14(4):1261–1295.  Publisher: Institute of Mathematical Statistics.
  * Zhang and Schoeps, (1997) Zhang, Z. and Schoeps, N. (1997).  On robust estimation of effect size under semiparametric models.  Psychometrika, 62(2):201–214.

## Appendix

### ANOES R code

In this section, we illustrate how to perform an ANOES in the relational
memory datasets using the functions from the RESI R package, and the complete
code for our data analyses are also included. The two datasets are described
in section 5 and have already been loaded into R. Each of these two datasets
contains 4 variables “accuracy” (the outcome), “group”, “age” and “gender”. We
first fit a glm object by specifying the variable accuracy as the outcome, it
can be a logistic regression model if we specify family = "binomial" or a
linear model if family = "gaussian" (default). Note, for both models, the
assumptions are violated because, for the logistic regression the data are not
Bernoulli distributed (they are proportions) and for the Gaussian model, the
errors are not normal. Then the ANOES based on the RESI can be implemented
with anoes(model.full = glm.fit), where glm.fit is the glm object we just
created.

The complete codes for our analyses are:

⬇

# \------------------------------------------- #

# ANOES in the relational memory data #

# \------------------------------------------- #

# Install the package RESI from Github

devtools::install_github("https://github.com/statimagcoll/RESI")

# Load in package(s)

library(RESI)

## The two datasets have already been loaded, called "SZ" and "EP"

## Each of them contains 3 variables called "accuracy", "age" and "gender"

## Set seed number (the ‘anoes()‘ function uses bootstraps to build CIs by
default)

set.seed(1213)

# Schizophrenia (SZ) patients vs Healthy control (HC) subjects

## fit a logistic model using ‘glm()‘

SZ_log_mod <- glm(accuracy ~ group + age + gender, data = SZ, family =
"binomial")

## conduct ANOES based on the logistic model

anoes(model.full = SZ_log_mod)

## fit a mutiple linear model

SZ_lin_mod <- glm(accuracy ~ group + age + gender, data = SZ, family =
"gaussian")

## conduct ANOES based on the linear model

anoes(model.full = SZ_lin_mod)

# Early psychosis (EP) patients vs healthy control (HC) subjects

## fit a logistic model using ‘glm()‘

EP_log_mod <- glm(accuracy ~ group + age + gender, data = EP, family =
"binomial")

## conduct ANOES based on the logistic model

anoes(model.full = EP_log_mod)

## fit a mutiple linear model

EP_lin_mod <- glm(accuracy ~ group + age + gender, data = EP, family =
"gaussian")

## conduct ANOES based on the linear model

anoes(model.full = EP_lin_mod)

By default, anoes uses the robust (sandwich) estimator (robust.var = TRUE,
default) to estimate the RESI, if robust.var = FALSE, it uses (4) to calculate
the parametric test statistics (2.2) and corresponding parametric version of
RESI estimates. anoes uses 1,000 bootstraps (nboot = 1000) to construct the
corresponding CIs. alpha is the significance level at which the CIs are
estimated, by default, alpha = 0.05;

If we are interested in estimating the overall effect size of a subset of
factors (e.g., age and gender in our example), we can specify a reduced model
to compare with (e.g., model.reduced = glm(accuracy $\sim$ group, data = data,
family = "binomial")) instead. By default, model.reduced = NULL. The function
will then output a table containing the overall effect size of age and gender
after controlling for the group. Here is an example:

⬇

> mod1 = glm(accuracy ~ group + age + gender, data = SZ, family = "gaussian")

> mod0 = glm(accuracy ~ group, data = SZ, family = "gaussian")

> anoes(model.full = mod1, model.reduced = mod0)

Chi-squared df p-val RESI LL UL

Tested 16.99 2 0 0.399 0.207 0.72

Residual NA 94 NA NA NA NA

where the RESI on the “tested” row is the RESI estimate for age and gender
together after controlling for group.

### Mathematical details

To show (8): the distribution is a non-central F distribution.
Under the homoskedasticity and Normality assumption, we have $Y\sim
N(X\beta,\sigma^{2}\mathbf{I})$.

$\displaystyle T_{(p)}^{2}/m_{1}$
$\displaystyle=n\hat{\sigma}^{-2}(\hat{\beta}_{OLS}-\beta_{0})^{T}X^{T}X(\hat{\beta}_{OLS}-\beta_{0})/m_{1}$
$\displaystyle=\frac{\sigma^{-2}n(\hat{\beta}_{OLS}-\beta_{0})^{T}X^{T}X(\hat{\beta}_{OLS}-\beta_{0})/m_{1}}{\sigma^{-2}\hat{\sigma^{2}}}$
$\displaystyle=\frac{\sigma^{-2}n(\hat{\beta}_{OLS}-\beta_{0})^{T}X^{T}X(\hat{\beta}_{OLS}-\beta_{0})}{\sigma^{-2}Y^{T}(I-H)Y/(n-m)}$

  1. 1.

In the numerator, under mild conditions, the quadratic form
$\sigma^{-2}n(\hat{\beta}_{OLS}-\beta_{0})^{T}X^{T}X(\hat{\beta}_{OLS}-\beta_{0})\sim\chi^{2}(m_{1};\lambda_{1})$,
where $\lambda_{1}$ is the non-centrality parameter:

$\displaystyle\lambda_{1}$
$\displaystyle=\sigma^{-2}n(\mathbb{E}\hat{\beta}_{OLS}-\beta_{0})^{T}X^{T}X(\mathbb{E}\hat{\beta}_{OLS}-\beta_{0})$
$\displaystyle=\sigma^{-2}n(\beta-\beta_{0})^{T}X^{T}X(\beta-\beta_{0})$
$\displaystyle=nS_{\beta}^{2}$

This quadratic form in the numerator is essentially the regression sum of
squares (SSR) in linear regression models.

  2. 2.

In the denominator, the quadratic from
$\sigma^{-2}Y^{T}(I-H)Y\sim\chi^{2}(n-m;\lambda_{2})$, where $\lambda_{2}$ is
the non-centrality parameter. It can be shown that $\lambda_{2}$ is 0:

$\displaystyle\lambda_{2}$
$\displaystyle=\sigma^{-2}\left[(X\beta)^{T}(I-H)X\beta\right]$
$\displaystyle=\sigma^{-2}\left[\beta^{T}X^{T}X\beta-\beta^{T}X^{T}HX\beta\right]$
$\displaystyle=\sigma^{-2}\left[\beta^{T}X^{T}X\beta-\beta^{T}X^{T}X(X^{T}X)^{-1}X^{T}X\beta\right]$
$\displaystyle=\sigma^{-2}\left[\beta^{T}X^{T}X\beta-\beta^{T}X^{T}X\beta\right]$
$\displaystyle=0$

Therefore, in the denominator, $\sigma^{-2}Y^{T}(I-H)Y$ follows a central Chi-
squared distribution $\chi^{2}_{n-m}$. This quadratic form is essentially the
error sum of squares (SSE) in linear regression models.

  3. 3.

Under homoskedasticity and Normality, the SSR and SSE are independent.
Therefore, we can prove that the distribution of $T_{(p)}^{2}/m_{1}$ is the
non-central F distribution $F(m_{1},n-m;nS^{2}_{\beta})$