Chapter 6 Work on your own data

6.1 Yuheng Sun

6.1.1 Load-in data

# install.packages("readxl")
library(readxl)
# annoying code to get data from website
temp = tempfile(fileext = ".xlsx")
dataURL <- "https://stulp.gmw.rug.nl/schier/data/LHT_data.xlsx"
download.file(dataURL, destfile = temp, mode = 'wb')

sun <- readxl::read_xlsx(temp)

6.1.2 Visualisation I: variation in fitness

library(dplyr)
library(ggplot2)

fit_per_ind <- sun |> # "|>" means "then do this"
  # for each bird
  group_by(BirdID) |>
  # summarise total number of recruits
  summarise(total_recruits = sum(NumberOfAnnualRecruit, na.rm = TRUE)) 

ggplot(fit_per_ind, aes(x = total_recruits)) + 
  geom_histogram(binwidth = 1, fill = "lightblue") +
  labs(x = "Total number of recruits", y = "Frequency") +
  theme_minimal() 

6.1.3 Visualisation II: age and recruitment

library(ggplot2)

ggplot(sun, aes(x = Age, y = NumberOfAnnualRecruit, colour = RearingEnv, fill = RearingEnv)) + 
  geom_jitter(alpha = 0.8) +
  geom_smooth(method = "lm") +
  labs(x = "Age", y = "Number of recruits", 
       colour = "rearing\nenvironment", fill = "rearing\nenvironment") +
  scale_x_continuous(breaks = c(0:9)) +
  scale_colour_brewer(palette = "Set2", labels = c("N" = "Noisy", "Q" = "Quiet"),
                      aesthetics = c("colour", "fill")) +
  theme_classic() +
  theme(legend.position = c(0.9, 0.9))

6.1.4 Visualisation III: senescence trajectories

“What I would like to visualise: Reproductive senescence trajectories of the birds reared in different early-life environments More specifically, the y-axis would be ARO, the x-axis would be DeltaAge, and the groups would be RearingEnv = Q and RearingEnv = N.(in other word, a curve for Q and a curve for N, respectively)”

ggplot(sun, aes(x = DeltaAge, y = NumberOfAnnualRecruit, 
                colour = RearingEnv, fill = RearingEnv)) +
  geom_jitter(size = 1, colour = "lightgrey") +
  # Add points for averages
  geom_point(stat = "summary", fun = "mean", size = 3, alpha = 0.8) +
  # Add errorbars
  geom_errorbar(stat = "summary", fun.data = "mean_se", 
                size = 1, alpha = 0.8, width = 0) +
  # Add cubic function
  geom_smooth(method = "lm", formula = "y ~ x + I(x^2) + I(x^3)") +
  facet_wrap(~factor(RearingEnv, levels = c("N", "Q"), labels = c("Noisy", "Quiet"))) +
  labs(x = "Delta age (years)", y = "Number of recruits", 
       colour = "rearing\nenvironment", fill = "rearing\nenvironment") +
  #scale_x_continuous(breaks = c(0:9)) +
  scale_colour_brewer(palette = "Set2", labels = c("N" = "Noisy", "Q" = "Quiet"),
                      aesthetics = c("colour" , "fill")) +
  theme_minimal() +
  theme(
    legend.position = "none",
    panel.grid.minor = element_blank()
  )

6.1.5 Visualisation IV: model estimates

You’d like to control for this: Year: the year when the bird bred NatalEnv: the environment where the bird was born (not necessarily the same as the rearing environment), also with two levels (Q/N) BroodSizeHatchling: size of the brood where the bird was born MeanAge: mean of the age that the bird presented in the dataset. This is to control for the between-individual effect NatalBroodRef: reference of the natal brood FosterBrood: reference of the rearing brood BirdID: bird ID

Below model leaves out NatalBroodRef and FosterBrood are not included because that requires more sophisticated models.

This is a simple intercept only model:

# install.packages("lme4")
library(lme4)

mod <- lmer(NumberOfAnnualRecruit ~ DeltaAge + Year + NatalEnv + BroodSizeHatchling +
              MeanAge + (1|BirdID), data = sun)
