Chapter 28 Exercise Solutions

These are the solutions to the exercises used in the course.

The solutions use these packages (others are loaded where they are needed):

library(dplyr)
library(tidyr)
library(ggplot2)

28.1 Californian bird diversity

1. Import the data as an object called birds, then check the columns with str().

birds <- read.csv("CourseData/suburbanBirds.csv")
str(birds)

## 'data.frame':    17 obs. of  5 variables:
##  $ Name        : chr  "Alamotos" "Ramona" "Verona" "Valle Vista" ...
##  $ Year        : int  1946 1946 1947 1950 1955 1956 1959 1962 1962 1966 ...
##  $ HabitatIndex: num  10 9.5 9.5 9.5 11 9 8 8 10 5 ...
##  $ nIndividuals: int  48 30 38 42 44 27 31 26 32 26 ...
##  $ nSpecies    : int  12 13 15 11 13 14 13 12 14 8 ...

There are 17 suburbs and the five columns Name, Year, HabitatIndex, nIndividuals and nSpecies. (HabitatIndex has one missing value, NA.)

2. Calculate the mean, minimum, and maximum of nSpecies.

mean(birds$nSpecies)

## [1] 9.647059

min(birds$nSpecies)

## [1] 3

max(birds$nSpecies)

## [1] 15

There is more than one way to do this. For example, range(birds$nSpecies) gives the minimum and maximum together, and summary(birds$nSpecies) gives all three (and more).

3. Create a new column with birds$Age <- 1975 - birds$Year and check that it looks sensible.

birds$Age <- 1975 - birds$Year
range(birds$Age)

## [1]  1 29

The suburbs were between 1 and 29 years old when they were surveyed.

4. Plot Age against nSpecies.

plot(birds$Age, birds$nSpecies)

5. Add a smoother.

plot(birds$Age, birds$nSpecies)
lines(lowess(birds$Age, birds$nSpecies))

Species richness increases with the age of the suburb: the newest suburbs have only 3 to 6 species, while those built 20 or more years before the survey have 11 to 15. The smoother suggests that richness rises steeply over the first 15 to 20 years and then levels off. A possible explanation is that the gardens, parks and trees that come with development take years to mature, and as they do they provide habitat for more species.

28.2 Wrangling the Amniote Life History Database

1. When you have imported the data, use dim to check the dimensions of the whole data frame (you should see that there are 36 columns and 21322 rows). Use names to look at the names of all columns in the data in amniote.

amniote <- read.csv("CourseData/Amniote_Database_Aug_2015.csv",
  na.strings = "-999"
)

dim(amniote)

## [1] 21322    36

names(amniote)

##  [1] "class"                                
##  [2] "order"                                
##  [3] "family"                               
##  [4] "genus"                                
##  [5] "species"                              
##  [6] "subspecies"                           
##  [7] "common_name"                          
##  [8] "female_maturity_d"                    
##  [9] "litter_or_clutch_size_n"              
## [10] "litters_or_clutches_per_y"            
## [11] "adult_body_mass_g"                    
## [12] "maximum_longevity_y"                  
## [13] "gestation_d"                          
## [14] "weaning_d"                            
## [15] "birth_or_hatching_weight_g"           
## [16] "weaning_weight_g"                     
## [17] "egg_mass_g"                           
## [18] "incubation_d"                         
## [19] "fledging_age_d"                       
## [20] "longevity_y"                          
## [21] "male_maturity_d"                      
## [22] "inter_litter_or_interbirth_interval_y"
## [23] "female_body_mass_g"                   
## [24] "male_body_mass_g"                     
## [25] "no_sex_body_mass_g"                   
## [26] "egg_width_mm"                         
## [27] "egg_length_mm"                        
## [28] "fledging_mass_g"                      
## [29] "adult_svl_cm"                         
## [30] "male_svl_cm"                          
## [31] "female_svl_cm"                        
## [32] "birth_or_hatching_svl_cm"             
## [33] "female_svl_at_maturity_cm"            
## [34] "female_body_mass_at_maturity_g"       
## [35] "no_sex_svl_cm"                        
## [36] "no_sex_maturity_d"

2. We are interested in longevity (lifespan), body size and reproductive effort, and how this might vary depending on the taxonomy (specifically, with Class). Use select to pick relevant columns of the dataset and discard the others. Call the new data frame x. The relevant columns are the taxonomic variables (class, genus & species) and longevity_y, litter_or_clutch_size_n, litters_or_clutches_per_y, and adult_body_mass_g.

x <- amniote %>%
  select(
    class, genus, species,
    longevity_y, adult_body_mass_g,
    litter_or_clutch_size_n, litters_or_clutches_per_y
  )

3. Take a look at the first few entries in the species column. You will see that it is only the epithet, the second part of the Genus_species name, that is given.
   Use mutate and paste to convert the species column to a Genus_species by pasting the data in genus and species together. To see how this works, try out the following command, paste(1:3, 4:6). After you have created the new column, remove the genus column (using select and -genus).

x <- x %>%
  mutate(species = paste(genus, species)) %>%
  select(-genus)
head(x)

##   class               species longevity_y adult_body_mass_g
## 1  Aves Accipiter albogularis          NA           251.500
## 2  Aves      Accipiter badius          NA           140.000
## 3  Aves     Accipiter bicolor          NA           345.000
## 4  Aves  Accipiter brachyurus          NA           142.000
## 5  Aves    Accipiter brevipes          NA           203.500
## 6  Aves Accipiter castanilius          NA           159.375
##   litter_or_clutch_size_n litters_or_clutches_per_y
## 1                      NA                        NA
## 2                    3.25                         1
## 3                    2.70                        NA
## 4                      NA                        NA
## 5                    4.00                         1
## 6                      NA                        NA

4. What is the longest living species in the record? Use arrange to sort the data from longest to shortest longevity (longevity_y), and then look at the top of the file using head to find out. (hint: you will need to use reverse sort (-)). Cut and paste the species name into Google to find out more!

x <- x %>% arrange(-longevity_y)
head(x)

