Overview and evaluation of various frequentist test statistics using constrained statistical inference in the context of linear regression

Within the framework of constrained statistical inference, we can test informative hypotheses, in which, for example, regression coefficients are constrained to have a certain direction or be in a specific order. A large amount of frequentist informative test statistics exist that each come with different versions, strengths and weaknesses. This paper gives an overview about these statistics, including the Wald, the LRT, the Score, the $\overline{F}$- and the D-statistic. Simulation studies are presented that clarify their performance in terms of type I and type II error rates under different conditions. Based on the results, it is recommended to use the Wald and $\overline{F}$-test rather than the LRT and Score test as the former need less computing time. Furthermore, it is favorable to use the degrees of freedom corrected rather than the naive mean squared error when calculating the test statistics as well as using the $\overline{F}$- rather than the ${\overline{\chi }}^2$-distribution when calculating the p-values.

Introduction

Imagine a researcher wants to examine a novel psychotherapy program. A randomized experiment is set up with three treatment groups. One is a control group (X = 0), one participates in an established, standard psychotherapy program (X = 1) and one participates in the novel psychotherapy program (X = 2). No covariates are considered. The researcher is interested in the group means of the dependent variable Y, which denotes the score on a mental health questionnaire. Studies like this are usually conducted to show the superiority of the novel treatment over the standard treatment, as well as the superiority of the standard treatment over the control group. Thus, the researcher assumes that μ₂ > μ₁ > μ₀. However, following classical null hypothesis testing procedures, we usually first test a hypothesis like H₀ : μ₂ = μ₁ = μ₀ against H_a: not H₀, that is, not all three means are equal. If we can reject H₀ in favor of H_a, a second step often follows, in which we execute pairwise comparisons to determine which means are equal and which means are not equal. This implies multiple testing, which brings along the risk of an inflated type I error rate. The framework of constrained statistical inference (Silvapulle and Sen, 2005; Hoijtink, 2012) allows us to test so-called informative hypotheses, meaning that we can test the null hypothesis H₀ : μ₂ = μ₁ = μ₀ against the ordered hypothesis H_a : μ₂ > μ₁ > μ₀ in a single step. Thus, in contrast to classical null hypothesis testing, researchers have the advantages that they can formulate their hypotheses of interest directly, instead of making a detour via another hypothesis, while additionally avoiding to increase the risk for inflated type I error rates.

Informative hypothesis testing can be conducted by means of the Bayesian (see, e.g., Hoijtink et al., 2008; Hoijtink, 2012) as well as the frequentist (see, e.g., Barlow et al., 1972; Robertson et al., 1988; Silvapulle and Sen, 2005) approach, where the latter is the focus of this paper. The Bayesian approach is implemented in the R (R Core Team, 2020) package bain (Gu et al., 2020). The frequentist approach is implemented in SAS/STAT^® by means of the PLM procedure (for instructions, see Chapter 87 of SAS Institute Inc., 2015) as well as in several R packages including restriktor (Vanbrabant, 2020) and ic.infer (Grömping, 2010). Recent work of Keck et al. (2021) also demonstrated how to integrate informative hypothesis testing into the EffectLiteR (Mayer and Dietzfelbinger, 2019) package.

Restriktor and ic.infer use a broad range of test statistics, which are presented in Silvapulle and Sen (2005). However, research in the field of constrained statistical inference often uses the famous $\overline{F}$-statistic (see, e.g., Kuiper and Hoijtink, 2010; Vanbrabant et al., 2015) and neglects the distance statistic (D-statistic). Furthermore, each test statistic comes in various versions, for example depending on which estimate is used for the mean squared error or the variance-covariance matrix, and oftentimes, it is not obvious which software program uses which test statistic. There are also different options regarding the distributions that can be used to compute the p-values (${\overline{\chi }}^2,\overline{F}$). At the same time, small sample properties of informative test statistics are mostly unknown. Finally, simulation studies that examine the performance of informative test statistics are lacking in the constrained statistical inference literature.

The aim of this paper is twofold. First, we want to give an overview of a broad range of different informative test statistics, including the Wald test, the likelihood-ratio test (LRT), the Score test, the $\overline{F}$- and the D-statistic as well as their different versions. Second, we want to clarify how those test statistics perform when sample and effect sizes, hypotheses and the distribution used for calculating the p-values vary. Note that we only consider the regression setting, where all variables are observed. The paper is structured as follows: We start by presenting the univariate linear regression model to explain all necessary terminology that is used in the following section, where we define the test statistics. These test statistics include “regular” as well as informative test statistics to illustrate the link between them. We also discuss different versions of these test statistics. Subsequently, we report about simulation studies that we conducted. We introduce the design of the studies, that included a broad range of sample sizes as well as effect sizes, and we outline type I and type II error rates. We conclude with a short discussion. Supplementary materials are provided and will be referenced throughout the paper.

Univariate linear regression model

The univariate linear regression model for an observation i can be defined as:

$\begin{array}{ll}{y}_i={\beta }_0{x}_{i0}+{\beta }_1{x}_{i1}+{\beta }_2{x}_{i2}+...+{\beta }_p{x}_{ip}+{\epsilon }_i=x_i^{\prime }\beta +{\epsilon }_i,& (1)\end{array}$

where y_i is the value of the response variable for observation i = 1, 2, …, n, x_i0 is 1 and x_i1, …, x_ip are the values of the p regressors for observation i, which are assumed to be fixed (in terms of repeated sampling). β₀, …, β_p are the regression coefficients and ε_i is a residual for observation i. In matrix notation, the model can be written as y = Xβ + ε, where X is called the design matrix.

This regression model relies on several assumptions. First, we assume that the expected value of ε_i is zero. That is, E(ε_i) = 0 for all i. In matrix notation, this is expressed as E(ε) = 0, which implies that E(y) = Xβ, meaning that there is a linear relationship between E(y) and the columns of X. Second, we assume that x_i is non-stochastic and X is of full column rank. Third, we assume that the error term has a constant variance: $Var({\epsilon }_i)={\sigma }_{\epsilon }^2$ for all i. This implies that $Var(y_i)={\sigma }_{\epsilon }^2$ for all i. Fourth, we assume that the covariance of any two error terms is zero, that is Cov(ε_i, ε_j) = 0 for all (i, j), where i ≠ j.

