Overview

We wanted to look at whether experiencing cold snaps or generally cold weather influenced the timing of breeding in the subsequent year. For convenience, I’m using exactly the same dataset that was presented in Winkler et al (2020). That paper includes records from 2002-2016; we could of course expand those years, but these are already filtered based on Kelly’s work to exclude treatments that were deemed to invasive, so it at least seems like a good place to start.

The dataset includes basic info on the dates of each nesting attempt for each year that met the inclusion criteria. I’m combining those data with the Ithaca Airport weather data and then wrangling to get paired observations of lay date + weather in one year with lay date in the subsequent year. To account for overall weather/laying times in the subsequent year, I’m just using a random effect for that next year. I’m accounting for individual consistency in lay date with a fixed effect of when they started the previous year. So the models look something like this (with various iterations of temperature metrics):

ci_doy ~ y-1_ci_doi + y-1_temperature + (1|year) + (1|band_number)

Consistency of lay date

The Ecology paper spends a lot of time showing that lay date is consistent between years, but just to include that here, the basic plot looks like this.

pacman::p_load("lmer", "tidyverse")

d <- read.delim(here::here("1_raw_data/tres_lay_dates.txt"))

ggplot(data = d, mapping = aes(x = prev_lay, y = doy_lay)) + geom_point(size = 0.8, alpha = 0.5, col = "slateblue") +
  theme_classic() + geom_smooth(col = "coral3", method = "gam") + ylab("Lay DOY year n (May 1 = 1)") +
  xlab("Lay DOY year n-1 (May 1 = 1)")

So there really is a pretty strong correlation at least through those first ~25 days of the egg laying period. The Ecology paper includes those fans of outliers to the bottom right and upper left, but I actually wonder if it might be better to remove them (though it isn’t clear how to do that in an objective way). I suspect a lot of those are cases where a very early abandonment was not detected in one year or the other. I’m fitting the smoothing line here as a GAM rather than a linear model so it isn’t dragged down as much by those outliers. In any case, it’s a good start for asking how much lay dates change over-and-above what is expected from the individual identity based on temperature experienced.

Add in weather

I’m adding in several different weather metrics here for first year of each pair of observations. We could also add in weather in the second year, but I think for now it’s easier to try to capture that variation with a random effect of year. I also don’t have exact hatch or fledge dates in here because they weren’t stored in the Ecology paper, so I’m just using the estimated ranges. For a final analysis we should add those back in from the database and FoxPro. This code chunk is just looping through and adding in weather summarized in each of the time periods:

And in these two ways:

dw <- read.delim(here::here("1_raw_data/Input_Weather_Data.txt"))
dw <- subset(dw, dw$station == "Airport")

d2 <- subset(d, is.na(d$prev_lay) == FALSE)

for(i in 1:nrow(d2)){
  sub <- subset(dw, dw$year == d2$current_year[i] - 1)
  sub_lay <- subset(sub, sub$doy > d2$prev_lay[i] + 120 - 1 & sub$doy < d2$prev_lay[i] + d2$clutch[i] + 1 + 120 & sub$hour > 6 & sub$hour < 20)
  sub_inc <- subset(sub, sub$doy > d2$prev_lay[i] + 120 + d2$clutch[i] + 1 & sub$doy < d2$prev_lay[i] + d2$clutch[i] + 15 + 120 & sub$hour > 6 & sub$hour < 20)
  sub_d8 <- subset(sub, sub$doy > d2$prev_lay[i] + 120 + d2$clutch[i] + 15 & sub$doy < d2$prev_lay[i] + d2$clutch[i] + 24 + 120 & sub$hour > 6 & sub$hour < 20)
  sub_d15 <- subset(sub, sub$doy > d2$prev_lay[i] + 120 + d2$clutch[i] + 24 & sub$doy < d2$prev_lay[i] + d2$clutch[i] + 33 + 120 & sub$hour > 6 & sub$hour < 20)
  
  lay_his <- as.data.frame(unique(sub_lay$doy))
  inc_his <- as.data.frame(unique(sub_inc$doy))
  d8_his <- as.data.frame(unique(sub_d8$doy))
  d15_his <- as.data.frame(unique(sub_d15$doy))
  
  for(k in 1:nrow(lay_his)){
    lay_his$day_high[k] <- max(na.omit(subset(sub_lay$avg_temp_C, sub_lay$doy == lay_his[k, 1])))
  }
  for(k in 1:nrow(inc_his)){
    inc_his$day_high[k] <- max(na.omit(subset(sub_inc$avg_temp_C, sub_inc$doy == inc_his[k, 1])))
  }
  for(k in 1:nrow(d8_his)){
    d8_his$day_high[k] <- max(na.omit(subset(sub_d8$avg_temp_C, sub_d8$doy == d8_his[k, 1])))
  }
  for(k in 1:nrow(d15_his)){
    d15_his$day_high[k] <- max(na.omit(subset(sub_d15$avg_temp_C, sub_d15$doy == d15_his[k, 1])))
  }
  
  d2$avg_lay[i] <- mean(na.omit(sub_lay$avg_temp_C))
  d2$cs_lay[i] <- nrow(subset(lay_his, lay_his$day_high < 18.5))
  
  d2$avg_inc[i] <- mean(na.omit(sub_inc$avg_temp_C))
  d2$cs_inc[i] <- nrow(subset(inc_his, inc_his$day_high < 18.5))
  
  d2$avg_d8[i] <- mean(na.omit(sub_d8$avg_temp_C))
  d2$cs_d8[i] <- nrow(subset(d8_his, d8_his$day_high < 18.5))
  
  d2$avg_d15[i] <- mean(na.omit(sub_d15$avg_temp_C))
  d2$cs_d15[i] <- nrow(subset(d15_his, d15_his$day_high < 18.5))
}

