Controlling false discovery rates in factorial experiments with between-subjects and within-subjects tests

Background The False Discovery Rate (FDR) controls the expected number of false positives among the positive test results. It is not straightforward how to conduct a FDR controlling procedure in experiments with a factorial structure, while at the same time there are between-subjects and within-subjects factors. This is because there are P-values for different tests in one and the same response along with P-values for the same test and different responses. Findings We propose a procedure resulting in a single P-value per response, calculated over the tests of all the factorial effects. FDR control can then be based on the set of single P-values. Conclusions The proposed procedure is very easy to apply and is recommended for all designs with factors applied at different levels of the randomization, such as cross-over designs with added between-subjects factors. Trial registration NCT00959790


Findings
The control of false positive test results has enjoined considerable attention in the statistical literature. For an overview of methods in case there are many comparisons among treatments, we refer to [1]. More recently, Benjamini and Hochberg [2] and Storey and Tibshirani [3] proposed procedures that control the False discovery Rate (FDR). This is the expected fraction of false positive results among all positive results. The procedures are particularly suited for the analysis of multiple response variables. However, they do not address explicitly the case that there are several tests for one and the same response variable, let alone the presence of several sources of random variation that are to be used for the tests. The purpose of the present communication is to develop an explicit procedure for this case. http://www.biomedcentral.com/1756-0500/6/204 On the last day of each period, subjects completed a physical exercise test. At three time points, blood samples were taken. This defines a within-period factor 'time' , which is a repeated measurement factor.
Levels of 21 oxylipids were determined in the blood samples; the 168 samples were processed in a completely randomized order.

Statistical model
The data were studied using the following statistical model.
with μ pqr = μ 0 + βB p + δD pq + γ B p D pq + τ r T r + η r T r B p + θ r T r D pq + κ r T r B p D pq .
In formula (1), y pqr is the level of an oxylipid from subject p (p = 1 . . . 28), period q (q = 1, 2) and time r (r = 0, 1, 2). The measurement is the sum of an expected value modeled with μ pqr and random contributions modeled with the terms e 2p , e 1pq , and e 0pqr . The expected value of the measurement y ijk is detailed in formula (2). We make a distinction between parameters, which are to be estimated from the data, and experimental variables, which indicate the BMI group, the diet, and the time point relevant to the observation. There are 11 parameters, given in Greek alphabet, and four experimental variables, given in Latin alphabet. First, the average difference between the lean and obese groups is modeled with parameter β and experimental variable B p . This variable takes the value 1 if subject p is obese and 0 otherwise. The parameter β thus models the increase in oxylipid level for an obese subject relative to a lean subject.
The average difference between diet 1 and diet 2 is modeled with parameter δ and experimental variable D pq . This variable takes the value 1 if subject p is given diet 2 and 0 otherwise. The parameter δ thus models the increase in oxylipid level for diet 2 relative to diet 1.
Next, the parameter γ models the interaction between diet and BMI group. If γ = 0, the difference between the diets does not depend on the BMI group. If γ = 0, the difference between the diets depends on the BMI group. The parameters that model the average change over time are τ 1 and τ 2 , respectively (τ 0 is taken to be zero). The corresponding experimental variables are T 1 and T 2 . The first of these takes the value of 1 at time point 1 and 0 otherwise; the second experimental variable takes the value of 1 at time point 2 and 0 otherwise. So the time changes are modeled relative to time point 0.
Further, the parameters η r , θ r and κ r model the interaction between BMI group and time, the interaction between diet and time and the three-factor interaction between BMI group, diet and time, respectively.
The three random terms in formula (1) model the random error between subjects, the random error within subjects and the random error within periods, respectively. We assume that the three random terms are independent of each other and normally distributed with variances σ 2 2 , σ 2 1 and σ 2 0 , respectively. The subjects can be considered as random samples from two specific populations. Therefore, the 28 e 2p are independent and we can validly carry out an F test to assess the difference in BMI level between the two populations.
Further, the subjects were randomly allocated to a treatment order. Therefore, the 28 differences e 1p1 − e 1p2 are independent and we can validly carry out F test to assess the effect of diet and its interaction with BMI group.
There could not be a random allocation of the time points to the blood samples. For this reason, the correlations between e pq0 and e pq1 , between e pq0 and e pq2 , and between e pq1 and e pq2 might not be equal. This would invalidate the analysis of variance F tests for the main effect of time and the interactions involving time. Fortunately, the problem posed by unequal correlations can be solved by applying a correction factor to the degrees of freedom for the F-tests due to Greenhouse and Geisser [4].
Sometimes, other assumptions on the random terms are reasonable, which may lead to other denominators of the F tests being appropriate. We refer to [5] for an extensive discussion of this issue.

