High-dimensional repeated measures.

Recently, new tests for main and simple treatment effects, time effects, and treatment by time interactions in possibly high-dimensional multigroup repeated-measures designs with unequal covariance matrices have been proposed. Technical details for using more than one between-subject and more than one within-subject factor are presented in this article. Furthermore, application to electroencephalography (EEG) data of a neurological study with two whole-plot factors (diagnosis and sex) and two subplot factors (variable and region) is shown with the R package HRM (high-dimensional repeated measures).

Introduction

Classical methods such as the Huynh–Feldt (Huynhand Felat 1976) or the Greenhouse–Geisser (Greenhouse and Geisser 1959; Geisser and Greenhouse 1958) approximation do not maintain the type I error rate in general in high-dimensional repeated measures. “High-dimensional” means that the number of repeated measurements d is larger than the number of subjects per group. In an article by Happ et al. (2016), various simulations have been performed, demonstrating the problems that these classical methods exhibit, particularly in unbalanced, heteroscedastic, and high-dimensional designs. A new approximate test had been proposed by Brunner et al. (2012) for two treatment groups. This approach was extended in Happ et al. (2016) to multiple treatment groups, along with asymptotic results for the situation such that the dimension d of repeated measurements and the total number of subjects N tend to infinity.

In the present work, the test proposed in Greenhouse and Geisser (1959) is generalized to multiple whole and sub-plot factors. Furthermore, the package HRM for the software environment R (R core Team 2016) is presented and application is shown on high-dimensional electroencephalography (EEG) data from Staffen et al. (2014).

Data example

From 160 subjects, electroencephalography (EEG) data have been measured in a study by Staffen et al. (2014). The patients were classified as either having Alzheimer’s disease (AD), mild cognitive impairment (MCI), or subjective cognitive complaints without any significant clinical deficits (SCC+ or SCC–). These groups are compared using four EEG variables (activity, complexity, mobility, brain rate) measured at 10 different brain regions (frontal left, frontal right, parietal left, parietal right, temporal left, temporal right, occipital left, occipital right, central left, central right). In total, measurements have been taken per subject. The numbers of subjects per diagnostic group and sex are given in Table 1. Note that the number of subjects in each group is smaller than the number of repeated measurements per subject. This high-dimensional data set is used as an example in section 6 in order to illustrate application of the new method proposed in this article.

Table 1.

Number of subjects per diagnostic group and sex.

	AD	MCI	SCC+	SCC-
Male	12	27	14	6
Female	24	30	31	16

Statistical model

We focus on the case of two crossed whole-plot factors A and S and two subplot factors B and C, to keep the notation simple. Also, interpretation of the interaction effects might get difficult if too many factors are used in the model.

We refer to factor A with a factor levels as the group factor and to S with s factor levels as the subgroup factor. Factors B and C are referred to as variable and region with b and c factor levels, respectively, where C is crossed with B. The terms used here for the factors are motivated by the EEG data from section 2, which are analyzed in section 6.

For the statistical model, we consider the independent random vectors

which are the dimensional data vectors of subject k in the ith group and jth subgroup where is a c dimensional random vector for , , and .

In each group, there are subjects, and overall there are subjects. Therefore, the random variables are the observations from the kth subject in the ith group, , and the jth subgroup, , for the variable and region . The expectations are vectors, and the variance–covariance matrices are matrices. The sample mean of the groups is denoted by the vector ; the sample variance is given by the matrix

We consider hypothesis tests for testing main effects (e.g., main effect of the region or group), as well as interaction effects of two, three, or four factors. Thus, in total, there are 15 different hypotheses we are interested in. We denote the respective null hypotheses by , where

In order to give two examples, the null hypothesis is the hypothesis of no main effect of the factor A. As we are interpreting the factor A as the group factor, we can describe as the null hypothesis of no main effect of the factor group. On the other hand, the null hypothesis is the hypothesis of no interaction between the factors A, S, B, and C. A more detailed description of these hypotheses is provided in subsections 3.1 and 3.2.

It is well known that each of these linear hypotheses can be formulated as a quadratic form with a projection matrix . For example, is the matrix used for testing the interaction effect of the factors A, S, and B. To be more precise, the matrix is of the form

where denotes the Kronecker product. The part to the left of the Kronecker product sign, namely, the matrix , describes the aspects of the hypotheses that are concerning the whole-plot factors, while the matrix to the right, , corresponds to the subplot factors. For stating the hypotheses, let denote the d-dimensional identity matrix. Furthermore, we define , the d-dimensional vector consisting only of ones, and the matrix , as well as the centering matrix . The representation of as a Kronecker product of four matrices is useful as it allows us to easily construct matrices used for testing no main and interaction effects. For example, let us consider the hypothesis of no interaction effect between the sub-plot factors B and C. In this case , , and . The dimensions of these four matrices correspond to the respective levels of the factors A, S, B, and C. Since we only want to test the interaction effect between B and C in this example, we can simply average over the levels of the factors A and S. This is achieved by using . For the matrix , we can use a Kronecker product of centering matrices with dimensions according to the levels of factors B and C. More details on the construction of the matrices are given in section 3.1.

