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')
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 NA
s, 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 ]
}