Rough Draft - 9915 assignment
Ashley Fang
2026-11-02
Introduction
Anchor article
“Light pollution is causing birds like the Australian magpie-lark to sing for longer” (ABC News, August 22, 2025).
https://www.abc.net.au/news/science/2025-08-22/light-pollution-extending-birdsong-hours/105673490
The ABC article reports on a study of 583 daytime bird species. On average, birds in the brightest areas were vocally active about 50 minutes longer each day. The article highlights the Australian magpie-lark (Grallina cyanoleuca) as one species that began calling earlier in bright areas.
The study behind the anchor article
Pease, B. S., & Gilbert, N. A. (2025). Light pollution prolongs avian activity. Science, 389(6762), 818-821. https://doi.org/10.1126/science.adv9472
The researchers used recordings from volunteer-run microphones around the world (BirdWeather) and compared when 583 daytime bird species sang with how much light pollution each location had, using satellite measurements.
Additional supporting studies/papers
My anchor article reports on a study of how light pollution affects bird behaviour. To understand it, I looked at five peer-reviewed papers on light at night and birds. This section summarizes what each one measured and which statistical method it used.
| # | Paper | What was measured | Main statistical method | Relevance |
|---|---|---|---|---|
| 1 | Pease & Gilbert (2025). Light pollution prolongs avian activity. Science, 389(6762), 818-821. | The time of first and last calls each day (relative to sunrise and sunset) for 583 bird species, compared with satellite measures of light pollution at each microphone | Regression models that test how call timing changes with light level, and whether this depends on traits like eye size, nest type and migration. Models account for repeated measures by species, place and week | The main study behind my anchor article |
| 2 | Dominoni et al. (2013). Individual-based measurements of light intensity provide new insights into the effects of artificial light at night on daily rhythms of urban-dwelling songbirds. Journal of Animal Ecology. | How much light each bird was exposed to at night, and when its daily activity started and ended | Linear Mixed-Effects Model - Because the study tracked individual wild birds over consecutive days, multiple data points were collected from the same individuals, creating nested/repeated measures data. They used LMMs to account for this non-independence by including individual bird identity as a random effect | Measures light at the level of the individual bird, and its data and R code are open |
| 3 | Da Silva et al. (2017). Experimental illumination of a forest: no effects of lights of different colours on the onset of the dawn chorus in songbirds. Royal Society Open Science. | When dawn singing started at eight forest edges lit with white, green, red or no light | Linear mixed-effect models (LMMs), which they fitted using maximum likelihood (ML) - By utilizing this specific mixed-modeling approach on a relative response variable, the researchers were able to perform direct day-to-day comparisons within each site. This effectively controlled for potentially confounding local variables like weather, changing seasons, or specific geographic layout differences across the eight forest sites | A controlled experiment that found no effect, which shows that results can differ from correlational studies |
| 4 | Osterhaus et al. (2025). Nocturnal flight call monitoring reveals in-flight behavioral alteration by avian migrants in response to artificial light at night. Biological Conservation. | Flight calls of migrating birds at night, in relation to weather, timing and lighting | Generalized linear mixed modeling (GLMM). Their code is shared as an R Markdown file | Very recent, and it extends the topic to migrating birds |
| 5 | Da Silva, Valcu & Kempenaers (2015). Light pollution alters the phenology of dawn and dusk singing in common European songbirds. Philosophical Transactions of the Royal Society B, 370(1667), 20140126. | Whether six songbird species sang at dawn and at dusk on each day at 12 woodland sites, half of them lit by street lighting | Generalized linear mixed models (a regression for yes/no outcomes that also accounts for repeated measures at the same sites), fit in R | Shows light affects when in the season birds start singing, not only the time of day |
What these papers have in common
Need to do this part still…what are the studies experiments or observations? What kinds of models do they use? Which methods will I demonstrate in my own analysis?
Data
The data come from the open data and code package shared by Pease and Gilbert on Zenodo (https://doi.org/10.5281/zenodo.15635107). I used the morning onset file, which has 2,640,731 rows. Each row is one bird species at one location on one date. I focused on one species, the Northern Cardinal, because the paper found it responds strongly to light (see following code which shows how I came to this decision). The key columns are the time of the first call (onset_min, in minutes from sunrise, where negative means before sunrise) and the light pollution at that location (light, satellite radiance).
Organizing the Raw Data from the Pease & Gilbert (2025) Study
library(tidyverse)
load("data/onset_data_conf_0.75_det_100_grid_10.RData")
ls()## [1] "final"glimpse(final)## Rows: 2,640,731
## Columns: 10
## $ lat <dbl> -25.768, -25.768, -25.768, -25.768, -25.768, -25.768, …
## $ lon <dbl> 28.2704, 28.2704, 28.2704, 28.2704, 28.2704, 28.2704, …
## $ date <date> 2023-03-01, 2023-03-04, 2023-03-05, 2023-03-06, 2023-…
## $ grid_ID_cell_5 <dbl> 906, 906, 906, 906, 906, 906, 906, 906, 906, 906, 906,…
## $ grid_ID_cell_10 <dbl> 237, 237, 237, 237, 237, 237, 237, 237, 237, 237, 237,…
## $ grid_ID_cell_15 <dbl> 110, 110, 110, 110, 110, 110, 110, 110, 110, 110, 110,…
## $ week <dbl> 9, 9, 10, 10, 10, 10, 11, 11, 11, 11, 11, 11, 12, 12, …
## $ sci_name <chr> "Cossypha caffra", "Cossypha caffra", "Cossypha caffra…
## $ value <dbl> 6.133333, 43.766667, 204.550000, -38.966667, 29.966667…
## $ avg_rad <dbl> 21.89, 21.89, 21.89, 21.89, 21.89, 21.89, 21.89, 21.89…# Seeing how much data is in this whole package
nrow(final)## [1] 2640731onset <- final %>%
select(sci_name, lat, lon, date, week, value, avg_rad) %>%
rename(onset_min = value, light = avg_rad) %>%
drop_na() %>%
mutate(log_light = log10(light + 1)) # light is very skewed
summary(onset)## sci_name lat lon date
## Length :2640731 Min. :-42.90 Min. :-158.123 Min. :2023-03-01
## N.unique : 583 1st Qu.: 38.87 1st Qu.: -85.735 1st Qu.:2023-08-14
## N.blank : 0 Median : 44.39 Median : -69.981 Median :2024-01-01
## Min.nchar: 9 Mean : 41.84 Mean : -37.698 Mean :2023-11-18
## Max.nchar: 31 3rd Qu.: 51.50 3rd Qu.: 5.927 3rd Qu.:2024-02-28
## Max. : 63.44 Max. : 175.282 Max. :2024-03-31
## week onset_min light log_light
## Min. : 1.00 Min. :-551.35 Min. : 0.00 Min. :0.0000
## 1st Qu.: 9.00 1st Qu.: -2.35 1st Qu.: 1.55 1st Qu.:0.4065
## Median :13.00 Median : 53.63 Median : 5.23 Median :0.7945
## Mean :20.25 Mean : 76.91 Mean : 11.45 Mean :0.8236
## 3rd Qu.:32.00 3rd Qu.: 154.95 3rd Qu.: 14.04 3rd Qu.:1.1772
## Max. :53.00 Max. : 616.85 Max. :432.49 Max. :2.6370Figuring out which species to use from the raw data (narrowing down stuff)
onset %>% count(sci_name, sort = TRUE) %>% head(10)## # A tibble: 10 × 2
## sci_name n
## <chr> <int>
## 1 Cyanistes caeruleus 77901
## 2 Erithacus rubecula 76525
## 3 Haemorhous mexicanus 75196
## 4 Pica pica 69312
## 5 Parus major 65426
## 6 Passer domesticus 65326
## 7 Cyanocitta cristata 59405
## 8 Cardinalis cardinalis 55372
## 9 Turdus merula 53598
## 10 Corvus brachyrhynchos 49421Picking a species from the data to work with
my_species <- "Cardinalis cardinalis"
one_sp <- onset %>% filter(sci_name == my_species)
my_species## [1] "Cardinalis cardinalis"nrow(one_sp)## [1] 55372Stats
Variables and Hypothesis
Dependent variable (y): onset
(onset_min), the time of a cardinal’s first call of the
day, in minutes from sunrise. It is a continuous numeric variable.
Negative values mean the bird called before sunrise.
Independent variable (x): light pollution
(log_light), the satellite-measured brightness at the
recording location, put on a log10 scale because the raw values are very
skewed. It is also a continuous numeric variable.
Other variables:
Week and month (
week,date): when the observation was made. Season may be a confounding variable, because birds call at different times across the year and the amount of light may differ by season too.Location (
lat,lon,site): where the microphone is. Many observations come from the same location, so they are not fully independent.Species (
sci_name): I kept only one species, the Northern Cardinal, so species does not vary in my analysis.
Hypothesis and prediction: birds may treat artificial light as daylight, so I predict that cardinals start calling earlier where there is more light pollution (onset ~ light). The null hypothesis is that light has no effect on onset (slope = 0), and the alternative is that onset changes with light (slope not equal to 0).
Plotting data
set.seed(42)
one_sp %>%
slice_sample(n = min(5000, nrow(one_sp))) %>%
mutate(month = month(date, label = TRUE)) %>%
ggplot(aes(x = log_light, y = onset_min, colour = month)) +
geom_point(alpha = 0.4) +
geom_smooth(method = "lm", colour = "black") +
labs(x = "Light pollution (log10 scale)",
y = "First call (minutes from sunrise)",
colour = "Month",
title = my_species) +
theme_light()
Each dot is the time of the first cardinal call at one location on one day. To keep the plot readable I only show a random sample of 5,000 of the 55,372 observations (the line in this plot is fitted to those 5,000 points, while the model below uses all of them).
Linear model
My question: do cardinals start calling earlier where there is more light at night?
m1 <- lm(onset_min ~ log_light, data = one_sp)
summary(m1)##
## Call:
## lm(formula = onset_min ~ log_light, data = one_sp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -467.43 -80.27 -43.34 54.99 403.80
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 64.349 1.076 59.814 < 2e-16 ***
## log_light -6.046 0.992 -6.095 1.1e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 108.1 on 55370 degrees of freedom
## Multiple R-squared: 0.0006705, Adjusted R-squared: 0.0006525
## F-statistic: 37.15 on 1 and 55370 DF, p-value: 1.101e-09Checking assumptions
library(ggfortify)
autoplot(m1, smooth.colour = NA)
The linear model estimates that each 10-fold increase in light is linked to first calls about 6 minutes earlier (linear regression: t55370 = -6.1, p < 0.001, R² = 0.0007). However, light explains less than 0.1% of the variation in call time. With over 55,000 observations, even a very small effect gets a tiny p-value, so the result is statistically clear but small. Call timing also depends on season, weather and location, which this simple model leaves out. Many observations also come from the same locations, so they are not fully independent (this is called pseudoreplication).
Fixing pseudoreplication: one row per site and season
Many rows come from the same microphone, so I averaged first-call time within each site and season. Now each site counts once per season. I defined the breeding season as weeks 13 to 28 (about April to mid-July), because the paper found the effects were strongest then, and kept only averages built from at least 5 days.
one_sp <- one_sp %>%
mutate(site = paste(lat, lon),
season = if_else(week >= 13 & week <= 28,
"Breeding (Apr-Jul)", "Non-breeding"))
site_means <- one_sp %>%
group_by(site, season) %>%
summarise(mean_onset = mean(onset_min),
log_light = mean(log_light),
n_days = n(), .groups = "drop") %>%
filter(n_days >= 5)
nrow(one_sp) # rows before## [1] 55372nrow(site_means) # rows after## [1] 1763ANCOVA: does the effect of light depend on season?
m2 <- lm(mean_onset ~ season * log_light, data = site_means)
autoplot(m2, smooth.colour = NA)
anova(m2)## Analysis of Variance Table
##
## Response: mean_onset
## Df Sum Sq Mean Sq F value Pr(>F)
## season 1 58312 58312 17.9126 2.432e-05 ***
## log_light 1 12911 12911 3.9660 0.04658 *
## season:log_light 1 9513 9513 2.9222 0.08755 .
## Residuals 1759 5726205 3255
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(m2)##
## Call:
## lm(formula = mean_onset ~ season * log_light, data = site_means)
##
## Residuals:
## Min 1Q Median 3Q Max
## -154.535 -40.877 -2.487 36.661 220.719
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 78.462 6.621 11.851 <2e-16 ***
## seasonNon-breeding 2.124 7.629 0.278 0.7808
## log_light -15.619 6.374 -2.450 0.0144 *
## seasonNon-breeding:log_light 12.383 7.244 1.709 0.0875 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 57.06 on 1759 degrees of freedom
## Multiple R-squared: 0.0139, Adjusted R-squared: 0.01222
## F-statistic: 8.267 on 3 and 1759 DF, p-value: 1.841e-05ggplot(site_means, aes(x = log_light, y = mean_onset, colour = season)) +
geom_point(alpha = 0.3) +
geom_smooth(method = "lm", se = TRUE) +
labs(x = "Light pollution (log10 scale)",
y = "Average first call (minutes from sunrise)",
colour = NULL, title = my_species) +
theme_light()
one_sp %>%
select(light, log_light) %>%
pivot_longer(everything(), names_to = "scale", values_to = "value") %>%
mutate(scale = recode(scale, light = "Raw radiance",
log_light = "log10(radiance + 1)")) %>%
ggplot(aes(x = value)) +
geom_histogram(bins = 40) +
facet_wrap(~ scale, scales = "free") +
labs(x = NULL, y = "Number of observations") +
theme_light()
Glossary (work in progress)
| Term | Plain meaning |
|---|---|
| Artificial light at night | Light from streetlights, buildings and signs that makes the night brighter than it would naturally be. |
| Radiance | A satellite measure of how much light is coming up from the ground at a location. |
| Onset | The time of a bird’s first call of the day, here measured in minutes from sunrise. |
| Random effect | A group, such as an individual bird or a site, that the model allows to differ from the others. |
| Citizen science | Research that uses data collected by volunteers. |
| Linear model | A straight line fitted through data to show how one thing changes as another changes. |
| Slope | How much the outcome changes for each one-unit increase in the predictor. |
| p-value | How surprising a result would be if there were really no effect. Small values suggest the effect is probably real. |
| R-squared (R²) | The share of the variation in the outcome that the model explains. 0 means none, 1 means all of it. |
| Residual | The gap between a data point and the line the model drew. |
| Pseudoreplication | Treating repeated measurements from the same place as if they were independent, which can make an effect look stronger than it is. |