General hypothesis

  • Building on Ryan’s 2020 PNAS paper about the timing of ‘cold snaps’ relative to breeding for tree swallows. He showed that in our Ithaca TRES population:

    • Springs are warmer over last 100 years
    • Despite that, the date of the last ‘cold snap’ (1/2/3 days with high < 18.5 C) has not changed
    • Birds breed earlier by ~9 days
    • Because of breeding earlier they are more often exposed to cold snaps while nestlings are at vulnerable age
    • Year to year, the date of last cold snap relative to breeding is strong predictor of breeding success

  • Idea for this paper is to essentially repeat that study but with many sites/years/species and to explicitly compare breeding performance for aerial insectivores (swallows, martins, flycatchers) to non aerial-insectivores (blubirds, chickadees, nuthatches, titmice, etc) that should not be as impacted by short cold snaps.

  • Predictions are the same as in the PNAS paper with the additions that:

    • Pattern should be present for aerial insectivores but not (or weaker) in non aerial insectivores
    • Geographic variation in timing of cold snaps relative to lay date/warming might modify where relationship is strongest

Cleaning Up

Cleaning breeding data for inclusion

  • Clean input data to these criteria
    • Has location.
    • Has clear date for at least egg laying/incubation start (ideally for hatch). When hatch date isn’t available, infer hatch based on lay date, clutch size, and species average incubation.
    • Has fate information fledge/failed (ideally with number fledged).
    • At least 500 records for the species in total.
    • For analyses by grid, has at least 10 records for the grid/year for that species.
    • I took out a few species that fit these criteria but weren’t worth including (e.g., starlings & house sparrows have tons of nest failures that I think must be intentional destructions?).
  • Filtering biologically implausible values
    • I plotted clutch size, lay date, and fledged number for each species to look at the distribution. Based on those plots, I filtered out values that were implausible for each species. I guess this was kind of a judgment call, but most of these are clearly implausible (e.g., breeding mid-January or a clutch size of 80). Probably this cleaning isn’t perfect. The exact ranges used are included in the main R code. We could make this into a supplemental table or something or just leave it in the code.

Cleaning and processing temperature data

  • Import temperature from NOAA using rnoaa package and choose stations that:
    • Have records starting in 1940 or earlier with max daily temperature (this leaves ~2,600 stations). Many of these stations go back <1900 but the number gets smaller that early.
  • Processing these temperature records:
    • Bin stations into the hexagonal grids (see below).
    • Average temperature for each day/year for all stations that appear in the same grid.
    • For each year/grid, identify the latest day that 1/2/3 day cold snap occurred in each year. Threshold = 18.5 C maximum daily temperature.
    • For counting cold snaps, limited to days of year 50 to 240 (otherwise you get fall cold snaps). Should we limit this more or by different grids?
    • Counted these

Standardizing breeding performance

  • For each species, calculate change in breeding phenology over time? For TRES we had records back to the 60/70s, but most of these only start in mid 90s so not sure how much we can detect here.

  • For each breeding record, calculate relative reproductive success. Will need to account for lat/long variation in average clutch/fledge size somehow. Also need to standardize across all species so that RS is a z-score of RS relative to species mean and relative to the expected RS for that particular location. This is actually going to be tricky because of course to the extent that temperature changes the clutch size pattern with latitude it might change some over the 25 years of data.

  • For each breeding attempt, calculate the offset in days from hatching to the last cold snap of the year using above criteria. Could get into more weather data specifically experienced for each nest but maybe that is best to skip at first.

  • Ask does date of last cold snap predict fledging RS (or fledge yes/no) and does the pattern differ for aerial insectivores vs. non aerial insectivores. Maybe need to also run this as phylogenetic correction for relationships?

Defining space used

Breeding grid

I made a hexagonal grid with 3 degree width hexagons. Because these are equal degree size, the actual land area varies a bit as you go north. I don’t think that is a problem since the total area doesn’t really factor into any analyses. I also plotted 1 or 2 degree grids, but those seem too small to work with in terms of having decent numbers of breeding records, etc. Of course the grid size is somewhat arbitrary here.


NOAA Stations

I downloaded daily weather data from a little over 2600 NOAA weather stations that started recording <1940 and had records from present day (some do have gaps in time). Here they are. You can see that these really only cover the US with a few in Canada, but the breeding records are also almost exclusively US so I think it’s fine to not worry about getting more records from farther north unless we find access to breeding records to go along with them.


Basic plots of breeding data