Distribution of temperature

Just to get an idea of how much variation there is, I’m plotting histograms for the average temperature and number of ‘cold snap’ days in each of the periods defined above. Here is the distribution for average daytime temperatures:

ggplot(data = d2, mapping = aes(x = avg_lay)) + geom_histogram(fill = "goldenrod") + theme_classic() + xlab("Avg. Laying")
ggplot(data = d2, mapping = aes(x = avg_inc)) + geom_histogram(fill = "goldenrod") + theme_classic() + xlab("Avg. Incubation")
ggplot(data = d2, mapping = aes(x = avg_d8)) + geom_histogram(fill = "goldenrod") + theme_classic() + xlab("Avg. Day 0-8")
ggplot(data = d2, mapping = aes(x = avg_d15)) + geom_histogram(fill = "goldenrod") + theme_classic() + xlab("Avg. Day 8-15")

And for the number of ‘cold snap’ days:

ggplot(data = d2, mapping = aes(x = cs_lay)) + geom_histogram(fill = "goldenrod", bins = length(unique(d2$cs_lay))) + theme_classic() + xlab("CS Days Laying")
ggplot(data = d2, mapping = aes(x = cs_inc)) + geom_histogram(fill = "goldenrod", bins = length(unique(d2$cs_inc))) + theme_classic() + xlab("CS Days Incubation")
ggplot(data = d2, mapping = aes(x = cs_d8)) + geom_histogram(fill = "goldenrod", bins = length(unique(d2$cs_d8))) + theme_classic() + xlab("CS Days 0-8")
ggplot(data = d2, mapping = aes(x = cs_d15)) + geom_histogram(fill = "goldenrod", bins = length(unique(d2$cs_d15))) + theme_classic() + xlab("CS Days 8-15")

So there is a fair amount of variation in temperature experienced during these periods and in the number of cold snap days, which isn’t surprising.

Temperature vs. future laying

OK, now we’ve gotten around to being able to plot relationships between temperature and the time of laying in the next year. First I’m plotting this just for the raw values and then I’ll plot values that take into account change in laying date relative to each bird’s own expectation. First plotting for average temperatures then for number of cold snap days.

d2$delta_lay <- d2$doy_lay - d2$prev_lay

ggplot(data = d2, mapping = aes(x = avg_lay, y = delta_lay)) + geom_point(col = "slateblue", size = 0.7, alpha = 0.5) + 
  geom_smooth(method = "lm", col = "coral3") + theme_classic() + xlab("Avg n-1 C Laying") + ylab("Delta Lay Date")
ggplot(data = d2, mapping = aes(x = avg_inc, y = delta_lay)) + geom_point(col = "slateblue", size = 0.7, alpha = 0.5) + 
  geom_smooth(method = "lm", col = "coral3") + theme_classic() + xlab("Avg n-1 C Incubation") + ylab("Delta Lay Date")
ggplot(data = d2, mapping = aes(x = avg_d8, y = delta_lay)) + geom_point(col = "slateblue", size = 0.7, alpha = 0.5) + 
  geom_smooth(method = "lm", col = "coral3") + theme_classic() + xlab("Avg n-1 C Day 0-8") + ylab("Delta Lay Date")
ggplot(data = d2, mapping = aes(x = avg_d15, y = delta_lay)) + geom_point(col = "slateblue", size = 0.7, alpha = 0.5) + 
  geom_smooth(method = "lm", col = "coral3") + theme_classic() + xlab("Avg n-1 C Day 8-15") + ylab("Delta Lay Date")

ggplot(data = d2, mapping = aes(x = cs_lay, y = delta_lay)) + geom_point(col = "slateblue", size = 0.7, alpha = 0.5) + 
  geom_smooth(method = "lm", col = "coral3") + theme_classic() + xlab("Cold Days Laying") + ylab("Delta Lay Date")
ggplot(data = d2, mapping = aes(x = cs_inc, y = delta_lay)) + geom_point(col = "slateblue", size = 0.7, alpha = 0.5) + 
  geom_smooth(method = "lm", col = "coral3") + theme_classic() + xlab("Cold Days Incubation") + ylab("Delta Lay Date")