The model can be estimated by means of different approaches such as ordinary least squares (OLS) or maximum likelihood (ML). It can be shown that under the presented assumptions, the OLS estimates of β are BLUE (best linear unbiased estimators, see, e.g., Seber and Lee, 2012). Using an example including four predictors, the following model is fitted:

$\begin{array}{ll}{y}_i={\beta }_0{x}_{i0}+{\beta }_1{x}_{i1}+{\beta }_2{x}_{i2}+{\beta }_3{x}_{i3}+{\beta }_4{x}_{i4}+{\epsilon }_i,& (2)\end{array}$

and ${\widehat{\beta }}_0,{\widehat{\beta }}_1,{\widehat{\beta }}_2,{\widehat{\beta }}_3$, and ${\widehat{\beta }}_4$ are obtained via OLS estimation. We may be interested in hypotheses concerning a single parameter like H₀ : β₁ = 0 vs. H_a : β₁ ≠ 0 or we might be interested in hypotheses about nested model comparisons like H₀ : β₁ = 0 ∧ β₂ = 0 vs. H_a : β₁ ≠ 0 ∨ β₂ ≠ 0 ∨ β₃ ≠ 0 ∨ β₄ ≠ 0. We can compute various important quantities that are used in hypothesis testing and that are characterized by a hat on top of it. Note that the hat indicates that estimation of the model parameters takes place in an unrestricted way, which will change once we test informative hypotheses. First, an unbiased estimator for the mean squared error is:

$\begin{array}{ll}{\widehat{\sigma }}_{\epsilon }^2={Ŝ}_{corrected}^2=\frac{\widehat{RSS}}{n-k},& (3)\end{array}$

where k is the column rank of X and $\widehat{RSS}$ is the estimated residual sum of squares $\sum _{i=1}^n{ê}_i^2$, where ê_i = y_i − ŷ_i and ŷ_i are the model predicted values of the response variable. Note that by considering k, we yield a small-sample correction for the mean squared error, as opposed to simply using:

$\begin{array}{ll}{Ŝ}_{naive}^2=\frac{\widehat{RSS}}{n},& (4)\end{array}$

which corresponds to the maximum likelihood estimator of ${\sigma }_{\epsilon }^2$.

The variance-covariance matrix of the estimated regression coefficients $\widehat{\beta }$ is usually computed as:

$\begin{array}{ll}VCOV(\widehat{\beta })=\frac{1}{n}{\widehat{I}}_1^{-1},& (5)\end{array}$

where $\widehat{I}$₁ is the unit information matrix:

$\begin{array}{ll}{\widehat{I}}_1=\frac{1}{n{Ŝ}_{corrected}^2}{X}^{\prime }X.& (6)\end{array}$

Note that if certain model assumptions are violated, for example if the error term does not have a constant variance, robust versions of the standard errors (Huber, 1967; White, 1980) and the variance-covariance matrix (Zeileis, 2006) can be used.

We can also test hypotheses about linear or non-linear combinations of regression parameters, like H₀ : β₁ + β₂ = 0 ∧ β₃ + β₄ = 0 vs. H_a : β₁ + β₂ ≠ 0 ∨ β₃ + β₄ ≠ 0. Note that in this paper, we will focus only on hypotheses containing linear combinations of regression coefficients. These combinations are specified by means of the R-matrix and each part of the hypothesis can be expressed as a row in R:

$\begin{array}{ll}{r}_1^{\prime }=\left(\begin{array}{ccccc}0& 1& 1& 0& 0\end{array}\right),& (7)\end{array}$

$\begin{array}{ll}{r}_2^{\prime }=\left(\begin{array}{ccccc}0& 0& 0& 1& 1\end{array}\right),& (8)\end{array}$

leading to the full constraint matrix:

$\begin{array}{ll}R=\left(\begin{array}{ccccc}0& 1& 1& 0& 0\\ 0& 0& 0& 1& 1\end{array}\right).& (9)\end{array}$

Then the hypothesis of interest can be expressed as H₀ : Rβ = 0 vs. H_a : Rβ ≠ 0. Note that all kinds of hypotheses, including the single parameter case and comparisons of nested models, as discussed before, can be expressed by means of the R-matrix.

In case our hypothesis of interest contains inequality constraints, like H_a : β₁ + β₂ > 0∨ β₃ + β₄ > 0, R still looks the same, but we need to fit a model where we enforce the inequality constraints on the regression coefficients. This can be done by means of quadratic programming, for example using the subroutine solve.QP() of the R package quadprog (Turlach and Weingessel, 2019). It implements the dual method of Goldfarb and Idnani (1982, 1983). If we apply this method in the linear regression context, it has the following form (see “Data Sheet 1” in the Supplementary materials for further explanations):

$\begin{array}{ll}min(-y^{\prime }X\beta +\frac{1}{2}{\beta }^{\prime }{X}^{\prime }X\beta )\text{ with the constraints }R\beta \ge {\beta }_0.& (10)\end{array}$

Note that all quantities based on an inequality constrained model are denoted by a tilde on top of them. Assume that the unconstrained estimates ${\widehat{\beta }}^{\prime }$ are (0.100 −0.1300.100 −0.2400.250), but the inequality constrained estimates ${\tilde{\beta }}^{\prime }$ may be (0.110 −0.1100.120 −0.2300.240), where the estimates of β₀, β₃ and β₄ may also change slightly, even though they already satisfied the constraints in the unrestricted estimation. The restricted estimation will also lead to different residuals than the unrestricted estimation.

If our hypothesis of interest contains equality constraints, for example H_a : β₁ + β₂ = 0 ∨ β₃ + β₄ = 0, the equality constrained estimates $\overline{\beta }$ can also be found via quadratic programming. Note that here, H_a from informative hypothesis testing equals H₀ from classical null hypothesis testing. Similarly, all estimated quantities with a bar on top are both the quantities from the equality constrained fit in informative hypothesis testing and the quantities obtained based on H₀ in classical null hypothesis testing, which are in fact equality constrained estimates as well. The corresponding mean squared error terms for the inequality and equality constrained case are defined as follows:

$\begin{array}{ll}{\tilde{S}}_{corrected}^2=\frac{\widetilde{RSS}}{n-k},& (11)\end{array}$

$\begin{array}{ll}{\tilde{S}}_{naive}^2=\frac{\widetilde{RSS}}{n},& (12)\end{array}$

$\begin{array}{ll}{\overline{S}}_{corrected}^2=\frac{\overline{RSS}}{n-(k-h)},& (13)\end{array}$

$\begin{array}{ll}{\overline{S}}_{naive}^2=\frac{\overline{RSS}}{n},& (14)\end{array}$

where $\widetilde{RSS}$ is the residual sum of squares of the inequality constrained fit $\sum _{i=1}^n{ẽ}_i^2$, where ẽ_i = y_i − ỹ_i and ỹ_i are the model predicted values of the response variable. Furthermore, $\overline{RSS}$ is the residual sum of squares under the equality constrained fit $\sum _{i=1}^n{ē}_i^2$, where ē_i = y_i − ȳ_i and ȳ_i are the model predicted values of the response variable. Finally, h is the row rank of R.

Similarly, we can define the unit information matrices of the inequality and equality constrained fits:

$\begin{array}{ll}{\tilde{I}}_1=\frac{1}{n{\tilde{S}}_{corrected}^2}{X}^{\prime }X,& (15)\end{array}$

$\begin{array}{ll}{\overline{I}}_1=\frac{1}{n{\overline{S}}_{corrected}^2}{X}^{\prime }X.& (16)\end{array}$

Note that X from the inequality constrained fit equals X from the unconstrained fit. The estimates $\widehat{\beta },\tilde{\beta }$ and $\overline{\beta }$ as well as the corresponding mean squared error terms and unit information matrices are used in the test statistics that are presented in the subsequent section.

Hypothesis testing

In order to give a broad overview about different test statistics, we present regular test statistics used in classical null hypothesis testing, as well as informative test statistics used in informative hypothesis testing. Note that an overview table containing all test statistics is provided at the end of each section. All test statistics can be applied in the setting of linear regression. “Data Sheet 2” in the Supplementary materials shows how these test statistics are implemented in R code.

Classical null hypothesis testing

The test statistics from classical null hypothesis testing that we will explain include the Wald test, the LRT, the Score test, the F-test as well as the t-test. The large sample test statistics, that is the Wald test, the LRT and the Score test, can be defined as follows Buse (1982):

$\begin{array}{ll}Wald=n(R\widehat{\beta })\prime (R{\widehat{I}}_1^{-1}{R}^{\prime }{)}^{-1}(R\widehat{\beta }),& (17)\end{array}$

$\begin{array}{ll}LRT=-2·[\ell (\overline{\beta })-\ell (\widehat{\beta })],& (18)\end{array}$

$\begin{array}{ll}Score=\frac{1}{n}S{(\overline{\beta })}^{\prime }{\overline{I}}_1^{-1}S(\overline{\beta }),& (19)\end{array}$

where $\ell (\overline{\beta })$ is the log-likelihood evaluated at $\overline{\beta }$, $\ell (\widehat{\beta })$ is the log-likelihood evaluated at $\widehat{\beta }$ and $S(\overline{\beta })=\frac{\partial }{\partial \overline{\beta }}\ell (\overline{\beta })$ is the score function evaluated at $\overline{\beta }$. All three test statistics follow asymptotically a χ²-distribution under the null hypothesis with df = h, if the model is correct.

Note that all three test statistics implicitly depend on S² in the information matrices (see Equation 6) and in the log-likelihoods. In the regression setting, since we always know what the residual degrees of freedom are, we can use ${Ŝ}_{corrected}^2$ instead of ${Ŝ}_{naive}^2$ to obtain the corrected instead of naive test statistic versions. That way, we can use the F-distribution with df₁ = h, df₂ = n − p to obtain the p-values, which is more precise in small samples compared to the χ²-distribution.

Note that the LRT, the Wald and the Score test are asymptotically equivalent. However, it has been shown that the values of the Wald test are always slightly larger than the values of the LRT, which in turn are always slightly larger than the values of the Score test (Buse, 1982, p. 157). Thus, using the same critical χ² value, the tests may have different power properties, which can be one aspect guiding the choice between them. Another aspect may be the time it takes to compute the three tests. For the Wald test, we need to fit the unconstrained model, whereas for the Score test, we need to fit the equality constrained model and for the LRT, we need to fit both the unconstrained and equality constrained model. In many cases, fitting the unconstrained model takes the least amount of time, which is why the Wald test is chosen often. However, in some cases, for example if the equality constrained model has a lot less parameters than the unconstrained model, it may be faster to fit the equality constrained model compared to the unconstrained model.

The F-test can be calculated as Seber and Lee (2012, p. 100):

$\begin{array}{ll}{F}_{corrected}=\frac{\frac{1}{h}[\overline{RSS}-\widehat{RSS}]}{{Ŝ}_{corrected}^2}.& (20)\end{array}$

Another test statistic version results from using ${Ŝ}_{naive}^2$ instead of ${Ŝ}_{corrected}^2$, which we denote as F_naive. Seber and Lee (2012, p. 100) show that F_corrected can be re-written to contain the unit information matrix:

$\begin{array}{ll}{F}_{corrected}^{info}=\frac{n}{h}(R\widehat{\beta })\prime (R{\widehat{I}}_1^{-1}{R}^{\prime }{)}^{-1}(R\widehat{\beta }),& (21)\end{array}$

where the superscript “info” refers to the information matrix. When ${Ŝ}_{naive}^2$ instead of ${Ŝ}_{corrected}^2$ is used in constructing the unit information matrix, we call this test statistic $F_{naive}^{info}$. If the model is specified correctly, F_corrected follows an F-distribution with df₁ = h, df₂ = n − k under the null hypothesis.

