Analysis

Author

Steffi LaZerte

Published

September 23, 2024

Setup

library(tidyverse)
library(gt)
library(DHARMa)
library(patchwork)
library(broom)

source("XX_functions.R")

# Metrics
v <- read_csv("Data/Datasets/vultures_final.csv") |>
  # Round non-integer values of population counts
  mutate(across(c(contains("pop"), contains("raw")), round)) 

# Raw counts
raw <- read_csv("Data/Datasets/vultures_clean_2023.csv")

# Predicted GAM models
pred <- read_csv("Data/Datasets/vultures_gams_pred.csv")

Set plotting defaults

y_breaks <- seq(200, 300, by = 5)

Using This Report

This report contains all the stats and analyses for the final data set.

Each section (e.g., “Timing of kettle formation”) contains the analysis for a particular type of data.

These analyses have several parts

  1. Descriptive statistics - min, max, median, etc.
  2. Model results - results of linear or general linear models, depending on the measures
  3. Figures - visualizing the models on top of the data
  4. Model Checks - DHARMa plots assessing model fit and assumptions
  5. Sensitivity - Checking influential observations
  6. Full Model Results - the raw summary output of the models, for completeness

Models

Most models are standard linear regressions.

Some, involving the count data (number of birds) are either Negative binomial or Poisson generalized linear models.

The Patterns of migration (skew and kurtosis) are linear models without predictors (intercept only), which effectively makes them t-tests (they produce the exact same estimates and statistics), but allow us to use DHARMA for model checking, etc.

Model Results

Model results are presented in a tabulated form, which is also saved to CSV as Data/Tables/table_MEASURE.csv.

They are also presented as the “Full Model Results”. This is the output of using summary(model) in R and is just a complete version of the output in case anything needs to be inspected at a later date. In short, I don’t think you’ll need it, but I’ve included it just in case.

The tabulated results contain both estimates and P values, as well as overall model R2 values for linear models (but not for general linear models, as they don’t really apply). I could also add confidence intervals as well if that’s useful.

R2 values are presented as regular and adjusted. The adjusted values take into account samples sizes and the number of predictors.

We can interpret the Estimates directly, especially in linear models.

For example, if there is a significant effect of year on mig_start_doy, and the estimate is 0.176 we can say:

“There was a significant effect of year on the start date of kettle formations, such that the start date increased by 0.176 days per year (t = XXX; P-value = XXX). Thus kettle formations started ~4.4 days later at the end of the study compared to the start.”

Considering that there are 25 years in this study.

For Negative Binomial or Poisson models, we will transform the estimates to “Incident Rate Ratios” by exponentiating the estimate.

For example, if there is a significant effect of year on mig_pop_max in a generalized linear model with the negative binomial family, and the estimate is 0.03625499, we can first calculate the exponential:

exp(0.03625499)
[1] 1.03692

Then we can say:

“There was a significant effect of year on the maximum number of vultures predicted in a kettle (Est = 0.0363; t = XXX; P-value = XXX), such that the number of birds increased by 3% per year.”

If our exponentiated result was 3, we would say a “300% increase”, or a “3-fold” increase.

These exponentiated estimates are included in the outputs as well as “Estimate Exp”

Variables

All variables are generally referred to by their name in the data set, e.g., mig_start_doy.

See the glossary of variables names in the Calculate Metrics > Data Details section.

Figures

You can click on any figure expand it out.

These figures are here to demostrate the patterns we’re seeing in the data and the stats. I imagine you would want different figures for publication, so just let me know what they should look like!

DHARMa plots

DHARMa is a package for simulating residuals to allow model checking for all types of models (details).

Generally these plots can be interpreted similarly to interpreting QQ Normal and Residual plots for linear models (although technically they’re using simulated residuals, not model residuals).

For context…

For example, we would interpret the following plot as showing decent model fit and assumptions.

The QQ plot on the left is nearly straight, there are no tests that flag problems in the distribution (KS test), Dispersion or Outliers. We also see no pattern in the residuals and no problems have been highlighted.

Similarly, I wouldn’t be unduly worried about the scary red line in this plot, as I don’t really find the residuals to show much pattern and there are no other problems highlighted, except that one outlier (the red asterisk).

This one starts to look a bit problematic…

This one is awful!

Sensitivity

The DHARMa model checks will usually catch outlier problems, but out of an abundance of caution, we’ll do a quick review of any high Cook’s D values and will rerun the models without certain data points to ensure they’re still acceptable.

These sections first calculate and show the Cook’s D for all observations in each model. In these tables, values highlighted in red are more than 4/25 (a standard Cook’s D cutoff 4/sample size). These indicate potential influential data points for a particular model.

Then I run a comparison of the model before and the after removing that point.

In every case here the pattern was either strengthened or weakened by the removal, but the overall pattern did not change (no Estimate signs changed, and significance changed from sig to non-sig or vice versa).

  • Strengthened patterns showed greater Estimates and smaller P values
  • Weakened patterns showed smaller Estimates and greater P values

Tables with coefficients for both original and new model are shown, with colour intensity showing strength in the pattern.

Timing of kettle formation

When do kettles form? Has this timing changed over the years?

  • Look for changes in the DOY of the 5%, 25%, 50%, 75%, 95%, of passage as well as the DOY with the highest predicted count.
  • mig_start_doy, peak_start_doy, p50_doy, peak_end_doy, mig_end_doy, max_doy

Descriptive stats

Code
v |> 
  select(contains("start_doy"), contains("50_doy"), contains("end_doy"), "max_doy") |>
  desc_stats()
measure mean sd min median max n
mig_start_doy 259.48 2.90 254 259 265 25
peak_start_doy 266.96 2.21 263 267 272 25
p50_doy 271.64 2.25 268 271 277 25
mig_end_doy 284.00 3.11 277 284 289 25
peak_end_doy 276.48 2.14 273 276 281 25
max_doy 271.60 2.80 266 271 278 25

Here we look at linear regressions of year by kettle timing dates.

m1 <- lm(mig_start_doy ~ year, data = v)
m2 <- lm(peak_start_doy ~ year, data = v)
m3 <- lm(p50_doy ~ year, data = v)
m4 <- lm(peak_end_doy ~ year, data = v)
m5 <- lm(mig_end_doy ~ year, data = v)
m6 <- lm(max_doy ~ year, data = v)

models <- list(m1, m2, m3, m4, m5, m6)

Tabulated Model Results

All show a significant increase in doy except mig_end_doy. Data saved to Data/Datasets/table_timing.csv

Code
t1 <- get_table(models)
write_csv(t1, d)
fmt_table(t1)
term estimate std error statistic p value n
mig_start_doy ~ year (R2 = 0.22; R2-adj = 0.19)
(Intercept) −94.618 138.088 −0.685 0.500 25
year 0.176 0.069 2.564 0.017 25
peak_start_doy ~ year (R2 = 0.3; R2-adj = 0.27)
(Intercept) −44.588 99.801 −0.447 0.659 25
year 0.155 0.050 3.122 0.005 25
p50_doy ~ year (R2 = 0.26; R2-adj = 0.23)
(Intercept) −24.819 104.600 −0.237 0.815 25
year 0.147 0.052 2.834 0.009 25
peak_end_doy ~ year (R2 = 0.23; R2-adj = 0.2)
(Intercept) 11.143 101.511 0.110 0.914 25
year 0.132 0.050 2.614 0.016 25
mig_end_doy ~ year (R2 = 0.05; R2-adj = 0.01)
(Intercept) 96.771 163.107 0.593 0.559 25
year 0.093 0.081 1.148 0.263 25
max_doy ~ year (R2 = 0.25; R2-adj = 0.22)
(Intercept) −92.317 130.519 −0.707 0.486 25
year 0.181 0.065 2.788 0.010 25

Here we look at a linear regression of year by kettle timing dates. The idea is to explore whether the interaction is significant, which would indicate that the rate of change in migration phenology differs among the seasonal points (start, start of peak, peak, end of peak, end).

We look at the effect of year and type of measurement (migration start/end, peak start/end and time of 50% passage) on date (day of year).

Code
v1 <- select(v, ends_with("doy"), year) |>
  pivot_longer(-year, names_to = "measure", values_to = "doy")

First we look at the interaction term to see if the slopes of change among measurements are significantly different from one another.

The type III ANOVA shows no significant interaction.

Code
m0 <- lm(doy ~ measure * year, data = v1)
car::Anova(m0, type = "III")
Anova Table (Type III tests)

Response: doy
             Sum Sq  Df F value  Pr(>F)   
(Intercept)    3.06   1  0.5445 0.46183   
measure        9.29   5  0.3309 0.89360   
year          47.49   1  8.4613 0.00423 **
measure:year   7.53   5  0.2683 0.92974   
Residuals    774.60 138                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

We’ll look at a post-hoc comparison of the slopes to interpret them individually.

This is nearly identical to running the models separately, above

As with the individual models, we can see that each measure shows a significant positive slope with the exception of the End of Migration (mig_end_day).

Code
e <- emmeans::emtrends(m0, ~ measure, var = "year") |> 
  emmeans::test()

All show a significant increase in doy except mig_end_doy. Data saved to Data/Datasets/table_timing_combined.csv

Code
t2 <- e |> 
  mutate(
    model = "doy ~ measure * year",
    n = nrow(v)) |>
  rename(
    estimate = year.trend, 
    std_error = SE, 
    
    statistic = t.ratio, 
    p_value = p.value)

write_csv(t2, d)

fmt_table(t2) |>
  tab_header(title = "Post-hoc results") |>
  tab_footnote("Number of years in the data.", locations = cells_column_labels(n)) |>
  tab_footnote("T statistic.", locations = cells_column_labels("statistic"))