##      class                 species longevity_y adult_body_mass_g
## 1 Reptilia Chelonoidis duncanensis       177.0                NA
## 2 Reptilia  Aldabrachelys gigantea       152.0            117200
## 3 Reptilia          Testudo graeca       127.0              1430
## 4 Mammalia            Homo sapiens       122.5             62035
## 5 Mammalia   Balaenoptera physalus        95.0          38800000
## 6 Mammalia            Orcinus orca        90.0           4300000
##   litter_or_clutch_size_n litters_or_clutches_per_y
## 1                      NA                        NA
## 2                    13.5                  2.000000
## 3                     5.0                  3.200993
## 4                     1.0                  0.485000
## 5                     1.0                  0.400000
## 6                     1.0                  0.210000

5. Do the same thing but this time find the shortest lived species. Is this result plausible? What might have gone wrong?

x <- x %>% arrange(longevity_y)
head(x)

##      class              species longevity_y adult_body_mass_g
## 1 Mammalia    Lepus nigricollis  0.08333333          2196.875
## 2 Mammalia  Notoryctes caurinus  0.08333333            34.000
## 3 Mammalia Allactaga balikunica  0.08333333                NA
## 4 Mammalia    Allactaga bullata  0.08333333                NA
## 5 Mammalia       Geomys pinetis  0.08333333           195.750
## 6 Mammalia          Mus sorella  0.08333333            12.535
##   litter_or_clutch_size_n litters_or_clutches_per_y
## 1                    1.59                     7.150
## 2                    1.50                        NA
## 3                    2.52                        NA
## 4                    2.52                        NA
## 5                    1.77                     1.865
## 6                    5.20                     3.000

The two shortest-lived species, the northern marsupial mole (Notoryctes caurinus) and the Indian hare (Lepus nigricollis), are recorded as living 1 month (0.083 years). That is not plausible for mammals of this size (the hare weighs over 2 kg). These values are almost certainly errors in the database, for example a value in months entered as years, or a juvenile life stage recorded as the life span. Real data sets often contain errors like this, which is why it is worth checking extreme values before analysing them.

6. Use summarise and group_by to make a table summarising min, median and max life spans (longevity_y) for the three taxonomic classes in the database. Remember that you need to tell R to remove the NA values using a na.rm = TRUE argument.

x %>%
  group_by(class) %>%
  summarise(
    min = min(longevity_y, na.rm = TRUE),
    median = median(longevity_y, na.rm = TRUE),
    max = max(longevity_y, na.rm = TRUE)
  )

## # A tibble: 3 × 4
##   class       min median   max
##   <chr>     <dbl>  <dbl> <dbl>
## 1 Aves     0.75    12.5    75 
## 2 Mammalia 0.0833   8.67  122.
## 3 Reptilia 0.2     11.0   177

7. Body size is thought to be associated with life span. Let’s treat that as a hypothesis and test it graphically. Sketch what the graph would look like if the hypothesis were true, and if it was false. Plot adult_body_mass_g vs. longevity_y (using base R graphics). You should notice that this looks a bit messy.

plot(x$adult_body_mass_g, x$longevity_y)

8. Use mutate to create a new log-transformed variables, logMass and logLongevity. Use these to make a “log-log” plot. You should see that makes the relationship more linear, and easier to “read”.

x <- x %>%
  mutate(
    logMass = log(adult_body_mass_g),
    logLongevity = log(longevity_y)
  )

plot(x$logMass, x$logLongevity)

9. Is there a trade-off between reproductive effort and life span? Think about this as a hypothesis - sketch what the graph would look like if that were true, and if it was false. Now use the data to test that hypothesis: Use mutate to create a variable called logOffspring which is the logarithm of (number of litters/clutches per year × number of babies in each litter/clutch), i.e. log(litters_or_clutches_per_y * litter_or_clutch_size_n). Then plot logOffspring vs. logLongevity.

x <- x %>%
  mutate(logOffspring = log(
    litter_or_clutch_size_n * litters_or_clutches_per_y
  ))

plot(x$logOffspring, x$logLongevity)

10. To answer the final question (differences between taxonomic classes) you could now use filter to subset to particular classes and repeat the plot to see whether the relationship holds universally.

aves <- x %>%
  filter(class == "Aves")

plot(aves$logOffspring, aves$logLongevity)
title("Aves")

mammalia <- x %>%
  filter(class == "Mammalia")

plot(mammalia$logOffspring, mammalia$logLongevity)
title("Mammalia")

reptilia <- x %>%
  filter(class == "Reptilia")

plot(reptilia$logOffspring, reptilia$logLongevity)
title("Reptilia")

28.3 Temperature effects on egg laying dates

1. Import the egg data and take a look at it with head or str.

eggDates <- read.csv("CourseData/eggDates.csv")
head(eggDates)

##   boxNumber y2013 y2014 y2016 y2017 y2018 y2019
## 1         1   116   103    NA   107   111    NA
## 2         2    NA    NA    NA   114   118    NA
## 3         3    NA   102   108    NA    NA    NA
## 4         4   121   103   121   155   111   110
## 5         5   135   100   108   102   106   108
## 6         6   122   113   122    NA   124   149

2. Use pivot_longer to reformat the data, and make sure Year is numeric.

eggDates <- eggDates %>%
  pivot_longer(
    cols = starts_with("y"),
    names_to = "Year", values_to = "dayNumber",
    names_prefix = "y"
  ) %>%
  mutate(Year = as.numeric(Year))

head(eggDates)

## # A tibble: 6 × 3
##   boxNumber  Year dayNumber
##       <int> <dbl>     <int>
## 1         1  2013       116
## 2         1  2014       103
## 3         1  2016        NA
## 4         1  2017       107
## 5         1  2018       111
## 6         1  2019        NA

(Without names_prefix, you could instead remove the “y” afterwards with mutate(Year = gsub(pattern = "y", replacement = "", Year)).)

3. Calculate the mean egg date per year using summarise.

meanEgg <- eggDates %>%
  group_by(Year) %>%
  summarise(meanEggDate = mean(dayNumber, na.rm = TRUE))

meanEgg

## # A tibble: 6 × 2
##    Year meanEggDate
##   <dbl>       <dbl>
## 1  2013        125.
## 2  2014        108.
## 3  2016        115.
## 4  2017        112.
## 5  2018        117.
## 6  2019        111.

4. Import the weather data and take a look at it with head or str.