summary(mod)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## NumberOfAnnualRecruit ~ DeltaAge + Year + NatalEnv + BroodSizeHatchling +  
##     MeanAge + (1 | BirdID)
##    Data: sun
## 
## REML criterion at convergence: 851.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1309 -0.6277 -0.2805  0.3611  5.2639 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  BirdID   (Intercept) 0.1871   0.4325  
##  Residual             1.3813   1.1753  
## Number of obs: 256, groups:  BirdID, 115
## 
## Fixed effects:
##                    Estimate Std. Error t value
## (Intercept)        78.53942   45.27273   1.735
## DeltaAge           -0.19072    0.07713  -2.473
## Year               -0.03959    0.02259  -1.752
## NatalEnvQ           0.48756    0.23413   2.082
## BroodSizeHatchling  0.19558    0.09210   2.123
## MeanAge             0.32972    0.11037   2.987
## 
## Correlation of Fixed Effects:
##             (Intr) DeltAg Year   NtlEnQ BrdSzH
## DeltaAge     0.293                            
## Year        -1.000 -0.293                     
## NatalEnvQ   -0.125 -0.035  0.120              
## BrdSzHtchln  0.040  0.014 -0.049  0.130       
## MeanAge      0.316  0.094 -0.320 -0.067  0.029
sun <- sun |>
  mutate(predictions = predict(mod))

Below graph has “raw” datapoints (so not controlled for any variables), whereas the fitted lines go through the prediction of mixed model (model above).

ggplot(sun, aes(x = DeltaAge, y = NumberOfAnnualRecruit, 
                colour = RearingEnv, fill = RearingEnv)) +
  geom_jitter(size = 1, colour = "lightgrey") +
  # Add points for averages
  geom_point(stat = "summary", fun = "mean", size = 3, alpha = 0.8) +
  # Add errorbars
  geom_errorbar(stat = "summary", fun.data = "mean_se", 
                size = 1, alpha = 0.8, width = 0) +
  # Add cubic function
  # THIS IS DIFFERENT NOW
  geom_smooth(
    aes(y = predictions),
    method = "lm", formula = "y ~ x + I(x^2) + I(x^3)"
  ) +
  facet_wrap(~factor(RearingEnv, levels = c("N", "Q"), labels = c("Noisy", "Quiet"))) +
  labs(x = "Delta age (years)", y = "Number of recruits", 
       colour = "rearing\nenvironment", fill = "rearing\nenvironment") +
  #scale_x_continuous(breaks = c(0:9)) +
  scale_colour_brewer(palette = "Set2", labels = c("N" = "Noisy", "Q" = "Quiet"),
                      aesthetics = c("colour" , "fill")) +
  theme_minimal() +
  theme(
    legend.position = "none",
    panel.grid.minor = element_blank()
  )

6.2 Life history theory data

These data come from this beautiful resource: https://esajournals.onlinelibrary.wiley.com/doi/10.1890/15-0846R.1

See here for explanation variables.

library(readr)
library(dplyr)
library(ggplot2)
lht <- read_csv("https://stulp.gmw.rug.nl/schier/data/Amniote_Database_Aug_2015.csv", 
                na = c("", "NA", -999))

6.2.1 Select sample and create new variables

lht <- lht |>
  # only select mammals
  filter(class == "Mammalia") |>
  # create new variable (named colour) translating mammals into three classes
  mutate(
    colour = case_when(
      species == "sapiens" ~ "Humans",
      order == "Primates" ~ "Primates",
      TRUE ~ "Other mammals"
    )
  )

6.2.2 Visualisation

ggplot(lht, aes(x = longevity_y, y = female_maturity_d)) +
  # first draw other mammals
  geom_point(data = filter(lht, colour == "Other mammals"), 
             aes(colour = "Other mammals"), size = 0.5, alpha = 0.8) +
  # then draw primates
  geom_point(data = filter(lht, colour == "Primates"), 
             aes(colour = "Primates"), size = 1, alpha = 0.8) +
  # then draw humans
  geom_point(data = filter(lht, colour == "Humans"), 
             aes(colour = "Humans"), size = 2, alpha = 0.8) +
  # log scale for x axis
  scale_x_continuous(trans = "log10", breaks = c(0.1, 1, 10, 100), 
                     labels = c(0.1, 1, 10, 100)) +
  # log scale for y axis
  scale_y_continuous(trans = "log10") +
  scale_color_manual(
    values = c("Other mammals" = "#CECED6",
               "Primates" = "#16D9F7",
               "Humans" = "#C71EB9"), 
    breaks = c("Other mammals", "Primates", "Humans")) +
  labs(x = "longevity (years)", y = "maturity female (days)", colour = NULL,
       caption = "SOURCE: Myhrvold et al., 2015") +
   theme(
     legend.position = "top",
     panel.grid.minor = element_blank(),
     axis.text = element_text(size = 20),
     axis.title = element_text(size = 25),
     strip.background = element_rect(),
     strip.text = element_text(size = 20),
     plot.title = element_text(size = 22, hjust = 0),
     legend.text = element_text(size = 20),
     plot.background = element_rect(colour = "lightgrey", fill = "white", linewidth = 1)
    )