Robert I. Kabacoff - R in action

Next, consider the Wilcoxon–Mann–Whitney U test. In chapter 7, we examined the difference in the probability of imprisonment in Southern versus non-Southern US states using the wilcox.test() function. Using an exact Wilcoxon rank sum test, we’d get

>library(MASS)

>UScrime <- transform(UScrime, So = factor(So))

>wilcox_test(Prob ~ So, data=UScrime, distribution="exact")

Exact Wilcoxon Mann-Whitney Rank Sum Test

data: Prob by So (0, 1)

Z = -3.7, p-value = 8.488e-05

alternative hypothesis: true mu is not equal to 0

suggesting that incarceration is more likely in Southern states. Note that in the previous code, the numeric variable So was transformed into a factor. This is because the coin package requires that all categorical variables be coded as factors. Additionally, the astute reader may have noted that these results agree exactly with the results of the wilcox.test() in chapter 7. This is because the wilcox.test() also computes an exact distribution by default.

Finally, consider a k-sample test. In chapter 9, we used a one-way ANOVA to evaluate the impact of five drug regimens on cholesterol reduction in a sample of 50 patients. An approximate k-sample permutation test can be performed instead, using this code:

>library(multcomp)

>set.seed(1234)

>oneway_test(response~trt, data=cholesterol, distribution=approximate(B=9999))

Approximative K-Sample Permutation Test

data: response by

trt (1time, 2times, 4times, drugD, drugE) maxT = 4.7623, p-value < 2.2e-16

Here, the reference distribution is based on 9,999 permutations of the data. The random number seed was set so that your results would be the same as mine. There’s clearly a difference in response among patients in the various groups.

12.2.2Independence in contingency tables

We can use permutation tests to assess the independence of two categorical variables using either the chisq_test() or the cmh_test() function. The latter function is used when the data is stratified on a third categorical variable. If both variables are ordinal, we can use the lbl_test() function to test for a linear trend.

In chapter 7, we applied a chi-square test to assess the relationship between Arthritis treatment and improvement. Treatment had two levels (Placebo, Treated), and Improved had three levels (None, Some, Marked). The Improved variable was encoded as an ordered factor.

If you want to perform a permutation version of the chi-square test, you could use the following code:

>library(coin)

>library(vcd)

>Arthritis <- transform(Arthritis, Improved=as.factor(as.numeric(Improved)))

>set.seed(1234)

>chisq_test(Treatment~Improved, data=Arthritis,

distribution=approximate(B=9999))

Approximative Pearson’s Chi-Squared Test

data: Treatment by Improved (1, 2, 3) chi-squared = 13.055, p-value = 0.0018

This gives you an approximate chi-square test based on 9,999 replications. You might ask why you transformed the variable Improved from an ordered factor to a categorical factor. (Good question!) If you’d left it an ordered factor, coin() would have generated a linear x linear trend test instead of a chi-square test. Although a trend test would be a good choice in this situation, keeping it a chi-square test allows you to compare the results with those reported in chapter 7.

12.2.3Independence between numeric variables

The spearman_test() function provides a permutation test of the independence of two numeric variables. In chapter 7, we examined the correlation between illiteracy rates and murder rates for US states. You can test the association via permutation, using the following code:

>states <- as.data.frame(state.x77)

>set.seed(1234)

>spearman_test(Illiteracy~Murder, data=states,

distribution=approximate(B=9999))

Approximative Spearman Correlation Test

data: Illiteracy by Murder

Z = 4.7065, p-value < 2.2e-16

alternative hypothesis: true mu is not equal to 0

Based on an approximate permutation test with 9,999 replications, the hypothesis of independence can be rejected. Note that state.x77 is a matrix. It had to be converted into a data frame for use in the coin package.

12.2.4Dependent two-sample and k-sample tests

Dependent sample tests are used when observations in different groups have been matched, or when repeated measures are used. For permutation tests with two paired groups, the wilcoxsign_test() function can be used. For more than two groups, use the friedman_test() function.

In chapter 7, we compared the unemployment rate for urban males age 14–24 (U1) with urban males age 35–39 (U2). Because the two variables are reported for each of the 50 US states, you have a two-dependent groups design (state is the matching variable). We can use an exact Wilcoxon Signed Rank Test to see if unemployment rates for the two age groups are equal:

>library(coin)

>library(MASS)

>wilcoxsign_test(U1~U2, data=UScrime, distribution="exact")

Exact Wilcoxon-Signed-Rank Test

data: y by x (neg, pos) stratified by block

Z = 5.9691, p-value = 1.421e-14

alternative hypothesis: true mu is not equal to 0

Based on the results, you’d conclude that the unemployment rates differ.

12.2.5Going further

The coin package provides a general framework for testing that one group of variables is independent of a second group of variables (with optional stratification on a blocking variable) against arbitrary alternatives, via approximate permutation tests. In particular, the independence_test() function allows the user to approach most traditional tests from a permutation perspective, and to create new and novel statistical tests for situations not covered by traditional methods. This flexibility comes at a price: a high level of statistical knowledge is required to use the function appropriately. See the vignettes that accompany the package (accessed via vignette("coin")) for further details.