Post-hoc results
measure estimate std error df statistic1 p value n2
doy ~ measure * year
max_doy 0.181 0.062 138 2.909 0.004 25
mig_end_doy 0.093 0.062 138 1.497 0.137 25
mig_start_doy 0.176 0.062 138 2.830 0.005 25
p50_doy 0.147 0.062 138 2.370 0.019 25
peak_end_doy 0.132 0.062 138 2.121 0.036 25
peak_start_doy 0.155 0.062 138 2.490 0.014 25
1 T statistic.
2 Number of years in the data.

Note: Lines without CI ribbons are not significant

Code
doy_figs <- v |>
  select(year, contains("doy")) |>
  pivot_longer(-year, names_to = "measure", values_to = "doy") |>
  left_join(mutate(t1, measure = str_trim(str_extract(model, "^[^~]*")), sig = p_value <= 0.05) |> 
              filter(term == "year") |>
              select(measure, sig), by = "measure")

ggplot(doy_figs, aes(x = year, y = doy, group = measure, colour = measure, fill = sig)) +
  theme_bw() +
  geom_point() +
  geom_smooth(method = "lm", se = TRUE) +
  scale_y_continuous(breaks = y_breaks) +
  scale_colour_viridis_d() +
  scale_fill_manual(values = c("TRUE" = "grey60", "FALSE" = "#FFFFFF00"), guide = "none") +
  labs(caption = "Lines without CI Error ribbons represent non-significant models")
`geom_smooth()` using formula = 'y ~ x'

The DHARMa model checks don’t highlight any particular outlier problems, but out of an abundance of caution, we’ll do a quick review of any high Cook’s D values.

Values highlighted in red are less than 4/25, a standard Cook’s D cutoff (4/sample size). These indicate potential influential data points for a particular model.

Code
get_cooks(models) |>
  gt_cooks(width = "80%")
Cook's Distances
year mig_start_doy peak_start_doy p50_doy peak_end_doy mig_end_doy max_doy
1998 0.05 0.03 0.04 0.04 0.00 0.05
1999 0.17 0.08 0.03 0.00 0.01 0.03
2000 0.14 0.11 0.08 0.09 0.03 0.01
2001 0.14 0.00 0.00 0.03 0.07 0.00
2002 0.01 0.04 0.08 0.00 0.02 0.15
2003 0.01 0.00 0.00 0.00 0.00 0.01
2004 0.08 0.05 0.06 0.06 0.09 0.02
2005 0.03 0.01 0.01 0.01 0.00 0.00
2006 0.01 0.01 0.00 0.01 0.06 0.02
2008 0.03 0.02 0.01 0.01 0.02 0.00
2009 0.01 0.00 0.03 0.07 0.11 0.02
2010 0.02 0.01 0.00 0.00 0.01 0.01
2011 0.04 0.10 0.10 0.07 0.01 0.14
2012 0.00 0.00 0.00 0.00 0.00 0.00
2013 0.02 0.02 0.02 0.07 0.02 0.03
2014 0.00 0.02 0.02 0.03 0.04 0.01
2015 0.01 0.00 0.01 0.01 0.07 0.03
2016 0.06 0.07 0.05 0.01 0.00 0.04
2017 0.00 0.01 0.02 0.06 0.09 0.00
2018 0.00 0.02 0.01 0.00 0.00 0.06
2019 0.01 0.02 0.05 0.04 0.04 0.04
2020 0.07 0.22 0.25 0.18 0.06 0.14
2021 0.12 0.13 0.09 0.16 0.12 0.07
2022 0.15 0.12 0.06 0.02 0.00 0.08
2023 0.10 0.02 0.01 0.00 0.01 0.00

Check influential points

Now let’s see what happens if we were to omit these years from the respective analyses.

Tables with coefficients for both original and new model are shown, with colour intensity showing strength in the pattern.

For the start of migration, if we omit 1999, the pattern is stronger, but the same.

Code
compare(m1, 1999)
model parameter Estimate Std. Error t value Pr(>|t|)
Original - mig_start_doy (Intercept) -94.618 138.088 -0.685 0.500
Drop 1999 - mig_start_doy (Intercept) -161.289 141.525 -1.140 0.267
Original - mig_start_doy year 0.176 0.069 2.564 0.017
Drop 1999 - mig_start_doy year 0.209 0.070 2.973 0.007

For the start of peak migration, the date of 50% passage and the end of peak migration, if we omit 2020, the pattern is weaker, but the same.

Code
compare(m2, 2020)
model parameter Estimate Std. Error t value Pr(>|t|)
Original - peak_start_doy (Intercept) -44.588 99.801 -0.447 0.659
Drop 2020 - peak_start_doy (Intercept) 7.053 95.767 0.074 0.942
Original - peak_start_doy year 0.155 0.050 3.122 0.005
Drop 2020 - peak_start_doy year 0.129 0.048 2.712 0.013
Code
compare(m3, 2020)
model parameter Estimate Std. Error t value Pr(>|t|)
Original - p50_doy (Intercept) -24.819 104.600 -0.237 0.815
Drop 2020 - p50_doy (Intercept) 32.436 99.102 0.327 0.747
Original - p50_doy year 0.147 0.052 2.834 0.009
Drop 2020 - p50_doy year 0.119 0.049 2.411 0.025
Code
compare(m4, 2020)
model parameter Estimate Std. Error t value Pr(>|t|)
Original - peak_end_doy (Intercept) 11.143 101.511 0.110 0.914
Drop 2020 - peak_end_doy (Intercept) 58.397 99.350 0.588 0.563
Original - peak_end_doy year 0.132 0.050 2.614 0.016
Drop 2020 - peak_end_doy year 0.108 0.049 2.193 0.039

Therefore I wouldn’t be concerned

In addition to the sensitivity check using cook’s distances, here we will conduct a more targeted exploration of some specific years.

Gaps around the peak

First, we consider removing 1998, 2011, 2013 as they have sampling gaps around the peak of migration

y <- c(1998, 2011, 2013)
m1a <- update(m1, data = filter(v, !year %in% y))
m2a <- update(m2, data = filter(v, !year %in% y))
m3a <- update(m3, data = filter(v, !year %in% y))
m4a <- update(m4, data = filter(v, !year %in% y))
m5a <- update(m5, data = filter(v, !year %in% y))
m6a <- update(m6, data = filter(v, !year %in% y))

models_a <- list(m1a, m2a, m3a, m4a, m5a, m6a)

Comparing the summary tables shows only superficial differences. If anything, it looks like removing these ‘gappy’ years only strengthens the existing patterns. Those which are significant are slightly more so, and the slopes are slightly greater.

Data saved to Data/Datasets/table_supplemental.csv

Original

Code
fmt_summary(models, intercept = FALSE)
model Parameter Estimate Std. Error T P
mig_start_doy ~ year year 0.176 0.069 2.564 0.017
peak_start_doy ~ year year 0.155 0.050 3.122 0.005
p50_doy ~ year year 0.147 0.052 2.834 0.009
peak_end_doy ~ year year 0.132 0.050 2.614 0.016
mig_end_doy ~ year year 0.093 0.081 1.148 0.263
max_doy ~ year year 0.181 0.065 2.788 0.010
Code
to <- get_table(models) |>
  select(model, term, estimate, std_error, statistic, p_value) |>
  mutate(type = "original")
write_csv(to, d)

Years removed

Code
fmt_summary(models_a, intercept = FALSE)
model Parameter Estimate Std. Error T P
mig_start_doy ~ year year 0.191 0.073 2.630 0.016
peak_start_doy ~ year year 0.163 0.049 3.355 0.003
p50_doy ~ year year 0.157 0.050 3.159 0.005
peak_end_doy ~ year year 0.140 0.048 2.924 0.008
mig_end_doy ~ year year 0.091 0.090 1.012 0.324
max_doy ~ year year 0.194 0.057 3.412 0.003
Code
ta <- get_table(models_a) |>
  select(model, term, estimate, std_error, statistic, p_value) |>
  mutate(type = "gaps_removed")
write_csv(ta, d, append = TRUE)

Missing the end of migration

Second, we consider removing 2001, 2006, 2023 as they were three years where it looks like sampling ended earlier than the end of migration.

y <- c(2001, 2006, 2023)
m1b <- update(m1, data = filter(v, !year %in% y))
m2b <- update(m2, data = filter(v, !year %in% y))
m3b <- update(m3, data = filter(v, !year %in% y))
m4b <- update(m4, data = filter(v, !year %in% y))
m5b <- update(m5, data = filter(v, !year %in% y))
m6b <- update(m6, data = filter(v, !year %in% y))

models_b <- list(m1b, m2b, m3b, m4b, m5b, m6b)

Removing these years from the analyses, reduced the strength of the patterns, but not by much. Further, there were few changes to the magnitude of the slopes.

Data saved to Data/Datasets/table_supplemental.csv

Original

Code
fmt_summary(models, intercept = FALSE)
model Parameter Estimate Std. Error T P
mig_start_doy ~ year year 0.176 0.069 2.564 0.017
peak_start_doy ~ year year 0.155 0.050 3.122 0.005
p50_doy ~ year year 0.147 0.052 2.834 0.009
peak_end_doy ~ year year 0.132 0.050 2.614 0.016
mig_end_doy ~ year year 0.093 0.081 1.148 0.263
max_doy ~ year year 0.181 0.065 2.788 0.010

Years removed

Code
fmt_summary(models_b, intercept = FALSE)
model Parameter Estimate Std. Error T P
mig_start_doy ~ year year 0.166 0.074 2.237 0.037
peak_start_doy ~ year year 0.155 0.058 2.680 0.014
p50_doy ~ year year 0.151 0.062 2.442 0.024
peak_end_doy ~ year year 0.128 0.059 2.166 0.043
mig_end_doy ~ year year 0.077 0.089 0.856 0.402
max_doy ~ year year 0.173 0.076 2.277 0.034
Code
tb <- get_table(models_b) |>
  select(model, term, estimate, std_error, statistic, p_value) |>
  mutate(type = "end_removed")
