In statlab/permuter: permuter: Permutation Tests

Vignette Info

This example comes from the textbook Section 7.11.4 (p 361-364).

This illustrative example, with a few changes made, comes from Pesarin (1991b). In it $n_1 = 20$ seeds of a given plant are sown in a standard (untreated) plot of soil and $n_2 = 20$ more seeds are sown in soil treated with a fertilizer. By assumption, seeds are randomized to treatments. Also imagine that the expected treatment effects of interest are:

a. improved germination probability, indicated by variable $O$ b. improved index of production, indicated by the bivariate quantitative response $(X, Y)$, where $X$ is an index measured by the weight of adult plants and $Y$ is an index measured by the total surface area of their leaves.

Note that $O$ is typically binary categorical, taking value $1$ if the seeds germinate and $0$ otherwise. $X$ and $Y$ are positive quantities when $O = 1$. The underlying non-degenerate distribution $P$ is presumed to be unknown. $O$ also plays the role of inclusion indicator of missing values with respect to quantitative responses $(X, Y)$; in fact, the pair $(X_{ij}, Y_{ij})$ is actually observed on the $ji$th unit if $O = 1$, whereas if $O = 0$ (i.e. the $ji$th unit did not germinate), $(X_{ij}, Y_{ij})$ cannot be observed and so appears to be missing or censored, $i = 1,\dots,n, j = 1,2$. In a way, the data here appears to be jointly censored by a random mechanism, depending on treatment.

In such a context, we may write the test hypotheses as $$ \begin{aligned} H_0 &: { P_1 = P_2} = {(O_1, X_1, Y_1) \,{\buildrel d \over =}\, (O_2, X_2, Y_2)} \ &= \left\lbrace \left[ O_1 \,{\buildrel d \over =}\, O_2\right] \bigcap \left[ (X_1, Y_1) \,{\buildrel d \over =}\, (X_2, Y_2) \mid \mathbf{O} \right]\right\rbrace \end{aligned} $$ against $$H_1 : \left\lbrace \left[ O_1 \,{\buildrel d \over <}\, O_2\right] \bigcup \left[ X_1 \,{\buildrel d \over <}\, X_2 \mid \mathbf{O} \right]\bigcup \left[ Y_1 \,{\buildrel d \over <}\, Y_2 \mid \mathbf{O} \right]\right\rbrace$$

because, conditionally on germination, the expected treatment effects are assumed to produce higher values in both component variables. Moreover, let us assume that partial expected effects on the first component variable $\left[X \mid \mathbf{O}\right]$ are marginally homoscedastic and additive on location, whereas effects on $\left[Y \mid \mathbf{O}\right]$ are presumed to act on both location and second moment (as in multi-aspect testing problems).

These side assumptions are, on the one hand, consistent with the impression that observed second-variable data show increments in both mean and variance; on the other hand, they are consistent with the fact that most of the weight in this particular kind of plant comes from its leaves. In addition, these side assumptions contribute towards showing the versatility of NPC methods. Hence, the set of sub-alternatives may be written as:

$$ \begin{aligned} H_{11}^\mathbf{O} &: \left\lbrace \mathbf{E}(O_1) < \mathbf{E}(O_2)\right\rbrace \ H_{12}^{(\mathbf{X},\mathbf{Y})\mid\mathbf{O}} &: \left\lbrace \mathbf{E}(X_1) < \mathbf{E}(X_2)\right\rbrace \ H_{13}^{(\mathbf{X},\mathbf{Y})\mid\mathbf{O}} &: \left\lbrace \mathbf{E}(Y_1) < \mathbf{E}(Y_2)\right\rbrace \ H_{14}^{(\mathbf{X},\mathbf{Y})\mid\mathbf{O}} &: \left\lbrace \mathbf{E}(Y_1^2) < \mathbf{E}(Y_2^2)\right\rbrace \end{aligned} $$

library(knitr)
library(permuter)
data(germina)
set.seed(50)

B <- 1000
n <- nrow(germina)
p <- ncol(germina) - 1
tst <- rep(0, 4)
distr <- matrix(0, nrow = B, ncol = 4)

The data from two groups are contained in the dataframe germina. The number of samples in each group is $n = 20$. The number of non-missing data (germinated seeds) in the no fertilizer group is $\nu_1 = 12$ and the number of non-missing data in the fertilizer group is $\nu_2 =17$. $\nu_j$ for $j=1,2$ represents the actual sample size of observed responses in each group. Under the null hypothesis, the group assignments are exchangeable; thus, to obtain null distributions for our test statistics, we must permute group assignments.

The hypotheses are decomposed into four sub-hypotheses in which the alternatives are all one-sided. Thus, according to the previous theory, partial tests may be based on the following set of statistics: $$ \begin{aligned} T_1 &= \nu_2, \ T_2 &= \sum_i X_{2i}\gamma_2 - \sum_i X_{1i}\gamma_1, \ T_3 &= \sum_i Y_{2i}\gamma_2 - \sum_i Y_{1i}\gamma_1, \ T_4 &= \sum_i Y_{2i}^2\gamma_2 - \sum_i Y_{1i}^2\gamma_1, \ \end{aligned} $$