The one-sample t-test is defined as Allen(1997, p. 67):

$\begin{array}{ll}t=\frac{\widehat{\beta }-\overline{\beta }}{S{E}_{\widehat{\beta }}},& (22)\end{array}$

where $\overline{\beta }$ is the value of β under the null hypothesis and $S{E}_{\widehat{\beta }}$ is the standard error of $\widehat{\beta }$. Under the null hypothesis, t is t-distributed with df = n − k, if the model is correct. Note that if h = 1 the t- and F-statistic are related in a certain way, which is t² = F.

It is widely known that the one-sample t-test can be used for testing both two-sided hypotheses like H₀ : β = 0 against H_a : β ≠ 0 as well as one-sided hypotheses like H₀ : β = 0 against H_a : β > 0 or H_a : β < 0. The test statistic stays the same in both cases, but the p-value is computed differently. That is, when testing a two-sided hypothesis, half of the significance level is allocated to each side of the t-distribution, whereas when testing a one-sided hypothesis, all of it is allocated to one side of the t-distribution. That means that the cut-off levels, denoting from which point on the t-statistic can be considered to be significant, change. The two-sided p-value, which is the default output of most statistical software, simply adds up the probabilities of the negative and positive version of the observed t-value (t_obs), independently of whether it was in fact positive or negative:

$\begin{array}{l}{p}_{two-sided}=2·P(t>|t_{obs}|)\\ =P(t>t_{obs})+P(t<-t_{obs}).& (23)\end{array}$

Since the t-distribution is symmetric, P(t > t_obs) is the same as P(t < −t_obs). When we are interested in the one-sided p-value and H_a : β > 0, the p-value is obtained as:

$\begin{array}{ll}{p}_{one-sided}=P(t>t_{obs}),& (24)\end{array}$

whereas if H_a : β < 0, the p-value is obtained as:

$\begin{array}{ll}{p}_{one-sided}=P(t<t_{obs}).& (25)\end{array}$

Note that in case the obtained t-value is a positive number and we are interested in H_a : β > 0 or in case t is a negative number and we are interested in H_a : β < 0, the one-sided p-value can be obtained by dividing the two-sided p-value by 2.

In summary, the t-statistic is a special case, since this statistic from the classical null hypothesis testing framework can be used for testing an informative hypothesis, as long as the hypothesis only contains one parameter. If we are interested in more than one parameter, we can no longer use the t-statistic, but have to use an informative test statistic. Table 1 shows an overview about all presented regular test statistics.

TABLE 1

Table 1. Overview of all presented regular test statistics.

Informative hypothesis testing

Informative test statistics are often a modified version of the regular test statistics. In case the model is correct, the large sample informative test statistics, including the LRT, the Wald test, the Score test and the D-statistic, asymptotically follow a ${\overline{\chi }}^2$-distribution under the null hypothesis, which is a mixture of χ²-distributions. The small sample informative test statistic, that is the $\overline{F}$-statistic, follows an $\overline{F}$-distribution under the null hypothesis, if the model is correctly specified. The $\overline{F}$-distribution is a mixture of F-distributions. Note that similar to classical null hypothesis testing, we can use the corrected instead of naive mean squared error to obtain the large sample test statistics. In that way, we can calculate the p-values by means of the $\overline{F}$-distribution instead of the ${\overline{\chi }}^2$-distribution to obtain more precise results in small sample sizes.

The LRT_corrected test statistic can be calculated as follows Silvapulle and Sen (2005, p. 157):

$\begin{array}{ll}LR{T}_{corrected}=-2·[\ell (\overline{\beta })-\ell (\tilde{\beta })],& (26)\end{array}$

where $\ell (\overline{\beta })$ is the log-likelihood evaluated at $\overline{\beta }$ and $\ell (\tilde{\beta })$ is the log-likelihood evaluated at $\tilde{\beta }$. $\ell (\overline{\beta })$ has been calculated using ${\overline{S}}_{corrected}^2$ and $\ell (\tilde{\beta })$ has been calculated using ${\tilde{S}}_{corrected}^2$. If ${\overline{S}}_{naive}^2$ and ${\tilde{S}}_{naive}^2$ were used instead, we would obtain LRT_naive.

The Wald statistic can be found in Silvapulle and Sen (2005, p. 154):

$\begin{array}{ll}Wal{d}_{corrected}^{info}=\frac{n}{{Ŝ}_{corrected}^2}{(R\tilde{\beta })}^{\prime }{(R{W}^{-1}{R}^{\prime })}^{-1}(R\tilde{\beta }),& (27)\end{array}$

where $W=\frac{1}{n}{X}^{\prime }X$. The Wald version where we use ${Ŝ}_{naive}^2$ instead of ${Ŝ}_{corrected}^2$ is called $Wal{d}_{naive}^{info}$. Both versions implicitly contain $\widehat{I}$₁ (see Equation 6), which can also be replaced by $\tilde{I}$₁. Note that $Wal{d}_{naive}^{info}$ will give different results, especially in small sample sizes, due to the missing correction. Assuming $VCOV(\widehat{\beta })$ is defined as in Equation 5, we can re-write the Wald statistic as:

$\begin{array}{lll}Wal{d}^{VCOV}& ={[R\tilde{\beta }]}^{\prime }{[RVCOV(\widehat{\beta })R^{\prime }]}^{-1}[R\tilde{\beta }],& (28)\end{array}$

which is identical to $Wal{d}_{corrected}^{info}$. Note that we can also replace $VCOV(\widehat{\beta })$ by a more robust sandwich-estimator, which is not commonly done in the applied literature.

The D-statistic is calculated as follows (Silvapulle and Sen, 2005, p. 159):

$\begin{array}{ll}{D}_{corrected}=\frac{2·n}{{Ŝ}_{corrected}^2}[d(\overline{\beta })-d(\tilde{\beta })],& (29)\end{array}$

where $d(\overline{\beta })$ and $d(\tilde{\beta })$ are the values of the following two functions at their solutions (see “Data Sheet 3” in the Supplementary materials for further information):