ggplot(data = d2, mapping = aes(x = cs_d8, y = delta_lay)) + geom_point(col = "slateblue", size = 0.7, alpha = 0.5) + 
  geom_smooth(method = "lm", col = "coral3") + theme_classic() + xlab("Cold Days 0-8") + ylab("Delta Lay Date")
ggplot(data = d2, mapping = aes(x = cs_d15, y = delta_lay)) + geom_point(col = "slateblue", size = 0.7, alpha = 0.5) + 
  geom_smooth(method = "lm", col = "coral3") + theme_classic() + xlab("Cold Days 8-15") + ylab("Delta Lay Date")

So it looks like there ar some interesting relationships here where birds that experience warmer conditions in one year advance their lay date by more days in the following year, but like we talked about this could all be a result of the fact that late laying birds both experience warmer temperatures and also move earlier in the following year (because of regression to the mean or some other reason) rather than anything directly to do with the temperature they actually experienced. So we need to account for laying date in the previous year and ask whether temperature adds any additional explanatory value over-and-above date.

Models accounting for lay date

Here I’m fitting some models that control for previous year lay date. I tried these a couple different ways. I think that it makes sense to have the lay date predictor fit as a GAM or spline like response since that plot at the top shows how the tails of the distribution curve away. I initially fit GAMs with a curved response for lay date but linear response for previous temperature. But I also wanted random effects for year and identity. I added those in to make GAMMs and those kind of worked, but the documentation is more complex for those GAM models with random effects and I’m not sure I was doing it right and the tables were kind of confusing. So, what I’m showing here is LMMs but with a spline fit for the previous lay date effect so that it is smoothed. I also fit just regular LMMs. The good thing is that all four of these approaches and the choice of number of knots in the spline seem to yield the same results as far as significance of effects in the model tables. So this should be worked out more as far as what makes the most sense, but I don’t think the choice really changes anything qualitatively. I’m showing a couple of those choice below but not all of them.

GAM for temperature

These ones are simple GAMs with smoothed previous lay date plus previous year temperature as predictors. No random effects.

library(mgcv)

d2$temp <- d2$avg_lay
m1 <- gam(doy_lay ~ s(prev_lay) + temp, data = d2)
d2$temp <- d2$avg_inc
m2 <- gam(doy_lay ~ s(prev_lay) + temp, data = d2)
d2$temp <- d2$avg_d8
m3 <- gam(doy_lay ~ s(prev_lay) + temp, data = d2)
d2$temp <- d2$avg_d15
m4 <- gam(doy_lay ~ s(prev_lay) + temp, data = d2)

tab_model(m1, m2, m3, m4,
          pred.labels = c("Intercept", "Previous Temperature", "Previous Lay Date (spline)"),
          dv.labels = c("Avg. Lay", "Avg. Incubation", "Avg Day 0-8", "Avg Day 8-15"))
  Avg. Lay Avg. Incubation Avg Day 0-8 Avg Day 8-15
Predictors Estimates CI p Estimates CI p Estimates CI p Estimates CI p
Intercept 14.21 11.69 – 16.72 <0.001 15.65 12.19 – 19.12 <0.001 15.61 11.88 – 19.33 <0.001 23.83 19.21 – 28.46 <0.001
Previous Temperature -0.02 -0.16 – 0.13 0.823 -0.09 -0.27 – 0.09 0.324 -0.08 -0.26 – 0.10 0.372 -0.45 -0.66 – -0.24 <0.001
Previous Lay Date (spline) 3.55 <0.001 3.57 <0.001 3.49 <0.001 3.89 <0.001
Observations 612 612 612 612
R2 0.156 0.158 0.157 0.182

Same thing but for cold snap days.

library(mgcv)

d2$temp <- d2$cs_lay
m1 <- gam(doy_lay ~ s(prev_lay) + temp, data = d2)
d2$temp <- d2$cs_inc
m2 <- gam(doy_lay ~ s(prev_lay) + temp, data = d2)
d2$temp <- d2$cs_d8
m3 <- gam(doy_lay ~ s(prev_lay) + temp, data = d2)
d2$temp <- d2$cs_d15
m4 <- gam(doy_lay ~ s(prev_lay) + temp, data = d2)

tab_model(m1, m2, m3, m4,
          pred.labels = c("Intercept", "Previous Cold Days", "Previous Lay Date (spline)"),
          dv.labels = c("Cold Days Lay", "Cold Days Inc", "Cold Days 0-8", "Cold Days 8-15"))
  Cold Days Lay Cold Days Inc Cold Days 0-8 Cold Days 8-15
Predictors Estimates CI p Estimates CI p Estimates CI p Estimates CI p
Intercept 14.67 13.89 – 15.45 <0.001 13.91 13.10 – 14.71 <0.001 13.87 13.27 – 14.47 <0.001 13.76 13.22 – 14.30 <0.001
Previous Cold Days -0.30 -0.55 – -0.04 0.022 0.01 -0.20 – 0.21 0.958 0.04 -0.28 – 0.36 0.796 0.20 -0.18 – 0.59 0.300
Previous Lay Date (spline) 3.41 <0.001 3.52 <0.001 3.55 <0.001 3.53 <0.001
Observations 612 612 612 612
R2 0.163 0.156 0.157 0.158