Hypotheses on main effects

The hypothesis for testing the main effect of the factor A can be written as

or equivalently as

The projection matrix for testing this hypothesis is given by

and the vectors , , are given by

for all , and . The other hypotheses for testing main effects can be written in a similar manner using the following matrices

Derivations of these equivalent results are done in a similar manner as in Happ et al. (2016) and can be easily extended to more than four factors.

Hypotheses on interaction effects

Similar to the main effects, we can write the null hypotheses regarding interactions between two or more factors as quadratic forms. Specifically, the hypothesis of no interaction effect involving all four factors (A, S, B, and C) can be written as

for all , , , and . This formulation is equivalent to the quadratic form

where the matrix is given by

In general, the matrix , , can be written as a Kronecker product of four matrices, where each of these four matrices operates on the components of a vector representing the levels of a specific factor. If a particular factor is selected for the formulation of an interaction or main effect, its corresponding matrix will be the centering matrix P with dimension r. Otherwise it will be , where r is the number of levels of the factor , and hence in our case . In the peceding example, we are testing the hypothesis of no interaction effect of all four factors. Therefore, the centering matrix is used for each of the four factors .

The null hypothesis of no interaction effect between the three factors A,S, and C can be written as

for all , , and . This is equivalent to

Therefore, the matrix is given by

Here, only the factor B is not selected in formulating the interaction effect. Therefore, we need to average over the levels of factor B, and consequently its corresponding matrix is .

Hypotheses for all other interaction effects can be constructed in a similar manner. The equivalent formulations can be derived analogously to those for the main effects.

Test statistic

To test the hypotheses, already stated, we propose a generalization of the procedure introduced in Happ et al. (2016). Specifically, let

where for and . Then we consider the test statistic

where is one of the 15 projection matrices from before, and

where The numerator and denominator can each be approximated by scale multiples of random variables, where the degrees of freedom are allowed to be any positive real number. Thus, we approximate under

Calculating expectation and variance of the numerator and denominator by using results about the distribution of quadratic forms—see, for example, Searle et al. (1997) — we obtain expressions for the degrees of freedom and of these two distributions. Now noting that the numerator and denominator are independent,

under . Straightforward calculations show , and therefore

under . As this approximation goes back to Box (1954a; 1954b), it is referred to as Box approximation. Let and as defined in section 3 and let and denote the diagonal elements of and , respectively. Then the degrees of freedom are given by

More details regarding the decomposition of are provided in sections 3 and 3.1. In Eqs. (15) and (16), the matrices (respectively ) are in general unknown and have to be estimated. Instead of estimating the matrices directly, we simply estimate their traces. Appropriate trace estimators have already been used in Brunner et al. (2012) and Happ et al. (2016), and they are given by

for and , and . Those four estimators satisfy the following statements:

That is, each of these estimators is unbiased.

R package HRM

The test described in section 4 has been implemented for the statistical software environment R (R Core Team 2016) in the package HRM (high-dimensional repeated measures, Happ et al. 2016) and is available via CRAN. There are currently four functions implemented:

hrm.test.matrix (data, alpha),

hrm.test.dataframe (data, alpha, group, subgroup, factor1, factor2, subject, response),

hrm_test (formula, data, alpha = 0.05, subject),

hrm.plot (data, group, factor1, subject, response, xlab = “dimension”, ylab = “means”).

The function hrm.test.matrix can be used for testing main and simple treatment effects (whole-plot factor), the main effect of the factor time (subplot factor), and the interaction between treatment and time. A more detailed description and thorough simulation study for this particular function can be found in Happ et al. (2016).

Tests for models with up to two whole-plot factors or two subplot factors can be performed with hrm.test.dataframe. Here the variables group, subgroup, factor1, and factor2 are used to indicate which columns of the data matrix are used as factors. If, for example, factor2 is not specified, then only one subplot factor, factor1, is being used for testing. But at least one wholeplot factor and one subplot factor need to be specified. The variable subject is necessary to identify the subjects, and response indicates in which column the measured response variable of the data set is stored. hrm.test.dataframe tests main effects as well as all interaction effects. If four factors are used, 15 hypotheses are tested by this function. In this regard, it should be noted that the results are not adjusted for multiple testing.

For the function hrm_test, a formula object can be used instead of specifying in which columns the whole- and subplot factors can be found. A formula object is of the form

In this context

is an abbreviation for