weather <- read.csv("CourseData/AarslevTemperature.csv")
str(weather)

## 'data.frame':    2559 obs. of  4 variables:
##  $ YEAR : int  2012 2012 2012 2012 2012 2012 2012 2012 2012 2012 ...
##  $ MONTH: int  4 4 4 4 4 4 4 4 4 4 ...
##  $ DAY  : int  11 12 13 14 15 16 17 18 19 20 ...
##  $ TEMP : num  8.1 6.6 5.9 6.2 6.5 4.3 4.3 7 8.9 6.5 ...

5. Use filter to subset to the months of interest (February to April) and then summarise the data to calculate the mean temperature in this period.

weather <- weather %>%
  filter(MONTH %in% 2:4) %>%
  group_by(YEAR) %>%
  summarise(meanSpringTemp = mean(TEMP))

head(weather)

## # A tibble: 6 × 2
##    YEAR meanSpringTemp
##   <int>          <dbl>
## 1  2012           7.86
## 2  2013           1.60
## 3  2014           6.03
## 4  2015           3.53
## 5  2016           4.50
## 6  2017           4.77

6. Join the two yearly datasets together using left_join.

joinedData <- left_join(meanEgg, weather, by = c("Year" = "YEAR"))
joinedData

## # A tibble: 6 × 3
##    Year meanEggDate meanSpringTemp
##   <dbl>       <dbl>          <dbl>
## 1  2013        125.           1.60
## 2  2014        108.           6.03
## 3  2016        115.           4.50
## 4  2017        112.           4.77
## 5  2018        117.           3.07
## 6  2019        111.           6.12

7. Plot a graph of meanSpringTemp on the x-axis and meanEggDate on the y-axis.

plot(joinedData$meanSpringTemp, joinedData$meanEggDate,
  xlab = "Mean February-April temperature (°C)",
  ylab = "Mean first egg date (day of year)"
)

Is laying date associated with spring temperature? There looks to be a negative relationship: in warmer springs, the birds start laying earlier. But we should not be too confident. There are only six years of data, and much of the pattern comes from one year: 2013 had a very cold spring and the latest mean laying date. With more years we could test the relationship properly (you will learn how to do this with linear regression later in the course).

28.4 Virtual dice

Let’s try the same kind of thing with the roll of (virtual) dice.

Here’s how to do one roll of the dice:

diceRoll <- 1:6
sample(diceRoll, 1)

[1] 3

1. Simulate 10 rolls of the dice, and make a table of results.

result <- sample(diceRoll, 10, replace = TRUE)
table(result)

## result
## 2 3 4 5 6 
## 2 3 1 1 3

Your table will probably look different to this, because it is a random process. You may notice that some numbers in the table are missing if some numbers were never rolled by our virtual dice.

2. Now simulate 90 rolls of the dice, and plot the results as a bar plot using geom_bar in ggplot. Add a horizontal line using geom_hline to show the expected result based on what you know about probability.

n <- 90
result <- data.frame(result = sample(diceRoll, n, replace = TRUE))

ggplot(result, aes(x = result)) +
  geom_bar() +
  geom_hline(yintercept = n / 6)

Figure 28.1: Barplot of 90 simulated dice throws

3. Try adjusting the code to simulate dice rolls with small (say, 30) and large (say, 600, or 6000, or 9000) samples. Observe what happens to the proportions, and compare them to the expected value. What does this tell you about the importance of sample size when trying to estimate real phenomena?

You only need to edit the n <- 90 line in the code above

The main message here is that as sample size increases you are more likely to obtain a good estimate of the true value of the phenomenon. You may also notice that, what would be considered a good sample size for the coin flipping (i.e. it recovers the true probability of 0.5 reasonably well) is NOT adequate for getting a good estimate of the probabilities for the dice.

This is because of the different number of possibilities: as the range of possible outcomes increases, the sample size requirements increase. In other words, choosing a good sample size is context-dependent.

Figure 28.2: Barplots of 30 and 9000 simulated dice throws

28.5 Sexual selection in Hercules beetles

Task 1: Explore the data.

hercules <- read.csv("CourseData/herculesBeetle.csv")
str(hercules)

## 'data.frame':    20 obs. of  2 variables:
##  $ width: num  14.3 17 14.7 19.6 19.4 ...
##  $ sex  : chr  "Male" "Male" "Male" "Male" ...

# Set the factor levels so that diff() gives Male minus Female
hercules$sex <- factor(hercules$sex, levels = c("Female", "Male"))

ggplot(hercules, aes(x = sex, y = width)) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(width = 0.1, height = 0, colour = "grey40") +
  labs(x = "Sex", y = "Width (cm)")

meanVals <- hercules %>%
  group_by(sex) %>%
  summarise(n = n(), mean = mean(width), sd = sd(width), .groups = "drop")
meanVals

## # A tibble: 2 × 4
##   sex        n  mean    sd
##   <fct>  <int> <dbl> <dbl>
## 1 Female    10  15.0  1.97
## 2 Male      10  17.4  2.22

Males are wider on average (about 17.4 cm) than females (about 15 cm), but the two groups overlap a lot.

Task 2: Formulate hypotheses.

Null hypothesis: sex is unrelated to width, so the sex labels are exchangeable.

Alternative hypothesis: males are wider than females. The research idea (sexual selection on males) predicts a direction before we look at the data, so a one-sided test is appropriate. (If we had no prior expectation of which sex is larger, we would use a two-sided test.)

Task 3: Calculate the observed statistic. The statistic is the mean width of males minus the mean width of females, so positive values support the alternative.

observedDiff <- hercules %>%
  group_by(sex) %>%
  summarise(mean = mean(width), .groups = "drop") %>%
  pull(mean) %>%
  diff()
observedDiff

## [1] 2.337433

Task 4: Construct a null distribution. We shuffle the sex labels, keeping 10 males and 10 females. Under the null hypothesis, sex is unrelated to width, so any allocation of the labels to the 20 measurements is equally plausible.

# Check the code with a few randomisations first
replicate(
  5,
  hercules %>%
    mutate(sex = sample(sex)) %>%
    group_by(sex) %>%
    summarise(mean = mean(width), .groups = "drop") %>%
    pull(mean) %>%
    diff()
)