write_csv(tb, d, append = TRUE)

Remove all problematic years

Out of curiosity, let’s try removing all 6 years. This reduces our data from 25 to 19 samples, so is a substantial change and we may see non-significant results merely because we have removed too many samples.

y <- c(1998, 2011, 2013, 2001, 2006, 2023)
m1c <- update(m1, data = filter(v, !year %in% y))
m2c <- update(m2, data = filter(v, !year %in% y))
m3c <- update(m3, data = filter(v, !year %in% y))
m4c <- update(m4, data = filter(v, !year %in% y))
m5c <- update(m5, data = filter(v, !year %in% y))
m6c <- update(m6, data = filter(v, !year %in% y))

models_c <- list(m1c, m2c, m3c, m4c, m5c, m6c)

Removing these years from the analyses, reduced the strength of the patterns, but not by much. Further, there are few changes to the magnitude of the slopes.

Data saved to Data/Datasets/table_supplemental.csv

Original

Code
fmt_summary(models, intercept = FALSE)
model Parameter Estimate Std. Error T P
mig_start_doy ~ year year 0.176 0.069 2.564 0.017
peak_start_doy ~ year year 0.155 0.050 3.122 0.005
p50_doy ~ year year 0.147 0.052 2.834 0.009
peak_end_doy ~ year year 0.132 0.050 2.614 0.016
mig_end_doy ~ year year 0.093 0.081 1.148 0.263
max_doy ~ year year 0.181 0.065 2.788 0.010

Years removed

Code
fmt_summary(models_c, intercept = FALSE)
model Parameter Estimate Std. Error T P
mig_start_doy ~ year year 0.180 0.081 2.213 0.041
peak_start_doy ~ year year 0.165 0.059 2.818 0.012
p50_doy ~ year year 0.165 0.061 2.715 0.015
peak_end_doy ~ year year 0.138 0.057 2.428 0.027
mig_end_doy ~ year year 0.073 0.101 0.722 0.480
max_doy ~ year year 0.190 0.069 2.773 0.013
Code
tc <- get_table(models_c) |>
  select(model, term, estimate, std_error, statistic, p_value) |>
  mutate(type = "all_removed")
write_csv(tc, d, append = TRUE)

Conclusions

I do not believe that these potentially problematic sampling years are driving the patterns we see in this data.

At worst, by removing these years we increase the error around a slope. But the magnitudes of the slopes do not change and even the significance is only ever reduced only a little.

Code
models <- append(list(m0), models)
map(models, summary)
[[1]]

Call:
lm(formula = doy ~ measure * year, data = v1)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.8473 -1.6117 -0.1487  1.6871  6.3348 

Coefficients:
                             Estimate Std. Error t value Pr(>|t|)   
(Intercept)                -92.316742 125.108782  -0.738  0.46183   
measuremig_end_doy         189.088235 176.930536   1.069  0.28706   
measuremig_start_doy        -2.300905 176.930536  -0.013  0.98964   
measurep50_doy              67.497738 176.930536   0.381  0.70342   
measurepeak_end_doy        103.459276 176.930536   0.585  0.55967   
measurepeak_start_doy       47.728507 176.930536   0.270  0.78775   
year                         0.180995   0.062223   2.909  0.00423 **
measuremig_end_doy:year     -0.087877   0.087996  -0.999  0.31972   
measuremig_start_doy:year   -0.004884   0.087996  -0.055  0.95582   
measurep50_doy:year         -0.033550   0.087996  -0.381  0.70359   
measurepeak_end_doy:year    -0.049029   0.087996  -0.557  0.57832   
measurepeak_start_doy:year  -0.026046   0.087996  -0.296  0.76769   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.369 on 138 degrees of freedom
Multiple R-squared:  0.9195,    Adjusted R-squared:  0.9131 
F-statistic: 143.3 on 11 and 138 DF,  p-value: < 2.2e-16


[[2]]

Call:
lm(formula = mig_start_doy ~ year, data = v)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.7823 -2.4867 -0.0717  2.1044  3.6894 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept) -94.61765  138.08785  -0.685   0.5001  
year          0.17611    0.06868   2.564   0.0173 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.615 on 23 degrees of freedom
Multiple R-squared:  0.2223,    Adjusted R-squared:  0.1885 
F-statistic: 6.576 on 1 and 23 DF,  p-value: 0.01734


[[3]]

Call:
lm(formula = peak_start_doy ~ year, data = v)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.7905 -1.2410 -0.6356  1.5194  3.9842 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)   
(Intercept) -44.58824   99.80097  -0.447  0.65922   
year          0.15495    0.04964   3.122  0.00479 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.89 on 23 degrees of freedom
Multiple R-squared:  0.2976,    Adjusted R-squared:  0.2671 
F-statistic: 9.745 on 1 and 23 DF,  p-value: 0.004794


[[4]]

Call:
lm(formula = p50_doy ~ year, data = v)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.4303 -1.5778 -0.2186  1.2237  4.3069 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)   
(Intercept) -24.81900  104.59963  -0.237   0.8145   
year          0.14745    0.05202   2.834   0.0094 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.981 on 23 degrees of freedom
Multiple R-squared:  0.2589,    Adjusted R-squared:  0.2266 
F-statistic: 8.033 on 1 and 23 DF,  p-value: 0.009399


[[5]]

Call:
lm(formula = peak_end_doy ~ year, data = v)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.2636 -1.1873 -0.3398  1.1323  3.4725 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept)  11.14253  101.51140   0.110   0.9135  
year          0.13197    0.05049   2.614   0.0155 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.922 on 23 degrees of freedom
Multiple R-squared:  0.229, Adjusted R-squared:  0.1955 
F-statistic: 6.832 on 1 and 23 DF,  p-value: 0.01552


[[6]]

Call:
lm(formula = mig_end_doy ~ year, data = v)

Residuals:
   Min     1Q Median     3Q    Max 
-6.847 -2.009  0.177  1.966  4.618 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept)  96.77149  163.10693   0.593    0.559
year          0.09312    0.08112   1.148    0.263

Residual standard error: 3.089 on 23 degrees of freedom
Multiple R-squared:  0.05419,   Adjusted R-squared:  0.01306 
F-statistic: 1.318 on 1 and 23 DF,  p-value: 0.2628


[[7]]

Call:
lm(formula = max_doy ~ year, data = v)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.0362 -2.1131 -0.1222  1.6018  6.3348 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept) -92.31674  130.51911  -0.707   0.4865  
year          0.18100    0.06491   2.788   0.0104 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.472 on 23 degrees of freedom
Multiple R-squared:  0.2526,    Adjusted R-squared:  0.2201 
F-statistic: 7.774 on 1 and 23 DF,  p-value: 0.01045
Code
model_check_figs(models)

Duration of migration

How long is migration? Has it changed in length?

  • Look for changes in the number of days over which migration and peak migration occur
  • mig_dur_days, peak_dur_days

Descriptive stats

Code
v |> 
  select(contains("dur")) |>
  desc_stats()
measure mean sd min median max n
mig_dur_days 24.52 3.62 16 25 31 25
peak_dur_days 9.52 1.19 7 10 12 25
m1 <- lm(mig_dur_days ~ year, data = v)
m2 <- lm(peak_dur_days ~ year, data = v)

models <- list(m1, m2)

Tabulated Results

No significant results

Code
t <- get_table(models)
write_csv(t, "Data/Datasets/table_duration.csv")
fmt_table(t)
term estimate std error statistic p value n
mig_dur_days ~ year (R2 = 0.03; R2-adj = -0.01)
(Intercept) 191.389 192.063 0.996 0.329 25
year −0.083 0.096 −0.869 0.394 25
peak_dur_days ~ year (R2 = 0.02; R2-adj = -0.02)
(Intercept) 55.731 63.706 0.875 0.391 25
year −0.023 0.032 −0.725 0.476 25

Note: Lines without Std Error ribbons are not significant

Code
dur_figs <- v |>
  select(year, contains("dur")) |>
  pivot_longer(-year, names_to = "measure", values_to = "dur")

ggplot(dur_figs, aes(x = year, y = dur, group = measure, colour = measure)) +
  theme_bw() +
  geom_point() +
  stat_smooth(method = "lm", se = FALSE) +
  #scale_y_continuous(breaks = y_breaks) +
  scale_colour_viridis_d()+
  labs(caption = "Lines without Std Error ribbons represent non-significant models")
`geom_smooth()` using formula = 'y ~ x'

Code
model_check_figs(models)

The DHARMa model checks don’t highlight any particular outlier problems, but out of an abundance of caution, we’ll do a quick review of any high Cook’s D values.

Values highlighted in red are less than 4/25, a standard Cook’s D cutoff (4/sample size). These indicate potential influential data points for a particular model.

Code
get_cooks(models) |>
  gt_cooks()
Cook's Distances
year mig_dur_days peak_dur_days
1998 0.02 0.00
1999 0.14 0.20
2000 0.01 0.00
2001 0.00 0.02
2002 0.03 0.06
2003 0.00 0.00
2004 0.00 0.00
2005 0.01 0.00
2006 0.09 0.12
2008 0.00 0.00
2009 0.13 0.10
2010 0.03 0.04
2011 0.00 0.00
2012 0.00 0.00
2013 0.00 0.04
2014 0.03 0.01
2015 0.08 0.00
2016 0.04 0.06
2017 0.07 0.05
2018 0.00 0.08
2019 0.06 0.00
2020 0.00 0.00
2021 0.00 0.00
2022 0.08 0.09
2023 0.09 0.04

Check influential points

Now let’s see what happens if we were to omit these years from the respective analyses.