$\begin{array}{ll}f(\beta )={(\widehat{\beta }-\beta )}^{\prime }W(\widehat{\beta }-\beta )\text{ under the constraint }R\beta =0,& (30)\end{array}$

$\begin{array}{ll}f(\beta )={(\widehat{\beta }-\beta )}^{\prime }W(\widehat{\beta }-\beta )\text{ under the constraint }R\beta \ge 0.& (31)\end{array}$

When minimizing these functions, we treat $\widehat{\beta }$ and W as known constants. Note that in the regression case, D_corrected is identical to $Wal{d}_{corrected}^{info}$ and Wald^VCOV, as long as ${Ŝ}_{corrected}^2$ is used. In contrast, if we switch to using ${Ŝ}_{naive}^2$, we obtain D_naive, in which case $D_{naive}=Wal{d}_{naive}^{info}$.

The $\overline{F}$-statistic can be found in (Silvapulle and Sen, 2005, p. 29):

$\begin{array}{ll}{\overline{F}}_{corrected}=\frac{\overline{RSS}-\widetilde{RSS}}{{Ŝ}_{corrected}^2}.& (32)\end{array}$

According to Silvapulle and Sen (2005, p. 29), including the constant $\frac{1}{h}$ from the regular F-statistic in the $\overline{F}$-statistic is not necessary, as it does not affect the results. Again, when using ${Ŝ}_{naive}^2$ instead of ${Ŝ}_{corrected}^2$, we obtain ${\overline{F}}_{naive}$. We can re-write the $\overline{F}$-statistic similarly to how we re-wrote the F-statistic. Assuming that we use ${Ŝ}_{corrected}^2$ to compute the unit information matrix, we obtain:

$\begin{array}{ll}{\overline{F}}_{corrected}^{info}=n{(R\tilde{\beta })}^{\prime }(R{\widehat{I}}_1^{-1}{R}^{\prime }{)}^{-1}(R\tilde{\beta }).& (33)\end{array}$

Again, $\widehat{I}$₁ can be replaced by $\tilde{I}$₁.

There are various versions of the Score statistic. $Scor{e}_{corrected}^U$ can be found in Silvapulle and Sen (2005, p. 159):

$\begin{array}{ll}Scor{e}_{corrected}^U=\frac{1}{n·{Ŝ}_{corrected}^2}{U}^{\prime }{(R{W}^{-1}{R}^{\prime })}^{-1}U,& (34)\end{array}$

where $U=R{W}^{-1}[S(\tilde{\beta })-S(\overline{\beta })]$. When using ${Ŝ}_{naive}^2$ as compared to ${Ŝ}_{corrected}^2$, we obtain $Scor{e}_{naive}^U$. Another version of the Score statistic, $Scor{e}_{corrected}^{null-info}$, is defined as follows Silvapulle and Silvapulle (1995, p. 342):

$\begin{array}{ll}Scor{e}_{corrected}^{null-info}=\frac{1}{n}{[S(\overline{\beta })-S(\tilde{\beta })]}^{\prime }{\overline{I}}_1^{-1}[S(\overline{\beta })-S(\tilde{\beta })],& (35)\end{array}$

where $\overline{I}$₁ has been calculated by means of ${\overline{S}}_{corrected}^2$ (see Equation 13). In contrast, if we use ${\overline{S}}_{naive}^2$, we obtain $Scor{e}_{naive}^{null-info}$.

Furthermore, $Scor{e}_{corrected}^{info}$ can be calculated as Silvapulle and Sen (2005, p. 166):

$\begin{array}{ll}Scor{e}_{corrected}^{info}=\frac{1}{n}{P}^{\prime }(R{\widehat{I}}_1^{-1}{R}^{\prime }{)}^{-1}P,& (36)\end{array}$

where $P=R{\widehat{I}}_1^{-1}[S(\tilde{\beta })-S(\overline{\beta })]$ and $\widehat{I}$₁ is calculated using ${Ŝ}_{corrected}^2$ and can be replaced by either $\tilde{I}$₁ or $\overline{I}$₁. If we use ${Ŝ}_{naive}^2$ to calculate $\widehat{I}$₁, we obtain $Scor{e}_{naive}^{info}$. Silvapulle and Sen (2005, p. 166) mention another way to express $Scor{e}_{corrected}^{info}$:

$\begin{array}{ll}Scor{e}_{corrected}^{info,Robertson}=\frac{1}{n}{[S(\tilde{\beta })-S(\overline{\beta })]}^{\prime }{\widehat{I}}_1^{-1}[S(\tilde{\beta })-S(\overline{\beta })],& (37)\end{array}$

where the superscript “Robertson” indicates that this is the version defined by Robertson et al. (1988), $\widehat{I}$₁ is calculated using ${Ŝ}_{corrected}^2$ and can be replaced by either $\tilde{I}$₁ or $\overline{I}$₁. Assuming that $VCOV(\widehat{\beta })$ is defined as in Equation 5, $Scor{e}_{corrected}^{info}$ can be re-written as:

$\begin{array}{ll}Scor{e}^{VCOV}=V^{\prime }{[RVCOV(\widehat{\beta })R^{\prime }]}^{-1}V,& (38)\end{array}$

where $V=RVCOV(\widehat{\beta })[S(\tilde{\beta })-S(\overline{\beta })]$, again allowing for a more robust sandwich-estimator of $VCOV(\widehat{\beta })$ to be inserted. Table 2 gives an overview about all the informative test statistics that were presented.

TABLE 2

Table 2. Overview of all presented informative test statistics.

P-values

There are two approaches for calculating the p-value of informative test statistics (Silvapulle and Sen, 2005). In this paper, we use the approach where we first calculate the weights of the respective mixture distribution (${\overline{\chi }}^2,\overline{F}$). Note that the sum of the weights from 0 to q is one, where q is the rank of X under the null hypothesis.

