Practical B

Exercises

Again, just like last time it is wise to start with fixing the random seed.

set.seed(123)

Generate two random samples of 10 numbers from a normal distribution with the below specifications. Test the null hypothesis that the population mean is 0.

\(\mu = 0\) and \(\sigma = 2\)
\(\mu = 1.5\) and \(\sigma = 2\)

x <- rnorm(10, mean = 0, sd = 2)
t.test(x)

## 
##  One Sample t-test
## 
## data:  x
## t = 0.24742, df = 9, p-value = 0.8101
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -1.215341  1.513843
## sample estimates:
## mean of x 
## 0.1492513

x <- rnorm(10, 1.5, 2)
t.test(x)

## 
##  One Sample t-test
## 
## data:  x
## t = 2.9202, df = 9, p-value = 0.01703
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.432058 3.402430
## sample estimates:
## mean of x 
##  1.917244

Write a function that generates a random sample of n numbers from a normal distribution with a user defined mean (i.e. a mean that you can choose when running the function) and standard deviation 1, and returns the p.value for the test that the mean is 0.

p.value.t <- function (n, mu) {
  x <- rnorm(n, mu, 1)
  t.test(x)$p.value
}

p.value.t(n = 30, mu = 3)

## [1] 2.912617e-17

Use the function of Exercise 3 to generate 50 \(p\)-values with \(n=10,\mu=0\), and make a qqplot to compare distribution of the \(p\)-values with a uniform \([0,1]\) variable.

y <- numeric(50)
for (i in 1:50) {
  y[i] <- p.value.t(n = 10, mu = 0)
}

qqplot(x=qunif(ppoints(50)), y)

The p-values follow a uniform distribution.

In a study that examined the use of acupuncture to treat migraine headaches, consenting patients on a waiting list for treatment for migraine were randomly assigned in a 2:1:1 ratio to acupuncture treatment, a “sham” acupuncture treatment in which needles were inserted at non-acupuncture points, and waiting-list patients whose only treatment was self-administered (Linde et al., 2005). The “sham” acupuncture treatment was described to trial participants as an acupuncture treatment that did not follow the principles of Chinese medicine.

What is the conclusion when the outcome is classified according to numbers of patients who experienced a greater than 50% reduction in headaches over a four-week period, relative to a pre-randomization baseline?

Use the following data

data <- matrix(c(74, 71, 43, 38, 11, 65), nrow = 2, ncol = 3)
colnames(data) <- c("Acupuncture", "Sham", "Waiting list")
rownames(data) <- c("> 50% reduction", "< 50% reduction")
data

##                 Acupuncture Sham Waiting list
## > 50% reduction          74   43           11
## < 50% reduction          71   38           65

We start with calculating the \(X^2\)-test:

X2test <- 
  data %>%
  chisq.test()

X2test

## 
##  Pearson's Chi-squared test
## 
## data:  .
## X-squared = 32.486, df = 2, p-value = 8.825e-08

which is extremely significant. We can then calculate the expected cell frequencies

X2test$expected

##                 Acupuncture     Sham Waiting list
## > 50% reduction    61.45695 34.33113     32.21192
## < 50% reduction    83.54305 46.66887     43.78808

and the raw residual

X2test$observed - X2test$expected

##                 Acupuncture      Sham Waiting list
## > 50% reduction    12.54305  8.668874    -21.21192
## < 50% reduction   -12.54305 -8.668874     21.21192

as well as the Pearson residual

X2test$residuals

##                 Acupuncture      Sham Waiting list
## > 50% reduction    1.599991  1.479513    -3.737418
## < 50% reduction   -1.372296 -1.268963     3.205546

to infer the difference in observed and expected cell frequencies. Patients on the waiting list experience > 50% reduction much less than we would expect under independence of treatment and outcome.

Patients who received the acupuncture and sham acupuncture treatments were asked to guess their treatment at the end of their trial. What would you conclude from this data?

