# Should missing observations be included in the number of observations if correcting for multiple testing

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
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"]
par(mfrow=c(1,2))


# 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 ]
}

• 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). – Sextus Empiricus Dec 18 '17 at 17:50
• @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. – akraf - Reinstate Monica Dec 18 '17 at 18:46

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.
• 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? – akraf - Reinstate Monica Dec 18 '17 at 18:13