Tables with coefficients for both original and new model are shown, with colour intensity showing strength in the pattern.

For the number of days in peak migration, if we omit 1999, the pattern is stronger, but the same (i.e. not significant).

Code
compare(m2, 1999)
model parameter Estimate Std. Error t value Pr(>|t|)
Original - peak_dur_days (Intercept) 55.731 63.706 0.875 0.391
Drop 1999 - peak_dur_days (Intercept) 89.114 64.690 1.378 0.182
Original - peak_dur_days year -0.023 0.032 -0.725 0.476
Drop 1999 - peak_dur_days year -0.040 0.032 -1.229 0.232

Therefore I wouldn’t be concerned

In addition to the sensitivity check using cook’s distances, here we will conduct a more targeted exploration of some specific years.

Gaps around the peak

First, we consider removing 1998, 2011, 2013 as they have sampling gaps around the peak of migration

y <- c(1998, 2011, 2013)
m1a <- update(m1, data = filter(v, !year %in% y))
m2a <- update(m2, data = filter(v, !year %in% y))
models_a <- list(m1a, m2a)

Comparing the summary tables shows only superficial differences.

Data saved to Data/Datasets/table_supplemental.csv

Original

fmt_summary(models, intercept = FALSE)
to <- get_table(models) |>
  select(model, term, estimate, std_error, statistic, p_value) |>
  mutate(type = "original")
write_csv(to, d, append = TRUE)
model Parameter Estimate Std. Error T P
mig_dur_days ~ year year −0.083 0.096 −0.869 0.394
peak_dur_days ~ year year −0.023 0.032 −0.725 0.476

Years removed

fmt_summary(models_a, intercept = FALSE)
ta <- get_table(models_a) |>
  select(model, term, estimate, std_error, statistic, p_value) |>
  mutate(type = "gaps_removed")
write_csv(ta, d, append = TRUE)
model Parameter Estimate Std. Error T P
mig_dur_days ~ year year −0.100 0.108 −0.927 0.365
peak_dur_days ~ year year −0.023 0.035 −0.674 0.508

Missing the end of migration

Second, we consider removing 2001, 2006, 2023 as they were three years where it looks like sampling ended earlier than the end of migration.

y <- c(2001, 2006, 2023)
m1b <- update(m1, data = filter(v, !year %in% y))
m2b <- update(m2, data = filter(v, !year %in% y))
models_b <- list(m1b, m2b)

Comparing the summary tables shows only superficial differences.

Data saved to Data/Datasets/table_supplemental.csv

Original

fmt_summary(models, intercept = FALSE)
model Parameter Estimate Std. Error T P
mig_dur_days ~ year year −0.083 0.096 −0.869 0.394
peak_dur_days ~ year year −0.023 0.032 −0.725 0.476

Years removed

fmt_summary(models_b, intercept = FALSE)
tb <- get_table(models_b) |>
  select(model, term, estimate, std_error, statistic, p_value) |>
  mutate(type = "end_removed")
write_csv(tb, d, append = TRUE)
model Parameter Estimate Std. Error T P
mig_dur_days ~ year year −0.089 0.103 −0.871 0.394
peak_dur_days ~ year year −0.028 0.033 −0.823 0.420

Remove all problematic years

Out of curiosity, let’s try removing all 6 years. This reduces our data from 25 to 19 samples, so is a substantial change and we may see non-significant results merely because we have removed too many samples.

y <- c(1998, 2011, 2013, 2001, 2006, 2023)
m1c <- update(m1, data = filter(v, !year %in% y))
m2c <- update(m2, data = filter(v, !year %in% y))
models_c <- list(m1c, m2c)

Comparing the summary tables shows only superficial differences.

Data saved to Data/Datasets/table_supplemental.csv

Original

fmt_summary(models, intercept = FALSE)
model Parameter Estimate Std. Error T P
mig_dur_days ~ year year −0.083 0.096 −0.869 0.394
peak_dur_days ~ year year −0.023 0.032 −0.725 0.476

Years removed

fmt_summary(models_c, intercept = FALSE)
tc <- get_table(models_c) |>
  select(model, term, estimate, std_error, statistic, p_value) |>
  mutate(type = "all_removed")
write_csv(tc, d, append = TRUE)
model Parameter Estimate Std. Error T P
mig_dur_days ~ year year −0.107 0.120 −0.897 0.382
peak_dur_days ~ year year −0.027 0.037 −0.741 0.469

Conclusions

I do not believe that these potentially problematic sampling years are driving the patterns we see in this data.

Code
map(models, summary)
[[1]]

Call:
lm(formula = mig_dur_days ~ year, data = v)

Residuals:
    Min      1Q  Median      3Q     Max 
-8.6561 -1.5690  0.5929  2.0119  6.0949 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) 191.38914  192.06304   0.996    0.329
year         -0.08299    0.09552  -0.869    0.394

Residual standard error: 3.637 on 23 degrees of freedom
Multiple R-squared:  0.03178,   Adjusted R-squared:  -0.01032 
F-statistic: 0.7549 on 1 and 23 DF,  p-value: 0.3939


[[2]]

Call:
lm(formula = peak_dur_days ~ year, data = v)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.5577 -0.5117  0.1895  0.5572  2.3734 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) 55.73077   63.70604   0.875    0.391
year        -0.02298    0.03168  -0.725    0.476

Residual standard error: 1.206 on 23 degrees of freedom
Multiple R-squared:  0.02237,   Adjusted R-squared:  -0.02014 
F-statistic: 0.5262 on 1 and 23 DF,  p-value: 0.4755

Number of birds in kettles

What is the maximum kettle count? Has this changed?

  • Look for changes in the maximum predicted count and in the maximum observed count
  • variables mig_pop_max, mig_raw_max
  • In both cases, these values represent the maximum population counts for the migration period (25%-75%), but also correspond the the maximum population count overall.

Descriptive stats

Code
v |> 
  select("mig_pop_max", "mig_raw_max") |>
  desc_stats()
measure mean sd min median max n
mig_pop_max 428.68 236.76 146 342 1304 25
mig_raw_max 797.76 468.45 310 650 2375 25

Here we use a Negative Binomial Generalized Linear Model.

  • The linear regression was a bit iffy (especially for the raw counts which are more volatile), although a log transformation helped.
  • However, since the data is technically count data (i.e. kettle counts, predicted and raw), theoretically is it Poisson distributed data.
  • So I tried Poisson models, but they were overdispersed
  • Therefore, I use a Negative Binomial distribution to model the dispersion directly
  • This ends up with very similar results to log-transformed linear regression, but is a more appropriate given the type of data.