Interpretting the GAMs

So for average temperature the previous year, it seems that temperature during laying, incubation, and young nestlings doesn’t predict the next year’s lay date after accounting for the individual bird’s previous lay date. But, there is an association with average temperature during days 8-15, where experiencing warmer temperatures during that time results in breeding earlier the next year.

For the number of cold snap days, there seems to be a similar relationship with cold days during laying, though it isn’t as clear and it also isn’t supported in the LMMs that control for random effects (see below). Here it is saying that all else being equal, experiencing more cold days in year n is associated with initiating breeding later in year n + 1.

LMMs for temperature

These ones are LMMs that include a 4th degree basis spline for previous lay date plus one of the temperature metrics. Also include random effects for year and for female identity.

library(gamm4)
library(splines)

d2$temp <- d2$avg_lay
lmr1 <- lmer(doy_lay ~ bs(prev_lay, 4) + temp + (1|current_year) + (1|band), data = d2)
d2$temp <- d2$avg_inc
lmr2 <- lmer(doy_lay ~ bs(prev_lay, 4) + temp + (1|current_year) + (1|band), data = d2)
d2$temp <- d2$avg_d8
lmr3 <- lmer(doy_lay ~ bs(prev_lay, 4) + temp + (1|current_year) + (1|band), data = d2)
d2$temp <- d2$avg_d15
lmr4 <- lmer(doy_lay ~ bs(prev_lay, 4) + temp + (1|current_year) + (1|band), data = d2)

tab_model(lmr1, lmr2, lmr3, lmr4, 
          pred.labels = c("Intercept", "Previous Lay 1st", "Previous Lay 2nd", "Previous Lay 3rd", "Previous Lay 4th",
                          "Previous Temperature"),
          dv.labels = c("Avg. Lay", "Avg. Incubation", "Avg Day 0-8", "Avg Day 8-15"))
  Avg. Lay Avg. Incubation Avg Day 0-8 Avg Day 8-15
Predictors Estimates CI p Estimates CI p Estimates CI p Estimates CI p
Intercept 8.79 3.80 – 13.78 0.001 8.51 3.35 – 13.67 0.001 6.36 0.65 – 12.07 0.029 13.36 7.20 – 19.52 <0.001
Previous Lay 1st 0.13 -5.08 – 5.34 0.960 0.60 -4.64 – 5.85 0.822 0.29 -4.91 – 5.50 0.913 0.23 -4.95 – 5.41 0.930
Previous Lay 2nd 19.82 14.83 – 24.81 <0.001 20.26 15.06 – 25.46 <0.001 19.48 14.48 – 24.48 <0.001 20.24 15.27 – 25.21 <0.001
Previous Lay 3rd 1.35 -7.40 – 10.09 0.763 1.11 -7.64 – 9.87 0.803 0.22 -8.39 – 8.83 0.960 1.70 -6.90 – 10.29 0.699
Previous Lay 4th 9.51 1.97 – 17.04 0.013 10.38 2.44 – 18.31 0.010 8.92 1.21 – 16.63 0.023 10.23 2.71 – 17.75 0.008
Previous Temperature -0.10 -0.27 – 0.07 0.261 -0.09 -0.31 – 0.12 0.397 0.04 -0.16 – 0.25 0.692 -0.30 -0.52 – -0.07 0.010
Random Effects
σ2 27.44 27.61 27.56 27.57
τ00 2.27 band 2.16 band 2.22 band 2.10 band
2.17 current_year 1.97 current_year 2.05 current_year 1.31 current_year
ICC 0.14 0.13 0.13 0.11
N 14 current_year 14 current_year 14 current_year 14 current_year
341 band 341 band 341 band 341 band
Observations 612 612 612 612
Marginal R2 / Conditional R2 0.161 / 0.277 0.160 / 0.269 0.160 / 0.272 0.167 / 0.259

Same thing but for number of cold days.

library(gamm4)
library(splines)

d2$temp <- d2$cs_lay
lmr1 <- lmer(doy_lay ~ bs(prev_lay, 4) + temp + (1|current_year) + (1|band), data = d2)
d2$temp <- d2$cs_inc
lmr2 <- lmer(doy_lay ~ bs(prev_lay, 4) + temp + (1|current_year) + (1|band), data = d2)
d2$temp <- d2$cs_d8
lmr3 <- lmer(doy_lay ~ bs(prev_lay, 4) + temp + (1|current_year) + (1|band), data = d2)
d2$temp <- d2$cs_d15
lmr4 <- lmer(doy_lay ~ bs(prev_lay, 4) + temp + (1|current_year) + (1|band), data = d2)