## [1]  1.67250493 -0.13625029  0.06589415  1.40854982  0.68823537

B <- 10000

diffs <- data.frame(
  diffs = replicate(
    B,
    hercules %>%
      mutate(sex = sample(sex)) %>%
      group_by(sex) %>%
      summarise(mean = mean(width), .groups = "drop") %>%
      pull(mean) %>%
      diff()
  )
)

Task 5: Evaluate the evidence.

ggplot(diffs, aes(x = diffs)) +
  geom_histogram(bins = 30, fill = "steelblue", colour = "white") +
  geom_vline(xintercept = observedDiff, colour = "#B2182B", linewidth = 1.2) +
  labs(
    x = "Difference between means (Male - Female, cm)",
    y = "Number of randomised datasets"
  )

Because the alternative predicts a positive difference, we count randomised differences greater than or equal to the observed one, and use the +1 correction:

b <- sum(diffs$diffs >= observedDiff)
p_value <- (b + 1) / (B + 1)
b

## [1] 103

p_value

## [1] 0.01039896

For comparison, the exact one-sided p-value from all 184,756 possible allocations of 10 males and 10 females is 2,202/184,756 = 0.012, very close to the Monte Carlo estimate.

Task 6: Write a short results paragraph.

“Mean width was 17.37 cm in males and 15.03 cm in females (n = 10 per sex), a difference of 2.34 cm (Male minus Female). A one-sided Monte Carlo randomisation test, in which the sex labels were randomly reallocated 10,000 times while preserving group sizes, provided evidence that males are wider than females (p = 0.01). This is consistent with sexual selection on male size, although these data alone cannot show that sexual selection is the cause.”

Is the difference biologically important? A statistically clear difference is not automatically an important one. Here the difference (about 2.3 cm, or about 16% of mean female width) is fairly large relative to body size. To judge its importance we would want to know, for example, whether wider males win more fights or mate more often, how the difference compares with variation among populations, and whether width is a good measure of overall body size.

28.6 Sex differences in fine motor skills

Some people have suggested that there might be sex differences in fine motor skills in humans. Use the data collected from students in previous years to address this topic using t-tests. The relevant data set is called classData.csv, and the columns of interest are Gender and Precision.

Carry out a two-sample t-test.

1. Plot the data (e.g. with a box plot, or histogram)

classData <- read.csv("CourseData/classData.csv") %>%
  filter(Gender %in% c("Male", "Female")) %>% 
  filter(Year == 2019) #you may use a different filter

ggplot(classData, aes(x = Gender, y = Precision)) +
  geom_boxplot() +
  geom_jitter(width = 0.2)

2. Formulate null and alternative hypotheses.

Null hypothesis: the mean precision score is the same for males and females.

Alternative hypothesis: the mean precision score differs between males and females (a two-sided alternative, because we have no strong reason to expect a particular direction).

3. Use the t.test function to do the test.

t.test(Precision ~ Gender, data = classData)

## 
##  Welch Two Sample t-test
## 
## data:  Precision by Gender
## t = -1.1836, df = 22.601, p-value = 0.2489
## alternative hypothesis: true difference in means between group Female and group Male is not equal to 0
## 95 percent confidence interval:
##  -6.873689  1.873689
## sample estimates:
## mean in group Female   mean in group Male 
##             26.16667             28.66667

4. Write a sentence or two describing the results.

There was no clear difference in mean precision between females and males (Welch’s t-test: t = -1.18, df = 22.60, p = 0.249). The 95% confidence interval for the difference (female minus male) was -6.874 to 1.874, which includes zero. We therefore have no evidence of a difference, but this is not the same as showing that there is no difference: the interval shows which differences are still compatible with the data. (These are also observational data from a convenience sample, so we should be careful about generalising.)

28.7 Therapy for anorexia

A study was carried out looking at the effect of cognitive behavioural therapy on weight of people with anorexia. Weight was measured in week 1 and again in week 8. Use a paired t-test to assess whether weight changed between week 1 and week 8 during the treatment. (Think about this: the study has no control group, so can a change in weight on its own tell you that the treatment caused it?)

The data is called anorexiaCBT.csv

The data are in “wide format”. You may wish to convert it to “long format” depending on how you use the data. You can do that with the pivot_longer function, which rearranges the data:

1. Plot the data (e.g. with an interaction plot).

anorexiaCBT <- read.csv("CourseData/anorexiaCBT.csv",
  header = TRUE
)

anorexiaCBT_long <- anorexiaCBT %>%
  pivot_longer(
    cols = starts_with("Week"), names_to = "time",
    values_to = "weight"
  )

ggplot(anorexiaCBT_long, aes(
  x = time, y = weight,
  group = Subject
)) +
  geom_line() +
  geom_point() +
  xlab("Time") +
  ylab("Weight (lb)")

2. Formulate a null and alternative hypothesis.

Null hypothesis: the mean change in weight between week 1 and week 8 is zero.

Alternative hypothesis: the mean change in weight is not zero (two-sided). You could argue for a one-sided alternative (weight increases), but only if that was decided before looking at the data.

3. Use t.test to conduct a paired t-test.

Putting week 8 first means the estimate is week 8 minus week 1, so a weight gain is positive:

t.test(anorexiaCBT$Week08, anorexiaCBT$Week01, paired = TRUE)

## 
##  Paired t-test
## 
## data:  anorexiaCBT$Week08 and anorexiaCBT$Week01
## t = 2.2156, df = 28, p-value = 0.03502
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  0.2268902 5.7869029
## sample estimates:
## mean difference 
##        3.006897

4. Write a couple of sentences to report your result.

x <- t.test(anorexiaCBT$Week08, anorexiaCBT$Week01, paired = TRUE)

On average, weight increased by 3.0 lb (about 1.4 kg) between week 1 and week 8 (paired t-test: t = 2.22, df = 28.00, p = 0.035; 95% CI for the mean change 0.2 to 5.8 lb). This is evidence that weight changed during the treatment period. However, because there is no control group, it does not show that CBT caused the change: weight might have changed anyway over those weeks (for example because of other care the patients received, or because people tend to enter treatment when their weight is at its lowest).

28.8 Compare t-tests with randomisation tests

