Closed Testing {cherry} | R Documentation |
Performs the closed testing procedure for user-specified local test.
closed (test, hypotheses, alpha = 0.05, adjust=FALSE)
test |
A function that performs the local test. The function should accept a subvector of the |
hypotheses |
Identifiers of the collection of elementary hypotheses. |
alpha |
The significance level of the test procedure. If set to |
adjust |
Whether adjusted p-values should be calculated. |
The function closed
performs the closed testing procedure on the collection hypotheses
, testing all their intersection hypotheses and controlling the familywise error rate.
The function closed
returns an object of class closure
.
The number of intersection hypotheses is exponential in the number of elementary hypotheses. The number of elementary hypotheses is therefore limited to log2(.Machine$integer.max+1)
(typically 31) for computational reasons.
It is possible to set both adjust
to TRUE
and specify alpha
. In that case, adjusted p-values are calculated up to a value alpha
; all higher p-values are set to 1.
Jelle Goeman: j.j.goeman@lumc.nl
Goeman and Solari (2011) Statistical Science 26 (4) 584-597.
# Example: the birthwt data set from the MASS library
# We want to find variables associated with low birth weight
library(MASS)
fullfit <- glm(low~age+lwt+race+smoke+ptl+ht+ui+ftv, family=binomial, data=birthwt)
hypotheses <- c("age", "lwt", "race", "smoke", "ptl", "ht", "ui", "ftv")
# Define the local test to be used in the closed testing procedure
mytest <- function(hyps) {
others <- setdiff(hypotheses, hyps)
form <- formula(paste(c("low~", paste(c("1", others), collapse="+"))))
anov <- anova(glm(form, data=birthwt, family=binomial), fullfit, test="Chisq")
res <- anov$"Pr("[2] # for R >= 2.14.0
if (is.null(res)) res <- anov$"P("[2] # earlier versions
res
}
# perform the closed testing
cl <- closed(mytest, hypotheses)
cl
# how many variables among a chosen set are associated with the response?
pick(cl, c("ht", "lwt", "smoke", "ui"))
# adjusted p-values and a confidence distribution
cl <- closed(mytest, hypotheses, alpha=NA)
pick(cl, c("ht", "lwt", "smoke", "ui"))