tab_model(lmr1, lmr2, lmr3, lmr4, 
          pred.labels = c("Intercept", "Previous Lay 1st", "Previous Lay 2nd", "Previous Lay 3rd", "Previous Lay 4th",
                          "Previous Cold Days"),
          dv.labels = c("Cold Days Lay", "Cold Days Lay", "Cold Days 0-8", "Cold Days 8-15"))
  Cold Days Lay Cold Days Lay Cold Days 0-8 Cold Days 8-15
Predictors Estimates CI p Estimates CI p Estimates CI p Estimates CI p
Intercept 7.98 3.73 – 12.23 <0.001 6.70 2.10 – 11.30 0.004 7.14 3.02 – 11.26 0.001 6.43 2.21 – 10.66 0.003
Previous Lay 1st 0.34 -4.85 – 5.54 0.897 0.55 -4.76 – 5.86 0.840 0.20 -5.04 – 5.44 0.939 0.72 -4.51 – 5.94 0.788
Previous Lay 2nd 19.09 14.08 – 24.09 <0.001 19.93 14.72 – 25.15 <0.001 19.62 14.60 – 24.64 <0.001 19.96 14.96 – 24.96 <0.001
Previous Lay 3rd -0.74 -9.45 – 7.97 0.868 0.72 -8.02 – 9.46 0.872 0.28 -8.31 – 8.87 0.949 1.82 -6.98 – 10.62 0.685
Previous Lay 4th 8.82 1.27 – 16.36 0.022 9.77 1.92 – 17.62 0.015 9.35 1.78 – 16.92 0.015 9.63 2.10 – 17.16 0.012
Previous Cold Days -0.21 -0.53 – 0.10 0.179 0.05 -0.19 – 0.29 0.662 0.05 -0.30 – 0.41 0.762 0.32 -0.12 – 0.75 0.151
Random Effects
σ2 27.60 27.60 27.54 27.46
τ00 2.17 band 2.19 band 2.28 band 2.24 band
1.77 current_year 1.99 current_year 1.94 current_year 2.01 current_year
ICC 0.12 0.13 0.13 0.13
N 14 current_year 14 current_year 14 current_year 14 current_year
341 band 341 band 341 band 341 band
Observations 612 612 612 612
Marginal R2 / Conditional R2 0.164 / 0.268 0.159 / 0.270 0.158 / 0.270 0.161 / 0.273

Interpreting the LMMs

The relationship with day 8-15 temperature the previous year is supported in a very similar pattern in these LMMs even after adding in random effects for year and female identity. There isn’t really any support for any additional relationship with cold days.

We could of course summarize these temperature metrics differently and I’m not sure that the number of cold days is really that useful since it is a very rough metric… I’m also not sure that the age ranges I’ve used are the best and remember too that I’m here just assuming 14 day incubations rather than pulling out actual hatch dates from the database since they were not included in the data I downloaded from the Ecology paper. In any case, it seems like there might be something going on with day 8-15 temperatures, so I’m going to try to make some plots of that relationship that control for previous lay date by using residuals. This is also still including some of those outlier points in the first plot where they were super late layers in year n and it might be worth exploring thresholds to remove some of those points since they may have been renest attempts after an early initial failure.

Plot temperature controlling for lay date

Here I’m taking residuals from an LMM that includes only previous lay date and random effects and then plotting those residuals against previous year temperature. This should be a better indication of the relationship than the small plots above since it is accounting for variation already explained by previous late date independent of temperature.

m <- lmer(doy_lay ~ bs(prev_lay, 4) + (1|current_year) + (1|band), data = d2)
d2$m_resid <- residuals(m)

ggplot(data = d2, mapping = aes(x = avg_d15, y = m_resid)) + geom_point(col = "slateblue", size = 0.7, alpha = 0.5) +
  geom_smooth(col = "coral3", method = "lm") + theme_classic() + 
  xlab("Average C Day 8-15 Year n-1") + ylab("Residual Lay Date Year n") + 
  ylim(-10, 10)

Despite the significant p-value in the table above, this definitely isn’t a very impressive relationship. I’m also cropping out some observations here with very high positive residuals (they are included for fitting, but plotting them stretches the y-axis really high); it looks even less impressive with those included. So I don’t know…maybe there is a hint of an effect here, but I certainly wouldn’t stake a lot on it. This version above is not accounting for yearly variation and random individual effects, so I’ve added those adjustments in for the plot below.

ind <- as.data.frame(ranef(m)$band)
colnames(ind) <- "ind_ranef"
ind$band <- rownames(ind)
yr <- as.data.frame(ranef(m)$current_year)
colnames(yr) <- "yr_ranef"
yr$current_year <- rownames(yr)

d2 <- plyr::join(d2, ind, "band", type = "left", match = "first")
d2 <- plyr::join(d2, yr, "current_year", type = "left", match = "first")

d2$correct_lay <- d2$m_resid + d2$ind_ranef + d2$yr_ranef

