You don’t need to any R code to do this week’s homework, so we will devote today’s lab to doing all the test’s we covered this week in R. We covered a lot of tests this week!

You should know how to get the Wilcoxon Rank Sum, Sign and Signed Rank test statistics by hand (for very small datasets) but in practice (in your project and outside this class) it is more useful for you to have R do them for you. There’s also an easier way to do the randomization test we saw in Chapter 1, than by doing all the random groupings yourself.

There are some commonalities in most of the functions to do the tests. They all take a special R object called a formula of the form Response ~ Groups. We read this as “Response is modelled by Groups”. R estimates a different mean Response in each level of the Groups variable, and tests if they are different. “Response” and “Groups” should be the columns in the appropriate data.frame that contain the response and the grouping variable. All the functions also take a data argument to specify the data.frame in which the Response and Groups columns are stored.

Remember the Creativity dataset,

library(Sleuth2)
head(case0101)

We wanted to test for a difference in mean Score for each group in Treatment. Our formula would be Score ~ Treatment and we would also specify the data.frame with data = case0101.

Your turn: Practice writing formulas

For each of these datasets (we have seen in class or homework), identify the response and grouping variable (for example by using head) and write the formula for the relationship we are interested in.

ex0112    # fish oil and blood pressure
case0401  # O-ring failure
ex0211    # Guinea pig lifetimes

The coin package

The coin package does a much better job than base R at the non-parametric tests (Wilcoxon Rank Sum, Rank Sign test and the randomization test). You will need to install it once:

install.packages("coin")

Then load it:

library(coin)

Most of the tests illustrated here are two-sided, all the test functions take an argument alternative = "greater" if you need a one-sided test. Unlike the t-tests, a one-sided p-value in a non-parametric test isn’t always half the two-sided p-value.

Two independent groups

Randomization test

Remember Lab 2 and Homework 2 where I had you do a randomization test step-by-step. Well coin can do it for you with the function oneway_test. This is equivalent to the Creativity example in the book (pg 1.8), calculating the p-value relative to 1000 random groupings:

oneway_test(Score ~ Treatment, data = case0101, 
   distribution = "approximate")

But now its really easy to do many more groupings, here is a p-value based on 10000:

oneway_test(Score ~ Treatment, data = case0101, 
    distribution = approximate(B = 10000))

Run the above line again, you should find the p-value doesn’t vary as much as when it is based on 1000 random groupings.

Homework 2 (fish oil and blood pressure) can now be done this way too (two-sided test):

oneway_test(BP ~ Diet, data = ex0112, distribution = approximate(B = 10000))

In fact doing an exact test for the blood pressure data isn’t too burdensome, and oneway_test will do it for you too:

oneway_test(BP ~ Diet, data = ex0112, distribution = "exact")

You can think of this p-value as being the one you would have got if you had time to run all possible regroupings of the fish oil data in two groups of size 7.

Wilcoxon Rank Sum

The coin package does a much better job of the Wilcoxon Rank Sum test. Using the O-ring example from class and the normal approximation (coin calls this the asymptotic distribution)

wilcox_test(Incidents ~ Launch, data = case0401, 
    distribution = "asymptotic")

This is exactly 2 times the one we got in class, we did a one-sided test, this is two sided. We can also approximate the p-value by doing random regroupings

wilcox_test(Incidents ~ Launch, data = case0401, 
  distribution = "approximate")

Or find the exact p-value like we did in class:

wilcox_test(Incidents ~ Launch, data = case0401, distribution = "exact")

You can get the Wilcoxon rank sum statistic by saving the test result and using the statistic function

wt <- wilcox_test(Incidents ~ Launch, data = case0401, 
  distribution = "exact")
statistic(wt, type = "linear")

Contrary to wilcox.test, this is the statistic we calculated by hand.

You can also get confidence intervals by adding conf.int = TRUE

wt <- wilcox_test(Incidents ~ Launch, data = case0401, 
    distribution = "exact", conf.int = TRUE)
wt

Levene’s Test

Levene’s test as we did it doesn’t seem to be implemented anywhere in R. Here’s a relatively simple and flexible way to do it, using the guinea pig data from class. We’ll see more about lm in Chapter 7.

ex0211$deviations.sq <- residuals(lm(Lifetime ~ Group, data = ex0211))^2
# uses residuals + lm to get deviations and store them in a new column in
# ex0211 called deviations.sq

t.test(deviations.sq ~ Group, data = ex0211, var.equal = TRUE)

Welch’s t-test

Easy! For the guinea pig data:

t.test(Lifetime ~ Group, data = ex0211)

Paired groups

We’ll see these two tests on Friday.

Sign test

I haven’t figured out to do this in coin yet :( so back to base R. For the Schizophrenia study:

diffs <- case0202$Unaffect - case0202$Affected
# how many differences are postive
sum(diffs > 0)
# how many aren't zero
sum(diffs != 0)
binom.test(14, 15)

# or all in one go, gives exact p-value (as opposed to Normal approximation used in Sleuth)
binom.test(sum(diffs > 0), sum(diffs != 0))

Signed Rank test

On the schizophrenia study, I’m using a one-sided test here to match the example in Sleuth.

wilcoxsign_test(Unaffect ~ Affected, data = case0202, 
  alternative = "greater", distribution = "exact")

Note that this is a different use of the ~ to the two group tests. We aren’t getting a mean value of Unaffect for each value of Affected.

Again we can get the statistic out if we want:

wst <- wilcoxsign_test(Unaffect ~ Affected, data = case0202, 
        alternative = "greater", 
        distribution = exact(), zero.method = "Wilcoxon")          
statistic(wst, "linear")

There’s a weird issue with the Schizoprenia case study where it has some erroneous level of precision, so we get a slightly different answer.