In this example we reconsider the ZelazoKolb1972 dataset discussed
in the example about ordered means (Ex.1). As a reminder, the purpose of the
study was to test the claim that the walking exercises are associated with a
reduction in the mean age at which children start to walk. The results confirmed
strong evidence in favor of the order-constrained hypothesis $$ H_1: \mu_1 \leq \mu_2 \leq \mu_3 \leq \mu_4 $$.
The second question we like to answer is whether the differences between the means
are relevant. To answer this question we use the popular effect-size measure
Cohen's $d$ and is given by:
$$ d = \frac{\mu_{max} - \mu_{min}}{\sigma}, $$
where $\mu_{max}$ is the largest of the $k$ means and $\mu_{min}$ is the smallest of $k$ means and $\sigma$ is the pooled standard deviation within the populations.
Hypothesis $H_1$ can be reformulated such that the effect-sizes are included. According to Cohen, values of 0.2, 0.5 and 0.8 indicate a small, medium and large effect, respectively. The within group standard deviation ($\sigma$) equals 1.516:
\begin{equation*} \begin{array}{l} H_1 =\ \end{array} \begin{array}{l} (\mu_2 - \mu_1) / 1.516 \geq 0.2 \\ (\mu_3 - \mu_2) / 1.516 \geq 0.2 \\ (\mu_4 - \mu_3) / 1.516 \geq 0.2. \end{array} \end{equation*}
This hypothesis states that $\mu_2$ is at least 0.2 standard deviations larger than $\mu_1$. The other two constraints have analogue interpretations.
In what follows, we describe all steps to test this order constrained hypothesis.
library(restriktor)
The ZelazoKolb1972 dataset is available in restriktor and can be
called directly, see next step.
An ANOVA model is just a special case of the linear model. Therefore, we can make
use of the built-in linear model lm() function in R. The intercept (-1)
is removed, such that the regression coefficients reflect the group means. To fit
the unrestricted linear model type
fit.ANOVAd <- lm(Age ~ -1 + Group, data = ZelazoKolb1972)
The constraint syntax can be constructed by using the factor level names preceded by the factor name. To get the correct names, type
names(coef(fit.ANOVAd))
[1] "GroupActive" "GroupControl" "GroupNo" "GroupPassive"
Based on these factor level names, we can construct the constraint syntax:
myConstraints <- ' (GroupPassive - GroupActive) / 1.516 > 0.2
(GroupControl - GroupPassive) / 1.516 > 0.2
(GroupNo - GroupControl) / 1.516 > 0.2 '
# note that the constraint syntax is enclosed within single quotes.
The restricted ANOVA model is estimated using the restriktor() function. The first argument to restriktor() is the fitted linear model and the second argument is the constraint syntax.
restr.ANOVAd <- restriktor(fit.ANOVAd, constraints = myConstraints)
summary(restr.ANOVAd)
Call:
conLM.lm(object = fit.ANOVAd, constraints = myConstraints)
Restriktor: restricted linear model:
Residuals:
Min 1Q Median 3Q Max
-2.7083 -0.8500 -0.3500 0.6375 3.6250
Coefficients:
Estimate Std. Error t value Pr(>|t|)
GroupActive 10.12500 0.61907 16.355 1.191e-12 ***
GroupControl 11.70833 0.61907 18.913 8.774e-14 ***
GroupNo 12.35000 0.67815 18.211 1.736e-13 ***
GroupPassive 11.37500 0.61907 18.375 1.478e-13 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.5164 on 19 degrees of freedom
Standard errors: standard
Multiple R-squared remains 0.986
Generalized Order-Restricted Information Criterion:
Loglik Penalty Goric
-40.0142 3.0976 86.2236
## by default the iht() function prints a summary of all available
## hypothesis tests. A more detailed overview of the separate hypothesis
## tests can be requested by specifying the test = "" argument, e.g.,
## iht(restr.ANOVAd, test = "A").
iht(fit.ANOVAd, constraints = myConstraints)
Restriktor: restricted hypothesis tests ( 19 residual degrees of freedom ):
Multiple R-squared remains 0.986
Constraint matrix:
GroupActive GroupControl GroupNo GroupPassive op rhs active
1: -0.6596 0 0 0.6596 >= 0.2 no
2: 0 0.6596 0 -0.6596 >= 0.2 no
3: 0 -0.6596 0.6596 0 >= 0.2 no
Overview of all available hypothesis tests:
Global test: H0: all parameters are restricted to be equal (==)
vs. HA: at least one inequality restriction is strictly true (>)
Test statistic: 2.3840, p-value: 0.1689
Type A test: H0: all restrictions are equalities (==)
vs. HA: at least one inequality restriction is strictly true (>)
Test statistic: 2.3840, p-value: 0.1689
Type B test: H0: all restrictions hold in the population
vs. HA: at least one restriction is violated
Test statistic: 0.0000, p-value: 1
Type C test: H0: at least one restriction is false or active (==)
vs. HA: all restrictions are strictly true (>)
Test statistic: 0.0344, p-value: 0.4865
Note: Type C test is based on a t-distribution (one-sided),
all other tests are based on a mixture of F-distributions.
In Ex.1 we found a significant result for the informative hypothesis. However, if we include a small effect-size in the hypothesis test, the mean differences become irrelevant. Although, the results from hypothesis test Type B show that all order-constraints are in line with the data ($\bar{\text{F}}^{\text{B}}_{(0,1,2,3; 19)}$ = 0, p = 1), the results from hypothesis test Type A show that the null-hypothesis is not rejected in favor of the order-constrained hypothesis ($\bar{\text{F}}^{\text{A}}_{(0,1,2,3; 19)}$ = 2.384, p = .1689).