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:
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:
rnoaa
package and choose stations that:
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?
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.
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.
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.
Here are all the records pooled over years and grouped together by similar species (number of records per group in corner of each plot).
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).
This is the length of laying + incubation time for each species.
Cluch size distribution by species after cleaning. Note x and y axis scales differ.
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.
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.
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).
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.
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.
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.
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.
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…
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
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 |
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 |