Analysis of variance
An analysis of variance for one of the oxylipids, namely arachidonic acid, is given in Table 1.
The first two columns of the table lists the three error strata and the 10 sources of variation present in the data. An error stratum collects all effects that are tested against the same variance; see [6] for a formal definition of a stratum.
All the effects that are measured by contrasting subjects are in the between-subjects stratum. The difference between the groups, which constitutes the BMI main effect modeled with β in formula (2), is tested against the random error between subjects.
Each of the two diets was given to each of the subjects. For this reason, the main effect of diet (modeled with δ in formula (2)) and the interaction between BMI group and diet (modeled with γ ) are tested against the random error within subjects.
Finally, the three time points at which blood samples were taken define a third factor, time, whose main effect (modeled with τ 1 and τ 2 ) is to be tested against a random error within periods. The interactions between BMI category and time (modeled with η 1 and η 2 ), and between diet and time (modeled with θ 1 and θ 2 ) are also tested against this random error. http://www.biomedcentral.com/1756-0500/6/204 The same is the case for the three-factor interaction (κ 1 and κ 2 ). All these effects are in the within-periods stratum.
Further columns in the table give the degrees of freedom (df ) for each source of variation, the corresponding mean square (MS), the value of the individual F-ratio (F ij ), and the P-value (P ij ). The index i points to the error stratum, while the index j points to the F-test within a stratum.
The four F-tests in the within periods stratum were carried out using the Greenhouse-Geisser as a correction factor to the degrees of freedom. The calculation of this factor is implemented in most major statistical packages. Here, = 0.8103. Accordingly, the degrees of freedom needed for the calculation of the P-values for time and its interactions with the other two factors were 0.8103 × 2 = 1.6206 for the numerator and 0.8103 × 104 = 84.2712 for the denominator.
Under an individual false positive error rate of 0.05, the outcome for the main effects of BMI and time are highly significant. There is no evidence that the main effect of diet or any interaction effect is statistically significant.

FDR in factorial experiments with a single stratum
A factorial structure of the study design permits the evaluation of main effects and interactions. For two factors and m response variables there are thus 3 m tests to carry out.
The tests for main effects might not be needed once the interaction is declared statistically significant. This is an important notion, because the total number of the tests is a parameter for the FDR procedure. One could start with a procedure for the m tests on active interactions only. In a second step, the variables with significant interactions, s 1 , say, are removed from further consideration, and we are left with m − s 1 variables not having a proven interaction among the factors. We could then consider applying the FDR procedure on 2(m − s 1 ) main effect tests. However, it is unclear what the performance criteria of the joint first and second step are.
To circumvent the above problem, we propose to replace the three tests with one omnibus F-test to see whether the treatments differ. So we initially forget about the factorial structure of the treatments and just check whether there are differences between the treatment groups. For the responses where this is indeed the case, we suggest a follow up that does use the factorial structure, and assess the main effects and interactions using the corresponding P-values.
The proposed replacement of individual statistical tests can be carried out easily if all the comparisons between the experimental groups are tested against one and the same error. This is the case if there is just one error stratum, but also if there are several strata while the effect tests involve only one stratum. However, the proposed replacement is not straightforward to apply when effects are tested in several strata. For example, in the motivating study, the error used to test the contrast between lean and obese is different from the error used to test the contrast between the diets. This issue is discussed next.

FDR in factorial experiments with several strata
We propose calculating a combined P-value over all the F tests of a response variable as follows: Under the null hypothesis, this is an F statistic with n i = n ij degrees of freedom for the numerator and d i degrees of freedom for the denominator. 3. Suppose that the combined F -test in stratum i has an associated P -value of U(a, b) denotes a uniform distribution with minimum a and maximum b. 4. Combine the P -values by calculating is a random variable following a χ 2 [2E] distribution. 6. Apply an FDR control method to the list of overall P -values. 7. For variables selected in step (6), study all P ij to see which factors or interactions contributed to the significance of T E .
The procedure to combine P-values is due to Fisher [7]. See [8] for other options to combine P-values. The crucial condition for a correct application of Fisher's procedure is the independence of the P-values. This condition is satisfied if the tests involve different error strata.
Step 6 in our procedure results in a set of variables with an expected fraction of at most α of false positive results among all positive results, where α is the desired level of protection. So the FDR procedure selects variables that show factorial effects. However, the FDR procedure does not operate on the overall list of decisions based on the individual P ij studied in Step 7. In this aspect, our procedure is analogous to Fisher's protected least significance difference procedure [1] in one-way analysis of variance, because, in the latter procedure, differences between treatment groups are tested only if the overall F-test is statistically significant.