In the next section, you’ll learn about the lmPerm package. This package provides a permutation approach to linear models, including regression and analysis of variance.

12.3 Permutation tests with the lmPerm package

The lmPerm package provides support for a permutation approach to linear models. In particular, the lmp() and aovp() functions are the lm() and aov() functions modified to perform permutation tests rather than normal theory tests.

The parameters within the lmp() and aovp() functions are similar to those in the lm() and aov()functions, with the addition of a perm= parameter. The perm= option can take on the values "Exact", "Prob", or "SPR". Exact produces an exact test, based on all possible permutations. Prob samples from all possible permutations. Sampling continues until the estimated standard deviation falls below 0.1 of the estimated p-value. The stopping rule is controlled by an optional Ca parameter. Finally, SPR uses a sequential probability ratio test to decide when to stop sampling. Note that if the number of observations is greater than 10, perm="Exact" will automatically default to perm="Prob"; exact tests are only available for small problems.

To see how this works, we’ll apply a permutation approach to simple regression, polynomial regression, multiple regression, one-way analysis of variance, one-way analysis of covariance, and a two-way factorial design.

12.3.1Simple and polynomial regression

In chapter 8, we used linear regression to study the relationship between weight and height for a group of 15 women. Using lmp() instead of lm() generates the permutation test results shown in the following listing.

Listing 12.2 Permutation tests for simple linear regression

>library(lmPerm)

>set.seed(1234)

>fit <- lmp(weight~height, data=women, perm="Prob") [1] "Settings: unique SS : numeric variables centered"

>summary(fit)

Call:

lmp(formula = weight ~ height, data = women, perm = "Prob")

Residual standard error: 1.5 on 13 degrees of freedom

Multiple R-Squared: 0.991, Adjusted R-squared: 0.99

F-statistic: 1.43e+03 on 1 and 13 DF, p-value: 1.09e-14

To fit a quadratic equation, you could use the code in this next listing.

Listing 12.3 Permutation tests for polynomial regression

>library(lmPerm)

>set.seed(1234)

>fit <- lmp(weight~height + I(height^2), data=women, perm="Prob") [1] "Settings: unique SS : numeric variables centered"

>summary(fit)

Call:

lmp(formula = weight ~ height + I(height^2), data = women, perm = "Prob")

Residuals:

Min 1Q Median 3Q Max -0.5094 -0.2961 -0.0094 0.2862 0.5971

Coefficients:

Estimate Iter Pr(Prob)

Residual standard error: 0.38 on 12 degrees of freedom

Multiple R-Squared: 0.999, Adjusted R-squared: 0.999

F-statistic: 1.14e+04 on 2 and 12 DF, p-value: <2e-16

As you can see, it’s a simple matter to test these regressions using permutation tests and requires little change in the underlying code. The output is also similar to that produced by the lm() function. Note that an Iter column is added indicating how many iterations were required to reach the stopping rule.

12.3.2Multiple regression

In chapter 8, multiple regression was used to predict the murder rate from population, illiteracy, income, and frost for 50 US states. Applying the lmp() function to this problem, results in the following output.

Listing 12.4 Permutation tests for multiple regression

>library(lmPerm)

>set.seed(1234)

>states <- as.data.frame(state.x77)

>fit <- lmp(Murder~Population + Illiteracy+Income+Frost,

data=states, perm="Prob")

[1] "Settings: unique SS : numeric variables centered" > summary(fit)

Call:

lmp(formula = Murder ~ Population + Illiteracy + Income + Frost, data = states, perm = "Prob")

Residual standard error: 2.535 on 45 degrees of freedom

Multiple R-Squared: 0.567, Adjusted R-squared: 0.5285

F-statistic: 14.73 on 4 and 45 DF, p-value: 9.133e-08

Looking back to chapter 8, both Population and Illiteracy are significant (p < 0.05) when normal theory is used. Based on the permutation tests, the Population

variable is no longer significant. When the two approaches don’t agree, you should look at your data more carefully. It may be that the assumption of normality is untenable or that outliers are present.

12.3.3One-way ANOVA and ANCOVA

Each of the analysis of variance designs discussed in chapter 9 can be performed via permutation tests. First, let’s look at the one-way ANOVA problem considered in sections 9.1 on the impact of treatment regimens on cholesterol reduction. The code and results are given in the next listing.

Listing 12.5 Permutation test for One-Way ANOVA

>library(lmPerm)

>library(multcomp)

>set.seed(1234)

>fit <- aovp(response~trt, data=cholesterol, perm="Prob") [1] "Settings: unique SS "

>summary(fit)

The results suggest that the treatment effects are not all equal.

This second example in this section applies a permutation test to a one-way analysis of covariance. The problem is from chapter 9, where you investigated the impact of four drug doses on the litter weights of rats, controlling for gestation times. The next listing shows the permutation test and results.

Listing 12.6 Permutation test for one-way ANCOVA

>library(lmPerm)

>set.seed(1234)