ggplot(data = d2, mapping = aes(x = avg_d15, y = correct_lay)) + geom_point(col = "slateblue", size = 0.7, alpha = 0.5) +
  geom_smooth(col = "coral3", method = "lm") + theme_classic() + 
  xlab("Average C Day 8-15 Year n-1") + ylab("Band + Year Corrected Lay Date Year n") + 
  ylim(-10, 10)

So the relationship actually does look a little stronger to my eye after accounting for those other random effects and I guess this is why it comes out as significant in the model above.

Any thoughts on changes or different approaches that should be looked at here?

Temperature by clutch size plots

We also talked about looking at temperature effects on the next year’s clutch size. This has many of the same caveats described above since later nests tend to have smaller clutches so lay date needs to be accounted for in a similar way.

Here are plots of clutch size similar to those above for lay date. Points are slightly jittered for easier viewing.

d2$delta_lay <- d2$doy_lay - d2$prev_lay

ggplot(data = d2, mapping = aes(x = avg_lay, y = clutch)) + 
  geom_jitter(col = "slateblue", size = 0.7, alpha = 0.5, width = 0, height = 0.1) + 
  geom_smooth(method = "lm", col = "coral3") + theme_classic() + xlab("Avg n-1 C Laying") + ylab("Year n Clutch") +
  ylim(1, 8)
ggplot(data = d2, mapping = aes(x = avg_inc, y = clutch)) + 
  geom_jitter(col = "slateblue", size = 0.7, alpha = 0.5, width = 0, height = 0.1) + 
  geom_smooth(method = "lm", col = "coral3") + theme_classic() + xlab("Avg n-1 C Incubation") + ylab("Year n Clutch") +
  ylim(1, 8)
ggplot(data = d2, mapping = aes(x = avg_d8, y = clutch)) + 
  geom_jitter(col = "slateblue", size = 0.7, alpha = 0.5, width = 0, height = 0.1) + 
  geom_smooth(method = "lm", col = "coral3") + theme_classic() + xlab("Avg n-1 C Day 0-8") + ylab("Year n Clutch") +
  ylim(1, 8)
ggplot(data = d2, mapping = aes(x = avg_d15, y = clutch)) + 
  geom_jitter(col = "slateblue", size = 0.7, alpha = 0.5, width = 0, height = 0.1) + 
  geom_smooth(method = "lm", col = "coral3") + theme_classic() + xlab("Avg n-1 C Day 8-15") + ylab("Year n Clutch") +
  ylim(1, 8)

ggplot(data = d2, mapping = aes(x = cs_lay, y = clutch)) + 
  geom_jitter(col = "slateblue", size = 0.7, alpha = 0.5, width = 0.1, height = 0.1) + 
  geom_smooth(method = "lm", col = "coral3") + theme_classic() + xlab("Cold Days Laying") + ylab("Year n Clutch") +
  ylim(1, 8)
ggplot(data = d2, mapping = aes(x = cs_inc, y = clutch)) + 
  geom_jitter(col = "slateblue", size = 0.7, alpha = 0.5, width = 0.1, height = 0.1) + 
  geom_smooth(method = "lm", col = "coral3") + theme_classic() + xlab("Cold Days Incubation") + ylab("Year n Clutch") +
  ylim(1, 8)
ggplot(data = d2, mapping = aes(x = cs_d8, y = clutch)) + 
  geom_jitter(col = "slateblue", size = 0.7, alpha = 0.5, width = 0.1, height = 0.1) + 
  geom_smooth(method = "lm", col = "coral3") + theme_classic() + xlab("Cold Days 0-8") + ylab("Year n Clutch") +
  ylim(1, 8)
ggplot(data = d2, mapping = aes(x = cs_d15, y = clutch)) + 
  geom_jitter(col = "slateblue", size = 0.7, alpha = 0.5, width = 0.1, height = 0.1) + 
  geom_smooth(method = "lm", col = "coral3") + theme_classic() + xlab("Cold Days 8-15") + ylab("Year n Clutch") +
  ylim(1, 8)

It sure doesn’t look like there is much going on here. I’ll fit models like above just to check.

Temperature and clutch size models

LMMs with year + individual ID as random effects. I first fit these just like above with previous lay date as a cubic spline predictor, but now I’m thinking that doesn’t actually make much sense since it is modeling clutch size. For now I’m instead fitting current year lay date as a cubic spline and then asking whether previous year temperature explains any additional variation in clutch size. It looked like the results are basically the same either way. I could also include previous clutch size as a predictor, I suppose, but I’m leaving that out for now. I haven’t thought through these analyses as clearly so maybe there is something better/different that I should be doing…

Average temperatures previous year

library(gamm4)
library(splines)

d2$temp <- d2$avg_lay
lmr1 <- lmer(clutch ~ bs(doy_lay, 4) + temp + (1|current_year) + (1|band), data = d2)
d2$temp <- d2$avg_inc
lmr2 <- lmer(clutch ~ bs(doy_lay, 4) + temp + (1|current_year) + (1|band), data = d2)
d2$temp <- d2$avg_d8
lmr3 <- lmer(clutch ~ bs(doy_lay, 4) + temp + (1|current_year) + (1|band), data = d2)
d2$temp <- d2$avg_d15
lmr4 <- lmer(clutch ~ bs(doy_lay, 4) + temp + (1|current_year) + (1|band), data = d2)