Try re-fitting some of these tests as randomisation tests (or analyse the randomisation test data using t.test). Do they give approximately the same results?

Then try answering the question - “are people who prefer dogs taller than those who prefer cats?” using the classData.csv. Can you think of any problems with this analysis?

The problem with the analysis is that it is likely to be confounded by gender. Men are on average taller than women. If men and women also differ in whether they prefer cats or dogs, then a height difference between “cat people” and “dog people” could simply reflect the different proportions of men and women in the two groups, rather than anything to do with pet preference.

28.9 Apple tree crop yield

Import the data (apples.csv) and analyse it using the techniques you have learned in the ANOVA lecture, and the previous worksheet, to answer the question “What is the effect of tree spacing on apple yields?”

1. Import and look at the data (e.g. summary or str or head)

apples <- read.csv("CourseData/apples.csv")
summary(apples)

##     spacing       yield      
##  Min.   : 6   Min.   : 64.1  
##  1st Qu.: 6   1st Qu.:108.2  
##  Median :10   Median :147.1  
##  Mean   :10   Mean   :145.4  
##  3rd Qu.:14   3rd Qu.:176.5  
##  Max.   :14   Max.   :282.3  
##               NAs    :28

2. Ensure that the explanatory variable (spacing) is defined as a categorical variable (i.e. a “factor”, in R-speak). You can use mutate and as.factor functions for this.

apples <- apples %>%
  mutate(spacing = as.factor(spacing))

3. Plot the data using ggplot (a box plot with (optionally) added jittered points would be good).

ggplot(apples, aes(x = spacing, y = yield)) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(width = 0.1, height = 0) +
  xlab("Spacing (feet)") +
  ylab("Yield (kg per tree)")

4. Fit an ANOVA model using lm.

apple_mod <- lm(yield ~ spacing, data = apples)

5. Check the model using a diagnostic plot (i.e. using autoplot from the ggfortify package).

library(ggfortify)
autoplot(apple_mod)

Look for points in the Q-Q plot lying close to the line (normality) and for the three columns of residuals in the Residuals vs Fitted and Scale-Location plots having a broadly similar spread (equal variance).

6. Use anova to get the ANOVA summary.

anova(apple_mod)

## Analysis of Variance Table
## 
## Response: yield
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## spacing    2  35801 17900.3  9.3851 0.0002003 ***
## Residuals 89 169750  1907.3                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

7. You should see that there are differences among treatments. But where are those differences? Use summary on your model to find out.

summary(apple_mod)

## 
## Call:
## lm(formula = yield ~ spacing, data = apples)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -92.389 -30.577  -3.516  33.192 117.628 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  120.566      7.382  16.332  < 2e-16 ***
## spacing10     35.924     11.073   3.244 0.001659 ** 
## spacing14     44.107     10.966   4.022 0.000121 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 43.67 on 89 degrees of freedom
##   (28 observations deleted due to missingness)
## Multiple R-squared:  0.1742, Adjusted R-squared:  0.1556 
## F-statistic: 9.385 on 2 and 89 DF,  p-value: 0.0002003

8. Use a Tukey test to make all the pairwise comparisons among the treatments.

TukeyHSD gives the estimated differences with 95% confidence intervals (adjusted for multiple comparisons). Note that aov (like lm) drops the rows with missing yields:

apple_tukey <- TukeyHSD(aov(yield ~ spacing, data = apples))
apple_tukey

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = yield ~ spacing, data = apples)
## 
## $spacing
##            diff        lwr      upr     p adj
## 10-6  35.923571   9.530402 62.31674 0.0046832
## 14-6  44.106700  17.967561 70.24584 0.0003526
## 14-10  8.183128 -19.396838 35.76309 0.7598733

The letters version from agricolae gives the same conclusions:

library(agricolae)
HSD.test(apple_mod, "spacing", console = TRUE)

## 
## Study: apple_mod ~ "spacing"
## 
## HSD Test for yield 
## 
## Mean Square Error:  1907.304 
## 
## spacing,  means
## 
##       yield      std  r       se  Min   Max     Q25   Q50     Q75
## 10 156.4893 45.60411 28 8.253364 64.1 248.6 123.725 153.2 190.225
## 14 164.6724 50.41401 29 8.109816 72.3 282.3 134.100 168.2 202.700
## 6  120.5657 35.32755 35 7.382033 66.4 185.5  90.900 116.4 151.350
## 
## Alpha: 0.05 ; DF Error: 89 
## Critical Value of Studentized Range: 3.370849 
## 
## Groups according to probability of means differences and alpha level( 0.05 )
## 
## Treatments with the same letter are not significantly different.
## 
##       yield groups
## 14 164.6724      a
## 10 156.4893      a
## 6  120.5657      b

9. Write a few sentences that summarise your findings. What biological processes do you think drive the effects that you have detected?

Mean yield per tree differed among the three spacings (ANOVA: F = 9.39, d.f. = 2 and 89, p = 0.0002; see the box plot above). Trees planted 6 feet apart yielded the least (mean 120.6 kg per tree), compared with 156.5 kg at 10 feet and 164.7 kg at 14 feet. Compared with the 6-foot spacing, yield per tree was 35.9 kg higher at 10 feet (95% CI 9.5 to 62.3; Tukey HSD) and 44.1 kg higher at 14 feet (95% CI 18.0 to 70.2). The difference between 10 and 14 feet was smaller and not clear (8.2 kg, 95% CI -19.4 to 35.8).

What biological processes might drive this? Trees that are closer together compete more strongly for light, water and nutrients, so each tree has fewer resources to put into growth and fruit. Shading of the lower canopy by neighbouring trees may be particularly important for fruit production. Beyond about 10 feet, the extra space seems to make less difference, perhaps because competition between neighbours is already weak.

Note that these are yields per tree. A grower is usually more interested in yield per hectare: closely spaced trees produce less each, but there are many more of them, so the total yield per hectare could still be higher at the 6-foot spacing. Which spacing is “best” therefore depends on the question you ask.

10. Optional. Instead of using a Tukey test, use the alternative “relevel” approach to make the missing comparison.

apples2 <- apples %>%
  mutate(spacing = relevel(spacing, ref = "10"))
apple_mod2 <- lm(yield ~ spacing, data = apples2)
summary(apple_mod2)