However we can easily present log10 models if you think that’s simpler to present.

  • I also ran a third NB model (m2b), which omitted some of the really large raw kettle counts (see the Figure), because these really large counts resulted in some problems in the model residual figure (see Model Checks)
  • Removing the counts resulted in better model fit, and didn’t really change the interpretations (see the results under mig_raw_max where the sample size is only 23.
  • Therefore, I think I wouldn’t worry about the model checks and would keep model m2 (the results saved to file only include the first two models).
m1 <- MASS::glm.nb(mig_pop_max ~ year, data = v)
m2 <- MASS::glm.nb(mig_raw_max ~ year, data = v)
m2b <- MASS::glm.nb(mig_raw_max ~ year, data = filter(v, mig_raw_max < 1500))

models <- list(m1, m2, m2b)

Not used, but included for completeness

x1 <- lm(log10(mig_pop_max) ~ year, data = v)
x2 <- lm(log10(mig_raw_max) ~ year, data = v)

y1 <- glm(mig_pop_max ~ year, family = "poisson", data = v)
y2 <- glm(mig_raw_max ~ year, family = "poisson", data = v)

s <- simulateResiduals(x1, plot = TRUE)
title("log10 - mig_pop_max ~ year")

s <- simulateResiduals(x2, plot = TRUE)
title("log10 - mig_raw_max ~ year")

s <- simulateResiduals(y1, plot = TRUE)
title("Poisson - mig_pop_max ~ year")

s <- simulateResiduals(y2, plot = TRUE)
title("Poisson - mig_raw_max ~ year")

Tabulated Results

All the measures are significant. There was an increase in the maximum number of birds in kettles in both the raw counts as well as the predicted max.

Code
t <- get_table(models)
write_csv(filter(t, n != 23), "Data/Datasets/table_pop_mig.csv")
fmt_table(t)
term estimate estimate exp std error statistic p value n
mig_pop_max ~ year
(Intercept) −66.877 0.000 20.084 −3.330 0.001 25
year 0.036 1.037 0.010 3.630 0.000 25
mig_raw_max ~ year
(Intercept) −75.414 0.000 20.048 −3.762 0.000 25
year 0.041 1.042 0.010 4.092 0.000 25
(Intercept) −45.033 0.000 19.928 −2.260 0.024 23
year 0.026 1.026 0.010 2.587 0.010 23
Code
pop_figs <- v |>
  select(year, mig_pop_max, mig_raw_max) |>
  pivot_longer(-year, names_to = "measure", values_to = "pop")

ggplot(pop_figs, aes(x = year, y = pop, group = measure, colour = measure)) +
  theme_bw() +
  geom_point() +
  stat_smooth(method = MASS::glm.nb) +
  scale_colour_viridis_d(end = 0.8)
`geom_smooth()` using formula = 'y ~ x'

Code
model_check_figs(models)
qu = 0.25, log(sigma) = -3.825783 : outer Newton did not converge fully.

The DHARMa model checks don’t highlight any particular outlier problems, but out of an abundance of caution, we’ll do a quick review of any high Cook’s D values.

Values highlighted in red are less than 4/25, a standard Cook’s D cutoff (4/sample size). These indicate potential influential data points for a particular model.

Code
get_cooks(models[-3]) |>
  gt_cooks()
Cook's Distances
year mig_pop_max mig_raw_max
1998 0.07 0.01
1999 0.96 1.10
2000 0.12 0.04
2001 0.04 0.07
2002 0.01 0.00
2003 0.00 0.00
2004 0.14 0.03
2005 0.01 0.02
2006 0.00 0.02
2008 0.02 0.00
2009 0.00 0.00
2010 0.04 0.10
2011 0.00 0.00
2012 0.01 0.03
2013 0.01 0.01
2014 0.01 0.00
2015 0.02 0.00
2016 0.00 0.02
2017 0.00 0.04
2018 0.05 0.04
2019 0.04 0.02
2020 0.02 0.00
2021 0.00 0.00
2022 0.72 0.57
2023 0.00 0.09

Check influential points

Now let’s see what happens if we were to omit these years from the respective analyses.

Tables with coefficients for both original and new model are shown, with colour intensity showing strength in the pattern.

For the maximum population recorded during migration based on both predicted counts and raw counts, if we omit 1999, the patterns are stronger, but the same.

Code
compare(m1, 1999)
model parameter Estimate Std. Error z value Pr(>|z|)
Original - mig_pop_max (Intercept) -66.877 20.084 -3.330 0.001
Drop 1999 - mig_pop_max (Intercept) -91.053 18.059 -5.042 0.000
Original - mig_pop_max year 0.036 0.010 3.630 0.000
Drop 1999 - mig_pop_max year 0.048 0.009 5.373 0.000
Code
compare(m2, 1999)
model parameter Estimate Std. Error z value Pr(>|z|)
Original - mig_raw_max (Intercept) -75.414 20.048 -3.762 0
Drop 1999 - mig_raw_max (Intercept) -100.459 17.545 -5.726 0
Original - mig_raw_max year 0.041 0.010 4.092 0
Drop 1999 - mig_raw_max year 0.053 0.009 6.101 0

If we omit 2022, the patterns are weaker, but still the same.

Code
compare(m1, 2022)
model parameter Estimate Std. Error z value Pr(>|z|)
Original - mig_pop_max (Intercept) -66.877 20.084 -3.330 0.001
Drop 2022 - mig_pop_max (Intercept) -46.746 18.695 -2.500 0.012
Original - mig_pop_max year 0.036 0.010 3.630 0.000
Drop 2022 - mig_pop_max year 0.026 0.009 2.819 0.005
Code
compare(m2, 2022)
model parameter Estimate Std. Error z value Pr(>|z|)
Original - mig_raw_max (Intercept) -75.414 20.048 -3.762 0.000
Drop 2022 - mig_raw_max (Intercept) -58.542 19.243 -3.042 0.002
Original - mig_raw_max year 0.041 0.010 4.092 0.000
Drop 2022 - mig_raw_max year 0.032 0.010 3.383 0.001

Therefore I wouldn’t be concerned

In addition to the sensitivity check using cook’s distances, here we will conduct a more targeted exploration of some specific years.

Gaps around the peak

First, we consider removing 1998, 2011, 2013 as they have sampling gaps around the peak of migration

y <- c(1998, 2011, 2013)
m1a <- update(m1, data = filter(v, !year %in% y))
m2a <- update(m2, data = filter(v, !year %in% y))
models_a <- list(m1a, m2a)

Comparing the summary tables shows only superficial differences.

Data saved to Data/Datasets/table_supplemental.csv

Original

fmt_summary(models[1:2], intercept = FALSE)
to <- get_table(models[1:2]) |>
  select(model, term, estimate, std_error, statistic, p_value) |>
  mutate(type = "original")
write_csv(to, d, append = TRUE)
model Parameter Estimate Std. Error Z P
mig_pop_max ~ year year 0.036 0.010 3.630 0.000
mig_raw_max ~ year year 0.041 0.010 4.092 0.000

Years removed

fmt_summary(models_a, intercept = FALSE)
ta <- get_table(models_a) |>
  select(model, term, estimate, std_error, statistic, p_value) |>
  mutate(type = "gaps_removed")
write_csv(ta, d, append = TRUE)
model Parameter Estimate Std. Error Z P
mig_pop_max ~ year year 0.034 0.011 3.111 0.002
mig_raw_max ~ year year 0.040 0.011 3.641 0.000

Missing the end of migration

Second, we consider removing 2001, 2006, 2023 as they were three years where it looks like sampling ended earlier than the end of migration.

y <- c(2001, 2006, 2023)
m1b <- update(m1, data = filter(v, !year %in% y))
m2b <- update(m2, data = filter(v, !year %in% y))
models_b <- list(m1b, m2b)

Comparing the summary tables shows only superficial differences.

Data saved to Data/Datasets/table_supplemental.csv

Original

fmt_summary(models[1:2], intercept = FALSE)
model Parameter Estimate Std. Error Z P
mig_pop_max ~ year year 0.036 0.010 3.630 0.000
mig_raw_max ~ year year 0.041 0.010 4.092 0.000

Years removed

fmt_summary(models_b, intercept = FALSE)
tb <- get_table(models_b) |>
  select(model, term, estimate, std_error, statistic, p_value) |>
  mutate(type = "end_removed")
write_csv(tb, d, append = TRUE)
model Parameter Estimate Std. Error Z P
mig_pop_max ~ year year 0.034 0.012 2.934 0.003
mig_raw_max ~ year year 0.034 0.011 3.079 0.002

Remove all problematic years

Out of curiosity, let’s try removing all 6 years. This reduces our data from 25 to 19 samples, so is a substantial change and we may see non-significant results merely because we have removed too many samples.

y <- c(1998, 2011, 2013, 2001, 2006, 2023)
m1c <- update(m1, data = filter(v, !year %in% y))
m2c <- update(m2, data = filter(v, !year %in% y))
models_c <- list(m1c, m2c)

Comparing the summary tables shows only superficial differences.

Data saved to Data/Datasets/table_supplemental.csv

Original

fmt_summary(models[1:2], intercept = FALSE)
model Parameter Estimate Std. Error Z P
mig_pop_max ~ year year 0.036 0.010 3.630 0.000
mig_raw_max ~ year year 0.041 0.010 4.092 0.000

Years removed

fmt_summary(models_c, intercept = FALSE)
tc <- get_table(models_c) |>
  select(model, term, estimate, std_error, statistic, p_value) |>
  mutate(type = "all_removed")
write_csv(tc, d, append = TRUE)
model Parameter Estimate Std. Error Z P
mig_pop_max ~ year year 0.031 0.013 2.389 0.017
mig_raw_max ~ year year 0.032 0.012 2.595 0.009

Conclusions

I do not believe that these potentially problematic sampling years are driving the patterns we see in this data.

Code
map(models, summary)
[[1]]

Call:
MASS::glm.nb(formula = mig_pop_max ~ year, data = v, init.theta = 7.041467497, 
    link = log)

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept) -66.876552  20.083997   -3.33 0.000869 ***
year          0.036255   0.009989    3.63 0.000284 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Negative Binomial(7.0415) family taken to be 1)

    Null deviance: 39.927  on 24  degrees of freedom
Residual deviance: 25.555  on 23  degrees of freedom
AIC: 327.07

Number of Fisher Scoring iterations: 1

              Theta:  7.04 
          Std. Err.:  1.98 

 2 x log-likelihood:  -321.07 

[[2]]

Call:
MASS::glm.nb(formula = mig_raw_max ~ year, data = v, init.theta = 7.009118042, 
    link = log)

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept) -75.413580  20.048025  -3.762 0.000169 ***
year          0.040803   0.009971   4.092 4.27e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Negative Binomial(7.0091) family taken to be 1)

    Null deviance: 44.412  on 24  degrees of freedom
Residual deviance: 25.564  on 23  degrees of freedom
AIC: 357.4

Number of Fisher Scoring iterations: 1

              Theta:  7.01 
          Std. Err.:  1.96 

 2 x log-likelihood:  -351.398 

[[3]]

Call:
MASS::glm.nb(formula = mig_raw_max ~ year, data = filter(v, mig_raw_max < 
    1500), init.theta = 9.016076348, link = log)

Coefficients:
              Estimate Std. Error z value Pr(>|z|)   
(Intercept) -45.033115  19.928099  -2.260  0.02383 * 
year          0.025653   0.009916   2.587  0.00968 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Negative Binomial(9.0161) family taken to be 1)

    Null deviance: 30.328  on 22  degrees of freedom
Residual deviance: 23.409  on 21  degrees of freedom
AIC: 319.15

Number of Fisher Scoring iterations: 1

              Theta:  9.02 
          Std. Err.:  2.65 

 2 x log-likelihood:  -313.148 

Number of residents

How many resident vultures are there? Has this changed?

  • Look for changes in the median number of residents
  • res_pop_median, res_raw_median

Descriptive stats

Code
v |> 
  select("res_pop_median", "res_raw_median") |>
  desc_stats()
measure mean sd min median max n
res_pop_median 5.56 2.71 2 5 12 25
res_raw_median 4.84 2.17 1 4 8 25

Here we use a Poisson Generalized Linear Model.

  • This yields similar results to a linear regression with log transformed values.
  • However, since the data is technically count data (i.e. kettle counts, predicted and raw), theoretically is it Poisson distributed data.
  • It is not overdispersed so we stick with Poisson (rather than Negative Binomial as we did above)
m1 <- glm(res_pop_median ~ year, family = "poisson", data = v)
m2 <- glm(res_raw_median ~ year, family = "poisson", data = v)

models <- list(m1, m2)

Not used, but included for completeness

x1 <- lm(log10(res_pop_median) ~ year, data = v)
x2 <- lm(log10(res_raw_median) ~ year, data = v)

s <- simulateResiduals(x1, plot = TRUE)
title("log10 - res_pop_median ~ year")

s <- simulateResiduals(x2, plot = TRUE)
title("log10 - res_raw_median ~ year")