If the residuals of our data are normally distributed, we can use the multivariate normal probability function as well as the ic.weight() function of the R package ic.infer (Grömping, 2010) to compute the weights. These calculations are also implemented in the R package restriktor (Vanbrabant, 2020). Once we have computed the weights, the p-values of the observed ${\overline{\chi }}^2$-value (${\overline{\chi }}_{obs}^2$) and of the observed $\overline{F}$-value (${\overline{F}}_{obs}$) are obtained as follows Silvapulle and Sen (2005, pp. 86 and 99):

$\begin{array}{ll}Pr({\overline{\chi }}^2\ge {\overline{\chi }}_{obs}^2)={\displaystyle \sum _{i=0}^q}{w}_i(H_0,H_a)\text{Pr}[(h-q+i){\chi }_{h-q+i}^2\ge {\overline{\chi }}_{obs}^2],& (39)\end{array}$

$\begin{array}{ll}Pr(\overline{F}\ge {\overline{F}}_{obs})={\displaystyle \sum _{i=0}^q}{w}_i(H_0,H_a)\text{Pr}[(h-q+i)F_{h-q+i,n-p}\ge {\overline{F}}_{obs}].& (40)\end{array}$

It can be expected that the p-values are very similar, irrespective of whether they are calculated based on the ${\overline{\chi }}^2$- or $\overline{F}$-distribution, as long as sample sizes are large. However, for small sample sizes, the $\overline{F}$-distribution should yield more accurate results.

Simulation studies

We conducted several simulation studies to examine the impact of different conditions on the performance of the presented test statistics in terms of type I and type II error rates. We were interested in the effects of sample and effect sizes, the number of regression parameters considered in H_a as well as the distribution used for calculating the p-values. Our main motivation was to provide a reference framework for applied researchers who wish to test informative hypotheses, helping them to chose the optimal test statistic(s) in the present situation.

Design

We generated a design matrix X, including data for five regression coefficients ${\beta }^{\prime }=({\beta }_1{\beta }_2{\beta }_3{\beta }_4{\beta }_5)$ and considered effect sizes of f² = 0.02 (small), f² = 0.10 (medium) and f² = 0.35 (large) and sample sizes of 10, 25, 50, 100, 500, 1000, 2000, and 10000. For examining the type I error rate, we generated a random outcome Y, whereas for examining the type II error rate, we fixed all βs to 0.1 and generated y with a random error term that was specific for the effect size used. Since $f^2=\frac{R^2}{1-R^2}$, where R² is the determination coefficient, we can calculate the error terms of y by plugging in the f²-specific value of R² in

$\begin{array}{ll}{S}_y^2=[\beta Cor(X)\beta ]\times \frac{1-R^2}{R^2},& (41)\end{array}$

where Cor(X) is the correlation matrix of the design matrix X. The number of replications was 1000.

We considered two different kinds of R matrices, where the first one was defined as follows:

$\begin{array}{ll}{R}_1=\left(\begin{array}{cccccc}0& 1& 0& 0& 0& 0\end{array}\right).& (42)\end{array}$

This represents the hypothesis that only β₁ is greater than zero: H_a : β₁ > 0. The second R matrix was defined as:

$\begin{array}{ll}{R}_2=\left(\begin{array}{cccccc}0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right),& (43)\end{array}$

stating that at least one of the regression coefficients, except the intercept, are greater than zero: H_a : β₁ > 0 ∨ β₂ > 0 ∨ β₃ > 0 ∨ β₄ > 0 ∨ β₅ > 0.

To compute the test statistics, we used ${Ŝ}_{naive}^2$ and ${Ŝ}_{corrected}^2$ as well as ${\tilde{S}}_{naive}^2,{\tilde{S}}_{corrected}^2,{\overline{S}}_{naive}^2$ and ${\overline{S}}_{corrected}^2$ and to compute the p-values, we used the ${\overline{\chi }}^2$- as well as the $\overline{F}$-distribution. In addition to the manual calculation of the test statistics, we also included the test statistics as reported by the R package restriktor.

Type I results

Test statistics were first applied the way they are presented in the referenced literature. That is, $Wal{d}_{naive}^{info}$ makes use of ${Ŝ}_{naive}^2$, whereas all other test statistics make use of ${Ŝ}_{corrected}^2$ (or ${\tilde{S}}_{corrected}^2,{\overline{S}}_{corrected}^2$). For calculating the p-values, the ${\overline{\chi }}^2$-distribution is used for LRT_corrected, $Wal{d}_{naive}^{info}$, $Wal{d}_{corrected}^{info}$, Wald^VCOV, D_corrected, $Scor{e}_{corrected}^U$, $Scor{e}_{corrected}^{null-info}$, $Scor{e}_{corrected}^{info}$ and Score^VCOV. The $\overline{F}$-distribution is used for calculating the p-values for the $\overline{F}$-statistic, the F-distribution is used for calculating the p-values for the F-statistic and the t-distribution is used for calculating the p-values for the t-statistic. Note that restriktor always uses ${Ŝ}_{corrected}^2$ (or ${\tilde{S}}_{corrected}^2,{\overline{S}}_{corrected}^2$) for all available test statistics and always calculates the p-value based on the $\overline{F}$-distribution. Tables 3, 4 show the results.

TABLE 3

Table 3. Type I error rates when using R₁ and applying the test statistics as outlined in the referenced books.

TABLE 4

Table 4. Type I error rates when using R₂ and applying the test statistics as outlined in the referenced books.

We can observe that when using R₁ (see Table 3), that is when testing a hypothesis concerning only one regression parameter, type I error rates are identical between F and t_two−sided as well as between $\overline{F}$ and t_one−sided, showing the link between classical null hypothesis testing and informative hypothesis testing. When using R₂ (see Table 4), that is when testing a hypothesis concerning multiple regression parameters, problems with type I error rates seem to occur earlier as compared to when using R₁. More specifically, problematic type I error rates occur as early as with n = 500 or n = 100 when using R₂, but only start occurring with n = 50 or n = 25 when using R₁. Apart from that, $Scor{e}_{corrected}^U$ and $Wal{d}_{naive}^{info}$ show the highest type I error rates for both R matrices, whereas $\overline{F}$ and ${\overline{F}}_{restriktor}$ show the most appropriate type I error rates for both R matrices. This is because the $\overline{F}$-distribution is more precise in small sample sizes as compared to the ${\overline{\chi }}^2$-distribution.

