# rcode19

## Computer-intensive methods: R code for Chapter 19 examples

Click on a function argument for a short description of its meaning. The variable names are plucked from the examples further below.

Calculate a P-value for a χ2 goodness-of-fit test using simulation:
`chisq.test(, , )`

Calculate a bootstrap standard error:
`boot(, , )`

Other new methods:
Bootstrap 95% confidence interval.

### Example 19.1 Two-digit “psychic”

Hypothesis testing using simulation to determine whether all two-digit numbers occur with equal probability when people choose two-digit numbers haphazardly.

Read and examine data. Each row is the two-digit number chosen by a different individual.

```haphazard <- read.csv(url("http://www.zoology.ubc.ca/~schluter/WhitlockSchluter/wp-content/data/chapter19/chap19e1TwoDigitNumbers.csv"))
```

Plot the frequency distribution of two-digit numbers chosen by volunteers.

```hist(haphazard\$numberSelected, right = FALSE, breaks = seq(10,100, by = 1),
las = 1, col = "firebrick", xlab = "Number thought of by volunteer",
ylim=c(0,50))
```

Frequency table of the two-digit numbers chosen by volunteers. Notice that many numbers between 0 and 99 are not represented in the data. To make sure that numbers not chosen are nevertheless included in the calculations, we must convert the variable to a factor and specify all possible categories.

```table(haphazard\$numberSelected)
factorData <- factor(haphazard\$numberSelected, levels = 10:99)
observedFreq <- table(factorData)
observedFreq
```

Calculate the observed value of the test statistic. Here, we are using the χ2 goodness-of-fit statistic. We are lazy and use the `chisq.test` function to do the calculations. R will give you a warning because of the low expected frequencies. This is exactly why we need to use simulation instead of the χ2 goodness-of-fit-test.

```chiData <- chisq.test(observedFreq)\$statistic
chiData
```

We'll need to know the sample size (number of rows, in this case) to run the simulation.

```n <- nrow(haphazard)
n
```

Simulate a single random sample of n numbers from 10 to 99 (with replacement), where n is the number of individuals in the sample.

```randomTwoDigit <- sample(10:99, size = n, replace = TRUE)
```

Calculate the test statistic for this single simulated sample2). Remember to convert to a factor to make sure that all categories are included in the calculations. Under the null hypothesis, the expected frequency is the same for all numbers.

```randomTwoDigit <- factor(randomTwoDigit, levels = 10:99)
observedFreq <- table(randomTwoDigit)
expectedFreq <- n/length(10:99)
chiSim <- sum( (observedFreq - expectedFreq)^2 / expectedFreq )
chiSim
```

Repeat this process many times. The following is a loop repeated `nSim` times. In each iteration `i`, a random sample of numbers is simulated, a frequency table is calculated, and the χ2 statistic is computed. The statistic is saved in the `i`th element of the vector `results`.

```nSim <- 10000
results <- vector()
for(i in 1:nSim){
randomTwoDigit <- sample(10:99, size = n, replace = TRUE)
randomTwoDigit <- factor(randomTwoDigit, levels = 10:99)
simFreq <- table(randomTwoDigit)
expectedFreq <- n/length(10:99)
results[i] <- sum( (simFreq - expectedFreq)^2 / expectedFreq )
}
```

Plot the frequency distribution of χ2 values from the simulation. This is the simulated null distribution for the test statistic.

```hist(results, right = FALSE, breaks = 50)
```

The P-value is the fraction of similated chisquare values equalling or exceededing the observed value of the test statistic in the data.

```P <- sum(results >= chiData)
P
```

R also has a built-in method to simulate the null distribution and calculate the P-value. Using it simply involves providing a couple of arguments to the `chisq.test` function. `B` is the number of iterations desired. If you don't provide the probabilities of each outcome under the null hypothesis as an argument, R assumes that all the categories are equiprobable under the nuill hypothesis.

```factorData <- factor(haphazard\$numberSelected, levels = 10:99)
observedFreq <- table(factorData)
chisq.test(observedFreq, simulate.p.value = TRUE, B = 10000)
```

### Example 19.2 Chimp language centers

Bootstrap estimation of the median asymmetry score for Brodmann’s area 44 in chimpanzees.

```chimp <- read.csv(url("http://www.zoology.ubc.ca/~schluter/WhitlockSchluter/wp-content/data/chapter19/chap19e2ChimpBrains.csv"))
```

Histogram of asymmetry scores of the 20 chimps.