## 
## Call:
## lm(formula = yield ~ spacing, data = apples2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -92.389 -30.577  -3.516  33.192 117.628 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  156.489      8.253  18.961  < 2e-16 ***
## spacing6     -35.924     11.073  -3.244  0.00166 ** 
## spacing14      8.183     11.571   0.707  0.48128    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 43.67 on 89 degrees of freedom
##   (28 observations deleted due to missingness)
## Multiple R-squared:  0.1742, Adjusted R-squared:  0.1556 
## F-statistic: 9.385 on 2 and 89 DF,  p-value: 0.0002003

---

If you get this far, try using the ANOVA approach on one of the previous two-sample t-test examples (remember that ANOVA can be used when your single explanatory variable has TWO or more levels). You should notice that the results are the same whether you use the t.test function with var.equal = TRUE or the ANOVA approach with lm (for two groups, F = t²). The default Welch t-test gives slightly different results, because it does not assume equal variances.

28.10 Chirping crickets

1. Import the data

chirps <- read.csv("CourseData/chirps.csv")

2. Use mutate to convert Fahrenheit to Celsius (Google it)

chirps <- chirps %>%
  mutate(temperatureC = (temperature - 32) * (5 / 9))
head(chirps)

##   chirps temperature temperatureC
## 1   20.0        88.6     31.44444
## 2   16.0        71.6     22.00000
## 3   19.8        93.3     34.05556
## 4   18.4        84.3     29.05556
## 5   17.1        80.6     27.00000
## 6   15.5        75.2     24.00000

3. Plot the data

(A <- ggplot(chirps, aes(x = temperatureC, y = chirps)) +
  geom_point())

4. Fit a linear regression model with lm

chirp_mod <- lm(chirps ~ temperatureC, data = chirps)

5. Look at diagnostic plots to evaluate the model

library(ggfortify)
autoplot(chirp_mod)

6. Use anova to figure out if the effect of temperature is statistically significant.

anova(chirp_mod)

## Analysis of Variance Table
## 
## Response: chirps
##              Df Sum Sq Mean Sq F value    Pr(>F)    
## temperatureC  1 28.435 28.4348  29.248 0.0001195 ***
## Residuals    13 12.639  0.9722                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

7. Use summary to obtain information about the coefficients and \(R^2\)-value.

summary(chirp_mod)

## 
## Call:
## lm(formula = chirps ~ temperatureC, data = chirps)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.6181 -0.6154  0.0916  0.7669  1.5549 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   6.95531    1.79534   3.874 0.001918 ** 
## temperatureC  0.36540    0.06756   5.408 0.000119 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.986 on 13 degrees of freedom
## Multiple R-squared:  0.6923, Adjusted R-squared:  0.6686 
## F-statistic: 29.25 on 1 and 13 DF,  p-value: 0.0001195

8. Summarise the model in words.

There is a statistically significant association between temperature and chirp frequency (F = 29.2482, d.f. = 1,13, p < 0.001). The equation of the fitted model is: Chirp Freq = 0.37(\(\pm\) 0.07) \(\times\) Temp + 6.96(\(\pm\) 1.80). The model explains 69% of the variation in chirp frequency (\(R^2\) = 0.692).

9. Add model fit line to the plot.

A + geom_smooth(method = "lm")

10. Can I use cricket chirp frequency as a kind of thermometer?

Yes, using the equation from the model. Ask an instructor if you can’t figure this out.

28.11 Fish behaviour

1. Import the data set, fishPersonality.csv

fishPersonality <- read.csv("CourseData/fishPersonality.csv")

2. Plot the data (e.g. as a box plot)

ggplot(fishPersonality, aes(
  x = colouration, y = spawning,
  fill = personality
)) +
  geom_boxplot()

3. Fit an ANOVA model using lm.

mod_A <- lm(spawning ~ personality + colouration +
  personality:colouration, data = fishPersonality)

4. Look at diagnostic plots from this model (autoplot)

library(ggfortify)
autoplot(mod_A)

The diagnostic plots look acceptable: no strong pattern in the residuals and no extreme points.

5. Use anova to get an Analysis of Variance summary table. Start with the interaction row.

anova(mod_A)

## Analysis of Variance Table
## 
## Response: spawning
##                         Df Sum Sq Mean Sq F value    Pr(>F)    
## personality              1   7781    7781  5.5032   0.02629 *  
## colouration              1  55862   55862 39.5074 8.519e-07 ***
## personality:colouration  1   4118    4118  2.9123   0.09898 .  
## Residuals               28  39591    1414                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The interaction is not significant (F = 2.91, d.f. = 1 and 28, p = 0.099), so there is no clear evidence that the effect of personality depends on colour pattern. This is not the same as showing that there is no interaction. The group means in the plot are clearly not parallel (proactive fish spawn about 54 more eggs than reactive fish in heterogeneous fish, but only about 9 more in homogeneous fish), and with only 8 fish per group the test has little power to detect an interaction of this size. The honest conclusion is “no clear evidence”, and a larger study would be needed to settle it.

6. The interaction is not needed, so drop it and refit the model with the main effects only.

mod_B <- lm(spawning ~ personality + colouration, data = fishPersonality)
anova(mod_B)

## Analysis of Variance Table
## 
## Response: spawning
##             Df Sum Sq Mean Sq F value    Pr(>F)    
## personality  1   7781    7781  5.1628   0.03067 *  
## colouration  1  55862   55862 37.0635 1.252e-06 ***
## Residuals   29  43708    1507                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Both main effects are significant: personality (F = 5.16, d.f. = 1 and 29, p = 0.031) and colouration (F = 37.1, d.f. = 1 and 29, p < 0.001). Colouration explains much more of the variation than personality (sum of squares 55862 vs 7781).

7. Get the coefficient summary and confidence intervals of the final model, and interpret them.

summary(mod_B)

## 
## Call:
## lm(formula = spawning ~ personality + colouration, data = fishPersonality)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -72.969 -22.781  -1.344  24.891  71.469 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              207.16      11.89  17.427  < 2e-16 ***
## personalityReactive      -31.19      13.73  -2.272   0.0307 *  
## colourationHomogeneous    83.56      13.73   6.088 1.25e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 38.82 on 29 degrees of freedom
## Multiple R-squared:  0.5928, Adjusted R-squared:  0.5648 
## F-statistic: 21.11 on 2 and 29 DF,  p-value: 2.195e-06

