# Multiple testing adjustment in linear models

I have troubles with multiple testing correction (For example, bonferroni adjustment. I'm aware that it's not most correct one). I have 1000 linear models like this:

Call:
lm(formula = count ~ type, data = table_lm_3)

Residuals:
Min       1Q   Median       3Q      Max
-1.96774  0.03226  0.03226  0.03226  0.66667

Coefficients:
Estimate Std. Error  t value     Pr(>|t|)
(Intercept)          1.96774    0.04983  39.488     < 2e-16 ***
typeA               -0.63441    0.16776  -3.782    0.000299 ***
typeE               -1.44393    0.11610  -12.437    < 2e-16 ***
typeN               -1.30108    0.28189  -4.616    1.47e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3924 on 80 degrees of freedom
Multiple R-squared:  0.6831,    Adjusted R-squared:  0.6712
F-statistic: 57.49 on 3 and 80 DF,  p-value: < 2.2e-16


In each model i have one overall p-value of Fisher test and three p-values from pairwise t-tests. So, i really confused which p-values i must adjust. There are three possibilities:

1) Correct overall p-value of Fisher test for 1000 observations (p-value*1000)

2)Adjust each p-values of pairwise t-tests for 1000 observations (p-valueA*1000,p-valueE*1000,p-valueN*1000)

3)Adjust each p-values of pairwise t-tests for 3000 observations (p-valueA*3000,p-valueE*3000,p-valueN*3000). Because i have 1000 models and inside each one i have 3 pairwise comparisons.

4)Adjust each p-values of pairwise t-tests for 3 observations (p-valueA*3,p-valueE*3,p-valueN*3). Because i have three pairwise comparisons inside each model.

Which option will be correct in my case?

I have a four levels of the factor type. For now we are curious just in differences of typeA, typeE and typeN against the reference factor. I mentioned we work with bioinformatic data. Lets imagine we have a gene A that consists from 3 parts (part 1, part 2,part 3). Each part could have a count from 0 to 8 in different samples. So it's defenetely were the integers, there are not any floats (no continuous scale). Initially we have a table with columns name and count like

Gene A part 1                  8
Gene A part 2                   1
Gene A part 3                   2
Gene B part 1                  0
Gene B part 2                   1
Gene B part 3                   2
Gene B part 4                   2


But because we was interesting to count genes but not a separate parts we did the follows:

Gene A count= mean(Gene A count part 1  + Gene A count part 2  + Gene A count part 3)


In our example it's Gene A count= mean (8+1+2)=3.67 In result we turned to floats. So we tested our hypotheses on huge matrix with floats values (continuous scale). But in our case any float value should be >=0 and <=8. There are not values greater than 8 and less than 0. We have a four groups of samples (Control and groups A,E,N) with different numbers of samples in each particular group. Later numbers of groups will be increased up to 6-8. And we are curious to find contrasts of Gene counts between Control group and other groups under investigation. We have not only one gene A but lets say a 1000 of different ones.

• Why does this need to be done as 1000 separate models? Is there some way to incorporate all your data into a single model that takes into account whatever is presently distinguishing among the 1000 models? – EdM Jul 17 '15 at 16:00
• Unfortunately not. It's a common practice of contemporary bioinformatics studies. Usually we perform numerous hypotheses testing to obtain the results. For example after exact Fisher test , we will have a thousand of independent p-values from each one of them. In this case quite clear on which observation numbers we should correct the p-values. But inside each lm model we have a several p-values and i'm not sure which one i must correct and on which total number of comparisons. – Denis Jul 20 '15 at 13:22
• The way the code and example are presented, it seems that you have 4 levels of the factor type. Are you only interested in differences of typeA, typeE and typeN against the reference factor, or also the differences among the 3 listed types? Are these actual count data (values like 0, 1, 2, 3...) or are the values on a continuous scale? Is the number of cases in each of the 1000 models the same? Do you have a specific hypothesis you are testing? What are your relative risks between false positives and false negatives? It would help to edit the question based on your answers to this comment. – EdM Jul 22 '15 at 22:11

## 2 Answers

A danger in analyzing this as 1000 separate models, one for each gene, is that you lose the ability to pool information from different comparisons.

As an example of a combined approach, consider the MulCom package in R. This package performs the types of comparisons you want, in the context of examining differential expression of many genes among multiple groups with respect to a reference, as an integrated analysis. This has the advantage of pooling information on error terms for statistical tests, rather than using the highly variable error estimates obtained in individual comparisons.

A combined approach also may simplify the correction for multiple testing. In the type of study you are doing, it is often more useful to control the false discovery rate (FDR, the fraction of "significant" results that aren't true) rather than correct individual p-values for multiple testing. A combined analysis like that done by MulCom allows for permutation tests to estimate the FDR directly from the data, and allows for optimization of test parameters. FDR control is well accepted in bioinformatics.

For FDR analysis the answer to your original question, about how many comparisons need to be taken into account when correcting p-values, becomes perhaps a bit simpler. You count the number of positive results of all comparisons (against reference) at your choice of parameter values, then determine the median number of positive results in the same analysis applied to multiple randomly permuted sample groups. The ratio of median positives in the permuted groups to the positives in the original data is the estimate of the FDR. See the paper for further details, or examine the R code.

With your count data rather than standard expression data, you may need to consider whether MulCom or a similar expression-focused analysis method is appropriate. Since you are willing to average counts across "parts" of genes to get a single value for each gene, however, your data then look like expression data appropriate for this type of analysis. If you decide that you need to do a different type of analysis because you have count data, you will be better off following an approach like that used by MulCom rather than the gene-by-gene correction of p-values that your question anticipated.

After discussion the problem with expert in biostatistics in my institute, i definitely can say that the direct answer in my case is number one. "Correct overall p-value of Fisher test for 1000 observations (p-value*1000)".