Here we demonstrate how to carry out a small simulation study of the robustness of the Type I error rate of a permutation test under two scenarios that violate exchangeability by having different shapes under the null hypothesis. One might encounter such a situation in practice in a study of growth in which unequal follow-up time is allotted to the groups. Our study is analogous to assessing the t-test assuming equal variances, when the variances are known to be unequal and/or the distributions non-normal under the null.

Scenario I: group 1 is distributed as a Normal random variable (section 1.10.5) with mean 0 and standard deviation 1, while group 2 is distributed as a Normal random variable with mean 0 and standard deviation 6. Here the distributions are both symmetric, but have markedly different variances.

Scenario II: group 1 is distributed as a Normal random variable with mean 6 and standard deviation 1, while group 2 is distributed as an exponential random variable (section 1.10.7) with mean 6 and standard deviation 6. Here the distributions have markedly different skewness, as well as variance.

We use an alpha level of 0.05, set the sample size in each group to 10, and generate 10,000 simulated data sets. Instead of the exact permutation test, we use the asymptotically equivalent Monte Carlo Hypothesis Test (Dwass, 1957), generating 1,000 samples from the permutation distribution for each simulation (thank heavens for fast computers!)

**R**

In R, we set up some constants and a vector for results that we will use. The group vector can be generated once, with the outcome generated inside the loop. The

`oneway_test()`function from the

`coin`library is used to run the permutation test, and return the p-value (the first element in the list provided by

`pvalue()`). Finally, the

`prop.test()`function (section 2.1.9) is used to calculate the 95% confidence interval for the level of the test.

For scenario I, the code is:

numsim <- 10000

numperm <- 1000

n <- 10

res <- numeric(numsim)

x <- c(rep(0, n), rep(1, n))

for (i in 1:numsim) {

y1 <- rnorm(n, 0, 1)

y2 <- rnorm(n, 0, 6)

y <- c(y1, y2)

res[i] <- pvalue(oneway_test(y ~ as.factor(x),

distribution=approximate(numperm)))[[1]]

}

The desired results can be generated as in section 2.1.9 by summarizing the res vector:

prop.test(sum(res<=0.05), numsim)

The result follows:

1-sample proportions test with continuity correction

data: sum(res <= 0.05) out of numsim, null probability 0.5

X-squared = 7574.221, df = 1, p-value < 2.2e-16

alternative hypothesis: true p is not equal to 0.5

95 percent confidence interval:

0.06009200 0.06984569

sample estimates:

p

0.0648

We see that the permutation test is somewhat anti-conservative given such an extreme difference in variances.

For scenario II, the code is the same with the exception of:

y1 <- rnorm(n, 6, 1)

y2 <- rexp(n, 1/6)

The yields the following results:

> prop.test(sum(res<=0.05), numsim)

1-sample proportions test with continuity correction

data: sum(res <= 0.05) out of numsim, null probability 0.5

X-squared = 5973.744, df = 1, p-value < 2.2e-16

alternative hypothesis: true p is not equal to 0.5

95 percent confidence interval:

0.1073820 0.1199172

sample estimates:

p

0.1135

The nominal Type I error rate of 0.05 is not preserved in this situation.

**SAS**

In SAS, we begin by making data for scenario 1, using the looping tools described in section 1.11.1.

data test;

do numsim = 1 to 10000;

do group = 0,1;

do i = 1 to 10;

if group = 0 then y = rand('NORMAL',0,1);

else y = rand('NORMAL',0,6);

output;

end;

end;

end;

run;

Next, we perform the permutation test for equal means as shown in section 2.4.3. The

`by`statement repeats the test for each simnum, while the

`ods`statements, as demonstrated throughout the examples in the book, prevent copious output and save relevant results. This was found using the

`ods trace on`statement, as discussed in section A.7.1.

ods select none;

ods output datascoresmc = kkout;

proc npar1way data = test;

by numsim;

class group;

var y;

exact scores = data / mc n = 1000;

run;

ods select all;

Next, we make a new data set from the saved output, saving just the lines we need, and generating a new variable which records whether the test was rejected or not. The names and values of the variables were found using a

`proc print`statement, which is not shown.

data props;

set kkout;

where (name1 = "MCP2_DATA");

reject = (nvalue1 lt .05);

run;

The estimated rejection probability in this setting is found as discussed in section 2.1.9:

proc freq data = props;

tables reject/binomial(level='1');

run;

The resulting output follows:

Cumulative Cumulative

reject Frequency Percent Frequency Percent

___________________________________________________________

0 9339 93.39 9339 93.39

1 661 6.61 10000 100.00

Proportion 0.0661

ASE 0.0025

95% Lower Conf Limit 0.0612

95% Upper Conf Limit 0.0710

Exact Conf Limits

95% Lower Conf Limit 0.0613

95% Upper Conf Limit 0.0711

ASE under H0 0.0050

Z -86.7800

One-sided Pr < Z <.0001

Two-sided Pr > |Z| <.0001

Sample Size = 10000

These results agree substantially with the results from R.

To produce scenario 2, we replace

if group = 0 then y = rand('NORMAL',0,1);

else y = rand('NORMAL',0,6);

with

if group = 0 then y = rand('NORMAL',6,1);

else y = rand('EXPONENTIAL') * 6;

with the following results:

Cumulative Cumulative

reject Frequency Percent Frequency Percent

___________________________________________________________

0 8877 88.77 8877 88.77

1 1123 11.23 10000 100.00

Proportion 0.1123

ASE 0.0032

95% Lower Conf Limit 0.1061

95% Upper Conf Limit 0.1185

Exact Conf Limits

95% Lower Conf Limit 0.1062

95% Upper Conf Limit 0.1187

ASE under H0 0.0050

Z -77.5400

One-sided Pr < Z <.0001

Two-sided Pr > |Z| <.0001

Sample Size = 10000

Again, these results are comfortably similar to those found in R and demonstrate that this extreme violation of the null results in a noticeable anti-conservative bias.

However, note that with additional observations, the anti-conservative bias diminishes considerably. With 100 observations per group, (replacing

`do i = 1 to 10`with

`do i = 1 to 100`) the following results are generated:

Cumulative Cumulative

reject Frequency Percent Frequency Percent

___________________________________________________________

0 9452 94.52 9452 94.52

1 548 5.48 10000 100.00

Proportion 0.0548

ASE 0.0023

95% Lower Conf Limit 0.0503

95% Upper Conf Limit 0.0593

Exact Conf Limits

95% Lower Conf Limit 0.0504

95% Upper Conf Limit 0.0594

ASE under H0 0.0050

Z -89.0400

One-sided Pr < Z <.0001

Two-sided Pr > |Z| <.0001

Sample Size = 10000

revealing a negligible bias.

Our conclusion: while the permutation test is somewhat robust to violations of exchangeability, when shapes and variances differ radically the results may be anti-conservatively biased.

It bears noting here that if there were no reason to expect different distributional shapes under the null, but we just happened to observed the large differences seen here, we would interpret these results as demonstrating the poor power of a permutation test based on the means to detect differences in shape alone.

Dwass M. Modified randomization tests for nonparametric hypotheses.

*Annals of Mathematical Statistics*, 28:181-187, 1957.

## No comments:

Post a Comment