confint(mod_B)

##                            2.5 %     97.5 %
## (Intercept)            182.84478 231.467719
## personalityReactive    -59.25997  -3.115034
## colourationHomogeneous  55.49003 111.634966

• The intercept (207.2) is the estimated spawning of the reference group: proactive, heterogeneous fish.
• personalityReactive (-31.2, 95% CI -59.3 to -3.1): reactive fish spawn about 31 fewer eggs than proactive fish of the same colour morph.
• colourationHomogeneous (83.6, 95% CI 55.5 to 111.6): homogeneous fish spawn about 84 more eggs than heterogeneous fish of the same personality.
• Each variable has only two levels, so these coefficients already compare the two groups and no post-hoc test is needed.

8. Sum up your findings with reference to the initial research questions.

• Both colour pattern and personality are associated with spawning. Colour pattern has the larger effect: homogeneous fish spawned about 84 more eggs than heterogeneous fish (95% CI 55 to 112; F = 37.1, d.f. = 1 and 29, p < 0.001), while proactive fish spawned about 31 more than reactive fish (95% CI 3 to 59; F = 5.16, d.f. = 1 and 29, p = 0.031).
• There was no clear evidence that the effect of personality depends on colour pattern (interaction F = 2.91, d.f. = 1 and 28, p = 0.099). The plot hints that the personality difference may be larger in heterogeneous fish, but with 8 fish per group this study cannot tell; it would need a larger sample.
• Do not conclude that “personality matters only in heterogeneous fish” from separate comparisons within each colour morph (e.g. Tukey tests on the four groups). A significant difference in one group and a non-significant one in the other does not show that the two differences are themselves different: that is exactly what the interaction test assesses.

28.12 Maze runner

1. Import the data and graph it (geom_boxplot). Try adding the points to the ggplot using the new (to you) function geom_dotplot(binaxis = "y", stackdir = "center").

maze <- read.csv("CourseData/maze.csv", stringsAsFactors = TRUE)
head(maze)

##     Age nErrors
## 1 Child       2
## 2 Child       4
## 3 Child       2
## 4 Child       5
## 5 Child       6
## 6 Child       0

(A <- ggplot(maze, aes(x = Age, y = nErrors)) +
  geom_boxplot() +
  geom_dotplot(binaxis = "y", stackdir = "center"))

2. Fit an appropriate GLM.

mod_A <- glm(nErrors ~ Age, data = maze, family = poisson)

3. Examine the diagnostic plots of the model (autoplot).

library(ggfortify)
autoplot(mod_A)

4. Get the analysis of variance (deviance) table (anova). What does this tell you?

anova(mod_A, test = "Chi")

## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: nErrors
## 
## Terms added sequentially (first to last)
## 
## 
##      Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                    19     36.672              
## Age   1   11.511        18     25.161 0.0006917 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

5. Obtain the summary table. What does this tell you?

summary(mod_A)

## 
## Call:
## glm(formula = nErrors ~ Age, family = poisson, data = maze)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  0.09531    0.30151   0.316   0.7519   
## AgeChild     1.09861    0.34815   3.156   0.0016 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 36.672  on 19  degrees of freedom
## Residual deviance: 25.161  on 18  degrees of freedom
## AIC: 71.845
## 
## Number of Fisher Scoring iterations: 5

6. Use the model to get the predicted average number of mistakes for each age group (plus/minus 95% confidence interval). Start from a data frame of the groups to predict for, newData <- data.frame(Age = c("Adult", "Child")). Remember that (i) the model summary is on the scale of the linear predictor, so exp(AgeChild) is the ratio of the Child mean to the Adult mean, not the Child mean itself; and (ii) the standard error of the Child mean is not in the summary table (it combines both coefficients), so use predict(..., se.fit = TRUE) as in the fox example: calculate the 95% CI as the fit plus/minus 1.96 times the standard error on the link scale, then back-transform. Ask for help if you get stuck.

The AgeChild coefficient, back-transformed, is the ratio of the Child mean to the Adult mean: children make about three times as many mistakes. (confint() gives a profile-likelihood interval, which can differ slightly from estimate plus/minus 1.96 SE.)

exp(coef(mod_A))

## (Intercept)    AgeChild 
##         1.1         3.0

exp(confint(mod_A))

##                 2.5 %   97.5 %
## (Intercept) 0.5711495 1.883736
## AgeChild    1.5656311 6.218153

The group means and their 95% confidence intervals come from predict():

newData <- data.frame(Age = c("Child", "Adult"))
pv <- predict(mod_A, newData, se.fit = TRUE)

newData <- newData %>%
  mutate(nErrors_LP = pv$fit) %>%
  mutate(lowerCI_LP = pv$fit - 1.96 * pv$se.fit) %>%
  mutate(upperCI_LP = pv$fit + 1.96 * pv$se.fit)

inverseFunction <- family(mod_A)$linkinv

newData <- newData %>%
  mutate(nErrors = inverseFunction(nErrors_LP)) %>%
  mutate(lowerCI = inverseFunction(lowerCI_LP)) %>%
  mutate(upperCI = inverseFunction(upperCI_LP))


A + geom_point(
  data = newData, aes(x = Age, y = nErrors),
  colour = "red"
) +
  geom_segment(data = newData, aes(
    x = Age, xend = Age,
    y = lowerCI, yend = upperCI
  ), colour = "red")

28.13 Forensic footprints

1. Import the reference data and create a binary (0/1) response for the probability of being male.

morph <- read.csv("CourseData/morphometry.csv") %>%
  mutate(male = ifelse(Sex == "Male", 1, 0))
head(morph)

##    Sex Height HandLength FootLength male
## 1 Male   1652        212        258    1
## 2 Male   1766        203        268    1
## 3 Male   1660        206        260    1
## 4 Male   1795        203        251    1
## 5 Male   1707        202        265    1
## 6 Male   1780        206        254    1

2. Plot the data.

(A <- ggplot(morph, aes(x = FootLength, y = male)) +
  geom_jitter(height = 0.05, alpha = 0.5))