where $\gamma_1 = \sqrt{\frac{\gamma_2}{\gamma_1}}$ and $\gamma_2 = \sqrt{\frac{\gamma_1}{\gamma_2}}$. Instead of using the number of missing values in the germinated group as a test statistic, we may wish to consider the odds ratio of missingness. This statistic is equivalent to the product of the number non-missing, non-fertilized seeds and the number of missing, fertilized seeds. With regards to the second, third, and fourth test statistics, the sums do not include values for seeds which did not germinate. Recognizing this, we can rewrite the test statistics as

$$ \begin{aligned} T_1 &= (20-\nu_1)\nu_2, \ T_2 &= \sqrt{\nu_1\nu_2}\left(\bar{X}_2 - \bar{X}_1\right), \ T_3 &= \sqrt{\nu_1\nu_2}\left(\bar{Y}_2 - \bar{Y}_1\right), \ T_4 &= \sqrt{\nu_1\nu_2}\left(\bar{Y}^2_2 - \bar{Y}^2_1\right) \ \end{aligned} $$

Note that in this form, the second, third, and fourth test statistics can be obtained by a rescaling of the test statistic returned by two_sample. However, we cannot use two_sample to obtain the null distributions of these test statistics, as they depend on the prefactor $\sqrt{\nu_1\nu_2}$ which changes with each permutation.

mnar_two_sample <- function(dat) {
    tst_stat <- rep(0, 4)
    o <- table(dat$Fertilizer, dat$Germinated)
    tst_stat[1] <- prod(diag(o))
    nu <- o[, 2]
    prefactor <- sqrt(prod(nu))
    fert <- dat$Fertilizer == 1
    tst_stat[2] <- prefactor * (mean(dat$Weight[fert], na.rm = TRUE) - mean(dat$Weight[!fert], 
        na.rm = TRUE))
    tst_stat[3] <- prefactor * (mean(dat$Surface[fert], na.rm = TRUE) - mean(dat$Surface[!fert], 
        na.rm = TRUE))
    tst_stat[4] <- prefactor * (mean(dat$Surface2[fert], na.rm = TRUE) - mean(dat$Surface2[!fert], 
        na.rm = TRUE))
    return(tst_stat)
}

tst <- mnar_two_sample(germina)
distr <- replicate(B, {
    perm <- germina$Fertilizer[sample(1:n)]
    germina_star <- cbind(Fertilizer = perm, germina[, -1])
    mnar_two_sample(germina_star)
})

Using the test statistics and null distributions, we call t2p to obtain the p-values for the four partial tests. Each test has a one-sided alternative.

partial_pvalues <- sapply(1:4, function(i) t2p(tst[i], distr[i, ], alternative = "greater"))
names(partial_pvalues) <- c("O", "X", "Y", "Y2")
partial_pvalues

As seen before, we may perform a NPC of four partial tests, as in $T'' = \psi_1(T_1, T_2, T_3, T_4)$, giving each partial test the same weight, implicitly weighing the Y-component twice and the others once. Of course, if this is compatible with the inferential objective of the experimenter, then this is the solution he or she is looking for. If the inferential objective considers three variables $(O, X, Y)$ as having the same importance, as seems reasonable in the present problem, then the overall solution becomes $T''' = \psi_3\left[T_1, T_2, T_Y^{''}\right]$, where $T_Y^{''} = \psi_2(T_3, T_4)$.

We proceed by combining the tests for the first and second moments of $Y$ using Fisher's combining function, then applying the min-P procedure (Bonferroni-Holms step-down method, using Tippett's combining function at each step) to adjust the partial p-values for multiple tests.

y_tst <- fisher(partial_pvalues[3:4])
obs_pvalue_distr <- sapply(1:4, function(i) pvalue_distr(distr[i, ], alternative = "greater"))
y_tst_distr <- apply(obs_pvalue_distr[, 3:4], 1, fisher)
y_pvalue_distr <- pvalue_distr(y_tst_distr, alternative = "greater")

new_partial_pvalues <- c(partial_pvalues[1:2], t2p(y_tst, y_tst_distr, alternative = "greater"))
names(new_partial_pvalues) <- c("O", "X", "Y")
cat("Partial p-values", "\n")
new_partial_pvalues

adjusted_pvalues <- fwe_minp(new_partial_pvalues, cbind(t(distr[1:2, ]), y_tst_distr), 
    combine = "tippett")
names(adjusted_pvalues) <- c("O", "X", "Y")
cat("Adjusted Partial p-values:", "\n")
adjusted_pvalues

Finally, to obtain a global p-value, we use Fisher's combining function to carry out NPC.

global_tst <- fisher(new_partial_pvalues)
global_distr <- apply(cbind(obs_pvalue_distr[, 1:2], y_pvalue_distr), 1, fisher)
global_pvalue <- t2p(global_tst, global_distr, alternative = "greater")
cat("Global p-value:\n")
global_pvalue

res <- cbind(c("Germinate", "Weight", "Surface+Surface^2", "Global"), c(new_partial_pvalues, 
    global_pvalue))
colnames(res) <- c("Variable", "P-value")
kable(res, row.names = F)

From these results we can conclude that there is a significant effect of the fertilizer on the tree growth. The greater effect seems to be related to the area of the leaves, rather than their weight. The fertilizer does not seem to be significantly associated with the germination of the seeds.

statlab/permuter documentation built on May 30, 2019, 9:45 a.m.