2. 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.

3. 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 (1)

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 (2)

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.

4. Test statistic

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

where for and . Then we consider the test statistic (10)

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

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 (12)

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, (13)

under . Straightforward calculations show , and therefore (14)

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 (15) (16)

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 (17) (18) (19) (20)

for and , and . Those four estimators satisfy the following statements: (21) (22) (23) (24)

That is, each of these estimators is unbiased.

5. 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:

High-dimensional repeated measures

All authors

Martin Happ, Solomon W. Harrar & Arne C. Bathke

https://doi.org/10.1080/15598608.2017.1307792

Published online:

21 April 2017

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

Display full size

• 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).

6. 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

• 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.

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.