Tabulated Results

All measures are significant, there is an increase in the number of resident birds over the years.

Code
t <- get_table(models)
write_csv(t, "Data/Datasets/table_pop_res.csv")
fmt_table(t)
term estimate estimate exp std error statistic p value n
res_pop_median ~ year
(Intercept) −77.157 0.000 23.209 −3.324 0.001 25
year 0.039 1.040 0.012 3.400 0.001 25
res_raw_median ~ year
(Intercept) −65.838 0.000 24.667 −2.669 0.008 25
year 0.034 1.034 0.012 2.734 0.006 25
Code
res_figs <- v |>
  select(year, res_pop_median, res_raw_median) |>
  pivot_longer(-year, names_to = "measure", values_to = "pop")

ggplot(res_figs, aes(x = year, y = pop, group = measure, colour = measure)) +
  theme_bw() +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "poisson")) +
  scale_colour_viridis_d(end = 0.8)
`geom_smooth()` using formula = 'y ~ x'

Code
model_check_figs(models)

The DHARMa model checks don’t highlight any particular outlier problems, but out of an abundance of caution, we’ll do a quick review of any high Cook’s D values.

Values highlighted in red are less than 4/25, a standard Cook’s D cutoff (4/sample size). These indicate potential influential data points for a particular model.

Code
get_cooks(models) |>
  gt_cooks()
Cook's Distances
year res_pop_median res_raw_median
1998 0.00 0.03
1999 0.05 0.08
2000 0.04 0.03
2001 0.30 0.37
2002 0.00 0.00
2003 0.01 0.00
2004 0.00 0.00
2005 0.00 0.00
2006 0.02 0.03
2008 0.04 0.07
2009 0.02 0.01
2010 0.00 0.00
2011 0.01 0.00
2012 0.01 0.00
2013 0.00 0.00
2014 0.04 0.02
2015 0.01 0.01
2016 0.03 0.03
2017 0.02 0.00
2018 0.10 0.03
2019 0.00 0.03
2020 0.16 0.03
2021 0.01 0.02
2022 0.02 0.01
2023 0.04 0.02

Check influential points

Now let’s see what happens if we were to omit these years from the respective analyses.

Tables with coefficients for both original and new model are shown, with colour intensity showing strength in the pattern.

For the daily median number of resident birds based on both predicted counts and raw counts, if we omit 2001, the patterns are stronger, but the same.

Code
compare(m1, 2001)
model parameter Estimate Std. Error z value Pr(>|z|)
Original - res_pop_median (Intercept) -77.157 23.209 -3.324 0.001
Drop 2001 - res_pop_median (Intercept) -93.177 24.827 -3.753 0.000
Original - res_pop_median year 0.039 0.012 3.400 0.001
Drop 2001 - res_pop_median year 0.047 0.012 3.824 0.000
Code
compare(m2, 2001)
model parameter Estimate Std. Error z value Pr(>|z|)
Original - res_raw_median (Intercept) -65.838 24.667 -2.669 0.008
Drop 2001 - res_raw_median (Intercept) -84.681 26.515 -3.194 0.001
Original - res_raw_median year 0.034 0.012 2.734 0.006
Drop 2001 - res_raw_median year 0.043 0.013 3.254 0.001

If we omit 2020, the patterns are weaker, but still the same.

Code
compare(m1, 2020)
model parameter Estimate Std. Error z value Pr(>|z|)
Original - res_pop_median (Intercept) -77.157 23.209 -3.324 0.001
Drop 2020 - res_pop_median (Intercept) -67.678 24.143 -2.803 0.005
Original - res_pop_median year 0.039 0.012 3.400 0.001
Drop 2020 - res_pop_median year 0.034 0.012 2.874 0.004

Therefore I wouldn’t be concerned

Code
map(models, summary)
[[1]]

Call:
glm(formula = res_pop_median ~ year, family = "poisson", data = v)

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -77.15702   23.20927  -3.324 0.000886 ***
year          0.03921    0.01153   3.400 0.000673 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 30.326  on 24  degrees of freedom
Residual deviance: 18.389  on 23  degrees of freedom
AIC: 109.34

Number of Fisher Scoring iterations: 4


[[2]]

Call:
glm(formula = res_raw_median ~ year, family = "poisson", data = v)

Coefficients:
             Estimate Std. Error z value Pr(>|z|)   
(Intercept) -65.83767   24.66741  -2.669  0.00761 **
year          0.03351    0.01226   2.734  0.00625 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 24.526  on 24  degrees of freedom
Residual deviance: 16.868  on 23  degrees of freedom
AIC: 104.42

Number of Fisher Scoring iterations: 4

Patterns of migration

Is the timing of migration skewed to earlier or later in the season? Is this distribution of migration flattened and wide (low kurtosis) or peaked (high kurtosis) Are these patterns consistent over the years?

Kurtosis can be used to indicates if all the birds pass through in a relatively quick ‘clump’ (thick-tailed, lower kurtosis), or whether the migration stretches out over a longer period of time (long-tailed, higher kurtosis). Possible implications for conservation or future survey designs?

  • Look for significant skew and kurtosis in migration (5-95%) and peak migration (25-75%)
  • mig_skew, peak_skew, mig_kurt, peak_kurt

Descriptive stats

Code
v |> 
  select(contains("skew"), contains("kurt")) |>
  desc_stats()
measure mean sd min median max n
mig_skew 0.01 0.21 -0.4117554 -0.0278960097 0.5067799 25
peak_skew 0.01 0.09 -0.1378556 -0.0006279026 0.2880288 25
all_skew −0.05 0.21 -0.4237418 -0.0403155341 0.3816406 25
mig_kurt −0.67 0.16 -0.9329652 -0.7095091193 -0.1280714 25
peak_kurt −1.12 0.03 -1.1883554 -1.1297710654 -1.0337500 25
all_kurt 0.26 0.56 -0.3467272 0.1957870557 2.4580910 25

Here we look at skew and excess kurtosis only against the intercept. This allows us to test for differences from 0.

  • Skew of 0 is normal, but below or above are considered left- and right-skewed
  • A normal distribution has an excess kurtosis of 0.

So here, we are not asking if skew or kurtosis changes over the years, just whether it is different from normal (0).

m1 <- lm(mig_skew ~ 1, data = v)
m2 <- lm(peak_skew ~ 1, data = v)
m3 <- lm(all_skew ~ 1, data = v)
m4 <- lm(mig_kurt ~ 1, data = v)
m5 <- lm(peak_kurt ~ 1, data = v)
m6 <- lm(all_kurt ~ 1, data = v)

models <- list(m1, m2, m3, m4, m5, m6)

Tabulated Results

Skew is not significantly different from zero in any period.

Kurtosis is significantly different from zero in all periods.

In the migration and peak migration periods, Kurtosis is negative, corresponding to thick-tailed distributions, squat distributions (closer to a uniform distribution).

However, in the ‘all’ period, kurtosis is actually significantly positive, corresponding to relatively more narrow, thin-tailed distributions.

Code
t <- get_table(models)
write_csv(t, "Data/Datasets/table_skew.csv")
fmt_table(t)
term estimate std error statistic p value n
mig_skew ~ 1 (R2 = 0; R2-adj = 0)
(Intercept) 0.010 0.042 0.236 0.815 25
peak_skew ~ 1 (R2 = 0; R2-adj = 0)
(Intercept) 0.007 0.018 0.377 0.709 25
all_skew ~ 1 (R2 = 0; R2-adj = 0)
(Intercept) −0.049 0.042 −1.184 0.248 25
mig_kurt ~ 1 (R2 = 0; R2-adj = 0)
(Intercept) −0.671 0.032 −20.956 0.000 25
peak_kurt ~ 1 (R2 = 0; R2-adj = 0)
(Intercept) −1.125 0.006 −188.017 0.000 25
all_kurt ~ 1 (R2 = 0; R2-adj = 0)
(Intercept) 0.256 0.112 2.287 0.031 25

The periods

The “all” period is defined individually as being the dates centred on the date at the 50% percentile passage is reached and includes all dates from this date to the end of the year, and that same number of dates prior to the p50 date. See the Skewness & Kurtosis section in the last report.

By definition, the migration period (5-95%) and the peak migration period (25-75%) are truncated distributions, so I don’t think we can really interpret kurtosis for these subsets. With the “all” period, I have attempted to subset the data to get a good measure of the migration period only, but honestly, kurtosis depends so much on where I define the cutoffs I’m not really sure that we can accurately measure it. Looking at the figures, I’d say they’re narrower than normal (so positive kurtosis, which does align with the stats for the “all” period), but that’s just a visual estimate.

This is figure showing the statistics.

Code
dist_figs <- v |>
  select(year, contains("skew"), contains("kurt")) |>
  pivot_longer(cols = -year, names_to = "measure", values_to = "statistic")

ggplot(filter(dist_figs, year != 2001), aes(x = statistic, y = measure, fill = measure)) +
  theme_bw() +
  geom_vline(xintercept = 0, colour = "grey20", linetype = "dotted") +
  geom_boxplot() +
  scale_fill_viridis_d()

This figure shows the distributions in each year over top a normal distribution, for illustration.

I think this shows positive kurtosis. The black lines are simulated normal distributions simulated using the mean, SD and counts from the data.

Code
normal <- v |>
  select(year, max_doy, mig_pop_total) |>
  left_join(summarize(filter(pred, doy > 240), sd = sd(rep(doy, round(count))), .by = "year"), 
            by = "year") |>
  mutate(doy = pmap(list(max_doy, sd, mig_pop_total),
                    \(x, y, z) as.integer(rnorm(mean = x, sd = y, n = z)))) |>
  unnest(doy) |>
  count(year, doy, name = "count")