We have records from NestWatch at Cornell for a bunch of species going back to early/mid 90s (some older, but sample sizes get much smaller). We also have records for Purple Martins from a recent MartinWatch dataset archived with a JAB paper that we could add. Ryan has also was thinking we might be able to add a similarly large database of European nest records if we want to bring on a couple more collaborators.

Between the two US datasets, there are >400,000 nesting records, though that will definitely be reduced somewhat after filtering for having the info we need, etc and what exactly we decide to include.

Included breeding records

Here are all the records pooled over years and grouped together by similar species (number of records per group in corner of each plot).


Breeding records

This is the distribution of lay dates by day of year for each species. Aerial insectivores vs. not and cavity nesters vs. not are split by color. Numbers are the total number of records for the species in the database (though at this moment that number includes records that will be filtered out by cleaning).

Lay date


Incubation length

This is the length of laying + incubation time for each species.


Clutch size

Cluch size distribution by species after cleaning. Note x and y axis scales differ.


Fledge number


Records by year


Breeding records by space

Now that the grid is established and space is linked to the cleaned records. I can plot breeding records across the entire grid. For each of the 24 species, we can make plots that show geographic variation in lay date, clutch size, number fledged, success/failure. We can also split those out by different years. Obviously that could create a million plots, so I’m not sure what will end up being the best use here. It would be kind of cool to make an interactive plot (shiny) that could toggle different species on and toggle between different breeding records.

I’m just producing a few here as examples.

Tree swallow


Purple martin


Eastern bluebird


Six species comp


To standardize??

I initially was thinking that it would make sense to standardize all of the breeding attributes (lay date, fledge y/n + number, clutch size) by each grid, so that the reported values would all be relative to the expectation for that grid. Other than lay date there actually aren’t huge differences for most of the species I looked at. I’m also now thinking that a better approach may be to include a random effect of grid id in our models. That would allow for varying means by grid for the outcome variables and I think would accomplish the same thing I was envisioning with using relative values, but would be a little easier to explain. I do have the code setup in the main R script to do normalization within each grid and then you can make plots like these with relative values instead.

Cold snaps

In the PNAS paper, Ryan used the the long term records from the one station next to our field site to look at change in the date of the last cold snap over the last 100+ years. We can do a similar thing using the NOAA daily weather records from each grid. What I’ve done here is to average all the daily records within a grid to get an overall daily high temperature for each grid for the last 100+ years (I also have minimum temperature downloaded if we want to look at that too, but haven’t done anything with it yet).

What is a cold snap?

This is where radar might come in. I’ve set a threshold of 18.5 C for now and cold snaps = 1/2/3 days in a row where the maximum temperature does not rise above that value. We can pretty easily change that to a different global threshold and run again as a sensitivity analysis. Another approach though would be to define regional thresholds based on insect biomass from the radar data if Adriaan can provide that. Then we could apply a particular number to each grid. We won’t be able to use insect biomass directly since it likely won’t cover every grid/date when we want records, but it could be a good way to validate that our threshold is reasonable.

Spatial variation in cold snaps

This is just averaging the date of the last 1, 2, or 3 day cold snap within each grid using records starting in 1900 (only grids with >80 years of total data are saved, some peripheral grids have a smaller number of years of data). My ending threshold for cold snaps is day 240 (~end of August). For very northern grids they sometimes have cold snaps that happen all the way at the end…not sure yet what the best way to deal with those is. Will want to pull out something about cold snaps in temporal proximity to actual nesting records.


Temporal change in cold snaps

One of the big points of the PNAS paper was that the date of the latest cold snap hasn’t really changed even though springs have gotten warmer overall. We can look at the same thing for each grid using the averaged weather station data. Here, each blue line is a loess smoothed fit of the time series for one grid. The red line is the loess smoothed curve + confidence interval for all the grids together. There is little–if any–apparent change over the last 100 years, though there is obviously a lot of geographic variation. I think some of the wiggliness of the uppermost lines is because some of those grids in some years have cold snaps all the way until the end of my cutoff (day 240). I could also color lines here by latitude or something, though I think it’s pretty self explanatory that the later cold snap lines are from farther north.


Temperature change anomaly

We could look at change in temperature at the stations in a similar way, but I’m not sure that actually makes the most sense. I think it may make more sense to use one of the established data products that quantifies temperature anomalies based on a comparison to historical periods because they have much fancier ways of comparing things. The NOAA version gives a temperature anomaly in a grid per month, so the resolution is a bit more rough, but that is probably more accurate anyway. The download is a raster, so we can clip out and do raster math to get an average anomaly for each of our hexagonal grids for spring months. Maybe just April/May/June?