3. Fit the binomial GLM - this is our “forensic tool”.

morphMod <- glm(male ~ FootLength, data = morph, family = binomial)

4. Check the model using DHARMa. The points fall on the QQ diagonal and the residual-versus-predicted plot is a random star-field with no patterns, so the binomial family is appropriate.

library(DHARMa)
sim_res <- simulateResiduals(fittedModel = morphMod)
plot(sim_res)

5. The analysis of deviance table and coefficient summary. Foot length is a highly significant predictor of sex, and the positive slope shows that the probability of being male increases with foot length.

anova(morphMod, test = "Chi")

## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: male
## 
## Terms added sequentially (first to last)
## 
## 
##            Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                         154    214.714              
## FootLength  1   117.29       153     97.423 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(morphMod)

## 
## Call:
## glm(formula = male ~ FootLength, family = binomial, data = morph)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -44.58985    7.12830  -6.255 3.97e-10 ***
## FootLength    0.17799    0.02841   6.265 3.74e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 214.714  on 154  degrees of freedom
## Residual deviance:  97.423  on 153  degrees of freedom
## AIC: 101.42
## 
## Number of Fisher Scoring iterations: 6

6. How good is the forensic tool? Predict the probability of being male for everyone in the reference data, classify as male when that probability exceeds 0.5, and compare with the truth. The confusion table and the overall accuracy show that the tool classifies about 90% of people correctly - good, but far from perfect because male and female foot lengths overlap.

morph <- morph %>%
  mutate(pMale = predict(morphMod, type = "response")) %>%
  mutate(prediction = ifelse(pMale > 0.5, "Male", "Female"))

# Confusion table: true sex versus predicted sex
table(trueSex = morph$Sex, prediction = morph$prediction)

##         prediction
## trueSex  Female Male
##   Female     68    7
##   Male        8   72

# Overall proportion correctly classified
mean(morph$prediction == morph$Sex)

## [1] 0.9032258

7. The crime scene. Predict the probability that the person with a 271 mm foot was male, with a 95% confidence interval, back-transforming to the probability scale with the inverse link function.

crimeScene <- data.frame(FootLength = 271)

pv <- predict(morphMod, newdata = crimeScene, se.fit = TRUE)
inverseFunction <- family(morphMod)$linkinv

crimeScene <- crimeScene %>%
  mutate(
    probMale = inverseFunction(pv$fit),
    lowerCI = inverseFunction(pv$fit - 1.96 * pv$se.fit),
    upperCI = inverseFunction(pv$fit + 1.96 * pv$se.fit)
  )
crimeScene

##   FootLength  probMale   lowerCI   upperCI
## 1        271 0.9745808 0.9176298 0.9924785

The probability of being male is very high (well above 0.9), so the evidence points strongly to a male culprit. Note, however, that “very likely” is not “certain”: one woman in the reference data has a foot this large, so on its own this would be suggestive rather than conclusive evidence.

8. The second print (249 mm). Repeat the prediction for a foot length of 249 mm.

crimeScene2 <- data.frame(FootLength = 249)
pv2 <- predict(morphMod, newdata = crimeScene2, se.fit = TRUE)

crimeScene2 <- crimeScene2 %>%
  mutate(
    probMale = inverseFunction(pv2$fit),
    lowerCI = inverseFunction(pv2$fit - 1.96 * pv2$se.fit),
    upperCI = inverseFunction(pv2$fit + 1.96 * pv2$se.fit)
  )
crimeScene2

##   FootLength  probMale   lowerCI   upperCI
## 1        249 0.4330632 0.3140979 0.5602808

Here the probability of being male is close to 0.5, so the tool is essentially guessing. A foot length of 249 mm sits almost exactly at the “tipping point” where males and females are equally likely (around 250 mm, roughly the midpoint between the average female and male foot lengths), so this print tells us very little about the culprit’s sex. This is an important lesson: a prediction is only as useful as the certainty attached to it.

28.14 Power analysis: questions in the chapter

These answers are approximate: because the simulations are random, your numbers will differ slightly.

1. What sample size would give you 80% power? About 26 per group. You can find this by increasing sampleSize until the power reaches 0.8:

set.seed(42)
sampleSize <- 26 # vary this until power > 0.8
result <- replicate(
  1000,
  t.test(
    rnorm(sampleSize, mean = 12, sd = 3),
    rnorm(sampleSize, mean = 14.4, sd = 3)
  )$p.value
)
mean(result < 0.05)

## [1] 0.804

2. What power would you have if the variation was larger (sd = 4)? With 10 per group, the power falls from about 40% to about 24%. More noise makes the signal harder to detect.

3. What power would you have if the difference was only 10%? The treatment mean would then be 12 \(\times\) 1.1 = 13.2 mm, and with 10 per group the power falls to about 13%. Small effects are much harder to detect.

28.15 Snails on the move

Simulate a t-test-based analysis in R to figure out what sample size would result in 80% power.

set.seed(42)
sampleSize <- 40 # vary this until power > 0.8
result <- replicate(
  1000,
  t.test(
    rnorm(sampleSize, mean = 7.4, sd = 2.76),
    rnorm(sampleSize, mean = 7.4 * 1.25, sd = 2.76)
  )$p.value
)
sum(result < 0.05) / 1000

## [1] 0.841

You should find that you need a sample size of about 40 per group.

28.16 Mouse lemur strength

What difference in strength could you reliably detect (with power >80%) with these sample sizes?

set.seed(42)
diff <- 185 # Vary this to give a difference to 2nd group t-test
result <- replicate(
  1000,
  t.test(
    rnorm(25, mean = 600, sd = 145),
    rnorm(8, mean = 600 - diff, sd = 145)
  )$p.value
)
sum(result < 0.05) / 1000

## [1] 0.828

You should find that you could detect a difference of about 185g or a mean for old animals of 415g which is 70% of the young animals i.e. a 30% reduction.

Arnold, Todd W. 2010. “Uninformative Parameters and Model Selection Using Akaike’s Information Criterion.” Journal of Wildlife Management 74 (6): 1175–78. https://doi.org/10.2193/2009-367.

Burnham, Kenneth P., and David R. Anderson. 2002. Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach. 2nd ed. New York: Springer.