ggplot(raw, aes(x = doy, y = count)) + 
  theme_bw() +
  geom_smooth(data = normal, method = "gam", colour = "black", linewidth = 1.5) +
  geom_point(size = 1, aes(colour = year), na.rm = TRUE) + 
  geom_line(data = pred, linewidth = 1, aes(group = year, colour = year), na.rm = TRUE) +
  scale_colour_viridis_c() +
  facet_wrap(~year, scales = "free")
`geom_smooth()` using formula = 'y ~ s(x, bs = "cs")'

Here is a scaled view (counts are divided by the yearly standard deviation) of the patterns of migration using our GAM models from the “Calculate Metrics” step.

Code
raw_scaled <- raw |>
  mutate(across(c(count), \(x) scale(x, center = FALSE)),
         .by = "year")

pred_scaled <- pred |>
  mutate(across(c(count, ci99_upper, ci99_lower), \(x) scale(x, center = FALSE)),
         .by = "year")
           
ggplot(raw_scaled, aes(x = doy, y = count)) + 
  theme_bw() +
  geom_point(aes(colour = year), na.rm = TRUE) +
  geom_line(data = pred_scaled, aes(colour = year, group = year), 
            size = 1, na.rm = TRUE) +
  scale_colour_viridis_c()
Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.

Code
model_check_figs(models)

The DHARMa model checks don’t highlight any particular outlier problems, but out of an abundance of caution, we’ll do a quick review of any high Cook’s D values.

Values highlighted in red are less than 4/25, a standard Cook’s D cutoff (4/sample size). These indicate potential influential data points for a particular model.

Code
get_cooks(models) |>
  gt_cooks()
Cook's Distances
year mig_skew ~ 1 peak_skew ~ 1 all_skew ~ 1 mig_kurt ~ 1 peak_kurt ~ 1 all_kurt ~ 1
1998 0.00 0.00 0.00 0.00 0.01 0.00
1999 0.00 0.00 0.00 0.00 0.00 0.04
2000 0.02 0.03 0.04 0.00 0.01 0.00
2001 0.18 0.04 0.14 0.14 0.00 0.00
2002 0.14 0.42 0.04 0.00 0.40 0.01
2003 0.01 0.04 0.02 0.03 0.03 0.01
2004 0.14 0.07 0.18 0.09 0.08 0.00
2005 0.00 0.00 0.00 0.00 0.01 0.00
2006 0.06 0.04 0.08 0.00 0.02 0.04
2008 0.06 0.03 0.09 0.00 0.00 0.00
2009 0.00 0.00 0.02 0.00 0.00 0.67
2010 0.01 0.04 0.01 0.01 0.01 0.01
2011 0.07 0.11 0.03 0.02 0.00 0.05
2012 0.00 0.00 0.01 0.01 0.00 0.00
2013 0.00 0.00 0.00 0.03 0.00 0.04
2014 0.00 0.00 0.01 0.00 0.00 0.02
2015 0.25 0.08 0.18 0.50 0.14 0.00
2016 0.02 0.02 0.04 0.00 0.00 0.00
2017 0.03 0.04 0.05 0.00 0.00 0.03
2018 0.00 0.01 0.00 0.12 0.20 0.03
2019 0.00 0.00 0.04 0.03 0.01 0.01
2020 0.01 0.02 0.00 0.00 0.01 0.03
2021 0.02 0.01 0.02 0.02 0.07 0.00
2022 0.01 0.04 0.02 0.02 0.04 0.00
2023 0.01 0.01 0.00 0.00 0.00 0.02

Check influential points

Now let’s see what happens if we were to omit these years from the respective analyses.

Tables with coefficients for both original and new model are shown, with colour intensity showing strength in the pattern.

For skew of the counts during the migration period, if we omit 2001, the patterns don’t really change (i.e., stay non-significant).

Code
compare(m1, 2002)
model parameter Estimate Std. Error t value Pr(>|t|)
Original - mig_skew ~ 1 (Intercept) 0.010 0.042 0.236 0.815
Drop 2002 - mig_skew ~ 1 (Intercept) -0.006 0.040 -0.141 0.889

For skew of the counts during the peak migration period, if we omit 2002, the patterns don’t really change (i.e., stay non-significant).

Code
compare(m2, 2002)
model parameter Estimate Std. Error t value Pr(>|t|)
Original - peak_skew ~ 1 (Intercept) 0.007 0.018 0.377 0.709
Drop 2002 - peak_skew ~ 1 (Intercept) -0.005 0.014 -0.337 0.739

For kurtosis during the peak migration period, if we omit 2002, the patterns don’t really change.

Code
compare(m5, 2002)
model parameter Estimate Std. Error t value Pr(>|t|)
Original - peak_kurt ~ 1 (Intercept) -1.125 0.006 -188.017 0
Drop 2002 - peak_kurt ~ 1 (Intercept) -1.128 0.005 -233.856 0

For skew of the counts during the “all” period, if we omit 2004, the patterns are stronger, and nearly show a trend, but are still non-significant.

Code
compare(m3, 2004)
model parameter Estimate Std. Error t value Pr(>|t|)
Original - all_skew ~ 1 (Intercept) -0.049 0.042 -1.184 0.248
Drop 2004 - all_skew ~ 1 (Intercept) -0.067 0.040 -1.696 0.103

For kurtosis during the “all” migration period, if we omit 2009, the patterns are less strong, but the same.

Code
compare(m6, 2009)
model parameter Estimate Std. Error t value Pr(>|t|)
Original - all_kurt ~ 1 (Intercept) 0.256 0.112 2.287 0.031
Drop 2009 - all_kurt ~ 1 (Intercept) 0.164 0.067 2.455 0.022

For skew of the counts during the migration period and during the “all” period if we omit 2015, the patterns are stronger, but the same (i.e., not significant).

Code
compare(m1, 2015)
model parameter Estimate Std. Error t value Pr(>|t|)
Original - mig_skew ~ 1 (Intercept) 0.010 0.042 0.236 0.815
Drop 2015 - mig_skew ~ 1 (Intercept) -0.011 0.038 -0.289 0.775
Code
compare(m3, 2015)
model parameter Estimate Std. Error t value Pr(>|t|)
Original - all_skew ~ 1 (Intercept) -0.049 0.042 -1.184 0.248
Drop 2015 - all_skew ~ 1 (Intercept) -0.067 0.039 -1.714 0.100

For kurtosis during the migration period, if we omit 2015, the patterns are stronger, but the same.

Code
compare(m4, 2015)
model parameter Estimate Std. Error t value Pr(>|t|)
Original - mig_kurt ~ 1 (Intercept) -0.671 0.032 -20.956 0
Drop 2015 - mig_kurt ~ 1 (Intercept) -0.694 0.024 -29.362 0

For kurtosis during the peak migration period, if we omit 2018, the patterns don’t really change.

Code
compare(m5, 2018)
model parameter Estimate Std. Error t value Pr(>|t|)
Original - peak_kurt ~ 1 (Intercept) -1.125 0.006 -188.017 0
Drop 2018 - peak_kurt ~ 1 (Intercept) -1.122 0.006 -200.726 0

For kurtosis during the migration period and during the peak migration period, if we omit 2018, the patterns are stronger.

Code
compare(m3, 2018)
model parameter Estimate Std. Error t value Pr(>|t|)
Original - all_skew ~ 1 (Intercept) -0.049 0.042 -1.184 0.248
Drop 2018 - all_skew ~ 1 (Intercept) -0.052 0.044 -1.192 0.246
Code
compare(m4, 2018)
model parameter Estimate Std. Error t value Pr(>|t|)
Original - mig_kurt ~ 1 (Intercept) -0.671 0.032 -20.956 0
Drop 2018 - mig_kurt ~ 1 (Intercept) -0.660 0.031 -21.030 0

Therefore I wouldn’t be concerned

Code
map(models, summary)
[[1]]

Call:
lm(formula = mig_skew ~ 1, data = v)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.42158 -0.12064 -0.03772  0.06249  0.49696 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.009823   0.041615   0.236    0.815

Residual standard error: 0.2081 on 24 degrees of freedom


[[2]]

Call:
lm(formula = peak_skew ~ 1, data = v)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.144701 -0.070028 -0.007474  0.035147  0.281183 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.006846   0.018135   0.377    0.709

Residual standard error: 0.09067 on 24 degrees of freedom


[[3]]

Call:
lm(formula = all_skew ~ 1, data = v)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.37426 -0.15526  0.00916  0.09960  0.43112 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.04948    0.04179  -1.184    0.248

Residual standard error: 0.209 on 24 degrees of freedom


[[4]]

Call:
lm(formula = mig_kurt ~ 1, data = v)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.26170 -0.09539 -0.03825  0.04156  0.54319 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.67126    0.03203  -20.96   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.1602 on 24 degrees of freedom


[[5]]

Call:
lm(formula = peak_kurt ~ 1, data = v)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.063659 -0.012237 -0.005074  0.007144  0.090947 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1.124697   0.005982    -188   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.02991 on 24 degrees of freedom


[[6]]

Call:
lm(formula = all_kurt ~ 1, data = v)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.60289 -0.39834 -0.06037  0.14526  2.20193 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)   0.2562     0.1120   2.287   0.0313 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.56 on 24 degrees of freedom

Reproducibility

Code
devtools::session_info()
─ Session info ───────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.4.0 (2024-04-24)
 os       Ubuntu 22.04.4 LTS
 system   x86_64, linux-gnu
 ui       X11
 language en_CA:en
 collate  en_CA.UTF-8
 ctype    en_CA.UTF-8
 tz       America/Winnipeg
 date     2024-09-23
 pandoc   3.1.1 @ /usr/lib/rstudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)