I’m downloading from here: http://berkeleyearth.org/archive/data/

NOAA also has a data product for this, but I found this Berkeley one easier to use and it is gridded at a smaller resolution (1 degree lat/long grids). The anomaly is expressed relative to a baseline of 1950-1980. I can pull out anomalies for every month/year, but here I’m showing the average anomaly for each month over the time range 1996-2020. Next I need to average this for each of the hexagonal bins to get an overall anomaly on our map and I could also pull that out for each year in our bins to match it specifically with breeding records if we want that.

Land temperature anomaly rasters



Anomaly by hexagon

To make the anomaly data comparable with everything else, I put it into the same hexagon grid used above. Here, the averag raster value is calculated for each polygon from the hex grid. This is showing averages for 1996-2020 for April/May/June, but it could also be calculated on a per year basis or for other months. It might make sense to link up the per year anomaly in the grid with all breeding records from that year. I think I wrote kind of a hacky way to do this that could probably be improved on, but it seems to work…


Anomaly vs. variability

One question we can ask at a grid level is whether there is any relationship between temperature anomaly and date of latest cold snap or change in latest cold snap (though it looks above like the cold snaps haven’t really changed). Probably for cold snaps we should actually calculate this in the same way by getting an average from 1950-1980 and then calculate a difference from that average for 1996-2020.

FILL IN HERE

Temperature + Breeding

Now to link temperature and breeding records. This is a very quick first pass using all records and looking at number fledged and fledged as a binary yes/no. Random effects are included for hexagonal grid ID and species ID. The effects here are exactly in the direction we expected, though the r-squared is very very low. These should probably be modified too though to include latitude itself (rather than just hex grouping) and maybe some other variables. This is also just a global last day of cold snap in that hex and not a specific link to the timing of the particular breeding attempt. I guess we probably also need to account for phylogeny in these models.

library(sjPlot)
mtab <- readRDS(here::here("3_r_scripts/initial_model.rds"))
mtab
  Number Fledged Fledged y/n
Predictors Estimates CI p Odds Ratios CI p
Intercept 3.47 3.08 – 3.86 <0.001 3.41 2.96 – 3.93 <0.001
Aerial Insectivore -0.09 -0.93 – 0.75 0.837 1.42 1.10 – 1.83 0.007
Last 3-day Cold Snap -0.02 -0.03 – -0.00 0.038 0.94 0.92 – 0.96 <0.001
Aerial * Cold Snap -0.11 -0.13 – -0.10 <0.001 0.96 0.94 – 0.98 <0.001
Random Effects
σ2 4.02 3.29
τ00 0.16 id.x 0.11 id.x
0.72 alpha4 0.10 alpha4
ICC 0.18 0.06
N 131 id.x 131 id.x
24 alpha4 24 alpha4
Observations 328223 328223
Marginal R2 / Conditional R2 0.002 / 0.181 0.010 / 0.069

Add in lay date

Here is the same thing but adding in lay date and an interaction between lay date and cold snap.

mtab2 <- readRDS(here::here("3_r_scripts/initial_model_add_lay.rds"))
mtab2
  Number Fledged Fledged y/n
Predictors Estimates CI p Odds Ratios CI p
Intercept 3.35 2.97 – 3.72 <0.001 3.23 2.76 – 3.79 <0.001
Aerial Insectivore 0.10 -0.69 – 0.90 0.795 1.50 1.12 – 2.02 0.007
Last 3-day Cold Snap -0.02 -0.03 – 0.00 0.061 0.93 0.91 – 0.96 <0.001
Lay Day of Year -0.32 -0.33 – -0.31 <0.001 0.90 0.89 – 0.91 <0.001
Aerial * Cold Snap -0.08 -0.10 – -0.06 <0.001 0.96 0.94 – 0.99 0.002
Cold Snap * Lay DOY 0.01 0.01 – 0.02 0.001 1.06 1.05 – 1.07 <0.001
Random Effects
σ2 3.87 3.29
τ00 0.26 id.x 0.17 id.x
0.64 alpha4 0.09 alpha4
ICC 0.19 0.07
N 129 id.x 129 id.x
24 alpha4 24 alpha4
Observations 287012 287012
Marginal R2 / Conditional R2 0.025 / 0.209 0.020 / 0.090