data <- matrix(c(82, 17, 30, 30, 26, 16), nrow = 3, ncol = 2)
colnames(data) <- c("Acupuncture", "Sham")
rownames(data) <- c("Chinese", "Other", "Don't know")
data

##            Acupuncture Sham
## Chinese             82   30
## Other               17   26
## Don't know          30   16

We again start with calculating the \(X^2\)-test:

X2test <- 
  data %>%
  chisq.test()

X2test

## 
##  Pearson's Chi-squared test
## 
## data:  .
## X-squared = 15.358, df = 2, p-value = 0.0004624

which is very significant. We can then calculate the expected cell frequencies

X2test$expected

##            Acupuncture     Sham
## Chinese       71.88060 40.11940
## Other         27.59701 15.40299
## Don't know    29.52239 16.47761

and the raw residual

X2test$observed - X2test$expected

##            Acupuncture        Sham
## Chinese     10.1194030 -10.1194030
## Other      -10.5970149  10.5970149
## Don't know   0.4776119  -0.4776119

as well as the Pearson residual

X2test$residuals

##            Acupuncture       Sham
## Chinese     1.19357319 -1.5976353
## Other      -2.01721641  2.7001078
## Don't know  0.08790214 -0.1176598

We find that people who are receiving true Acupuncture are more inclined to believe that they receive Chinese acupuncture than we would expect under independence, while people wo received Sham acupuncture are more inclined to believe that they receive Other type of acupuncture. Don't know is more or less similarly distributed over the observed and expected frequencies.

Write a function that chooses automatically whether to do the chisq.test() or the fisher.test(). Create the function such that it:

takes two vectors as input
creates a table from the two vectors
checks if there is any expected cell frequency that is smaller than 5
and that then performs and prints the results of the appropriate test.

contingencyTest <- function(x, y) {
  
  # Make a table out of the variables.
  tab <- table(x, y)
  
  #expected frequencies
  exp.freq <- (colSums(tab) %*% t(rowSums(tab))) / sum(tab)
  
  # Choose the correct test. 
  if (any(exp.freq < 5)) {
    results <- fisher.test(x, y)
    return(results)
  } else {
    results <- chisq.test(x, y)
    return(results)
  }
}

A more efficient function that tests whether the expected frequencies are smaller than a threshold \(t\) is:

contingencyTest2 <- function(x, y, t = 5){
  test <- suppressWarnings(chisq.test(x, y))
  if (any(test$expected < t)){
    test <- fisher.test(x, y)
  }
  return(test)
}

The suppressWarnings() function is used to surpress printing of the warnings that function chi-square puts out to the console when too low expected cell-frequencies are encountered. Since we decide to use Fisher’s exact test in those situations anyway, it is redundant to print the message.

Test the function with the dataset bacteria (from MASS) by testing independence between compliance (hilo) and the presence or absence of disease (y).

bacteria %$%
  contingencyTest(ap, hilo)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  x and y
## X-squared = 2.9352, df = 1, p-value = 0.08667

bacteria %$%
  contingencyTest2(ap, hilo)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  x and y
## X-squared = 2.9352, df = 1, p-value = 0.08667

Does your function work differently if we only put in the first 25 rows of the bacteria dataset?

bacteria[1:25, ] %$%
  contingencyTest(ap, hilo)

## 
##  Fisher's Exact Test for Count Data
## 
## data:  x and y
## p-value = 1
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##   0.1896616 10.6212859
## sample estimates:
## odds ratio 
##   1.408082

bacteria[1:25, ] %$%
  contingencyTest2(ap, hilo)

## 
##  Fisher's Exact Test for Count Data
## 
## data:  x and y
## p-value = 1
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##   0.1896616 10.6212859
## sample estimates:
## odds ratio 
##   1.408082

Yes, it performs differently: expected cell frequencies are smaller than 5 for this subset.

End of practical.

Practical B

Jolien Cremers (adapted from Gerko Vink)

Statistical Programming with R

Exercises