When using the $\overline{F}$-distribution instead of the ${\overline{\chi }}^2$-distribution when calculating the p-value for all test statistics, type I error rates are closer to the nominal level when sample sizes get smaller. This can be seen in Tables 5, 6 where a selection of test statistics are shown.

TABLE 5

Table 5. Type I error rates when using R₁, ${Ŝ}_{corrected}^2$ (or ${\tilde{S}}_{corrected}^2,{\overline{S}}_{corrected}^2$) and the $\overline{F}$-distribution for calculating the p-value.

TABLE 6

Table 6. Type I error rates when using R₂, ${Ŝ}_{corrected}^2$ (or ${\tilde{S}}_{corrected}^2,{\overline{S}}_{corrected}^2$) and the $\overline{F}$-distribution for calculating the p-value.

Furthermore, it can be observed that when using R₁, type I error rates increase when using LRT_corrected and $Scor{e}_{corrected}^U$ and n = 10 in contrast to n = 25. The same can only be observed for $Scor{e}_{corrected}^U$ when using R₂, but not for LRT_corrected, where the type I error rate decreases quite substantially instead. More results can be found in “Data Sheet 4” in the Supplementary materials.

Type II results

Figures 1, 2 show the type II error rates when applying the test statistics as in the referenced books.

FIGURE 1

Figure 1. Type II error rates when using R₁ and applying the test statistics as outlined in the referenced books.

FIGURE 2

Figure 2. Type II error rates when using R₂ and applying the test statistics as outlined in the referenced books.

Once more, we can observe that when using R₁ (see Figure 1), that is when testing a hypothesis concerning only one regression parameter, type II error rates are identical between F and t_two−sided as well as between $\overline{F}$ and t_one−sided, showing the link between classical null hypothesis testing and informative hypothesis testing. When using R₂ (see Figure 2), that is when testing a hypothesis concerning multiple regression parameters, problems with type II error rates seem to occur later (in terms of sample size) as compared to when using R₁. This was the other way around regarding the type I error rate and it demonstrates the nature of the relationship between type I and type II error rates: If one goes down, the other one goes up and vice versa.

The same mechanism can be observed when using the $\overline{F}$-distribution instead of the ${\overline{\chi }}^2$-distribution when calculating the p-value for all test statistics (Figures 3, 4): Type II error rates are increased in small sample sizes, since type I error rates had improved, that is, decreased. Again, further results can be found in “Data Sheet 5” in the Supplementary materials.

FIGURE 3

Figure 3. Type II error rates when using R₁, ${Ŝ}_{corrected}^2$ (or ${\tilde{S}}_{corrected}^2,{\overline{S}}_{corrected}^2$) and the $\overline{F}$-distribution for calculating the p-value.

FIGURE 4

Figure 4. Type II error rates when using R₂, ${Ŝ}_{corrected}^2$ (or ${\tilde{S}}_{corrected}^2,{\overline{S}}_{corrected}^2$) and the $\overline{F}$-distribution for calculating the p-value.

Discussion

In this paper, we gave an overview of a large number of different informative test statistics, including their different versions. Furthermore, we clarified how those test statistics perform in terms of type I and type II error rates under different conditions by means of simulation studies in the context of linear regression. We considered varying sample and effect sizes as well as two different constraint matrices, where one specified a hypothesis about one parameter and the other one specified a hypothesis about multiple parameters. Moreover, we considered the naive and corrected mean squared errors of the unconstrained, inequality and equality constrained models as part of the test statistics as well as the ${\overline{\chi }}^2$- and $\overline{F}$-distribution to calculate the p-values.

Based on our findings, the following recommendations can be made. Considering the time it takes to compute the informative test statistics, both the Wald and the $\overline{F}$-test versions are favorable, since they only need fitting of the inequality constrained model to obtain $\tilde{\beta }$ and $\tilde{I}$₁. Even if we do not use $\tilde{I}$₁ but use $\widehat{I}$₁ instead, the increase in time is small in the context of linear regression. The Score test and the LRT versions are less favorable, since they require fitting both the inequality constrained as well as the equality constrained model to obtain $\tilde{\beta }$ and $\overline{\beta }$, as well as the respective unit information matrices or log-likelihoods.

The D-statistic versions only require fitting the unconstrained model to obtain $\widehat{\beta }$. However, we then additionally need to compute the two functions $d(\overline{\beta })$ and $d(\tilde{\beta })$, which is as time-consuming as fitting the inequality constrained model. Thus, there is no advantage of using the D-statistic versions over the Wald and the $\overline{F}$-test versions in the context of linear regression. However, if the regression model was non-linear, computing the two functions would be significantly less computationally expensive than fitting the inequality constrained model.

Moreover, we recommend using the corrected mean squared error versions in the test statistics as well as using the $\overline{F}$-distribution for calculating the p-values, if sample sizes are small. This seems to keep type I error rates closer to the nominal level compared to using the naive mean squared error versions and using the ${\overline{\chi }}^2$-distribution for calculating the p-value. An additional interesting finding was that the relationship between LRT, Wald and Score test values that has been found in the unconstrained context also holds in the constrained context. That is, Wald test values are always slightly larger than LRT values, which in turn are always slightly larger than Score test values.

The limitations of our simulation studies include the following aspects. We treated all variables as manifest, even though variables of interest in the social and behavioral sciences are often latent in nature. Furthermore, we solely generated normal data despite the fact that violations against the normality assumption occur regularly. Moreover, we used orthogonal predictors without interactions albeit this is rarely the case in the social and behavioral sciences. And lastly, we only included the regular versions of the standard errors and the variance-covariance matrix. Future research should thus repeat the simulation studies in the context of Structural Equation Modeling (SEM) to take into account latent variables. Furthermore, the impact of non-normal data as well as correlated predictors with interactions and using the robust versions of the standard errors and the variance-covariance matrix should be examined. It may be that under these conditions, type I and type II error rates deviate from the results presented in this paper. Moreover, the properties of informative test statistics, especially concerning the D-statistic, should also be investigated in the context of non-linear models.

