5
$\begingroup$

I am confused whether one should exclude missing observations when adjusting p values for multiple testing. There seems to be no consensus among R function on whether to do this or not. stats::p.adjust(x) behaves different if you specify the default n = length(x) explicitly (NAs are counted) vs. if you do not specify the default explicitly (NAs are not counted). multtest::rawp2adjp(x) counts NAs. What is the correct behaviour?

EDIT: Clarifcation about what is meant by NA in this data was requested in the comments. The p values are calculated per residual after a mixed effects model fit to identify outliers in the data. The experimental procedure is complex and was undertaken by many experimentators in parallel, so errors are possible. Significance identifies data points, which are unexpectedly far from the fitted value, given the residual error is centered and normally distributed -> probably outlier [ref]. Some observations had to be removed prior to the fit, for example because there are too little observations for a certain group which causes problems with model fitting or were already missing prior to analysis because of experimental failure.

MWE

## Generate some p values and compare the three possibilities of 
## adjusting for multiple testing
n <- 10000
x <- pmin(rexp(n,rate =1/0.01), 1) # Generate some p values
x[sample(c(F,T), n/10, TRUE)] <- NA # delete some observations
# Three different methods of p value calculation
x1 <- p.adjust(x, method = 'holm') 
x2 <- p.adjust(x, method = 'holm', n = length(x))
x3 <- multtest::mt.rawp2adjp(x, proc = 'Holm')
x3 <- x3$adjp[order(x3$index),"Holm"]
# Compare p.adjust and mt.rawp2adjp
par(mfrow=c(1,2))
plot(x1, x3); title('p.adj(x) \n vs. mt.rawp2adjp')
plot(x2, x3); title('p.adj(x, len = length(x)) \n vs. mt.rawp2adjp')

Plot output showing x1 and x3 differ in their output while x1 and x2 do not

Appendix: Why does stats::p.adjust behave the way it does?

Beginning of p.adjust source code, R 3.4.3:

In the head, n is defined as length(p), however, R does not evaluate arguments until they are needed

function (p, method = p.adjust.methods, n = length(p)) 
{
    method <- match.arg(method)
    if (method == "fdr") 
        method <- "BH"
    nm <- names(p)
    p <- as.numeric(p)
    p0 <- setNames(p, nm)
    if (all(nna <- !is.na(p))) 
        nna <- TRUE

Here, p is stripped of all NAs, without n being needed up to this point

    p <- p[nna]
    lp <- length(p)

Now, n is used for the first time, which means length(p) is evaluated only now. Therefore if left to default setting, length(p[!is.na(p)]) is calculated.

    stopifnot(n >= lp)
[ remaining source code omitted ]
}
|cite|improve this question
$\endgroup$
  • $\begingroup$ I find a vector with p-values containing some NA's a bit strange. What does a missing value in the list mean? It was not measured, not evaluated, outside some range, the hypothesis test failed, etc? You can not correct for non-available p-values if you do not know the meaning of the non-availability. Note that the mt.rawp2adjp also has an option to take out the NA values (this parameter just happens to be false by default). $\endgroup$ – Sextus Empiricus Dec 18 '17 at 17:50
  • $\begingroup$ @MartijnWeterings The data vector stems from the design of the experiment. Some experiments failed -> NA. Some groups had too little data to have any intra-group variance -> NA. So NAs correspond to data being excluded from further consideration. $\endgroup$ – akraf Dec 18 '17 at 18:46
2
$\begingroup$

If the missing value causes the observation not to be included in the calculations of parameter estimates, it's contributing nothing to the end result (for better or worse) and should not be included in the p-value or adjusted p-value calculation. Its effect is the same as if it hadn't been included in the data set at all.

However, in some cases, missing values are not excluded from the calculations. They may be imputed, or, in the case of values that are censored (e.g., $x_1$ is not observed, but we know that $x_1 \geq 10$), included in the calculations but in a different way than if they had been observed. This is a murkier area. Clearly we wouldn't want either extreme - counted as if the observation was fully informative or counted as if the observation didn't exist at all - as the basis for p-value calculations, but it's not clear (and, indeed, problem-specific) how much "weight" between 0 and 1 the observation should get. Providing the ability to calculate adjusted p-values using the full observation count enables us to get a bound on the adjusted p-values that we'd have liked to calculate. If a particular value for a statistic isn't significant with a "sample size" = 100, it's not going to be significant with a "sample size" of less than 100 either, so the calculation with the sample size equal to the full number of observations does contain information useful for testing and evaluation.

To summarize: both calculations are useful, depending on the circumstances of the testing problem and how the estimation procedure treats missing values.

|cite|improve this answer
$\endgroup$
  • $\begingroup$ Thank you for your suggestion, however, what about the following counter example to your claim: "If a particular value for a statistic isn't significant with a "sample size" = 100, it's not going to be significant with a "sample size" of less than 100 either": The Bonferroni correction of the p values $\{0.100; 0.050; 0.020\}$ yields $\{0.300; 0.150; 0.060\}$, however, the same correction for $\{0.050; 0.020\}$ yields $\{0.100; 0.040\}$, hence, a p value of $0.020$ may be significant at the 5% level if the sample size is smaller, contrary to your claim? Or did I misunderstand you? $\endgroup$ – akraf Dec 18 '17 at 18:13
  • $\begingroup$ The sample size I'm referring to is that of the data, not of the number of p-values being corrected. It may be that I totally misread your question! If so, let me know, and I'll delete the answer as being non-responsive :( $\endgroup$ – jbowman Dec 18 '17 at 18:15
  • $\begingroup$ I think your answer makes sense even if I don't understand that specific claim. I have edited my original question. I will neither impute missing values nor are they censored. They are excluded from further consideration and are kept to facilitate the tracking of single data points through the analysis (their indices remain the same). Data points which are missing or excluded are set to NA downstream. Therefore, am I interpreting your answer correctly that you suggest to consider the length of the vector without the NAs for p value multiple testing correction? $\endgroup$ – akraf Dec 18 '17 at 19:24
  • $\begingroup$ Yes, you are correct, under those circumstances you would want to consider the length of the vector without the NAs. $\endgroup$ – jbowman Dec 18 '17 at 19:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.