where : refers to the interaction effect of and subplot factor. The function hrm_test supports up to two whole- and three subplot factors, but at least one of each has to be used and the maximum number of factors is four. The function automatically detects whether a column is a whole- or a subplot factor. For example, by using the notation from the EEG data from section 2 and the formula

the function hrm_test is testing for the main effects of the factors group, sex, variable, and region and the interaction effect between group and sex. By using

we only test for the interaction effect between group, sex, region, and variable.

The only difference between the function hrm.test.matrix and the functions hrm.test.dataframe and hrm_test is the way the data need to be organized. For the former function, one needs to use one data.frame for each group, the rows are considered to be the subjects, and the columns are the repeated measurements. For the functions hrm.test.dataframe and hrm_test, all data need to be stored in one data.frame, which contains columns as described before.

The function hrm.plot provides an easy and fast way to plot the profiles in case of exactly one whole- and one subplot factor. It is basically just a wrapper for the package ggplot2 (Wickham 2009). Similar to the function hrm.test.dataframe, the variables group, factor1 need to be specified in order to indicate which columns of the data.frame should be used as factors. An example of a profile plot generated by this function can be seen in Figure 1 for the EEG data. The corresponding R code is as follows:

Figure 1.

Profile plot of the diagnostic groups SCC+, SCC-, MCI and AD.

library(HRM)

hrm.plot(EEG, group=“group”, factor1=“dimension”, subject=“subject”,

response=“value”)

From Eq. (10), we can see that we have to calculate the trace of a matrix product. If we have many repeated measurements , say , then these matrices will have dimension . In this case, the calculation of matrix product and trace may take very long. A much faster method uses the identity

where are matrices for , and * is the Hadamard–Schur product.

Another problem is presented by calculating for and if the number of repeated measurements is large. We can use the fact that empirical variance–covariance matrices can be written as quadratic forms with an idempotent matrix :

where is a matrix containing the data of group i and subgroup j. Then using the fact that the trace of a matrix is invariant under cyclic permutations, we obtain

where is an matrix. These two improvements for working with high-dimensional data have already been used by Brunner et al. in a technical report associated with Brunner et al. (2012).

Data example revisited

We now apply the R package HRM from section 5 to the EEG data described in section 2. Here, the data are stored as a data.frame, and the repeated measurements corresponding to a subject are each given by a row. There is one column specifying each of group (diagnosis), subgroup (sex), the first subplot factor (variable), and the second subplot factor (region). The structure of the data.frame can be seen in Table 2. All 15 hypotheses are calculated with the R code

Table 2.

Overview of the data matrix.

Group	Value	Sex	Subject	Variable	Region
SCC+	2.80	F	1	1	1
SCC+	2.29	F	1	1	2
SCC+	1.47	F	1	1	3
SCC+	1.55	F	1	1	4
SCC+	2.16	F	1	1	5
SCC+	2.14	F	1	1	6

Note. The variables are coded with 1 to 4; the 10 regions are coded with 1 to 10. Sex is either female (F) or male (M) and the group is either SCC+, SCC–, MCI or AD.

library(HRM)

hrm_test(value group * sex * region * variable, data=EEG,

subject=“subject”)

where EEG is the data.frame from Table 2. The output of this code is given in Table 3.

Table 3.

Results of the tests for main and interaction effects without adjustments for multiple testing.

Hypothesis	df1	df2	crit	test	p value
group	3.04	116.29	2.67	1.66	0.18
sex	1.00	116.29	3.92	6.05	0.02
variable	1.51	178.72	3.36	4930.69	<0.001
region	5.59	246.82	2.18	196.60	<0.001
group: sex	3.04	116.29	2.67	0.18	0.91
group: variable	4.52	178.72	2.33	3.34	0.01
group: region	14.09	246.82	1.73	0.80	0.67
sex: variable	1.51	178.72	3.36	4.89	0.02
sex: region	5.59	246.82	2.18	0.99	0.43
variable: region	7.76	299.89	1.99	143.58	<0.001
group: sex: variable	4.52	178.72	2.33	0.26	0.92
group: sex: region	14.09	246.82	1.73	0.58	0.88
group: variable: region	18.75	299.89	1.63	1.04	0.41
sex: variable: region	7.76	299.89	1.99	1.26	0.27
group: sex: variable: region	18.75	299.89	1.63	0.74	0.77

It is not surprising to neurologists that the variables activity, complexity, and mobility yield rather different values. Our test confirms this by rejecting the hypothesis of no main variable effect. But also the factor region (e.g., parietal left, parietal right, frontal left, frontal right, etc.), and the interaction between variable and region are highly significant. From the whole-plot factors only sex is significant with a p value of about 0.0154. Furthermore, the interaction between group and variable, and the interaction between sex and variable are also significant with p values of 0.0085 and 0.0152.