tab_model(lmr1, lmr2, lmr3, lmr4, 
          pred.labels = c("Intercept", "Lay Date 1st", "Lay Date 2nd", "Lay Date 3rd", "Lay Date 4th",
                          "Previous Temperature"),
          dv.labels = c("Avg. Lay", "Avg. Incubation", "Avg Day 0-8", "Avg Day 8-15"))
  Avg. Lay Avg. Incubation Avg Day 0-8 Avg Day 8-15
Predictors Estimates CI p Estimates CI p Estimates CI p Estimates CI p
Intercept 6.20 5.45 – 6.94 <0.001 5.98 5.23 – 6.73 <0.001 5.58 4.77 – 6.38 <0.001 5.91 4.96 – 6.86 <0.001
Lay Date 1st -0.46 -1.26 – 0.34 0.259 -0.51 -1.31 – 0.29 0.214 -0.53 -1.32 – 0.27 0.192 -0.41 -1.22 – 0.40 0.320
Lay Date 2nd -0.27 -1.08 – 0.53 0.509 -0.35 -1.16 – 0.46 0.402 -0.30 -1.10 – 0.50 0.460 -0.27 -1.07 – 0.53 0.511
Lay Date 3rd -2.19 -3.64 – -0.75 0.003 -2.24 -3.68 – -0.79 0.002 -2.29 -3.73 – -0.84 0.002 -2.13 -3.59 – -0.68 0.004
Lay Date 4th -1.87 -2.79 – -0.96 <0.001 -1.93 -2.84 – -1.01 <0.001 -1.98 -2.89 – -1.07 <0.001 -1.83 -2.75 – -0.91 <0.001
Previous Temperature 0.00 -0.02 – 0.02 0.972 0.01 -0.01 – 0.04 0.229 0.03 0.01 – 0.06 0.008 0.01 -0.02 – 0.04 0.401
Random Effects
σ2 0.38 0.38 0.38 0.38
τ00 0.29 band 0.29 band 0.28 band 0.29 band
0.01 current_year 0.01 current_year 0.01 current_year 0.01 current_year
ICC 0.44 0.44 0.43 0.44
N 14 current_year 14 current_year 14 current_year 14 current_year
341 band 341 band 341 band 341 band
Observations 612 612 612 612
Marginal R2 / Conditional R2 0.066 / 0.476 0.068 / 0.475 0.077 / 0.478 0.067 / 0.474

So according to this model, warmer temperatures during nestling days 0-8 in the previous year are associated with larger than expected clutches in year n + 1 after controlling for lay date, year, and individual identity.

Same thing but for number of cold days.

library(gamm4)
library(splines)

d2$temp <- d2$cs_lay
lmr1 <- lmer(clutch ~ bs(doy_lay, 4) + doy_lay + temp + (1|current_year) + (1|band), data = d2)
d2$temp <- d2$cs_inc
lmr2 <- lmer(clutch ~ bs(doy_lay, 4) + temp + (1|current_year) + (1|band), data = d2)
d2$temp <- d2$cs_d8
lmr3 <- lmer(clutch ~ bs(doy_lay, 4) + temp + (1|current_year) + (1|band), data = d2)
d2$temp <- d2$cs_d15
lmr4 <- lmer(clutch ~ bs(doy_lay, 4) + temp + (1|current_year) + (1|band), data = d2)

tab_model(lmr1, lmr2, lmr3, lmr4,
          pred.labels = c("Intercept", "Lay Date 1st", "Lay Date 2nd", "Lay Date 3rd", "Lay Date 4th",
                          "Previous Cold Days"),
          dv.labels = c("Cold Days Lay", "Cold Days Lay", "Cold Days 0-8", "Cold Days 8-15"))
  Cold Days Lay Cold Days Lay Cold Days 0-8 Cold Days 8-15