Application
We apply the proposed procedure to the arachidonic acid response of the motivating example. In the betweensubjects stratum, there is nothing to combine, because there is just a single test carried out in this stratum. Recall that the P-value for the main effect of BMI is 0.004.
The two F-tests in the within-subjects stratum are combined by adding the mean squares of 0.0091 and 0.6465, dividing by 2, and dividing the result by the error mean square of 0.1887. The F-value for this stratum is 1.74, with two degrees of freedom for the numerator and 26 degrees of freedom for the denominator. The P-value is 0.20. This P-value suggests an absence of treatment effects.
For the within-periods stratum, we multiply the mean squares for time, BMI × time, diet × time, and the threefactor interaction BMI × diets × time with 2, add up and divide by 8. This results in a combined mean square of 1.2463. This mean square is tested against the error mean Finally, the three P-values are to be combined to one overall value. We take -2 times the natural logarithm, and add up. This gives X 2 = 87.945. The reference distribution for this statistic is the χ 2 [6] distribution. The statistic has a P-value of 8.09 × 10 −17 .
All the overall P-values according to the proposed procedure are shown in Figure 1. For 12 of the oxylipids, including arachidonic acid, P < 0.05. The application of the FDR-controlling procedures of [2] and [3], is visualized in Figure 2. The P-values are ordered and plotted against their order number. We restrict attention to the values below 0.2, and we use a boundary value of 0.05 for both procedures. The lower line gives the boundary values for the Benjamini-Hochberg procedure [2]. The largest Figure 2 Rejections for two FDR procedures. P-values below lower line: rejected by the Benjamini-Hochberg procedure [2]; P-values below upper line: rejected by the Storey-Tibshirani procedure [3]. http://www.biomedcentral.com/1756-0500/6/204 P-value below the line has order number 10. So the procedure reveals that 10 out of 21 oxylipids are affected by the experimental factors. The upper line in Figure 2 bears on the procedure of Storey and Tibshirani [3]. When compared with the Benjamini-Hochberg procedure, two more P-values are included in the set with q < 0.05. Note that the set now includes all oxylipids for which P < 0.05. This is not generally the case, however.
Some authors would favor error-control methods that are more conservative than FDR. For example, the wellknown Bonferroni correction would compare all 21 combined P-values with an error rate of 0.05/21. Clearly, the proposed FDR methods are more lenient than the Bonferroni correction in declaring that a variable is significantly affected by the study factors.
We like to point out that both FDR controlling procedures are sensitive to strong negative correlations between the P-values; see [9]. For the oxylipids, this is not really an issue because the average pairwise correlation among the oxylipids was +0.1. With three exceptions, all correlations were above -0.3; the smallest exceptional value was -0.5. We therefore think that our application of the FDR control is justified.
As a final issue, we had an equal interest in all the oxylipids and all the model parameters. In case of variables or parameters of primary interest, one option is to include only these variables or parameters. This will make the procedure more powerful, because non-significant values of the F statistic that are not of interest will tend to reduce the overall test statistic. Alternatively, there are options to introduce weights to the variables or parameters other than 1 for those of primary interest and 0 for those of secondary interest. However, a discussion of these options is beyond the scope of the present paper.

Availability of supporting data
The data set supporting the results of this article is included within the article and its additional file called FDR_overall_Pvalue_calculation.xlsx. The additional file shows for each of the 21 oxylipids, first the F ij values arranged in seven rows and 21 columns. The columns correspond to the oxylipids and the rows correspond to the seven statistical tests for each individual oxylipid. Next, the 21 values for the Greenhouse-Geisser epsilon statistic are given. Then we give the P-values for each of the three error strata arranged in three rows and 21 columns. The columns correspond to the oxylipids and the rows correspond to the between-subjects, withinsubjects and within-period strata, respectively. Finally, we give the value of the statistic T E , as calculated in step 4 of the proposed procedure, and the corresponding overall P-value for the factorial effects of the 21 oxylipids.
Abbreviations BMI: Body mass index; FDR: False discovery rate.