>fit <- aovp(weight ~ gesttime + dose, data=litter, perm="Prob") [1] "Settings: unique SS : numeric variables centered"

>summary(fit)

Based on the p-values, the four drug doses do not equally impact litter weights, controlling for gestation time.

12.3.4Two-way ANOVA

We’ll end this section by applying permutation tests to a factorial design. In chapter 9, we examined the impact of vitamin C on the tooth growth in guinea pigs. The two manipulated factors were dose (three levels) and delivery method (two levels). Ten guinea pigs were placed in each treatment combination, resulting in a balanced 3 x 2 factorial design. The permutation tests are provided in the next listing.

Listing 12.7 Permutation test for two-way ANOVA

>library(lmPerm)

>set.seed(1234)

>fit <- aovp(len~supp*dose, data=ToothGrowth, perm="Prob") [1] "Settings: unique SS : numeric variables centered"

>summary(fit)

At the .05 level of significance, all three effects are statistically different from zero. At the .01 level, only the main effects are significant.

It’s important to note that when aovp() is applied to ANOVA designs, it defaults to unique sums of squares (also called SAS Type III sums of squares). Each effect is adjusted for every other effect. The default for parametric ANOVA designs in R is sequential sums of squares (SAS Type I sums of squares). Each effect is adjusted for those that appear earlier in the model. For balanced designs, the two approaches will agree, but for unbalanced designs with unequal numbers of observations per cell, they won’t. The greater the imbalance, the greater the disagreement. If desired, specifying seqs=TRUE in the aovp() function will produce sequential sums of squares. For more on Type I and Type III sums of squares, see chapter 9,section 9.2.

12.4 Additional comments on permutation tests

R offers other permutation packages besides coin and lmPerm. The perm package provides some of the same functionality provided by the coin package and can act as an independent validation of that package. The corrperm package provides permutation tests of correlations with repeated measures. The logregperm package offers a permutation test for logistic regression. Perhaps most importantly, the glmperm package extends permutation tests to generalized linear models. Generalized linear models are described in the next chapter.

Permutation tests provide a powerful alternative to tests that rely on a knowledge of the underlying sampling distribution. In each of the permutation tests described, we were able to test statistical hypotheses without recourse to the normal, t, F, or chisquare distributions.

approach is advantageous. What if you wanted confidence intervals for the sample median, or the difference between two sample medians? There are no simple normal theory formulas here, and bootstrapping is the approach of choice. If the underlying distributions are unknown, if outliers are a problem, if sample sizes are small, or if parametric approaches don’t exist, bootstrapping can often provide a useful method of generating confidence intervals and testing hypotheses.

12.6 Bootstrapping with the boot package

The boot package provides extensive facilities for bootstrapping and related resampling methods. You can bootstrap a single statistic (for example, a median), or a vector of statistics (for example, a set of regression coefficients). Be sure to download and install the boot package before first use:

install.packages("boot")

The bootstrapping process will seem complicated, but once you review the examples it should make sense.

In general, bootstrapping involves three main steps:

1Write a function that returns the statistic or statistics of interest. If there is a single statistic (for example, a median), the function should return a number. If there is a set of statistics (for example, a set of regression coefficients), the function should return a vector.

2Process this function through the boot() function in order to generate R bootstrap replications of the statistic(s).

3Use the boot.ci() function to obtain confidence intervals for the statistic(s) generated in step 2.

Now to the specifics.

The main bootstrapping function is boot(). The boot() function has the format

bootobject <- boot(data=, statistic=, R=, ...)

The parameters are described in table 12.3.

Table 12.3 Parameters of the boot() function

The boot() function calls the statistic function R times. Each time, it generates a set of random indices, with replacement, from the integers 1:nrow(data). These indices are used within the statistic function to select a sample. The statistics are calculated on the sample and the results are accumulated in the bootobject. The bootobject structure is described in table 12.4.

You can access these elements as bootobject$t0 and bootobject$t.

Once you generate the bootstrap samples, you can use print() and plot() to examine the results. If the results look reasonable, you can use the boot.ci() function to obtain confidence intervals for the statistic(s). The format is

boot.ci(bootobject, conf=, type= )

The parameters are given in table 12.5.

Table 12.5 Parameters of the boot.ci() function

The type parameter specifies the method for obtaining the confidence limits. The perc method (percentile) was demonstrated in the sample mean example. The bca provides an interval that makes simple adjustments for bias. I find bca preferable in most circumstances. See Mooney and Duval (1993) for an introduction to these methods.

In the remaining sections, we’ll look at bootstrapping a single statistic and a vector of statistics.

12.6.1Bootstrapping a single statistic

The mtcars dataset contains information on 32 automobiles reported in the 1974 Motor Trend magazine. Suppose you’re using multiple regression to predict miles per gallon from a car’s weight (lb/1,000) and engine displacement (cu. in.). In addition to the standard regression statistics, you’d like to obtain a 95 percent confidence interval for the R-squared value (the percent of variance in the response variable explained by the predictors). The confidence interval can be obtained using nonparametric bootstrapping.