Predictors Estimates CI p Estimates CI p Estimates CI p Estimates CI p
Intercept 6.07 5.42 – 6.73 <0.001 6.29 5.61 – 6.97 <0.001 6.23 5.57 – 6.89 <0.001 6.20 5.54 – 6.86 <0.001
Lay Date 1st -0.57 -1.37 – 0.22 0.160 -0.50 -1.30 – 0.30 0.221 -0.45 -1.25 – 0.34 0.266 -0.44 -1.25 – 0.37 0.286
Lay Date 2nd -0.18 -0.98 – 0.62 0.664 -0.32 -1.13 – 0.49 0.436 -0.28 -1.08 – 0.53 0.498 -0.28 -1.08 – 0.53 0.498
Lay Date 3rd -2.14 -3.57 – -0.70 0.004 -2.23 -3.67 – -0.78 0.003 -2.21 -3.66 – -0.76 0.003 -2.17 -3.63 – -0.72 0.003
Lay Date 4th -1.97 -2.87 – -1.06 <0.001 -1.93 -2.85 – -1.01 <0.001 -1.87 -2.78 – -0.96 <0.001 -1.85 -2.77 – -0.93 <0.001
Previous Cold Days 0.06 0.03 – 0.10 0.001 -0.01 -0.04 – 0.01 0.304 -0.03 -0.07 – 0.02 0.254 -0.01 -0.06 – 0.04 0.714
Random Effects
σ2 0.38 0.38 0.38 0.38
τ00 0.28 band 0.29 band 0.29 band 0.29 band
0.01 current_year 0.01 current_year 0.01 current_year 0.01 current_year
ICC 0.44 0.44 0.44 0.44
N 14 current_year 14 current_year 14 current_year 14 current_year
341 band 341 band 341 band 341 band
Observations 612 612 612 612
Marginal R2 / Conditional R2 0.086 / 0.485 0.067 / 0.476 0.067 / 0.476 0.066 / 0.475

And according to this, experiencing more cold days during laying in year n is associated with larger clutches in year n + 1 after controlling for lay date, year, and individual identity.

Temperature by clutch residual plots

Despite the p-values for those two effects, it doesn’t seem like there is much going on with clutch size. I’m making residual corrected plots for those two significant p-value effects just to see. These are made following the same approach as above: residuals of lay date + random effect residuals plotted against clutch size. First for average temperature day 0-8

d2$temp <- d2$avg_d8
lmr3 <- lmer(clutch ~ bs(doy_lay, 4) + (1|current_year) + (1|band), data = d2)
d2$resid <- residuals(lmr3)

ind <- as.data.frame(ranef(lmr3)$band)
colnames(ind) <- "ind_ranef_c"
ind$band <- rownames(ind)
yr <- as.data.frame(ranef(lmr3)$current_year)
colnames(yr) <- "yr_ranef_c"
yr$current_year <- rownames(yr)

d2 <- plyr::join(d2, ind, "band", type = "left", match = "first")
d2 <- plyr::join(d2, yr, "current_year", type = "left", match = "first")

d2$correct_clutch <- d2$resid + d2$ind_ranef_c + d2$yr_ranef_c

ggplot(data = d2, mapping = aes(x = avg_d8, y = d2$correct_clutch)) + geom_point(col = "slateblue", size = 0.7, alpha = 0.6) + 
  geom_smooth(col = "coral3", method = "gam") + theme_classic() + ylim(-2.5, 2.5) + xlab("Average C Day 0-8 Year n-1") +
  ylab("Residual Clutch Size Year n")

d2$temp <- d2$cs_lay
lmr3 <- lmer(clutch ~ bs(doy_lay, 4) + (1|current_year) + (1|band), data = d2)
d2$resid2 <- residuals(lmr3)

ind <- as.data.frame(ranef(lmr3)$band)
colnames(ind) <- "ind_ranef_c2"
ind$band <- rownames(ind)
yr <- as.data.frame(ranef(lmr3)$current_year)
colnames(yr) <- "yr_ranef_c2"
yr$current_year <- rownames(yr)

d2 <- plyr::join(d2, ind, "band", type = "left", match = "first")
d2 <- plyr::join(d2, yr, "current_year", type = "left", match = "first")

d2$correct_clutch2 <- d2$resid2 + d2$ind_ranef_c2 + d2$yr_ranef_c2

ggplot(data = d2, mapping = aes(x = cs_lay, y = d2$correct_clutch2)) + geom_jitter(col = "slateblue", size = 0.7, alpha = 0.6, width = 0.1) + 
  geom_smooth(col = "coral3", method = "lm") + theme_classic() + xlab("Cold Days Laying Year n-1") + ylim(-2.5, 2.5) +
  ylab("Residual Clutch Size Year n")

I don’t know…these effects still look like basically nothing. I think we’re just seeing significant p-values because the sample size is so large for these analyses. Maybe there is a very small effect, but I certainly wouldn’t want to bank on it as support for our predictions in these experiments.

Conclusion

Lay date is strongly related to previous year late date, largely because individuals are consistent in when they lay, with year of laying also having a substantial effect. It looks like previous year temperature has a big impact on next year’s lay date, but most of this effect is actually driven by the fact that late breeding birds breed earlier the next year and also experience warmer temperatures.

That being said, controlling for previous lay date, year, and individual identity, there is some indication that females experiencing warmer temperatures during nestling days 8-15 breed earlier the following year, though it isn’t a very strong effect.

There isn’t much indication that clutch size is related to any attribute of previous year temperature except for some very week relationships.

Bibliography

Winkler, D. W., Hallinger, K. K., Pegan, T. M., Taff, C. C., Verhoeven, M. A., Chang Van Oordt, D., Stager, M., Uehling, J. J., Vitousek, M. N., Andersen, M. J., & others. (2020). Full lifetime perspectives on the costs and benefits of lay date variation in tree swallows. Ecology.