Finally, research in the social and behavioral sciences is often not only interested in inference concerning regression coefficients, but also regarding effects of interest. These effects may be average or conditional treatment effects, which are defined as a linear or non-linear combination of regression coefficients. The EffectLiteR approach (Mayer et al., 2016) provides a framework and R package for the estimation of average and conditional effects of a discrete treatment variable on a continuous outcome variable, conditioning on categorical and continuous covariates. Keck et al. (2021) already demonstrated how to integrate informative hypothesis testing into the EffectLiteR framework in the context of linear regression. The present paper provides interested readers who want to apply informative hypothesis testing concerning regression coefficients or effects of interest with practical information regarding test statistics as well as type I and type II error rates.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions

CK and YR contributed to conception and design of the study. CK performed the statistical analysis and wrote the first draft of the manuscript. CK, YR, and AM wrote sections of the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.

Funding

This work has been supported by the Research Foundation Flanders (FWO, Grant No. G002819N to YR and AM).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpsyg.2022.899165/full#supplementary-material

References

Allen, M. P. (1997). Understanding Regression Analysis. Boston, MA: Springer.

Google Scholar

Barlow, R. E., Bartholomew, D. J., Bremner, J. M., and Brunk, H. D. (1972). Statistical Inference under Order Restrictions. New York, NY: Wiley.

PubMed Abstract

Buse, A. (1982). The likelihood ratio, Wald and Lagrange multiplier tests: An expository note. Am. Stat. 36, 153–157. doi: 10.1080/00031305.1982.10482817

CrossRef Full Text | Google Scholar

Goldfarb, D., and Idnani, A. (1982). Dual and Primal-Dual Methods for Solving Strictly Convex Quadratic Programs, Berlin: Springer.

Google Scholar

Goldfarb, D., and Idnani, A. (1983). A numerically stable dual method for solving strictly convex quadratic programs. Math. Program. 27, 1–33. doi: 10.1007/BF02591962

CrossRef Full Text | Google Scholar

Grömping, U. (2010). Inference with linear equality and inequality constraints using R: The package ic.infer. J. Stat. Softw. 33, 1–31. doi: 10.18637/jss.v033.i10

CrossRef Full Text | Google Scholar

Gu, X., Hoijtink, H., Mulder, J., Van Lissa, C. J., Van Zundert, C., Jones, J., et al. (2020). Bain: Bayes factors for informative hypotheses. R package version 0.2.4.

PubMed Abstract

Hoijtink, H. (2012). Informative Hypotheses: Theory and Practice for Behavioral and Social Scientists. Boca Raton, FL: Chapman & Hall/CRC.

Google Scholar

Hoijtink, H., Klugkist, I., and Boelen, P. A. (2008). Bayesian Evaluation of Informative Hypotheses. New York, NY: Springer.

Google Scholar

Huber, P. (1967). “The behavior of maximum likelihood estimates under nonstandard conditions,” in Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability, Vol. 1, Statistics (Berkeley, CA: University of California Press), 221–233.

Google Scholar

Keck, C., Mayer, A., and Rosseel, Y. (2021). Integrating informative hypotheses into the EffectLiteR framework. Methodology 17, 307–325. doi: 10.5964/meth.7379

CrossRef Full Text | Google Scholar

Kuiper, R. M., and Hoijtink, H. (2010). Comparisons of means using exploratory and confirmatory approaches. Psychol. Methods 15, 69–86. doi: 10.1037/a0018720

PubMed Abstract | CrossRef Full Text | Google Scholar

Mayer, A., and Dietzfelbinger, L. (2019). EffectLiteR: Average and conditional effects. R package version 0.4–4.

PubMed Abstract

Mayer, A., Dietzfelbinger, L., Rosseel, Y., and Steyer, R. (2016). The EffectLiteR approach for analyzing average and conditional effects. Multivariate Behav. Res. 5, 374–391. doi: 10.1080/00273171.2016.1151334

PubMed Abstract | CrossRef Full Text | Google Scholar

R Core Team (2020). R: A Language and Environment for Statistical Computing. Vienna: R foundation for statistical computing.

Google Scholar

Robertson, T., Wright, F. T., and Dykstra, R. L. (1988). Order Restricted Statistical Inference. New York, NY: Wiley.

SAS Institute Inc. (2015). Sas/stat^® 14.1 User'S Guide. Cary, NC: SAS Institute Inc.

Seber, G. A. F., and Lee, A. J. (2012). Linear Regression Analysis. Hoboken, NJ: Wiley.

Google Scholar

Silvapulle, M. J., and Sen, P. K. (2005). Constrained Statistical Inference: Order, Inequality, and Shape Restrictions. Hoboken, NJ: Wiley.

Google Scholar

Silvapulle, M. J., and Silvapulle, P. (1995). A Score test against one-sided alternatives. J. Am. Stat. Assoc. 90, 342–349. doi: 10.1080/01621459.1995.10476518

CrossRef Full Text | Google Scholar

Turlach, B. A., and Weingessel, A. (2019). Quadprog: Functions to solve quadratic programming problems. R package version 1.5–8.

Vanbrabant, L. (2020). Restriktor: Constrained statistical inference. R package version 0.2–800.

Vanbrabant, L., Van de Schoot, R., and Rosseel, Y. (2015). Constrained statistical inference: Sample-size tables for ANOVA and regression. Front. Psychol. 5, 1565. doi: 10.3389/fpsyg.2014.01565

PubMed Abstract | CrossRef Full Text | Google Scholar

White, H. (1980). A heteroskedasticity-consistent covariance matrix and a direct test for heteroskedasticity. Econometrica 48, 817–838. doi: 10.2307/1912934

CrossRef Full Text | Google Scholar

Zeileis, A. (2006). Object-oriented computation of sandwich estimators. J. Stat. Softw. 16, 1–16. doi: 10.18637/jss.v016.i09

CrossRef Full Text | Google Scholar