```hist(chimp\$asymmetryScore, breaks = seq(-.5, 1.25, by = 0.25),
right = FALSE, col = "firebrick", las = 1)
```

Estimate the median asymmetry score from the data.

```median(chimp\$asymmetryScore)
```

Obtain a single bootstrap replicate. This involves sampling (with replacement) from the actual data. Calculate the median of this bootstrap replicate.

```bootSample <- sample(chimp\$asymmetryScore, size = nrow(chimp), replace = TRUE)
median(bootSample)
```

Repeat this process many times. The following is a loop repeated `B` times. In each iteration, the bootstrap replicate estimate (median) is recalculated and saved in the results vector `bootMedian`.

```B <- 10000
bootMedian <- vector()
for(i in 1:B){
bootSample <- sample(chimp\$asymmetryScore, size = nrow(chimp), replace = TRUE)
bootMedian[i] <- median(bootSample)
}
```

Draw a histogram of the bootstrap replicate estimates for median asymmetry. Note: your results won't be the identical to the one in Figure 19.2-2, because 10,000 random samples is not large enough to obtain the sampling distribution with extreme accuracy. Increase `B` in the above loop for greater accuracy.

```hist(bootMedian, breaks = 100, right = FALSE, col = "firebrick", las = 1)
```

Calculate the mean of the bootstrap replicate estimates.

```mean(bootMedian)
```

The bootstrap standard error is the standard deviation of the bootstrap replicate estimates.

```sd(bootMedian)
```

Use the percentiles of the bootstrap replicate estimates to obtain a bootstrap 95% confidence interval for the population median.

```quantile(bootMedian, probs = c(0.025, 0.975))
```

Use the `boot` package for bootstrap estimation.

```library(boot)
```

To begin, we need to define a new R function to tell the `boot` package what statistic we want to estimate, which is the median in this example. Let's call the new function `boot.median`. The syntax required by the package is not too complex, but our function must include a counter as an argument (as usual, we call it `i`).

```boot.median <- function(x, i){ median(x[i]) }
```

Use the boot function of the `boot` package to obtain many bootstrap replicate estimates and calculate the bootstrap standard error.

```bootResults <- boot(chimp\$asymmetryScore, boot.median, R = 10000)
bootResults
```

Draw a histogram of the bootstrap replicate estimates for median asymmetry.

```hist(bootResults\$t, breaks = 100, right = FALSE, col = "firebrick", las = 1)
```

Get the bootstrap 95% confidence interval for the population median. We show two methods. The percentile method is already familiar, but `boot` also provides an improved method, called "bias corrected and accelerated".

```boot.ci(bootResults, type = "perc")
boot.ci(bootResults, type = "bca")
```

### Figure 19.2-3 Sexual cannibalism in crickets revisited

Bootstrap estimation of the difference between medians of two groups, using times to mating (in hours) of female sagebrush crickets that were either starved or fed (data from Example 13.5). We show the method using R's `boot` package.

```cannibalism <- read.csv(url("http://www.zoology.ubc.ca/~schluter/WhitlockSchluter/wp-content/data/chapter13/chap13e5SagebrushCrickets.csv"))
```

Calculate the observed difference between the sample medians of the two groups.

```twoMedians <- tapply(cannibalism\$timeToMating, cannibalism\$feedingStatus, median)
diffMedian <- twoMedians -  twoMedians
diffMedian
```

Define a new function to calculate the difference between the medians of the bootstrap replicate estimates. Let's call it `boot.diffMedian`.

```boot.diffMedian <- function(x, i){
twoMedians <- tapply(x\$timeToMating[i], x\$feedingStatus[i], median)
diffMedian <- twoMedians -  twoMedians
}
```

Load the `boot` package, obtain the many bootstrap replicate estimates and get the bootstrap standard error.

```library(boot)
bootResults <- boot(cannibalism, boot.diffMedian, R = 10000)
bootResults
```

Histogram of the bootstrap replicate estimates of the difference between two medians (Figure 19.2-3).

```hist(bootResults\$t, breaks = 20, right = FALSE, col = "firebrick", las = 1)
```

Obtain the bootstrap 95% confidence interval for the difference between the population medians, using the percentile and the "bca" method.

```boot.ci(bootResults, type = "perc")
boot.ci(bootResults, type = "bca")
```
A table of observed frequencies.
Instructs R to use simulation instead of the theoretical χ2 distribution to approximate the P-value.
Number of simulated random samples to use.
A sample of data for a single variable.
Name of a function created by the user to calculate the statistic of interest.
Number of bootstrap replicates to use.