─ Packages ───────────────────────────────────────────────────────────────────
 ! package      * version  date (UTC) lib source
 P abind          1.4-5    2016-07-21 [?] CRAN (R 4.4.0)
 P backports      1.5.0    2024-05-23 [?] CRAN (R 4.4.0)
 P bit            4.0.5    2022-11-15 [?] CRAN (R 4.4.0)
 P bit64          4.0.5    2020-08-30 [?] CRAN (R 4.4.0)
 P boot           1.3-30   2024-02-26 [?] CRAN (R 4.3.3)
 P broom        * 1.0.6    2024-05-17 [?] CRAN (R 4.4.0)
 P cachem         1.1.0    2024-05-16 [?] CRAN (R 4.4.0)
 P car            3.1-2    2023-03-30 [?] CRAN (R 4.4.0)
 P carData        3.0-5    2022-01-06 [?] CRAN (R 4.4.0)
 P cli            3.6.2    2023-12-11 [?] CRAN (R 4.4.0)
 P codetools      0.2-19   2023-02-01 [?] CRAN (R 4.2.2)
 P colorspace     2.1-0    2023-01-23 [?] CRAN (R 4.4.0)
 P crayon         1.5.2    2022-09-29 [?] CRAN (R 4.4.0)
 P devtools       2.4.5    2022-10-11 [?] CRAN (R 4.4.0)
 P DHARMa       * 0.4.6    2022-09-08 [?] CRAN (R 4.4.0)
 P digest         0.6.35   2024-03-11 [?] CRAN (R 4.4.0)
 P doParallel     1.0.17   2022-02-07 [?] CRAN (R 4.4.0)
 P dplyr        * 1.1.4    2023-11-17 [?] CRAN (R 4.4.0)
 P ellipsis       0.3.2    2021-04-29 [?] CRAN (R 4.4.0)
 P emmeans        1.10.2   2024-05-20 [?] CRAN (R 4.4.0)
 P estimability   1.5.1    2024-05-12 [?] CRAN (R 4.4.0)
 P evaluate       0.23     2023-11-01 [?] CRAN (R 4.4.0)
 P fansi          1.0.6    2023-12-08 [?] CRAN (R 4.4.0)
 P farver         2.1.2    2024-05-13 [?] CRAN (R 4.4.0)
 P fastmap        1.2.0    2024-05-15 [?] CRAN (R 4.4.0)
 P forcats      * 1.0.0    2023-01-29 [?] CRAN (R 4.4.0)
 P foreach        1.5.2    2022-02-02 [?] CRAN (R 4.4.0)
 P fs             1.6.4    2024-04-25 [?] CRAN (R 4.4.0)
 P gap            1.5-3    2023-08-26 [?] CRAN (R 4.4.0)
 P gap.datasets   0.0.6    2023-08-25 [?] CRAN (R 4.4.0)
 P generics       0.1.3    2022-07-05 [?] CRAN (R 4.4.0)
 P ggplot2      * 3.5.1    2024-04-23 [?] CRAN (R 4.4.0)
 P glue           1.7.0    2024-01-09 [?] CRAN (R 4.4.0)
 P gt           * 0.10.1   2024-01-17 [?] CRAN (R 4.4.0)
 P gtable         0.3.5    2024-04-22 [?] CRAN (R 4.4.0)
 P hms            1.1.3    2023-03-21 [?] CRAN (R 4.4.0)
 P htmltools      0.5.8.1  2024-04-04 [?] CRAN (R 4.4.0)
 P htmlwidgets    1.6.4    2023-12-06 [?] CRAN (R 4.4.0)
 P httpuv         1.6.15   2024-03-26 [?] CRAN (R 4.4.0)
 P iterators      1.0.14   2022-02-05 [?] CRAN (R 4.4.0)
 P jsonlite       1.8.8    2023-12-04 [?] CRAN (R 4.4.0)
 P knitr          1.47     2024-05-29 [?] CRAN (R 4.4.0)
 P labeling       0.4.3    2023-08-29 [?] CRAN (R 4.4.0)
 P later          1.3.2    2023-12-06 [?] CRAN (R 4.4.0)
 P lattice        0.22-5   2023-10-24 [?] CRAN (R 4.3.1)
 P lifecycle      1.0.4    2023-11-07 [?] CRAN (R 4.4.0)
 P lme4           1.1-35.3 2024-04-16 [?] CRAN (R 4.4.0)
 P lubridate    * 1.9.3    2023-09-27 [?] CRAN (R 4.4.0)
 P magrittr       2.0.3    2022-03-30 [?] CRAN (R 4.4.0)
 P MASS           7.3-60   2023-05-04 [?] CRAN (R 4.3.1)
 P Matrix         1.6-5    2024-01-11 [?] CRAN (R 4.3.3)
 P memoise        2.0.1    2021-11-26 [?] CRAN (R 4.4.0)
 P mgcv           1.9-1    2023-12-21 [?] CRAN (R 4.3.2)
 P mime           0.12     2021-09-28 [?] CRAN (R 4.4.0)
 P miniUI         0.1.1.1  2018-05-18 [?] CRAN (R 4.4.0)
 P minqa          1.2.6    2023-09-11 [?] CRAN (R 4.4.0)
 P munsell        0.5.1    2024-04-01 [?] CRAN (R 4.4.0)
 P mvtnorm        1.2-5    2024-05-21 [?] CRAN (R 4.4.0)
 P nlme           3.1-165  2024-06-06 [?] CRAN (R 4.4.0)
 P nloptr         2.0.3    2022-05-26 [?] CRAN (R 4.4.0)
 P patchwork    * 1.2.0    2024-01-08 [?] CRAN (R 4.4.0)
 P pillar         1.9.0    2023-03-22 [?] CRAN (R 4.4.0)
 P pkgbuild       1.4.4    2024-03-17 [?] CRAN (R 4.4.0)
 P pkgconfig      2.0.3    2019-09-22 [?] CRAN (R 4.4.0)
 P pkgload        1.3.4    2024-01-16 [?] CRAN (R 4.4.0)
 P plyr           1.8.9    2023-10-02 [?] CRAN (R 4.4.0)
 P profvis        0.3.8    2023-05-02 [?] CRAN (R 4.4.0)
 P promises       1.3.0    2024-04-05 [?] CRAN (R 4.4.0)
 P purrr        * 1.0.2    2023-08-10 [?] CRAN (R 4.4.0)
 P qgam           1.3.4    2021-11-22 [?] CRAN (R 4.4.0)
 P R6             2.5.1    2021-08-19 [?] CRAN (R 4.4.0)
 P rbibutils      2.2.16   2023-10-25 [?] CRAN (R 4.4.0)
 P RColorBrewer   1.1-3    2022-04-03 [?] CRAN (R 4.4.0)
 P Rcpp           1.0.12   2024-01-09 [?] CRAN (R 4.4.0)
 P Rdpack         2.6      2023-11-08 [?] CRAN (R 4.4.0)
 P readr        * 2.1.5    2024-01-10 [?] CRAN (R 4.4.0)
 P remotes        2.5.0    2024-03-17 [?] CRAN (R 4.4.0)
   renv           1.0.7    2024-04-11 [1] CRAN (R 4.4.0)
 P rlang          1.1.3    2024-01-10 [?] CRAN (R 4.4.0)
 P rmarkdown      2.27     2024-05-17 [?] CRAN (R 4.4.0)
 P rstudioapi     0.16.0   2024-03-24 [?] CRAN (R 4.4.0)
 P sass           0.4.9    2024-03-15 [?] CRAN (R 4.4.0)
 P scales         1.3.0    2023-11-28 [?] CRAN (R 4.4.0)
 P sessioninfo    1.2.2    2021-12-06 [?] CRAN (R 4.4.0)
 P shiny          1.8.1.1  2024-04-02 [?] CRAN (R 4.4.0)
 P stringi        1.8.4    2024-05-06 [?] CRAN (R 4.4.0)
 P stringr      * 1.5.1    2023-11-14 [?] CRAN (R 4.4.0)
 P tibble       * 3.2.1    2023-03-20 [?] CRAN (R 4.4.0)
 P tidyr        * 1.3.1    2024-01-24 [?] CRAN (R 4.4.0)
 P tidyselect     1.2.1    2024-03-11 [?] CRAN (R 4.4.0)
 P tidyverse    * 2.0.0    2023-02-22 [?] CRAN (R 4.4.0)
 P timechange     0.3.0    2024-01-18 [?] CRAN (R 4.4.0)
 P tzdb           0.4.0    2023-05-12 [?] CRAN (R 4.4.0)
 P urlchecker     1.0.1    2021-11-30 [?] CRAN (R 4.4.0)
 P usethis        2.2.3    2024-02-19 [?] CRAN (R 4.4.0)
 P utf8           1.2.4    2023-10-22 [?] CRAN (R 4.4.0)
 P vctrs          0.6.5    2023-12-01 [?] CRAN (R 4.4.0)
 P viridisLite    0.4.2    2023-05-02 [?] CRAN (R 4.4.0)
 P vroom          1.6.5    2023-12-05 [?] CRAN (R 4.4.0)
 P withr          3.0.0    2024-01-16 [?] CRAN (R 4.4.0)
 P xfun           0.44     2024-05-15 [?] CRAN (R 4.4.0)
 P xml2           1.3.6    2023-12-04 [?] CRAN (R 4.4.0)
 P xtable         1.8-4    2019-04-21 [?] CRAN (R 4.4.0)
 P yaml           2.3.8    2023-12-11 [?] CRAN (R 4.4.0)

 [1] /home/steffi/Projects/vulture_migration/renv/library/linux-ubuntu-jammy/R-4.4/x86_64-pc-linux-gnu
 [2] /home/steffi/.cache/R/renv/sandbox/linux-ubuntu-jammy/R-4.4/x86_64-pc-linux-gnu/9a444a72

 P ── Loaded and on-disk path mismatch.

──────────────────────────────────────────────────────────────────────────────
Back to top