Supplementary code to ‘Fueling conditions at staging sites can mitigate Arctic warming effects in a migratory bird’ by Rakhimberidev et al., submitted
code by Eldar Rakhimberdiev
You might enjoy html version of this supplementary.
0. Set up output and run options
set.seed(123)0.1 To run or to load?
Parts of this code take time to run, so if you want to go through examples fast and do not want to install all the software
Run_everything=FALSE # Change to TRUE if you want run everything by your computer0.2 Online or local data?
This script can operate with online or local data. If Use_local_data set to TRUE working directory expected to have folders data. If Run_everything set to TRUE results directory will be created and all results will be saved there.
Use_local_data=FALSE # Change to TRUE if you want to use local data0.3 Checks and downloads
if (Use_local_data) {
if (!'results' %in% list.dirs(full.names=FALSE, recursive=FALSE)) {
if (Run_everything) {
dir.create('results')
} else {
stop('results directory not found!')
}
}
if (!'data' %in% list.dirs(full.names=FALSE, recursive=FALSE)) stop('data directory not found!')
} else {
if (!'results' %in% list.dirs(full.names=FALSE, recursive=FALSE)) dir.create('results')
download.file('https://git.io/vXAZ1', destfile='./results/Snow_remote.csv')
download.file('https://git.io/vXgZG', destfile='./results/migration_model.RDS', mode='wb')
download.file('https://git.io/vXED0', destfile='./results/Arrivals_migration_model.csv')
download.file('https://git.io/vFM79', destfile='./results/fueling_model.RDS',mode='wb')
download.file('https://git.io/v7D8z', destfile='./results/phenology_SEM.RDS',mode='wb')
download.file('https://git.io/vXEbX', destfile='./results/Arrival_weights.csv')
download.file('https://git.io/vFM72', destfile='./results/Fueling_slopes.csv')
download.file('https://git.io/vXueC', destfile='./results/survival_result.RDS',mode='wb')
download.file('https://git.io/vXuqH', destfile='./results/Mauritania_counts_model.RDS',mode='wb')
}
if (Use_local_data & !'data' %in% list.dirs(full.names=FALSE, recursive=FALSE)) stop('data directory not found') 0.4 Plot on screen or save plots as pdfs?
Parts of this code take time to run, so if you want to go through examples fast you can specify
Save_pdf=FALSE # Change to TRUE if you want to save pdf of figures.1. Reconstruction of snowmelt from remote sensing data
We used NOAA CDR snow cover extent based on remote sensed data Brown & Robinson 2011. Weekly snow cover data were downloaded on 24×24 km grid and subset by the Limosa lapponica taymyrensis breeding range Lappo et al. 2012. Following code provided in van Gils et al 2016, we extracted date 50% spring snow coverage as linearly interpolated time point at which mean snow cover first time in a year crossed 50%, results are presented in Supplementary Materials Table S1.
1.1 Download publically available data
Load required packages
library(ncdf4)
library(raster)
library(rgdal)
library(maptools)
library(maps)
library(RCurl)
library(fields)Download taimyrensis bar-trailed godwit (Limosa lapponica taymyrensis) range from this GitHub repo
if (Use_local_data) {
area.longlat<-readRDS('data/taimyrensis_godwit_area.RDS')
} else {
area.longlat<-readRDS(gzcon(url('https://git.io/vXRBz')))
}Download and read in R all snow data from NOAA
# first we get filename of the latest file
Page<-getURL('https://www.ncei.noaa.gov/data/snow-cover-extent/access/',
ssl.verifypeer = FALSE, dirlistonly = TRUE, ftp.use.epsv = FALSE)
filename<-paste0('https://www.ncei.noaa.gov/data/snow-cover-extent/access/nhsce_v01r01_19661004_',
substr(strsplit(Page, 'nhsce_v01r01_19661004_')[[1]][2], start=1, stop=11))
download.file(filename, destfile='snow.nc', mode="wb")
dat <- nc_open('snow.nc')
t <- ncvar_get(dat,'time') ## "days since 1966-10-03"
z <- ncvar_get(dat,'snow_cover_extent')
ylat <- ncvar_get(dat,'latitude')
xlon <- ncvar_get(dat,'longitude')
nc_close(dat)1.2 Projection and raster extent
lon <- raster(xlon)
lat <- raster(ylat)
snow <- brick(z)
start <- as.POSIXct("1966-10-03", "GMT")
date <- start + ((t-6)*24*60*60)
projSnow <- "+proj=stere +lat_0=90 +lon_0=10"
xy <- project(cbind(values(lon), values(lat)), projSnow)
area <- spTransform(area.longlat, CRS(projSnow))
Grid_centers<-SpatialPoints(xy, proj4string=CRS(projSnow))
Grid_centers_lon_lat<-spTransform(Grid_centers, CRS=CRS(proj4string(area.longlat)))1.3 Select and plot snow points
# select snow points within godwtis's area
Selected_Points<-which(!is.na(over(Grid_centers, area)))
# Ok now we want to take two points closest to the field site
Field_site<-cbind(106.0,72.8)
Selected_Points_near_field_station<-Selected_Points[
sort(order(spDists(coordinates(Grid_centers_lon_lat[Selected_Points,]),
Field_site, longlat=TRUE))[1:2])]
# Plot
Colors<-c('#1b9e77', '#d95f02', '#7570b3')
plot(area.longlat)
map('world', add=T, col=grey(0.9), fill=TRUE )
map('world', add=T, col=grey(0.5), lwd=2)
plot(area.longlat, border=grey(0.3), lwd=2, add=T, col='#bae4bc')
# plot all snow grid points around
plot(Grid_centers_lon_lat, col=Colors[1], add=TRUE, pch=19, cex=2)
# plot snow grid points withing godwits range
plot(Grid_centers_lon_lat[Selected_Points,], col=Colors[2], pch=19, cex=2, add=TRUE)
# field site
points(Field_site[1],Field_site[2], pch=17, cex=2, col=2)
# two closes points to the field site
plot(Grid_centers_lon_lat[Selected_Points_near_field_station,], col=Colors[3], pch=19, cex=2, add=TRUE)
1.4 Extract first date of 50% snowmelt
Res<-c()
for(i in 1982:2017) {
ind <- which(as.numeric(format(date, "%Y"))==i &
as.numeric(format(as.POSIXct(date), format='%m'))<8)
M<-snow@data@values[Selected_Points,ind]
M_Southern<-snow@data@values[Selected_Points_near_field_station,ind]
# we now want to get last point with non zero snow coverage
Mean_all<-mean(as.numeric(format(c(date[ind][
apply(M, 1, FUN=function(x) which(x==0)[1])-1]), '%j')))
Mean_Southern<-mean(as.numeric(format(date[ind][
apply(M_Southern, 1, FUN=function(x) which(x==0)[1])-1], '%j')))
Res<-rbind(Res, cbind(i, Mean_all, Mean_Southern))
}
Res<-as.data.frame(Res)
names(Res)<-c('Yr', 'jdate50', 'jdate50_south')
write.csv(Res, file='./results/Snow_remote.csv', row.names=FALSE)
plot(jdate50~Yr, data=Res, pch=19, col=Colors[2], ylim=range(c(Res[,2], Res[,3])), las=1, xlab='Year', cex=2)
summary(lm(jdate50~Yr, data=subset(Res, Yr>1995 & Yr<2017)))##
## Call:
## lm(formula = jdate50 ~ Yr, data = subset(Res, Yr > 1995 & Yr <
## 2017))
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.4612 -4.1142 0.1379 2.6605 9.8315
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1715.9698 375.7468 4.567 0.000211 ***
## Yr -0.7748 0.1873 -4.136 0.000561 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.198 on 19 degrees of freedom
## Multiple R-squared: 0.4738, Adjusted R-squared: 0.4461
## F-statistic: 17.11 on 1 and 19 DF, p-value: 0.00056121.5 Comparison of remote sensing data with ground observations
Snow_remote<-read.csv('./results/Snow_remote.csv')
if (Use_local_data) {
Snow_ground<-read.csv('./data/Snow_ground.csv')
} else {
Snow_ground<-read.csv('https://git.io/vXAnO')
}
# now we merge them
Snow_all<-merge(Snow_remote, Snow_ground)
plot(jdate50_south~ jdate50_ground, data=Snow_all, pch=19, cex=2, col=Colors[3], las=1,
xlab='Snow 50% date from ground observations (yday)',
ylab='Snow 50% date from remote data (yday)')
# draw regression line
abline(a=0, b=1, col=1, lwd=2)
abline(lm(jdate50_south~ jdate50_ground, data=Snow_all), lwd=2, col=Colors[3])
# now formally check the regression
print(summary(lm(jdate50_south~ jdate50_ground, data=Snow_all)))##
## Call:
## lm(formula = jdate50_south ~ jdate50_ground, data = Snow_all)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.7745 -1.8274 -0.8513 3.6726 6.4106
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.4266 20.7319 0.214 0.834
## jdate50_ground 0.9736 0.1288 7.557 1.72e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.454 on 15 degrees of freedom
## Multiple R-squared: 0.792, Adjusted R-squared: 0.7781
## F-statistic: 57.11 on 1 and 15 DF, p-value: 1.725e-06# slope does not differ from 1:
print(summary(lm(I(jdate50_south- jdate50_ground)~jdate50_ground, data=Snow_all)))##
## Call:
## lm(formula = I(jdate50_south - jdate50_ground) ~ jdate50_ground,
## data = Snow_all)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.7745 -1.8274 -0.8513 3.6726 6.4106
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.42659 20.73185 0.214 0.834
## jdate50_ground -0.02645 0.12883 -0.205 0.840
##
## Residual standard error: 4.454 on 15 degrees of freedom
## Multiple R-squared: 0.002802, Adjusted R-squared: -0.06368
## F-statistic: 0.04214 on 1 and 15 DF, p-value: 0.8401# and Intercept from 0
print(summary(lm(I(jdate50_south- jdate50_ground)~1, data=Snow_all)))##
## Call:
## lm(formula = I(jdate50_south - jdate50_ground) ~ 1, data = Snow_all)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.676 -1.677 -1.177 3.824 6.824
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.1765 1.0475 0.168 0.868
##
## Residual standard error: 4.319 on 16 degrees of freedom2. Arrival dates from citizen science data
We used data from citizen science project www.trektellen.nl and modelled them with approach similar to one used in Lindén & Mäntyniemi (2011). Model was estimated with JAGS sampler (Plummer 2003) via R2jags (Su & Yahima 2015) interface from R. To run code from this chapter you need to install JAGS software.
2.1 Load packages and data
library(R2jags)
if (Use_local_data) {
sites_years_dates<-read.csv('./data/godwit_trektellen_seven_major_sites.csv', as.is=TRUE)
} else {
sites_years_dates<-read.csv('https://git.io/vXEBH', as.is=TRUE)
}2.2 Formulate model in JAGS language
{
sink("migration_model.jags")
cat("
model {
# Model
for (Obs in 1:nObs) {
log(N[Obs])<-a0 + total_k[Year[Obs]] +
site_i[Site[Obs]]+Offset[Obs]-((Date[Obs]-
T_0k[Year[Obs]])^2)/(2*sigma[Year[Obs]]^2)-
0.9189385-log(sigma[Year[Obs]])
# Bayesian p-value
eval[Obs]<-p[Obs]*r/(1-p[Obs])
E[Obs]<-pow(( Y[Obs]-eval[Obs]), 2)/(eval[Obs]) # +0.5 following Kery & Schaub 2012
Y.new[Obs]~dnegbin( p[Obs] , r)
E.new[Obs]<-pow(( Y.new[Obs]-eval[Obs]), 2)/(eval[Obs])
}
#Priors
for (k in 1:nYears) {
total_k_src[k] ~ dnorm(0, 0.10)
T_0k[k]~dnorm(125, 0.1)
sigma[k]~dnorm(2, 0.1)T(0,)
}
for (i in 1:nSites) {
site_i_src[i] ~ dnorm(0, 0.10)
}
site_i <- site_i_src[] -mean(site_i_src[])
total_k <- total_k_src[]-mean(total_k_src[])
a0 ~ dnorm(6, 0.1)
r ~ dgamma(0.01,0.01)
for (Obs in 1:nObs) {
p[Obs]<- r/(r+N[Obs])
}
#Likelihood
for (Obs in 1:nObs) {
Y[Obs]~dnegbin( p[Obs] , r )
}
fit<-sum(E[])
fit.new<-sum(E.new[])
}",fill = TRUE)
sink()
}2.3 Prepare data and initial values
Data<-list('nObs'=nrow(sites_years_dates),
'nYears'=length(unique(sites_years_dates$Year)),
'nSites'=length(unique(sites_years_dates$Site)),
'Year'=as.numeric(as.factor(sites_years_dates$Year)),
'Site'=as.numeric(as.factor(sites_years_dates$Site)),
'Date'=sites_years_dates$Date,
'Offset'=sites_years_dates$Offset,
'Y'=sites_years_dates$Y)
nSites=length(unique(sites_years_dates$Site))
nYears=length(unique(sites_years_dates$Year))
inits=function() {list(
'a0' =runif(1,3,10),
'total_k_src'=rnorm(nYears, 0, 1),
'T_0k'=rnorm(nYears, 125, 2),
'site_i_src'=log(plogis(rnorm(nSites, 0, 0.5))),
'sigma'=rnorm(nYears, 8, 2))
}
jags.params=c("a0","total_k", "T_0k", "site_i", "sigma", "r", "fit", "fit.new")2.4 Run the model
if (Run_everything) {
# following run takes ~ 5 minutes brt
nchains=4
migration_model<-jags.parallel(data=Data,
inits,
parameters=jags.params,
"migration_model.jags",
n.chains = nchains,
n.thin = 3,
n.iter = 500,
n.burnin = 200,
working.directory = getwd())
saveRDS(migration_model, file='./results/migration_model.RDS')
} else {
if (Use_local_data) {
migration_model<-readRDS('./results/migration_model.RDS')
} else {
migration_model<-readRDS(gzcon(url('https://git.io/vXgZG')))
}
}2.5 Look at the output
print(migration_model)## Inference for Bugs model at "migration_model.jags", fit using jags,
## 4 chains, each with 500 iterations (first 200 discarded), n.thin = 3
## n.sims = 400 iterations saved
## mu.vect sd.vect 2.5% 25%
## T_0k[1] 1.243450e+02 5.920000e-01 1.231800e+02 1.239510e+02
## T_0k[2] 1.221750e+02 8.500000e-01 1.205680e+02 1.216060e+02
## T_0k[3] 1.233630e+02 6.160000e-01 1.221290e+02 1.229930e+02
## T_0k[4] 1.201340e+02 6.300000e-01 1.189310e+02 1.197350e+02
## T_0k[5] 1.188660e+02 8.900000e-01 1.167810e+02 1.184020e+02
## T_0k[6] 1.220490e+02 1.546000e+00 1.185970e+02 1.212800e+02
## T_0k[7] 1.201670e+02 9.100000e-01 1.185670e+02 1.195500e+02
## T_0k[8] 1.223600e+02 1.156000e+00 1.199330e+02 1.216450e+02
## T_0k[9] 1.231320e+02 1.075000e+00 1.210700e+02 1.224800e+02
## T_0k[10] 1.229310e+02 1.200000e+00 1.209720e+02 1.221060e+02
## T_0k[11] 1.225770e+02 1.761000e+00 1.194100e+02 1.213570e+02
## T_0k[12] 1.207730e+02 1.062000e+00 1.186510e+02 1.201110e+02
## T_0k[13] 1.250340e+02 1.515000e+00 1.225450e+02 1.239580e+02
## T_0k[14] 1.180880e+02 1.031000e+00 1.157200e+02 1.175100e+02
## T_0k[15] 1.215560e+02 7.450000e-01 1.202030e+02 1.210840e+02
## T_0k[16] 1.184420e+02 1.632000e+00 1.142740e+02 1.175340e+02
## T_0k[17] 1.217780e+02 1.259000e+00 1.194240e+02 1.210000e+02
## T_0k[18] 1.231100e+02 6.540000e-01 1.218060e+02 1.226490e+02
## T_0k[19] 1.225840e+02 8.890000e-01 1.210280e+02 1.220280e+02
## T_0k[20] 1.223690e+02 8.220000e-01 1.206900e+02 1.218020e+02
## T_0k[21] 1.243310e+02 1.386000e+00 1.216830e+02 1.234990e+02
## T_0k[22] 1.185030e+02 1.203000e+00 1.163990e+02 1.177410e+02
## T_0k[23] 1.232730e+02 9.180000e-01 1.214410e+02 1.227000e+02
## T_0k[24] 1.183180e+02 8.400000e-01 1.167500e+02 1.178220e+02
## a0 6.649000e+00 3.300000e-02 6.591000e+00 6.627000e+00
## fit 1.694172e+12 1.689301e+11 1.419677e+12 1.572122e+12
## fit.new 1.736452e+12 1.055523e+12 5.703880e+11 1.044406e+12
## r 5.580000e-01 1.500000e-02 5.290000e-01 5.490000e-01
## sigma[1] 6.535000e+00 2.920000e-01 6.040000e+00 6.324000e+00
## sigma[2] 8.487000e+00 5.730000e-01 7.533000e+00 8.070000e+00
## sigma[3] 7.159000e+00 3.830000e-01 6.490000e+00 6.903000e+00
## sigma[4] 7.036000e+00 3.350000e-01 6.515000e+00 6.784000e+00
## sigma[5] 7.190000e+00 4.830000e-01 6.413000e+00 6.871000e+00
## sigma[6] 8.204000e+00 8.140000e-01 6.867000e+00 7.590000e+00
## sigma[7] 6.773000e+00 4.150000e-01 6.084000e+00 6.483000e+00
## sigma[8] 7.809000e+00 5.660000e-01 6.863000e+00 7.435000e+00
## sigma[9] 8.928000e+00 7.300000e-01 7.781000e+00 8.380000e+00
## sigma[10] 8.220000e+00 8.290000e-01 6.896000e+00 7.608000e+00
## sigma[11] 1.196500e+01 1.417000e+00 9.768000e+00 1.094900e+01
## sigma[12] 9.458000e+00 6.220000e-01 8.446000e+00 9.016000e+00
## sigma[13] 1.126500e+01 1.381000e+00 9.168000e+00 1.028400e+01
## sigma[14] 6.790000e+00 6.310000e-01 5.827000e+00 6.331000e+00
## sigma[15] 6.866000e+00 3.010000e-01 6.345000e+00 6.669000e+00
## sigma[16] 1.261100e+01 1.181000e+00 1.046600e+01 1.178400e+01
## sigma[17] 1.021300e+01 1.118000e+00 8.522000e+00 9.426000e+00
## sigma[18] 7.441000e+00 3.310000e-01 6.910000e+00 7.199000e+00
## sigma[19] 6.927000e+00 3.550000e-01 6.378000e+00 6.670000e+00
## sigma[20] 6.903000e+00 3.680000e-01 6.262000e+00 6.635000e+00
## sigma[21] 9.236000e+00 8.820000e-01 7.865000e+00 8.625000e+00
## sigma[22] 7.912000e+00 6.980000e-01 6.754000e+00 7.459000e+00
## sigma[23] 7.076000e+00 4.930000e-01 6.299000e+00 6.716000e+00
## sigma[24] 6.661000e+00 4.190000e-01 5.963000e+00 6.406000e+00
## site_i[1] -1.228000e+00 6.700000e-02 -1.371000e+00 -1.270000e+00
## site_i[2] 2.930000e-01 8.800000e-02 1.200000e-01 2.350000e-01
## site_i[3] -6.700000e-01 7.300000e-02 -8.000000e-01 -7.190000e-01
## site_i[4] -6.100000e-02 1.000000e-01 -2.370000e-01 -1.320000e-01
## site_i[5] 1.005000e+00 8.500000e-02 8.520000e-01 9.440000e-01
## site_i[6] 1.290000e-01 8.700000e-02 -3.900000e-02 7.700000e-02
## site_i[7] 5.320000e-01 5.900000e-02 4.220000e-01 4.850000e-01
## total_k[1] 1.082000e+00 1.340000e-01 8.340000e-01 9.800000e-01
## total_k[2] -2.960000e-01 1.210000e-01 -5.330000e-01 -3.660000e-01
## total_k[3] 9.090000e-01 1.370000e-01 6.470000e-01 8.140000e-01
## total_k[4] 3.120000e-01 1.460000e-01 6.900000e-02 2.090000e-01
## total_k[5] 1.230000e-01 1.790000e-01 -2.060000e-01 -1.000000e-03
## total_k[6] -6.200000e-01 1.740000e-01 -9.160000e-01 -7.440000e-01
## total_k[7] 4.980000e-01 1.700000e-01 1.800000e-01 3.940000e-01
## total_k[8] 1.250000e-01 1.480000e-01 -1.710000e-01 2.700000e-02
## total_k[9] 3.490000e-01 1.340000e-01 1.020000e-01 2.650000e-01
## total_k[10] 2.470000e-01 1.590000e-01 -2.200000e-02 1.300000e-01
## total_k[11] -1.480000e+00 1.730000e-01 -1.824000e+00 -1.592000e+00
## total_k[12] 1.320000e-01 1.290000e-01 -1.060000e-01 4.700000e-02
## total_k[13] -3.400000e-01 1.230000e-01 -5.720000e-01 -4.220000e-01
## total_k[14] 1.350000e-01 1.470000e-01 -1.420000e-01 3.400000e-02
## total_k[15] 9.960000e-01 1.390000e-01 7.180000e-01 9.070000e-01
## total_k[16] -1.081000e+00 1.320000e-01 -1.340000e+00 -1.165000e+00
## total_k[17] -1.001000e+00 1.280000e-01 -1.225000e+00 -1.095000e+00
## total_k[18] 5.100000e-02 1.260000e-01 -1.940000e-01 -3.200000e-02
## total_k[19] 5.920000e-01 1.550000e-01 2.590000e-01 4.980000e-01
## total_k[20] 6.200000e-02 1.510000e-01 -2.400000e-01 -2.600000e-02
## total_k[21] -3.720000e-01 1.410000e-01 -6.340000e-01 -4.670000e-01
## total_k[22] -3.970000e-01 1.710000e-01 -7.090000e-01 -5.030000e-01
## total_k[23] -8.800000e-02 2.030000e-01 -4.730000e-01 -2.300000e-01
## total_k[24] 6.200000e-02 1.510000e-01 -2.510000e-01 -4.300000e-02
## deviance 2.310263e+04 1.266400e+01 2.307836e+04 2.309439e+04
## 50% 75% 97.5% Rhat n.eff
## T_0k[1] 1.243360e+02 1.247390e+02 1.254430e+02 1.015 150
## T_0k[2] 1.221310e+02 1.227260e+02 1.239190e+02 1.013 220
## T_0k[3] 1.233850e+02 1.238080e+02 1.245620e+02 1.019 120
## T_0k[4] 1.200960e+02 1.205220e+02 1.213200e+02 1.028 85
## T_0k[5] 1.189440e+02 1.194740e+02 1.204460e+02 1.007 240
## T_0k[6] 1.222400e+02 1.231070e+02 1.246110e+02 1.015 360
## T_0k[7] 1.201130e+02 1.207500e+02 1.220330e+02 1.011 400
## T_0k[8] 1.223360e+02 1.231010e+02 1.246740e+02 1.014 150
## T_0k[9] 1.230950e+02 1.238020e+02 1.255170e+02 0.999 400
## T_0k[10] 1.227220e+02 1.237110e+02 1.257710e+02 1.010 380
## T_0k[11] 1.225130e+02 1.238010e+02 1.259470e+02 1.002 400
## T_0k[12] 1.208520e+02 1.215420e+02 1.227650e+02 1.006 400
## T_0k[13] 1.248930e+02 1.258900e+02 1.287900e+02 0.999 400
## T_0k[14] 1.182320e+02 1.188120e+02 1.197100e+02 1.014 300
## T_0k[15] 1.215250e+02 1.220290e+02 1.230640e+02 1.005 270
## T_0k[16] 1.185680e+02 1.195750e+02 1.213290e+02 1.015 200
## T_0k[17] 1.217630e+02 1.226170e+02 1.242670e+02 1.013 160
## T_0k[18] 1.230720e+02 1.235580e+02 1.244680e+02 1.001 400
## T_0k[19] 1.226080e+02 1.231310e+02 1.242590e+02 1.004 400
## T_0k[20] 1.224330e+02 1.229370e+02 1.238590e+02 1.002 400
## T_0k[21] 1.243690e+02 1.252610e+02 1.269030e+02 1.001 400
## T_0k[22] 1.184210e+02 1.192160e+02 1.211960e+02 1.011 340
## T_0k[23] 1.233300e+02 1.238440e+02 1.251250e+02 1.028 120
## T_0k[24] 1.183710e+02 1.188670e+02 1.196790e+02 1.014 240
## a0 6.647000e+00 6.672000e+00 6.717000e+00 1.011 200
## fit 1.679602e+12 1.798219e+12 2.059171e+12 1.023 98
## fit.new 1.423153e+12 2.137605e+12 4.432844e+12 1.006 260
## r 5.570000e-01 5.670000e-01 5.890000e-01 1.005 330
## sigma[1] 6.512000e+00 6.711000e+00 7.187000e+00 1.003 370
## sigma[2] 8.448000e+00 8.879000e+00 9.776000e+00 0.999 400
## sigma[3] 7.126000e+00 7.397000e+00 8.011000e+00 1.002 400
## sigma[4] 7.011000e+00 7.260000e+00 7.710000e+00 1.008 230
## sigma[5] 7.127000e+00 7.502000e+00 8.281000e+00 1.014 290
## sigma[6] 8.088000e+00 8.679000e+00 9.812000e+00 1.003 400
## sigma[7] 6.742000e+00 7.019000e+00 7.697000e+00 1.016 230
## sigma[8] 7.731000e+00 8.122000e+00 9.025000e+00 1.000 400
## sigma[9] 8.822000e+00 9.382000e+00 1.059300e+01 0.997 400
## sigma[10] 8.166000e+00 8.741000e+00 1.015700e+01 1.015 140
## sigma[11] 1.182800e+01 1.278100e+01 1.541100e+01 1.001 400
## sigma[12] 9.424000e+00 9.829000e+00 1.095900e+01 1.000 400
## sigma[13] 1.109200e+01 1.195500e+01 1.429800e+01 1.002 400
## sigma[14] 6.722000e+00 7.131000e+00 8.269000e+00 1.054 61
## sigma[15] 6.838000e+00 7.042000e+00 7.534000e+00 1.004 300
## sigma[16] 1.250000e+01 1.331300e+01 1.521000e+01 1.003 390
## sigma[17] 1.002800e+01 1.078400e+01 1.292900e+01 1.004 400
## sigma[18] 7.412000e+00 7.649000e+00 8.143000e+00 1.002 400
## sigma[19] 6.882000e+00 7.125000e+00 7.714000e+00 0.999 400
## sigma[20] 6.899000e+00 7.149000e+00 7.666000e+00 0.998 400
## sigma[21] 9.082000e+00 9.735000e+00 1.139200e+01 1.015 250
## sigma[22] 7.811000e+00 8.321000e+00 9.533000e+00 1.030 140
## sigma[23] 7.023000e+00 7.362000e+00 8.161000e+00 1.003 400
## sigma[24] 6.603000e+00 6.865000e+00 7.646000e+00 1.019 140
## site_i[1] -1.224000e+00 -1.178000e+00 -1.116000e+00 1.005 340
## site_i[2] 2.950000e-01 3.580000e-01 4.560000e-01 1.029 400
## site_i[3] -6.680000e-01 -6.210000e-01 -5.350000e-01 0.998 400
## site_i[4] -7.000000e-02 7.000000e-03 1.310000e-01 1.012 180
## site_i[5] 9.990000e-01 1.063000e+00 1.177000e+00 1.012 180
## site_i[6] 1.300000e-01 1.880000e-01 2.870000e-01 1.024 97
## site_i[7] 5.330000e-01 5.730000e-01 6.450000e-01 0.999 400
## total_k[1] 1.080000e+00 1.171000e+00 1.359000e+00 1.003 400
## total_k[2] -3.090000e-01 -2.100000e-01 -6.300000e-02 1.008 240
## total_k[3] 9.020000e-01 9.930000e-01 1.173000e+00 0.998 400
## total_k[4] 3.000000e-01 4.030000e-01 6.260000e-01 1.002 400
## total_k[5] 1.170000e-01 2.500000e-01 4.860000e-01 0.998 400
## total_k[6] -6.360000e-01 -5.060000e-01 -2.380000e-01 1.017 160
## total_k[7] 4.780000e-01 6.060000e-01 8.540000e-01 1.000 400
## total_k[8] 1.330000e-01 2.170000e-01 3.990000e-01 0.998 400
## total_k[9] 3.450000e-01 4.280000e-01 6.500000e-01 1.008 220
## total_k[10] 2.420000e-01 3.510000e-01 5.740000e-01 1.000 400
## total_k[11] -1.480000e+00 -1.371000e+00 -1.136000e+00 1.010 190
## total_k[12] 1.280000e-01 2.130000e-01 4.040000e-01 0.999 400
## total_k[13] -3.420000e-01 -2.460000e-01 -1.040000e-01 0.997 400
## total_k[14] 1.270000e-01 2.290000e-01 4.330000e-01 1.006 280
## total_k[15] 9.980000e-01 1.092000e+00 1.260000e+00 1.008 210
## total_k[16] -1.081000e+00 -9.900000e-01 -8.010000e-01 1.002 400
## total_k[17] -1.007000e+00 -9.220000e-01 -7.320000e-01 1.003 400
## total_k[18] 4.500000e-02 1.420000e-01 2.780000e-01 1.001 400
## total_k[19] 5.810000e-01 6.970000e-01 8.700000e-01 1.015 400
## total_k[20] 6.300000e-02 1.630000e-01 3.430000e-01 1.000 400
## total_k[21] -3.760000e-01 -2.770000e-01 -9.800000e-02 0.998 400
## total_k[22] -4.100000e-01 -2.940000e-01 -3.700000e-02 1.001 400
## total_k[23] -9.600000e-02 4.400000e-02 3.240000e-01 1.002 400
## total_k[24] 5.700000e-02 1.560000e-01 3.670000e-01 1.010 240
## deviance 2.310210e+04 2.311026e+04 2.312809e+04 1.007 400
##
## For each parameter, n.eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
##
## DIC info (using the rule, pD = var(deviance)/2)
## pD = 80.2 and DIC = 23182.9
## DIC is an estimate of expected predictive error (lower deviance is better).cat('Bayesian p-value = ',
mean(migration_model$BUGSoutput$sims.list$fit.new<migration_model$BUGSoutput$sims.list$fit), '\n')## Bayesian p-value = 0.59cat('lack of fit ratio = ',
mean(migration_model$BUGSoutput$sims.list$fit)/mean(migration_model$BUGSoutput$sims.list$fit.new), '\n')## lack of fit ratio = 0.9756514Arrivals<-data.frame(Yr=unique(sites_years_dates$Year_real),
arrival=migration_model$BUGSoutput$mean$T_0k,
sigma=migration_model$BUGSoutput$mean$sigma)
write.csv(Arrivals, file="./results/Arrivals_migration_model.csv", row.names=F)3. Causality in phenological data with structural equations modelling
3.1 Load packages and data
Here we load arrival time estimated in 2.5 data on lugworm densities in the Dutch Wadden Sea.
library(R2jags)
Arrivals_WZ<-read.csv('./results/Arrivals_migration_model.csv')
Snowmelt<-read.csv('./results/Snow_remote.csv', as.is=TRUE)
if (Use_local_data) {
# read breeding data
Nests<-read.csv('./data/Nests.csv', as.is=TRUE)
# read crane flies
Tipula<-read.csv('./data/Tipula.csv', as.is=TRUE)
# read arrivals to Khatanga
Arr_Taimyr<-read.csv('./data/Arrivals_Khat.csv', as.is=TRUE)
} else {
# read breeding data
Nests<-read.csv('https://git.io/vXu3y', as.is=TRUE)
# read crane flies
Tipula<-read.csv('https://git.io/vXu3x', as.is=TRUE)
# read arrivals to Khatanga
Arr_Taimyr<-read.csv('https://git.io/vXusf', as.is=TRUE)
}Merge data
tmp<-merge(Arrivals_WZ[,1:2], Snowmelt[,1:2], all=TRUE)
names(tmp)[2:3]<-c('Arrivals_WZ', 'Snowmelt')
tmp<-merge(tmp, Arr_Taimyr[,c(1,3)], all=TRUE)
names(tmp)[4]<-'Arr_T'
tmp<-merge(tmp, subset(Tipula, Region=='South', select=c('Yr', 'yday')), all=TRUE)
names(tmp)[5]<-'Tipula'
Data_all<-subset(tmp, Yr>=1992)
Nests<-subset(Nests, Region=='South', select=c('Yr', 'yday'))
Nests<-Nests[order(Nests$Yr),]
Years2nests<-sapply(Nests$Yr, FUN=function(x) which(Data_all$Yr==x))3.2 Formulate model
{
sink("pathsimple.scaled.jags")
cat("
model {
# The model
for (i in 1:length(time)) {
Arr_WZ_exp[i] = b1.1*time[i]
Arr_WZ[i] ~ dnorm(Arr_WZ_exp[i], Arr_WZ_tau)
Snowmelt_exp[i] = b2.1*time[i]
Snowmelt[i] ~ dnorm(Snowmelt_exp[i], Snowmelt_tau)
Arr_T_exp[i] = b3.1*time[i] + b3.2*Snowmelt[i] + b3.3*Arr_WZ[i]
Arr_T[i] ~ dnorm(Arr_T_exp[i], Arr_T_tau)
Tipula_exp[i] = b4.1*time[i] + b4.2*Snowmelt[i]
Tipula[i] ~ dnorm(Tipula_exp[i], Tipula_tau)
Breeding_exp[i] = b5.1*time[i] + b5.2*Snowmelt[i] + b5.3 * Arr_T[i] + b5.4*Tipula[i]
Breeding[i] ~ dnorm(Breeding_exp[i], Breeding_tau)
# so I will just add another loop
}
for (nest in 1:length(Nests)) {
Nests[nest]~dnorm(Breeding[Years2nests[nest]], Nests_tau)
}
## Priors
b1.1 ~ dnorm(0,0.00001);
Arr_WZ_tau = pow(Arr_WZ_sigma, -2)
Arr_WZ_sigma ~ dnorm(0,1)T(0,)
b2.1 ~ dnorm(0,0.00001)
Snowmelt_tau = pow(Snowmelt_sigma, -2)
Snowmelt_sigma ~ dnorm(0,1)T(0,)
b3.1 ~ dnorm(0,0.00001); b3.2 ~ dnorm(0,0.00001);b3.3 ~ dnorm(0,0.00001)
Arr_T_tau = pow(Arr_T_sigma, -2)
Arr_T_sigma ~ dnorm(0,1)T(0,)
b4.1 ~ dnorm(0,0.00001); b4.2 ~ dnorm(0,0.00001)
Tipula_tau = pow(Tipula_sigma, -2)
Tipula_sigma ~ dnorm(0,1)T(0,)
b5.1 ~ dnorm(0,0.00001); b5.2 ~ dnorm(0,0.00001);
b5.3 ~ dnorm(0,0.00001); b5.4 ~ dnorm(0,0.00001);
Breeding_tau = pow(Breeding_sigma, -2)
Breeding_sigma ~ dnorm(0,1)T(0,)
Nests_tau = pow(Nests_sigma, -2)
Nests_sigma ~ dnorm(0,1)T(0,)
}
",fill=TRUE)
sink()
}3.3 Prepare data and initial values
Model does not have intercepts, so it is important to center data
Data_cent<-list(Arr_WZ=scale(subset(Data_all, select=Arrivals_WZ, drop=TRUE), scale=FALSE)[,1],
Snowmelt=scale(subset(Data_all, select=Snowmelt, drop=TRUE), scale=FALSE)[,1],
Arr_T=scale(subset(Data_all, select=Arr_T, drop=TRUE), scale=FALSE)[,1],
Tipula=scale(subset(Data_all, select=Tipula, drop=TRUE), scale=FALSE)[,1],
Nests=scale(Nests[,2], scale=FALSE)[,1],
Years2nests=Years2nests,
time=seq(-12, 12))
Inits<-function() {
list(
b1.1=runif(1,-1, 1),
b2.1=runif(1,-1, 1),
b3.1=runif(1,-1, 1),
b3.2=runif(1,-1, 1),
b3.3=runif(1,-1, 1),
b4.1=runif(1,-1, 1),
b4.2=runif(1,-1, 1),
b5.1=runif(1,-1, 1),
b5.2=runif(1,-1, 1),
b5.3=runif(1,-1, 1),
b5.4=runif(1,-1, 1),
Arr_WZ_sigma=1,
Snowmelt_sigma=1,
Arr_T_sigma=1,
Tipula_sigma=1,
Breeding_sigma=1,
Nests_sigma=1)
}
Parameters<-c('b1.1',
'b2.1',
'b3.1', 'b3.2', 'b3.3',
'b4.1', 'b4.2',
'b5.1', 'b5.2', 'b5.3', 'b5.4',
'Arr_WZ_sigma', 'Snowmelt_sigma', 'Arr_T_sigma', 'Breeding_sigma', 'Nests_sigma')3.4 Run the model
if (Run_everything) {
# following run takes ~ 4 minutes brt
phenology_SEM<-jags.parallel(data=Data_cent ,
Inits,
parameters=Parameters,
"pathsimple.scaled.jags",
n.chains = 5,
n.thin = 500,
n.iter = 1000000,
n.burnin = 100000)
saveRDS(phenology_SEM, file='./results/phenology_SEM.RDS')
}
phenology_SEM<-readRDS('./results/phenology_SEM.RDS')3.5 Explore and save model results
print(phenology_SEM)## Inference for Bugs model at "pathsimple.scaled.jags", fit using jags,
## 5 chains, each with 1e+06 iterations (first 1e+05 discarded), n.thin = 500
## n.sims = 9000 iterations saved
## mu.vect sd.vect 2.5% 25% 50% 75% 97.5%
## Arr_T_sigma 2.770 0.340 2.191 2.530 2.742 2.977 3.507
## Arr_WZ_sigma 1.984 0.269 1.534 1.796 1.960 2.145 2.589
## Breeding_sigma 0.759 0.572 0.026 0.301 0.648 1.098 2.100
## Nests_sigma 2.739 0.436 1.996 2.428 2.703 3.011 3.680
## Snowmelt_sigma 5.047 0.421 4.296 4.753 5.023 5.314 5.940
## b1.1 -0.049 0.059 -0.167 -0.088 -0.050 -0.010 0.066
## b2.1 -0.737 0.140 -1.011 -0.832 -0.738 -0.641 -0.463
## b3.1 -0.180 0.109 -0.396 -0.251 -0.179 -0.108 0.034
## b3.2 0.201 0.087 0.032 0.143 0.200 0.258 0.374
## b3.3 0.224 0.328 -0.439 0.005 0.231 0.447 0.839
## b4.1 -0.019 0.218 -0.447 -0.164 -0.016 0.127 0.403
## b4.2 0.365 0.158 0.051 0.258 0.367 0.472 0.671
## b5.1 0.156 0.291 -0.409 -0.040 0.151 0.349 0.735
## b5.2 0.322 0.200 -0.068 0.190 0.322 0.455 0.718
## b5.3 2.040 0.489 1.110 1.716 2.023 2.350 3.073
## b5.4 -0.337 0.252 -0.831 -0.502 -0.337 -0.171 0.159
## deviance 563.297 7.363 550.742 558.048 562.682 567.906 579.267
## Rhat n.eff
## Arr_T_sigma 1.001 9000
## Arr_WZ_sigma 1.001 9000
## Breeding_sigma 1.001 4900
## Nests_sigma 1.001 5500
## Snowmelt_sigma 1.001 9000
## b1.1 1.001 9000
## b2.1 1.001 9000
## b3.1 1.001 8700
## b3.2 1.001 7900
## b3.3 1.001 9000
## b4.1 1.001 9000
## b4.2 1.001 9000
## b5.1 1.001 9000
## b5.2 1.001 9000
## b5.3 1.001 9000
## b5.4 1.001 9000
## deviance 1.001 9000
##
## For each parameter, n.eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
##
## DIC info (using the rule, pD = var(deviance)/2)
## pD = 27.1 and DIC = 590.4
## DIC is an estimate of expected predictive error (lower deviance is better).Parameter.Index<-which(substr(names(phenology_SEM$BUGSoutput$sims.list),1,1)=='b')
SEM_coefs<-data.frame(
'mean'=sapply(phenology_SEM$BUGSoutput$sims.list[Parameter.Index],
FUN=function(x) mean(x)),
'sd'=sapply(phenology_SEM$BUGSoutput$sims.list[Parameter.Index],
FUN=function(x) sd(x)),
'median'=sapply(phenology_SEM$BUGSoutput$sims.list[Parameter.Index],
FUN=function(x) median(x)),
'2.5%'=sapply(phenology_SEM$BUGSoutput$sims.list[Parameter.Index],
FUN=function(x) quantile(x, 0.025)),
'97.5%'=sapply(phenology_SEM$BUGSoutput$sims.list[Parameter.Index],
FUN=function(x) quantile(x, 0.975)),
'P(|b|>0)'=sapply(phenology_SEM$BUGSoutput$sims.list[Parameter.Index],
FUN=function(x) max(mean(x>0), mean(x<0))), check.names=FALSE)
Vars<-data.frame(rbind(c('Time', 'Arrival Wadden Sea'),
c('Time','Snowmelt'),
c('Time', 'Arrival to Taimyr'),
c('Snowmelt', 'Arrival to Taimyr'),
c('Arrival Wadden Sea', 'Arrival to Taimyr'),
c('Time', 'Crane fly emergence'),
c('Snowmelt', 'Crane fly emergence'),
c('Time', 'Clutch Initiation'),
c('Snowmelt', 'Clutch Initiation'),
c('Arrival to Taimyr', 'Clutch Initiation'),
c('Crane fly emergence', 'Clutch Initiation')))
names(Vars)<-c('Effect of', 'on')
SEM_coefs<-cbind(Vars, SEM_coefs)
print(SEM_coefs, digits=2)## Effect of on mean sd median 2.5%
## b1.1 Time Arrival Wadden Sea -0.049 0.059 -0.050 -0.167
## b2.1 Time Snowmelt -0.737 0.140 -0.738 -1.011
## b3.1 Time Arrival to Taimyr -0.180 0.109 -0.179 -0.396
## b3.2 Snowmelt Arrival to Taimyr 0.201 0.087 0.200 0.032
## b3.3 Arrival Wadden Sea Arrival to Taimyr 0.224 0.328 0.231 -0.439
## b4.1 Time Crane fly emergence -0.019 0.218 -0.016 -0.447
## b4.2 Snowmelt Crane fly emergence 0.365 0.158 0.367 0.051
## b5.1 Time Clutch Initiation 0.156 0.291 0.151 -0.409
## b5.2 Snowmelt Clutch Initiation 0.322 0.200 0.322 -0.068
## b5.3 Arrival to Taimyr Clutch Initiation 2.040 0.489 2.023 1.110
## b5.4 Crane fly emergence Clutch Initiation -0.337 0.252 -0.337 -0.831
## 97.5% P(|b|>0)
## b1.1 0.066 0.80
## b2.1 -0.463 1.00
## b3.1 0.034 0.95
## b3.2 0.374 0.99
## b3.3 0.839 0.76
## b4.1 0.403 0.53
## b4.2 0.671 0.99
## b5.1 0.735 0.70
## b5.2 0.718 0.95
## b5.3 3.073 1.00
## b5.4 0.159 0.91write.csv(SEM_coefs, file='./results/SEM_coefs.csv')4. Migration speed and stopover duration from satellite tracking data
4.1 Download data
if (Use_local_data) {
# read breeding data
All_tracks<-read.csv('./data/barg_2016_tracks.csv', as.is=TRUE)
} else {
All_tracks<-read.csv('https://git.io/v7gXj', as.is=TRUE)
}
All_tracks$Time<-as.POSIXct(All_tracks$Time, tz='GMT')4.2 Estimate schedules
Individuals<-unique(All_tracks$ID)
Schedules<-data.frame()
for (i in 1:length(Individuals)) {
Track<-subset(All_tracks, ID==(Individuals[i]))
Year<-unique(as.numeric(format(Track$Time, format='%Y')))
New_time<-seq(min(Track$Time), max(Track$Time), by=60)
New_lon<-approx(x=Track$Time, y=Track$lon, xout=New_time)$y
New_lat<-approx(x=Track$Time, y=Track$lat, xout=New_time)$y
Schedules<-rbind(Schedules, data.frame(Bird_id= (Individuals[i]),
Dep_Mau=ifelse(New_lat[1]> 21.5, NA, New_time[which(New_lat>21.5)[1]]),
Arr_WZ=ifelse(New_lat[1]> 21.5, NA, New_time[which(New_lat>52.5)[1]]),
Dep_WZ=New_time[which(New_lon>10)[1]],
Arr_Tai= New_time[which(New_lon>80)[1]],
year=Year))
}
Schedules[,2]<-as.POSIXct(Schedules[,2], tz="GMT", origin='1970-01-01')
Schedules[,3]<-as.POSIXct(Schedules[,3], tz="GMT", origin='1970-01-01')
Schedules$Migration_to_WZ<-as.numeric((difftime(Schedules[,3],Schedules[,2], units='days')))
Schedules$Refueling<-as.numeric((difftime(Schedules[,4],Schedules[,3], units='days')))
Schedules$Migration_to_Tai<-as.numeric((difftime(Schedules[,5],Schedules[,4], units='days')))
print(mean(Schedules$Migration_to_Tai)) # 5.16## [1] 5.160503print(median(Schedules$Migration_to_Tai)) #5.72## [1] 5.722222print(mean(Schedules$Migration_to_WZ, na.rm=TRUE)) # 3.78## [1] 3.7767364.3 Estimate stopovers
Periods_all<-c()
for (i in 1:length(Individuals)) {
Track<-subset(All_tracks, ID==(Individuals[i]) & Time<as.POSIXct('2016/06/24'))
Year<-unique(as.numeric(format(Track$Time, format='%Y')))
Index_stationary<-rep(1, length=length(spDists(cbind(Track$lon, Track$lat), longlat=TRUE, segments=TRUE)))
Index_stationary[
which(spDists(cbind(Track$lon, Track$lat), longlat=TRUE, segments=TRUE)>25)]<-0#>25km=>movement
Starts<-unique(c(1, which(diff(Index_stationary)==1)+1))
Ends<-unique(c(which(diff(Index_stationary)==-1)+1, length(Index_stationary)))
Periods<-data.frame(Start_Index=Starts, End_Index=Ends,
Start_time=Track$Time[Starts], End_time=Track$Time[Ends])
Periods$Duration<- (as.numeric(Periods$End_time)- as.numeric(Periods$Start_time))/3600/24 # in days
Periods<-Periods[Periods$Duration>1,]
Periods$Mean_lon<-sapply(1:nrow(Periods),
FUN=function(x) mean(Track$lon[Periods$Start_Index[x]: Periods$End_Index[x]]))
Periods$Mean_lat<-sapply(1:nrow(Periods),
FUN=function(x) mean(Track$lat[Periods$Start_Index[x]: Periods$End_Index[x]]))
Periods$Year= Year
# and now I want to exclude ones in the wintering, Wadden Sea or breeding grounds
Periods<-Periods[!(Periods$Mean_lat<21.5 |
(Periods$Mean_lat>52.4 & Periods$Mean_lat<56 & Periods$Mean_lon > 4 & Periods$Mean_lon < 10) |
Periods$Mean_lon > 80),]
if (nrow(Periods)>=1) {
Periods_all<-rbind(Periods_all, Periods)
}
}
print(mean(Periods_all$Duration)) # 1.8 days on average## [1] 1.80428print(nrow(Periods_all)) # 6 stopovers for longer than 24 hours.## [1] 6print(length(which(Periods_all$Duration>2))) #only two stopover longer than 2 days.## [1] 25. Refueling rates and their relationship with lugworm abundance
5.1 Estimate yearly arrival mass from Castricum data
library(lme4)
if (Use_local_data) {
Capture_data<-read.csv('./data/Capture_data.csv', as.is=TRUE)
} else {
Capture_data<-read.csv("https://git.io/vXEDL", as.is=TRUE)
}
Capture_data$date<-as.Date(Capture_data$date)
# arrival mass females
aggr_data_females<-subset(Capture_data, Sex==2 & Site == "Castricum")
Arrival_weight_females<-aggregate(aggr_data_females$Mass,
by=list(aggr_data_females$Yr), FUN=mean)
names(Arrival_weight_females) <-c("Yr", "Arr_Weight")
aggr_data_females$fYr<-as.factor(aggr_data_females$Yr)
weight_model_females<-lmer(Mass~-1+(1|fYr), data=aggr_data_females)
mu_W0_females<-coef(weight_model_females)$fYr[,1]
sigma2_W0_females<-var(residuals(weight_model_females))
# arrival mass males
aggr_data_males<-subset(Capture_data, Sex==1 & Site == "Castricum")
Arrival_weight_males<-aggregate(aggr_data_males$Mass,
by=list(aggr_data_males$Yr), FUN=mean)
names(Arrival_weight_males) <-c("Yr", "Arr_Weight")
aggr_data_males$fYr<-as.factor(aggr_data_males$Yr)
weight_model_males<-lmer(Mass~-1+(1|fYr), data=aggr_data_males)
mu_W0_males<-coef(weight_model_males)$fYr[,1]
sigma2_W0_males<-var(residuals(weight_model_males))
Arrival_weights<-data.frame(Yr=as.numeric(levels(aggr_data_females$fYr)),
Male_arrival_weight=mu_W0_males,
Female_arrival_weight=mu_W0_females, Male_arrival_weight_sigma=sqrt(sigma2_W0_males),
Female_arrival_weight_sigma=sqrt(sigma2_W0_females))
write.csv(Arrival_weights, file='./results/Arrival_weights.csv', row.names=FALSE)5.2 Load packages and data
Here we load arrival time estimated in 2.5 data on lugworm densities in the Dutch Wadden Sea.
library(R2jags)
Arrivals<-read.csv('./results/Arrivals_migration_model.csv')
if (Use_local_data) {
Capture_data<-read.csv('./data/Capture_data.csv', as.is=TRUE)
Lugworms<-read.csv('./data/Lugworms.csv')
} else {
Capture_data<-read.csv("https://git.io/vXEDL", as.is=TRUE)
Lugworms<-read.csv('https://git.io/vXESR')
}
Lugworms<-subset(Lugworms, Yr>1995)
Capture_data$date<-as.Date(Capture_data$date)
names(Arrivals)[2]<-"arrival_day"
Data_to_model<-merge(subset(Capture_data, Site!='Castricum'), Arrivals)
Data_to_model<-subset(Data_to_model, Yr %in% Lugworms$Yr)
mu_T0<-Arrivals$arrival[Arrivals$Yr %in% unique(Data_to_model$Yr)]
sigma2_T0<-(Arrivals$sigma^2)[Arrivals$Yr %in% unique(Data_to_model$Yr)]
Lugworms_vec<-Lugworms[Lugworms$Yr %in% unique(Data_to_model$Yr),2]5.3 Lugworm data validation
While refueling in the Wadden Sea taymyrensis godwits concentrate near Terschelling Island (53.4°N, 5.34°E), that is 50 km away from long term lugworm sampling area at Balgazand (52.9°N, 4.87°E). In order to validate use of Balgzand data on lugworm abundance for godwits we used shorter, but covering both areas dataset - Synoptic Intertidal Benthic Survey (SIBES) - large scale grid sampling of benthic fauna in the Wadden Sea (Bijleveld et al. 2012, Compton et al. 2013). We used SIBES data from sampling stations within 500 meters from Terschelling coast (refueling area) and from Balgzand for all available years (2008 - 2014).
if (Use_local_data) {
Lugworms_SIBES<-read.csv('./data/lugworm_densities_two_sites.csv')
} else {
Lugworms_SIBES<-read.csv('https://git.io/vMei1')
}
print(cor.test(Lugworms_SIBES$Balgzand,Lugworms_SIBES$Terschelling))##
## Pearson's product-moment correlation
##
## data: Lugworms_SIBES$Balgzand and Lugworms_SIBES$Terschelling
## t = 3.7207, df = 5, p-value = 0.0137
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.2935283 0.9785580
## sample estimates:
## cor
## 0.8571243Data on both sites show correlated changes over years (Pearson’s r=0.86) allowing us to use Balgzand data as a proxy for refueling conditions for godwits.
5.4 Formulate model
{
sink("fueling_model.jags")
cat("
model {
# prior
for (i in 1:N) {
T0[i] ~ dnorm(mu_T0[i2k[i]], tau_T0[i2k[i], i2s[i]])
}
# Michaelins - Menten or functional response type II equation for fueling on lugworms
for(k in 1:K) {
Alpha_t[k,1] ~ dnorm(((11*Lugworms[k])/(K_sat_m+Lugworms[k])),
1/sigma2_alpha)
Alpha_t[k,2] ~ dnorm(((13.2*Lugworms[k])/(K_sat_f+Lugworms[k])),
1/sigma2_alpha_fem)
# 11 for males and 13.2 for females are fueling maximum rates measured by Prokosch
}
# Likelihood
for (i in 1:N) {
Wt_mu[i] <- mu_W0[i2k[i],i2s[i]]+(Time[i]-T0[i])*(Alpha_t[i2k[i], i2s[i]])
Wt[i] ~ dnorm(Wt_mu[i],Wt_tau[i2k[i], i2s[i]])
}
# params an hyper-priors
for (k in 1:K) {
for (s in 1:2) {
Wt_s2[k,s] <- sigma2_W0[s]
Wt_tau[k,s]<-1/Wt_s2[k,s]
}
}
sigma2_alpha~dunif(0, 5)
sigma2_alpha_fem~dunif(0, 5)
K_sat_m~dunif(1, 20)
K_sat_f~dunif(1, 20)
for(k in 1:K) {
tau_T0[k,1]<-1/(sigma2_T0[k]+sigma2_W0[1]/Alpha_t[k,1])
tau_T0[k,2]<-1/(sigma2_T0[k]+sigma2_W0[2]/Alpha_t[k,2])
}
}",fill = TRUE)
sink()
}5.5 Prepare data and initial values
Data=list(Wt=Data_to_model$Mass, Time=Data_to_model$yday,
i2k=as.numeric(as.factor(Data_to_model$Yr)),
i2s=Data_to_model$Sex, N=length(Data_to_model$Mass),
K=length(unique(Data_to_model$Yr)),
mu_W0=subset(Arrival_weights, Yr %in% Lugworms$Yr,
select=c(Male_arrival_weight, Female_arrival_weight)),
sigma2_W0=unique(subset(Arrival_weights, Yr %in% Lugworms$Yr,
select=c(Male_arrival_weight_sigma, Female_arrival_weight_sigma))^2)[1,]
,
mu_T0=mu_T0,
sigma2_T0=sigma2_T0, Lugworms=Lugworms_vec)
inits=function() {
list(sigma2_alpha=runif(1, 0, 5),
sigma2_alpha_fem=runif(1, 0, 5),
K_sat_m=runif(1, 1, 13),
K_sat_f=runif(1, 1, 13))
}
jags.params=c("Alpha_t", "sigma2_alpha", "sigma2_alpha_fem", "K_sat_m", "K_sat_f")5.6 Run the model
if (Run_everything) {
# following run takes ~ 2 minutes brt
fueling_model<- jags(data=Data, inits, parameters=jags.params,
"fueling_model.jags", n.chains = 2, n.thin = 10,
n.iter = 10000, n.burnin = 1000, working.directory = getwd())
saveRDS(fueling_model, file='./results/fueling_model.RDS')
}
fueling_model<-readRDS('./results/fueling_model.RDS')5.7 Explore and save model results
print(fueling_model)## Inference for Bugs model at "fuelling_model.jags", fit using jags,
## 2 chains, each with 10000 iterations (first 1000 discarded), n.thin = 10
## n.sims = 1800 iterations saved
## mu.vect sd.vect 2.5% 25% 50% 75%
## Alpha_t[1,1] 5.419 0.174 5.093 5.297 5.420 5.541
## Alpha_t[2,1] 6.060 0.170 5.731 5.942 6.056 6.179
## Alpha_t[3,1] 6.113 0.182 5.760 5.990 6.108 6.236
## Alpha_t[4,1] 5.199 0.304 4.675 4.982 5.187 5.390
## Alpha_t[5,1] 4.191 0.287 3.658 3.990 4.187 4.377
## Alpha_t[6,1] 3.972 0.462 3.130 3.654 3.955 4.270
## Alpha_t[7,1] 4.243 0.321 3.635 4.027 4.240 4.450
## Alpha_t[8,1] 6.675 0.244 6.225 6.510 6.668 6.833
## Alpha_t[9,1] 5.905 0.280 5.370 5.707 5.901 6.098
## Alpha_t[10,1] 4.647 0.324 4.032 4.415 4.641 4.868
## Alpha_t[11,1] 5.787 0.259 5.310 5.616 5.776 5.955
## Alpha_t[12,1] 6.345 0.276 5.826 6.152 6.336 6.536
## Alpha_t[13,1] 4.752 0.346 4.120 4.507 4.742 4.971
## Alpha_t[14,1] 5.236 0.167 4.921 5.127 5.233 5.341
## Alpha_t[15,1] 4.941 0.262 4.458 4.762 4.921 5.109
## Alpha_t[16,1] 5.833 0.293 5.287 5.630 5.822 6.041
## Alpha_t[17,1] 6.207 0.215 5.796 6.060 6.206 6.355
## Alpha_t[18,1] 5.775 0.320 5.160 5.555 5.765 5.982
## Alpha_t[19,1] 5.689 0.163 5.384 5.580 5.687 5.797
## Alpha_t[20,1] 7.997 0.249 7.551 7.825 7.979 8.158
## Alpha_t[21,1] 6.419 0.413 5.643 6.140 6.407 6.681
## Alpha_t[1,2] 6.335 0.256 5.875 6.153 6.322 6.509
## Alpha_t[2,2] 6.268 0.207 5.874 6.128 6.261 6.406
## Alpha_t[3,2] 6.816 0.233 6.379 6.649 6.812 6.970
## Alpha_t[4,2] 5.761 0.331 5.141 5.533 5.754 5.965
## Alpha_t[5,2] 4.898 0.349 4.243 4.658 4.887 5.114
## Alpha_t[6,2] 4.590 0.508 3.676 4.234 4.559 4.909
## Alpha_t[7,2] 4.969 0.402 4.225 4.693 4.945 5.234
## Alpha_t[8,2] 7.263 0.267 6.752 7.086 7.252 7.432
## Alpha_t[9,2] 6.845 0.295 6.292 6.639 6.837 7.029
## Alpha_t[10,2] 4.884 0.347 4.251 4.642 4.865 5.117
## Alpha_t[11,2] 6.949 0.322 6.330 6.739 6.948 7.150
## Alpha_t[12,2] 8.276 0.431 7.490 7.983 8.276 8.542
## Alpha_t[13,2] 5.654 0.404 4.930 5.366 5.630 5.918
## Alpha_t[14,2] 5.284 0.256 4.818 5.108 5.273 5.448
## Alpha_t[15,2] 5.843 0.280 5.346 5.644 5.825 6.032
## Alpha_t[16,2] 6.624 0.280 6.099 6.423 6.624 6.807
## Alpha_t[17,2] 7.193 0.235 6.762 7.033 7.184 7.351
## Alpha_t[18,2] 6.962 0.340 6.322 6.732 6.949 7.190
## Alpha_t[19,2] 6.706 0.217 6.299 6.559 6.693 6.844
## Alpha_t[20,2] 8.271 0.361 7.611 8.014 8.248 8.516
## Alpha_t[21,2] 7.050 0.553 6.019 6.662 7.036 7.402
## K_sat_f 15.482 0.999 13.681 14.823 15.422 16.110
## K_sat_m 13.886 0.957 12.055 13.246 13.857 14.498
## sigma2_alpha 0.642 0.271 0.277 0.446 0.587 0.773
## sigma2_alpha_fem 0.776 0.339 0.325 0.541 0.704 0.941
## deviance 46809.988 103.587 46607.906 46737.682 46810.388 46880.728
## 97.5% Rhat n.eff
## Alpha_t[1,1] 5.770 1.008 300
## Alpha_t[2,1] 6.396 1.000 1800
## Alpha_t[3,1] 6.482 1.000 1800
## Alpha_t[4,1] 5.825 1.003 680
## Alpha_t[5,1] 4.777 1.026 63
## Alpha_t[6,1] 4.927 1.000 1800
## Alpha_t[7,1] 4.907 1.001 1800
## Alpha_t[8,1] 7.174 1.031 92
## Alpha_t[9,1] 6.461 1.002 1800
## Alpha_t[10,1] 5.295 1.002 1000
## Alpha_t[11,1] 6.308 1.001 1700
## Alpha_t[12,1] 6.904 1.005 1800
## Alpha_t[13,1] 5.485 1.004 440
## Alpha_t[14,1] 5.580 1.001 1800
## Alpha_t[15,1] 5.494 1.004 550
## Alpha_t[16,1] 6.413 1.005 1800
## Alpha_t[17,1] 6.636 1.000 1800
## Alpha_t[18,1] 6.428 1.012 130
## Alpha_t[19,1] 6.024 1.011 160
## Alpha_t[20,1] 8.533 1.000 1800
## Alpha_t[21,1] 7.273 1.002 1400
## Alpha_t[1,2] 6.840 1.005 380
## Alpha_t[2,2] 6.696 1.003 570
## Alpha_t[3,2] 7.273 1.008 200
## Alpha_t[4,2] 6.456 1.002 790
## Alpha_t[5,2] 5.627 1.002 1800
## Alpha_t[6,2] 5.690 1.002 1800
## Alpha_t[7,2] 5.798 1.013 120
## Alpha_t[8,2] 7.812 1.003 1800
## Alpha_t[9,2] 7.463 1.013 140
## Alpha_t[10,2] 5.605 1.004 460
## Alpha_t[11,2] 7.600 1.002 1800
## Alpha_t[12,2] 9.175 1.002 1200
## Alpha_t[13,2] 6.511 1.003 590
## Alpha_t[14,2] 5.819 1.003 740
## Alpha_t[15,2] 6.411 1.003 760
## Alpha_t[16,2] 7.189 1.002 1000
## Alpha_t[17,2] 7.666 1.000 1800
## Alpha_t[18,2] 7.647 1.003 1800
## Alpha_t[19,2] 7.157 1.004 1200
## Alpha_t[20,2] 9.041 1.001 1800
## Alpha_t[21,2] 8.169 1.000 1800
## K_sat_f 17.657 1.001 1800
## K_sat_m 15.785 1.000 1800
## sigma2_alpha 1.310 1.002 1100
## sigma2_alpha_fem 1.585 1.003 640
## deviance 47016.197 1.002 1000
##
## For each parameter, n.eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
##
## DIC info (using the rule, pD = var(deviance)/2)
## pD = 5362.8 and DIC = 52172.8
## DIC is an estimate of expected predictive error (lower deviance is better).Years<-unique(as.numeric((Data_to_model$Yr)))
# save estimated slopes for survival model
Fueling_slopes<-data.frame(Yr=Years,
Fueling_slope_males=fueling_model$BUGSoutput$mean$Alpha_t[,1],
Fueling_slope_females=fueling_model$BUGSoutput$mean$Alpha_t[,2],
Fueling_slope_males_predicted=11*Lugworms_vec/
(fueling_model$BUGSoutput$mean$K_sat_m[1]+Lugworms_vec),
Fueling_slope_females_predicted=13.2*Lugworms_vec/
(fueling_model$BUGSoutput$mean$K_sat_f[1]+Lugworms_vec))
write.csv(Fueling_slopes, file='./results/Fueling_slopes.csv', row.names=FALSE)6. Departure fuel load effects on survival from capture-recapture data
We used Cormack–Jolly–Seber (CJS) model (Lebreton et al., 1992) to model apparent survival independently from resighting probability. Set of nested models was estimated in program MARK (White & Burnham, 1999) through RMark (Laake, 2013) interface. To run model estimation with MARK or RMark you should first install program MARK from here.
6.1 Load and combine data
library('RMark')
#get arrivals to the Wadden Sea
Arrivals_WS<-read.csv('./results/Arrivals_migration_model.csv')
#get arrival mass
Arrival_weights<-read.csv('./results/Arrival_weights.csv')
#get annual slopes
Fueling_slopes<-read.csv('./results/Fueling_slopes.csv')
if (Use_local_data) {
# load main CR data
Full_CH<-read.csv('./data/Full_CH.csv', colClasses='character', row.names='code')
#get departures from Wadden Sea through arrival to Taimyr
Arrival_Taimyr<-read.csv('./data/Arrivals_Khat.csv')
} else {
#load main CR data
Full_CH<-read.csv('https://git.io/vXEFF', colClasses='character', row.names='code')
#get departures from Wadden Sea through arrival to Taimyr
Arrival_Taimyr<-read.csv('https://git.io/vXEbj')
}
Full_CH$SEX<-as.factor(Full_CH$SEX)
# a. get arrivals to the Wadden Sea
Arrivals_WS<-read.csv('./results/Arrivals_migration_model.csv')
# b. get arrival mass
Arrival_weights<-read.csv('./results/Arrival_weights.csv')
# c. get annual slopes
Fueling_slopes<-read.csv('./results/Fueling_slopes.csv')
Departures_Wadden<-Arrival_Taimyr[,c(1,3)]
Departures_Wadden$yday<-Departures_Wadden$yday-5.5
names(Departures_Wadden)[2]<-'departure'
# now we combine all these measurements.
All_staging_data<-merge(merge(merge(Arrivals_WS[,1:2],
Arrival_weights), Fueling_slopes),Departures_Wadden)
All_staging_data$Stopover<-(All_staging_data$departure-All_staging_data$arrival)
All_staging_data$log_fueling_slope_females_st<-
scale(log(All_staging_data$Fueling_slope_females), center=TRUE, scale=FALSE)
All_staging_data$log_fueling_slope_males_st<-
scale(log(All_staging_data$Fueling_slope_males), center=TRUE, scale=FALSE)
All_staging_data$log_Stopover<-log((All_staging_data$departure-All_staging_data$arrival))
Departure_Weight_gain_Females<- cbind( group='f', subset(All_staging_data,
select=c(Yr, log_fueling_slope_females_st, log_Stopover)))
names(Departure_Weight_gain_Females)[3]<-'log_fueling_slope'
names(Departure_Weight_gain_Females)[4]<-'log_Stopover'
Departure_Weight_gain_Males<- cbind(group='m', subset(All_staging_data,
select=c(Yr, log_fueling_slope_males_st, log_Stopover)))
names(Departure_Weight_gain_Males)[3]<-'log_fueling_slope'
names(Departure_Weight_gain_Males)[4]<-'log_Stopover'
Departure_Weight_gain<-rbind(Departure_Weight_gain_Females, Departure_Weight_gain_Males)
names(Departure_Weight_gain)[2]<-'time'
Departure_Weight_gain<-subset(Departure_Weight_gain, time>2001 & time<2016)6.2 Prepare data in RMark
library('RMark')
Godwit.data<-process.data(Full_CH, model = "CJS", begin.time = 2002,
groups=c("SEX"))
Godwit.ddl<-make.design.data(Godwit.data, remove.unused = FALSE)
Godwit.ddl$Phi<-merge_design.covariates(Godwit.ddl$Phi, Departure_Weight_gain,
bygroup=TRUE, bytime=TRUE)
# specify first year after marking birds
Godwit.ddl$Phi$two_plus=0
Godwit.ddl$Phi$two_plus[Godwit.ddl$Phi$Age>=1]=1
Godwit.ddl$p$two_plus=0
Godwit.ddl$p$two_plus[Godwit.ddl$p$Age>=1]=1
# clean from previous runs in case there have been any
if (exists('cml')) {
suppressWarnings(rm(list=cml[,1]))
suppressWarnings(rm(list=cml[,2]))
}
# specify model
# Survival part
Phi.dot=list(formula=~1)
Phi.1=list(formula=~two_plus + log_fueling_slope*SEX+log_Stopover*SEX)
Phi.2=list(formula=~two_plus + log_fueling_slope*SEX+log_Stopover+SEX)
Phi.3=list(formula=~two_plus + log_fueling_slope+log_Stopover*SEX)
Phi.4=list(formula=~two_plus + log_Stopover*SEX)
Phi.5=list(formula=~two_plus + log_fueling_slope*SEX)
Phi.8=list(formula=~two_plus + Time*SEX)
Phi.9=list(formula=~two_plus + SEX)
# Resighting part
p.time=list(formula=~time)
cml=create.model.list("CJS")6.3 Run MARK model
if (Run_everything) {
Res<-mark.wrapper(cml, data=Godwit.data, ddl=Godwit.ddl, profile.int=F, delete=FALSE, brief=TRUE)
saveRDS(Res, file='./results/survival_result.RDS')
} else {
Res<-readRDS('./results/survival_result.RDS')
}6.4 Median c-hat test in MARK
We now have to estimate potential overdispersion and correct for it. This possibility currently does not exist in RMark, so we export data to MARK and do test there.
export.MARK(Godwit.data, 'godwit_data', Res, replace=TRUE)Estimated c-hat = 1.1186307 with sampling SE = 0.0045798 95% Conf. Interval c-hat = 1.0623081 to 1.1749534
So we can adjust our result
Res_adj<-adjust.chat(Res, chat=1.1186307)
print(Res_adj)## model
## 3 Phi(~two_plus + log_fuelling_slope + log_Stopover * SEX)p(~time)
## 1 Phi(~two_plus + log_fuelling_slope * SEX + log_Stopover * SEX)p(~time)
## 4 Phi(~two_plus + log_Stopover * SEX)p(~time)
## 2 Phi(~two_plus + log_fuelling_slope * SEX + log_Stopover + SEX)p(~time)
## 6 Phi(~two_plus + Time * SEX)p(~time)
## 7 Phi(~two_plus + SEX)p(~time)
## 5 Phi(~two_plus + log_fuelling_slope * SEX)p(~time)
## 8 Phi(~1)p(~time)
## npar QAICc DeltaQAICc weight QDeviance chat
## 3 20 19709.74 0.000000 5.501511e-01 3554.558 1.118631
## 1 21 19711.32 1.589415 2.485105e-01 3554.137 1.118631
## 4 19 19711.99 2.253282 1.783150e-01 3558.822 1.118631
## 2 20 19716.24 6.501699 2.131356e-02 3561.060 1.118631
## 6 19 19721.28 11.547681 1.709759e-03 3568.116 1.118631
## 7 17 19756.11 46.379426 4.670008e-11 3606.966 1.118631
## 5 19 19759.64 49.906185 8.007402e-12 3606.475 1.118631
## 8 15 19816.70 106.965190 0.000000e+00 3671.568 1.1186317. Population trends from winter counts in West Africa
We modelled change in total population of godwits in state space modelled suggested by Kéry & Shcaub (2012), but with negative binomial error distribution.
7.1 Load packages and data
library(R2jags)
if (Use_local_data) {
Counts<-read.csv('/data/Counts_Mauritania.csv')
} else {
Counts<-read.csv('https://git.io/vXuL5')
}
Areas<-sort(unique(Counts$Area))
Years<-seq( min(unique(Counts$Yr)), max(unique(Counts$Yr)), 1)
y.matrix<-matrix(ncol=length(Areas), nrow= length(Years))
for (i in 1:nrow(Counts)) {
y.matrix[which(Years==Counts$Yr[i]), which(Areas==Counts$Area[i])]<-Counts$Count[i]
}7.2 Formulate model
{
sink("winter_counts_model.jags")
cat("
model {
# priors
for (s in 1:S) {N.est[1,s] ~ dexp(exp.lam)}
mean.lambda ~ dunif(0.5,1.5)
for (s in 1:S) {r[s] ~ dgamma(0.01, 0.01)}
for (t in 1:T) {n_total[t]<-sum(N.est[t,])}
#likelihood
# state process
for (t in 1:(T-1)) {
for (s in 1:S) {
N.est[t+1,s] <- N.est[t,s]*mean.lambda
}
}
# Observation process
for (t in 1:T) {
for (s in 1:S) {
p[t,s]<-r[s]/(r[s]+N.est[t,s])
y[t, s] ~ dnegbin(p[t,s], r[s])
}
}
}
", fill=TRUE)
sink()
}7.3 Priors and initial values
exp.mean<-mean(as.double(y.matrix),na.rm=T)
exp.lam <-1/exp.mean
jags.data<-list(y=y.matrix, T=length(Years), S=length(Areas), exp.lam=exp.lam)
# initials
inits<-function() {list(
N.est=matrix(ncol=6, nrow=15, c(runif(6, 0, 15000), rep(NA, 6*15-6)), byrow=TRUE),
mean.lambda=runif(1, 0.5,1.5),
r=runif(6,0,10))}
jags.params=c("n_total", "mean.lambda", "r")7.4 Run the model
if (Run_everything) {
# this model runs for ~ 2 mins brt
Counts_model<- jags(data=jags.data, inits, parameters=jags.params,
"winter_counts_model.jags", n.chains = 2, n.thin = 10, n.iter = 20000,
n.burnin = 5000, working.directory = getwd())
saveRDS(Counts_model, file='./results/Mauritania_counts_model.RDS')
}
Counts_model<-readRDS('./results/Mauritania_counts_model.RDS')
print(Counts_model)8. Phenology plot (Fig. 1)
8.1 Load data
# read Snow
Snowmelt<-read.csv('./results/Snow_remote.csv', as.is=TRUE)
# read arrivals to Wadden Sea (predicted by arrival model)
Arrivals_WZ<-read.csv('./results/Arrivals_migration_model.csv')
# read Mauritania counts model resutls
Counts_model<-readRDS('./results/Mauritania_counts_model.RDS')
if (Use_local_data) {
# read breeding data
Nests<-read.csv('./data/Nests.csv', as.is=TRUE)
# read crane flies
Tipula<-read.csv('./data/Tipula.csv', as.is=TRUE)
# read arrivals to Khatanga
Arr_Taimyr<-read.csv('./data/Arrivals_Khat.csv', as.is=TRUE)
# read raw counts
Counts<-read.csv('./data/Counts_Mauritania.csv', as.is=TRUE)
} else {
# read breeding data
Nests<-read.csv('https://git.io/vXu3y', as.is=TRUE)
# read crane flies
Tipula<-read.csv('https://git.io/vXu3x', as.is=TRUE)
# read arrivals to Khatanga
Arr_Taimyr<-read.csv('https://git.io/vXusf', as.is=TRUE)
# read raw counts
Counts<-read.csv('https://git.io/vXuL5')
}8.2 Merge data
tmp<-merge(Arrivals_WZ[,1:2], Snowmelt[,c(1,2,3)], all=TRUE)
names(tmp)[2:4]<-c('Arrivals_WZ', 'Snowmelt_global', 'Snowmelt_local')
tmp<-merge(tmp, Arr_Taimyr[,c(1,3)], all=TRUE)
names(tmp)[5]<-'Arr_T'
tmp<-merge(tmp, subset(Tipula, Region=='South', select=c('Yr', 'yday')), all=TRUE)
names(tmp)[6]<-'Tipula'
Data_all<-subset(tmp, Yr>1992 & Yr<2017)
Data_all$Ref_time<-(Data_all$Arr_T - Data_all$Arrivals_WZ-5.5)
Nests_m<-merge(Nests[Nests$Region=='South',1:3], Snowmelt[,c(1:3)], all.x=TRUE)
Areas<-sort(unique(Counts$Area))
Years<-seq( min(unique(Counts$Yr)), max(unique(Counts$Yr)), 1)
y.matrix<-matrix(ncol=length(Areas), nrow= length(Years))
for (i in 1:nrow(Counts)) {
y.matrix[which(Years==Counts$Yr[i]), which(Areas==Counts$Area[i])]<-Counts$Count[i]
}
Predicted_densities<-as.data.frame(cbind(Years,
t(apply(Counts_model$BUGSoutput$sims.list$n_total, 2,
FUN=function(x) quantile(x, probs =c(0.025, 0.5, 0.975))))))
names(Predicted_densities)[2:4]<-c("lcl", "median","ucl")8.3 Linear models results (slopes and P-values)
summary(lm(Snowmelt_global~Yr, data=Data_all))##
## Call:
## lm(formula = Snowmelt_global ~ Yr, data = Data_all)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.1589 -4.4375 0.2873 3.4019 10.3533
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1627.8680 310.1606 5.248 2.89e-05 ***
## Yr -0.7309 0.1547 -4.724 0.000103 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.247 on 22 degrees of freedom
## Multiple R-squared: 0.5035, Adjusted R-squared: 0.481
## F-statistic: 22.31 on 1 and 22 DF, p-value: 0.0001031summary(lm(Snowmelt_local~Yr, data=Data_all))##
## Call:
## lm(formula = Snowmelt_local ~ Yr, data = Data_all)
##
## Residuals:
## Min 1Q Median 3Q Max
## -16.7736 -4.9896 -0.4366 4.2523 15.9928
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1632.4638 439.5342 3.714 0.00121 **
## Yr -0.7337 0.2193 -3.346 0.00292 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.436 on 22 degrees of freedom
## Multiple R-squared: 0.3373, Adjusted R-squared: 0.3071
## F-statistic: 11.2 on 1 and 22 DF, p-value: 0.002924summary(lm(Tipula~Snowmelt_local, data=Data_all))##
## Call:
## lm(formula = Tipula ~ Snowmelt_local, data = Data_all)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.989 -2.611 -1.482 3.470 8.287
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 118.9308 22.6589 5.249 0.000123 ***
## Snowmelt_local 0.3769 0.1401 2.691 0.017559 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.081 on 14 degrees of freedom
## (8 observations deleted due to missingness)
## Multiple R-squared: 0.3409, Adjusted R-squared: 0.2939
## F-statistic: 7.242 on 1 and 14 DF, p-value: 0.01756summary(lm(Arr_T ~Snowmelt_local, data=Data_all))##
## Call:
## lm(formula = Arr_T ~ Snowmelt_local, data = Data_all)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.1558 -1.7813 -0.4854 1.7680 4.4497
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 116.31850 11.43359 10.173 3.99e-09 ***
## Snowmelt_local 0.21972 0.07077 3.105 0.00583 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.93 on 19 degrees of freedom
## (3 observations deleted due to missingness)
## Multiple R-squared: 0.3366, Adjusted R-squared: 0.3017
## F-statistic: 9.64 on 1 and 19 DF, p-value: 0.005832summary(lm(yday~jdate50_south, data=Nests_m))##
## Call:
## lm(formula = yday ~ jdate50_south, data = Nests_m)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.949 -4.061 -1.184 2.490 10.612
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 77.8544 26.6718 2.919 0.01396 *
## jdate50_south 0.5561 0.1648 3.374 0.00621 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.674 on 11 degrees of freedom
## Multiple R-squared: 0.5086, Adjusted R-squared: 0.4639
## F-statistic: 11.38 on 1 and 11 DF, p-value: 0.006206summary(lm(Arrivals_WZ ~ Snowmelt_global, data=Data_all))##
## Call:
## lm(formula = Arrivals_WZ ~ Snowmelt_global, data = Data_all)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.7601 -0.7849 0.5056 1.0911 3.3503
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 132.36203 9.55779 13.849 2.42e-12 ***
## Snowmelt_global -0.06564 0.05866 -1.119 0.275
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.049 on 22 degrees of freedom
## Multiple R-squared: 0.05385, Adjusted R-squared: 0.01084
## F-statistic: 1.252 on 1 and 22 DF, p-value: 0.2752summary(lm(Ref_time ~ Snowmelt_global, data=Data_all))##
## Call:
## lm(formula = Ref_time ~ Snowmelt_global, data = Data_all)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.4723 -3.0721 0.4625 2.2868 5.7913
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -29.6800 17.4629 -1.700 0.10552
## Snowmelt_global 0.3355 0.1078 3.112 0.00574 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.37 on 19 degrees of freedom
## (3 observations deleted due to missingness)
## Multiple R-squared: 0.3376, Adjusted R-squared: 0.3028
## F-statistic: 9.685 on 1 and 19 DF, p-value: 0.005748.4 Plot the figure
if (Save_pdf) pdf('Fig. 1.pdf', width=6, height=9, useDingbats=FALSE)
Snow_global_color<-'#8BAB8D'
Tipula_color<-'#8783CE'
Khat_arr_color<-'#7B3575'
Nest_color<-'#CC79A7'
Arr_WZ_col<-'#7B3575'
Color_BA<-'#7B3575'
Ref_col<-'#FB7508'
point_lwd=0.5
point_cex=2
font_cex=0.6
axes.labels.cex=1
x_marg=1
par(mar=c(4,4,2,2))
m <- cbind(c(1, 1,1,1, 3,3, 5,5,5), c(2,2,2, 2,4,4, 6,6,6))
m## [,1] [,2]
## [1,] 1 2
## [2,] 1 2
## [3,] 1 2
## [4,] 1 2
## [5,] 3 4
## [6,] 3 4
## [7,] 5 6
## [8,] 5 6
## [9,] 5 6 layout(m)
#######################
# first will be panel with snow over years.
plot(Snowmelt_global~Yr, data=Data_all, type='n', ylim=c(140+3, 195-5),
las=1,axes=FALSE, ylab='', xlab='', xlim=c(1992, 2017))
box()
abline(h=121, col=grey(0.5))
abline(h=121-7.5, col=grey(0.5), lty=3)
abline(h=121+7.5, col=grey(0.5), lty=3)
abline(h=121, col=grey(0.5))
abline(h=152, col=grey(0.5))
abline(h=182, col=grey(0.5))
abline(h=152+15, col=grey(0.5), lty=2)
abline(h=182+7.5, col=grey(0.5), lty=3)
abline(h=152+15+7.5, col=grey(0.5), lty=3)
abline(h=152+7.5, col=grey(0.5), lty=3)
abline(h=121+15+7.5, col=grey(0.5), lty=3)
abline(v=1995, col=grey(0.5))
abline(v=2015, col=grey(0.5))
abline(v=2005, col=grey(0.5), lty=2)
abline(v=2000, col=grey(0.5), lty=3)
abline(v=2010, col=grey(0.5), lty=3)
axis(1)
mtext(side=1, 'year', line=2.5, cex=axes.labels.cex*0.8)
mtext(side=2, 'date of snowmelt (Taimyr)', line=1, cex=axes.labels.cex*0.8)
rect(1994, 152.15, 1996, 154.2, col='white', border=NA)
text(y=153.2, x=1995, '1 June', cex=axes.labels.cex)
rect(1994, 182.15, 1996, 184.2, col='white', border=NA)
text(y=183.2, x=1995, '1 July', cex=axes.labels.cex)
abline(h=152, col=grey(0.5))
abline(h=182, col=grey(0.5))
Snow_over_years_lm<-lm(Snowmelt_global~Yr, data=Data_all)
abline(Snow_over_years_lm, lwd=2)
# plot prediction lines..
XX_sn<-seq(1992, 2017, by=1)
Pred_sn<-predict(Snow_over_years_lm, newdata=data.frame(Yr=XX_sn), se.fit=TRUE)
Pred_sn_UCI<-Pred_sn$fit+1.96*Pred_sn$se.fit
Pred_sn_LCI<-Pred_sn$fit-1.96*Pred_sn$se.fit
lines(Pred_sn_UCI~XX_sn, lty=2)
lines(Pred_sn_LCI~XX_sn, lty=2)
for (i in 1:nrow(Data_all)) {
points(Snowmelt_global~Yr, data=Data_all[i,], pch=21,
cex=point_cex, bg=Snow_global_color, col=grey(0.4), lwd=point_lwd)
text(y=Data_all$Snowmelt_global[i], x=Data_all$Yr[i],
labels=substr(Data_all$Yr[i], 3,4) , col='white', cex=font_cex)
}
box()
#########################################
# Now panel with Taimyr phenology.
X_range<-range(Data_all$Snowmelt_local)
plot(Tipula~Snowmelt_local, data=Data_all, type='n', ylim=c(140+3, 195-5),
las=1, xlim=c(X_range[1]-x_marg, X_range[2]+x_marg),axes=FALSE, ylab='', xlab='')
box()
abline(h=121, col=grey(0.5))
abline(h=152+15, col=grey(0.5), lty=2)
abline(h=182+7.5, col=grey(0.5), lty=3)
abline(h=152+15+7.5, col=grey(0.5), lty=3)
abline(h=152+7.5, col=grey(0.5), lty=3)
abline(h=121+15+7.5, col=grey(0.5), lty=3)
abline(v=121, col=grey(0.5))
abline(v=152, col=grey(0.5))
abline(v=182, col=grey(0.5))
abline(v=121+15, col=grey(0.5), lty=2)
abline(v=152+15, col=grey(0.5), lty=2)
abline(v=182+7.5, col=grey(0.5), lty=3)
abline(v=152+15+7.5, col=grey(0.5), lty=3)
abline(v=152+7.5, col=grey(0.5), lty=3)
abline(v=121+15+7.5, col=grey(0.5), lty=3)
axis(1, at=c(152, 152+15, 182), labels=c('1 June', '15 June', '1 July'),
cex.axis=axes.labels.cex)
mtext(side=1, 'date of snowmelt (SE Taimyr)', line=2.5, cex=axes.labels.cex*0.8)
mtext(side=2, 'Taimyr phenology', line=1, cex=axes.labels.cex*0.8)
# and lables for them
rect(180, 152.15+15, 185, 154.2+15, col='white', border=NA)
text(y=153.2+15, x=182, '15 June', cex=axes.labels.cex)
rect(180, 152.15, 185, 154.2, col='white', border=NA)
text(y=153.2, x=182, '1 June', cex=axes.labels.cex)
rect(180, 182.15, 185, 184.4, col='white', border=NA)
text(y=183.2, x=182, '1 July', cex=axes.labels.cex)
abline(h=152, col=grey(0.5))
abline(h=182, col=grey(0.5))
lm_snow<-lm(Snowmelt_global~Snowmelt_local, data=Data_all)
abline(lm_snow, col=Snow_global_color, lwd=2)
# Tipula
lm_tipula<-lm(Tipula~Snowmelt_local, data=Data_all)
abline(lm_tipula, col=Tipula_color, lwd=2)
# Clutches
lm_nests<-lm(yday~jdate50_south, data=Nests_m)
abline(lm_nests, lwd=2, col=Nest_color)
# Arrivals
abline(lm(Arr_T~Snowmelt_local, data=subset(Data_all, Yr>1992)),
col=Khat_arr_color, lwd=2)
for (i in 1:nrow(Data_all)) {
points(Snowmelt_global~Snowmelt_local, data=Data_all[i,], pch=21,
bg=Snow_global_color, cex=point_cex, col=grey(0.5), lwd=point_lwd )
text(y=Data_all$Snowmelt_global[i], x=Data_all$Snowmelt_local[i],
labels=substr(Data_all$Yr[i], 3,4) , col='white', cex=font_cex)
}
for (i in 1:nrow(Data_all)) {
# Tipula
points(Tipula~Snowmelt_local, data=Data_all[i,], pch=21, bg=Tipula_color,
lwd=point_lwd, cex=point_cex, col=grey(0.5))
text(y=Data_all$Tipula[i], x=Data_all$Snowmelt_local[i],
labels=substr(Data_all$Yr[i], 3,4) , col='white', cex=font_cex)
# Arrivals
points(Arr_T~Snowmelt_local, data=subset(Data_all, Yr>=1992)[i,],
pch=21, cex=point_cex, col=grey(0.4), bg=Khat_arr_color, lwd=point_lwd)
text(y=Data_all$Arr_T[i], x=Data_all$Snowmelt_local[i],
labels=substr(Data_all$Yr[i], 3,4) , col='white', cex=font_cex)
}
# nests
for (i in 1:nrow(Nests_m)) {
points(yday~jdate50_south, data=Nests_m[i,], pch=21, lwd=point_lwd,
bg=Nest_color, cex=point_cex, col=grey(0.5))
text(y=Nests_m$yday[i], x=Nests_m$jdate50_south[i],
labels=substr(Nests_m$Yr[i], 3,4) , col='white',cex=font_cex)
}
box()
############################
# this will be panel B
plot(Arrivals_WZ~Snowmelt_global, data=Data_all, type='n', ylim=c(111, 135.5-5),
las=1, xlim=c(X_range[1]-x_marg, X_range[2]+x_marg),axes=FALSE, ylab='', xlab='')
box()
abline(h=121, col=grey(0.5))
abline(h=121-7.5, col=grey(0.5), lty=3)
abline(h=121+7.5, col=grey(0.5), lty=3)
abline(v=152, col=grey(0.5))
abline(v=182, col=grey(0.5))
abline(v=152+15, col=grey(0.5), lty=2)
abline(v=182+7.5, col=grey(0.5), lty=3)
abline(v=152+15+7.5, col=grey(0.5), lty=3)
abline(v=152+7.5, col=grey(0.5), lty=3)
abline(v=121+15+7.5, col=grey(0.5), lty=3)
axis(1, at=c(152, 152+15, 182), labels=c('1 June', '15 June', '1 July'),
cex.axis=axes.labels.cex)
mtext(side=1, 'date of snowmelt (Taimyr)', line=2.5, cex=axes.labels.cex*0.8)
mtext(side=2, 'arrival at Wadden Sea', line=1, cex=axes.labels.cex*0.8)
Arrivals_WZ_lm<-lm(Arrivals_WZ~Snowmelt_global, data=Data_all)
abline(Arrivals_WZ_lm, lwd=2)
XX<-seq(X_range[1]-x_marg, X_range[2]+x_marg, by=1)
Pred<-predict(Arrivals_WZ_lm, newdata=data.frame(Snowmelt_global=XX), se.fit=TRUE)
Pred_UCI<-Pred$fit+1.96*Pred$se.fit
Pred_LCI<-Pred$fit-1.96*Pred$se.fit
lines(Pred_UCI~XX, lty=2)
lines(Pred_LCI~XX, lty=2)
rect(178, 121.15, 185, 123.6, col='white', border=NA)
text(y=122.35, x=182, '1 May', cex=axes.labels.cex)
abline(h=121, col=grey(0.5))
for (i in 1:nrow(Data_all)) {
points(Arrivals_WZ~Snowmelt_global, data=Data_all[i,], pch=21, bg=Arr_WZ_col,
col=grey(0.4), lwd=point_lwd, cex=point_cex)
text(y=Data_all$Arrivals_WZ[i], x=Data_all$Snowmelt_globa[i],
labels=substr(Data_all$Yr[i], 3,4) , col='white', cex=font_cex)
}
box()
############################
# this will be panel D
plot(Ref_time~Snowmelt_global, data=Data_all, type='n', las=1, xlim=c(X_range[1]-x_marg,
X_range[2]+x_marg),axes=FALSE, xlab='', ylab='', ylim=c(13.5-0.5+2, 37-4.5+2))
box()
axis(2, las=1)
mtext(side=2, 'refueling time (d)', line=2.5, cex=axes.labels.cex*0.8)
axis(1, at=c(152, 152+15, 182), labels=c('1 June', '15 June', '1 July'),
cex.axis=axes.labels.cex)
mtext(side=1, 'date of snowmelt (Taimyr)', line=2.5, cex=axes.labels.cex*0.8)
abline(v=152, col=grey(0.5))
abline(v=182, col=grey(0.5))
abline(v=152+15, col=grey(0.5), lty=2)
abline(h=25, col=grey(0.5))
abline(h=25-7.5, col=grey(0.5), lty=3)
abline(h=25+7.5, col=grey(0.5), lty=3)
abline(v=182+7.5, col=grey(0.5), lty=3)
abline(v=152+15+7.5, col=grey(0.5), lty=3)
abline(v=152+7.5, col=grey(0.5), lty=3)
abline(v=121+15+7.5, col=grey(0.5), lty=3)
Ref_time_lm<-lm(Ref_time~Snowmelt_global, data=Data_all)
abline(Ref_time_lm, lwd=2)
# plot prediction lines..
XX<-seq(X_range[1]-x_marg, X_range[2]+x_marg, by=1)
Pred_ref<-predict(Ref_time_lm, newdata=data.frame(Snowmelt_global=XX), se.fit=TRUE)
Pred_ref_UCI<-Pred_ref$fit+1.96*Pred_ref$se.fit
Pred_ref_LCI<-Pred_ref$fit-1.96*Pred_ref$se.fit
lines(Pred_ref_UCI~XX, lty=2)
lines(Pred_ref_LCI~XX, lty=2)
for (i in 1:nrow(Data_all)) {
points(Ref_time~Snowmelt_global, data=Data_all[i,], pch=21, cex=point_cex,
bg=Ref_col, col=grey(0.4), lwd=point_lwd)
text(y=Data_all$Ref_time[i], x=Data_all$Snowmelt_globa[i],
labels=substr(Data_all$Yr[i], 3,4) , col='white', cex=font_cex)
}
# panel E
plot(median~Years, data=Predicted_densities, ylim=c(7000, 27000),
type='n', las=1, ylab='', axes=FALSE, xlab='', xlim=c(1992, 2017))
abline(h=15000, col=grey(0.5), lwd=1)
abline(h=25000, col=grey(0.5), lty=2)
abline(h=c(10000, 20000), col=grey(0.5), lty=3)
abline(h=15000, col=grey(0.5), lty=2)
abline(v=2005, lty=2, col=grey(0.5))
abline(v=1995, col=grey(0.5))
abline(v=2015, col=grey(0.5))
abline(v=2005, col=grey(0.5), lty=2)
abline(v=2000, col=grey(0.5), lty=3)
abline(v=2010, col=grey(0.5), lty=3)
axis(1)
mtext(side=1, 'year', line=2.5, cex=axes.labels.cex*0.8)
axis(2, at=seq(10000, 25000, by=5000), label=seq(10000, 25000, by=5000)/1000, las=1)
mtext(side=2, 'count (x1000)', line=2.5, cex=axes.labels.cex*0.8)
lines(x=Predicted_densities$Years, Predicted_densities$ucl, lty=2)
lines(x=Predicted_densities$Years, Predicted_densities$lcl, lty=2)
lines(median~Years, data=Predicted_densities, lwd=2)
points(rowSums(y.matrix)~ Years, pch=21, cex=point_cex, bg=Color_BA,
col=grey(0.4), lwd=point_lwd)
text(y=rowSums(y.matrix), x=Years, labels=substr(Years, 3,4), col='white', cex=font_cex)
box()
if (Save_pdf) dev.off()
9. Stayover duration - lugworms - godwit survival (Fig. 2)
9.1 Load data
fueling_model<-readRDS('./results/fueling_model.RDS')
# read arrivals to Wadden Sea (predicted by arrival model)
Arr_Wadden<-read.csv(file='./results/Arrivals_migration_model.csv')
Fueling_slopes<-read.csv('./results/Fueling_slopes.csv')
Res<-readRDS('./results/survival_result.RDS')
Res_adj<-adjust.chat(Res, chat=1.1186307)
Model<-Res_adj[[3]]
if (Use_local_data) {
Lugworms<-read.csv('./data/Lugworms.csv')
} else {
Lugworms<-read.csv('https://git.io/vXESR')
}
# No trend in lugworms:
summary(lm(lugworms_ad~Yr, data=subset(Lugworms, Yr>1995)))##
## Call:
## lm(formula = lugworms_ad ~ Yr, data = subset(Lugworms, Yr > 1995))
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.3739 -2.5079 0.4706 2.5608 6.5419
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -291.9543 280.4488 -1.041 0.311
## Yr 0.1529 0.1398 1.094 0.288
##
## Residual standard error: 3.879 on 19 degrees of freedom
## Multiple R-squared: 0.05922, Adjusted R-squared: 0.009708
## F-statistic: 1.196 on 1 and 19 DF, p-value: 0.2878# Average density:
summary(lm(lugworms_ad~1, data=subset(Lugworms, Yr>1995)))##
## Call:
## lm(formula = lugworms_ad ~ 1, data = subset(Lugworms, Yr > 1995))
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.9152 -3.4252 0.0448 2.0348 7.4248
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.7552 0.8507 17.34 1.61e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.898 on 20 degrees of freedom9.2 Prepare data
Arrivals_Khat$Date<-as.Date(Arrivals_Khat$Date, tz='GMT', format='%Y-%m-%d' )
Arrivals_Khat$Departure<-Arrivals_Khat$yday-5.5
names(Arr_Wadden)[2]<-'yday'
Stayovers<-merge(Arr_Wadden,
Arrivals_Khat[,-which(names(Arrivals_Khat)=='yday')], by='Yr')
Stayovers$Stayover<- Stayovers$Departure- Stayovers$yday
Stayovers<-merge(Stayovers, subset(Lugworms, Yr>1995), all.x=TRUE)
Stayovers<-merge(Stayovers, Fueling_slopes, all=TRUE)Create grid for prediction
Scale_females<-scale(log(All_staging_data$Fueling_slope_females))
Scale_males<-scale(log(All_staging_data$Fueling_slope_males))
Prediction.matrix<-expand.grid(log_fueling_slope=seq(1.25-1.717, 2.35-1.717, by=0.001),
log_Stopover=seq(2.75, 3.6, by=0.001))
Design_females<-matrix(nrow=nrow(Prediction.matrix), ncol=20, data=0)
Design_females[,1]<-1
Design_females[,2]<-1
Design_females[,3]<-Prediction.matrix$log_fueling_slope
Design_females[,4]<-Prediction.matrix$log_Stopover
Design_males<-matrix(nrow=nrow(Prediction.matrix), ncol=20, data=0)
Design_males[,1]<-1
Design_males[,2]<-1
Design_males[,3]<-Prediction.matrix$log_fueling_slope
Design_males[,4]<-Prediction.matrix$log_Stopover
Design_males[,5]<-1
Design_males[,6]<-Prediction.matrix$log_Stopover
Prediction.matrix$Prediction_females<-
plogis(Design_females%*%as.matrix(summary(Model)$beta[,1, drop=FALSE]))
Prediction.matrix$Prediction_males<-
plogis(Design_males%*%as.matrix(summary(Model)$beta[,1, drop=FALSE]))
Prediction.matrix$fueling_slope_females<-
exp(Prediction.matrix$log_fueling_slope+attr(Scale_females,"scaled:center"))
Prediction.matrix$Stopover<-exp(Prediction.matrix$log_Stopover)
Prediction.matrix$fueling_slope_males<-
exp(Prediction.matrix$log_fueling_slope+attr(Scale_males,"scaled:center"))
M1<-(lm(Stayover~Yr, data=Stayovers))
Stopover_before_now<-predict(M1, newdata=data.frame(Yr=c(1995, 2016)))
M3<-lm(Fueling_slope_males~Yr, data=Stayovers)
Slope_males_before_now<-predict(M3, newdata=data.frame(Yr=c(1995, 2016)))
M4<-lm(Fueling_slope_females~Yr, data=Stayovers)
Slope_females_before_now<-predict(M4, newdata=data.frame(Yr=c(1995, 2016)))
M5<-(lm(Stayover~Fueling_slope_females, data=Stayovers))
M6<-(lm(Stayover~Fueling_slope_males, data=Stayovers))
M7<-lm(lugworms_ad~Yr, data=Lugworms)
Design_females_s<-matrix(nrow=2, ncol=20, data=0)
Design_females_s[,1]<-1
Design_females_s[,2]<-1
Design_females_s[,3]<-log(Slope_females_before_now)-attr(Scale_females,"scaled:center")
Design_females_s[,4]<-log(Stopover_before_now)
Design_males_s<-matrix(nrow=2, ncol=20, data=0)
Design_males_s[,1]<-1
Design_males_s[,2]<-1
Design_males_s[,3]<-log(Slope_males_before_now)-attr(Scale_males,"scaled:center")
Design_males_s[,4]<-log(Stopover_before_now)
Design_males_s[,5]<-1
Design_males_s[,6]<-log(Stopover_before_now)9.3 Female plot
First I plot female figure (Figure 2) and then male figure for Supplementary Materials Fig S5.
if (Save_pdf) pdf(file='Figure 2.pdf', useDingbats=FALSE, width=6, height=9)
par(mfrow=c(3,2))
par(mar=c(4, 5,2,2 ))
# set up colors
Late_color<-'#CD6632'
Early_color<-'#6FC067'
Colors<-colorRampPalette(c("#b2182b",
"#d6604d",
"#f4a582",
"#fddbc7",
"#f7f7f7",
"#d1e5f0",
"#92c5de",
"#4393c3",
"#2166ac"), bias=0.5)
axes.labels.cex=1.2
# Panel 1 - stayover change over years
plot(Stayovers$Stayover~Stayovers$Yr, type='p', ylim=c(17, 35), ylab='',
las=1, xlim=c(1992, 2018), xlab='', cex=2, axes=F, yaxs='i', pch=21,
col='white', bg='black')
axis(1, cex.axis=axes.labels.cex)
axis(4, labels=FALSE)
par(xpd=TRUE)
mtext(side=1, 'year', line=3, cex=axes.labels.cex)
par(xpd=FALSE)
#axis(4, labels=FALSE)
axis(2, cex.axis=axes.labels.cex, las=1)
mtext(side=2, 'refueling time (d)', line=3,cex=axes.labels.cex,)
abline(lm(Stayover~Yr, data=Stayovers), lwd=2)
# 95%CI
XX<-seq(1992,2018, by=0.1)
Pred<-predict(M1, newdata=data.frame(Yr=XX), se.fit=TRUE)
UCI<-Pred$fit+1.96*Pred$se.fit
LCI<-Pred$fit-1.96*Pred$se.fit
lines(UCI~XX, lty=2)
lines(LCI~XX, lty=2)
box()
points(Stopover_before_now[1]~I(1995), pch=21, bg=Early_color, col='white', cex=4, lwd=3)
points(Stopover_before_now[2]~I(2016), pch=22, col=grey(0.45), bg=Late_color, cex=4, lwd=3)
#---------------------------------------
# panel 2
Prediction.matrix_cur<-subset(Prediction.matrix, Stopover>=17 &
Stopover<=35 &
fueling_slope_females<=10.5 &
fueling_slope_females>=4)
Female_Image<-as.image(x=Prediction.matrix_cur[,c(5,6)], Z=Prediction.matrix_cur$Prediction_females)
plot(0,0, type='n', xaxt='n', yaxt='n', xlab='', ylab='', bty='n') # make an empty plot
par(xpd=TRUE)
mtext("refueling rate (g/d)", side=1, line=3, cex=axes.labels.cex)
par(xpd=FALSE)
par(new=T)
par(mar=c(4, 5,2,2 ))
image(Female_Image,
col=Colors(256), xlab='', ylab='', las=1, ylim=c(17,35), xlim=c(4,10.5),
cex.axis=axes.labels.cex)
contour(Female_Image,add=T, levels=plogis(Design_females_s%*%as.matrix(summary(Model)$beta[,1, drop=FALSE]))[1],
lwd=3, col=grey(0.5),
labels=round(plogis(Design_females_s%*%as.matrix(summary(Model)$beta[,1, drop=FALSE]))[1],2),
labcex = 1)
contour(Female_Image,add=T, levels=c(0.8, 0.85, 0.9, 0.95), col=grey(0.5), labcex=1)
points(Stayover~Fueling_slope_females, data=Stayovers, pch=21, col='white', bg='black', cex=2)
box()
ScaleX<-seq(Slope_females_before_now[1], 9.65, length.out=11)
points(rep(Stopover_before_now[2], 11)~ScaleX, pch='|')
points(Stopover_before_now[1]~Slope_females_before_now[1], pch=21, bg=Early_color, col='white', cex=4, lwd=3)
points(Stopover_before_now[2]~Slope_females_before_now[1], pch=22, bg='white',col=Late_color, cex=4, lwd=3)
points(Stopover_before_now[2]~c(9.65), pch=22, bg='black', col=Late_color, cex=4, lwd=3)
points(Stopover_before_now[2]~Slope_females_before_now[2], pch=22, col=grey(0.45), bg=Late_color, cex=4, lwd=3)
##-------------------------------------
# Now panel 3
par(mar=c(4, 5,2,2 ))
plot(Stayovers$Fueling_slope_females~Stayovers$Yr, type='p', ylab='',
las=1, xlim=c(1992, 2018), xlab='', axes=F, yaxs='i', cex=2, pch=21,
bg='black', col='white', las=1, ylim=c(4,10.5))
axis(1, cex.axis=axes.labels.cex)
axis(3, labels=FALSE)
par(xpd=TRUE)
mtext(side=1, 'year', line=3, cex=axes.labels.cex)
par(xpd=FALSE)
axis(2, las=1, cex.axis=axes.labels.cex)
mtext(side=2, 'refueling rate (g/d)', line=3,cex=axes.labels.cex)
abline(lm(Fueling_slope_females~Yr, data=Stayovers), lwd=2)
# 95%CI
XX<-seq(1992,2018, by=0.1)
Pred<-predict(M4, newdata=data.frame(Yr=XX), se.fit=TRUE)
UCI<-Pred$fit+1.96*Pred$se.fit
LCI<-Pred$fit-1.96*Pred$se.fit
lines(UCI~XX, lty=2)
lines(LCI~XX, lty=2)
box()
points(Slope_females_before_now[1]~I(1995), pch=21, bg=Early_color, col='white', cex=4, lwd=3)
points(Slope_females_before_now[2]~I(2016), pch=22, col=grey(0.45), bg=Late_color, cex=4, lwd=3)
ScaleY<-seq(Slope_females_before_now[1], 9.65, length.out=11)
points(ScaleY~rep(2016, 11), pch=151)
points(Stopover_before_now[1]~Slope_females_before_now[1], pch=21, bg=Early_color, col='white', cex=4, lwd=3)
points(Slope_females_before_now[1]~I(2016), pch=22, bg='white',col=Late_color, cex=4, lwd=3)
points(I(9.65)~I(2016), pch=22, bg='black', col=Late_color, cex=4, lwd=3)
points(Stayovers$Fueling_slope_females~Stayovers$Yr, cex=2, pch=21, bg='black', col='white')
##-------------------------------------
# Now panel 4
par(mar=c(4, 5,2,2 ))
par(xaxs='i')
plot(Stayovers$Stayover~Stayovers$Fueling_slope_females, type='p', ylim=c(17,35), xlim=c(4,10.5) ,ylab='',
las=1, xlab='', cex=2, axes=F, yaxs='i', pch=21, col='white', bg='black')
axis(1, cex.axis=axes.labels.cex)
axis(3, labels=FALSE)
par(xpd=TRUE)
mtext(side=1, 'refueling rate (g/d)', line=3, cex=axes.labels.cex)
par(xpd=FALSE)
axis(2,cex.axis=axes.labels.cex, las=1)
mtext(side=2, 'refueling time (d)', line=3,cex=axes.labels.cex)
abline(lm(Stayover~Fueling_slope_females, data=Stayovers), lwd=2)
# 95%CI
XX<-seq(4,10.5, by=0.1)
Pred<-predict(M5, newdata=data.frame(Fueling_slope_females=XX), se.fit=TRUE)
UCI<-Pred$fit+1.96*Pred$se.fit
LCI<-Pred$fit-1.96*Pred$se.fit
lines(UCI~XX, lty=2)
lines(LCI~XX, lty=2)
box()
points(Stopover_before_now[1]~Slope_females_before_now[1], pch=21, bg=Early_color, col='white', cex=4, lwd=3)
points(Stopover_before_now[2]~Slope_females_before_now[2], pch=22, col=grey(0.45), bg=Late_color, cex=4, lwd=3)
##-------------------------------------
# Now panel 5
par(mar=c(4, 5,2,2 ))
par(xaxs='r')
plot(lugworms_ad~Yr, data=Lugworms, type='p', ylab='',
las=1, xlim=c(1992, 2018), xlab='', axes=F, yaxs='i',
cex=2, pch=21, bg='black', col='white',
ylim=c(min(Lugworms$lugworms_ad, na.rm=TRUE)-0.5, max(Lugworms$lugworms_ad, na.rm=TRUE)+0.5))
axis(1, cex.axis=axes.labels.cex)
axis(3, labels=FALSE)
par(xpd=TRUE)
mtext(side=1, 'year', line=3, cex=axes.labels.cex)
par(xpd=FALSE)
axis(2, las=1, cex.axis=axes.labels.cex)
mtext(side=2, 'lugworm density (ind./sq.m)', line=3,cex=axes.labels.cex)
abline(lm(lugworms_ad~Yr, data=Lugworms), lwd=2)
# 95%CI
XX<-seq(1992,2018, by=0.1)
Pred<-predict(M7, newdata=data.frame(Yr=XX), se.fit=TRUE)
UCI<-Pred$fit+1.96*Pred$se.fit
LCI<-Pred$fit-1.96*Pred$se.fit
lines(UCI~XX, lty=2)
lines(LCI~XX, lty=2)
box()
##-------------------------------------
# Now panel 6
par(mar=c(4, 5,2,2 ))
plot(fueling_model$BUGSoutput$mean$Alpha_t[,2]~Lugworms$lugworms_ad,
ylim=c(4,8.5),ylab='',
xlab='' ,las=1, type='n', axes=FALSE)
axis(2, cex.axis=axes.labels.cex, las=1)
mtext(side=2, 'refueling rate (g/d)', line=3, cex=axes.labels.cex)
axis(1, las=1, cex.axis=axes.labels.cex)
par(xpd=TRUE)
mtext(side=1, 'lugworm density (ind./sq.m)', line=3,cex=axes.labels.cex)
par(xpd=FALSE)
# 95%CI
XX<-seq(min(Lugworms$lugworms_ad), max(Lugworms$lugworms_ad), 0.1)
Y_fem<- t(sapply(XX, FUN=function(x)
quantile(x*13.2 /(x+fueling_model$BUGSoutput$sims.list$K_sat_f),
probs=c(0.025, 0.5,0.975))))
lines(Y_fem[,1]~XX, lty=2)
lines(Y_fem[,2]~XX, lwd=2)
lines(Y_fem[,3]~XX, lty=2)
points(fueling_model$BUGSoutput$mean$Alpha_t[,2]~Lugworms$lugworms_ad,
pch=21, cex=2, bg='black', col='white')
box()
#--------------------------------------------------------------------------
if (Save_pdf) dev.off()9.4 Check for temporal autocorrelation
# stayover duration
library(nlme)
M1<-lm(Stayover~Yr, data=Stayovers)
Stayovers$resid_stayover<-NA
Stayovers$resid_stayover[!is.na(Stayovers$Stayover)]<-residuals(M1)
acf(Stayovers$resid_stayover, na.action = na.pass)
M1_gls<-gls(Stayover~Yr, data=Stayovers, na.action = na.omit)
M1_gls_AR1<-gls(Stayover~Yr, data=Stayovers, na.action = na.omit, correlation = corAR1(form = ~ Yr))
AIC(M1_gls, M1_gls_AR1)## df AIC
## M1_gls 3 120.5769
## M1_gls_AR1 4 122.2908##########
# refueling rate
M2<-lm(Fueling_slope_females~Yr, data=Stayovers)
Stayovers$resid_fsf<-NA
Stayovers$resid_fsf[!is.na(Stayovers$Fueling_slope_females)]<-residuals(M2)
acf(Stayovers$resid_fsf, na.action = na.pass)
M2_gls<-gls(Fueling_slope_females~Yr, data=Stayovers, na.action = na.omit)
M2_gls_AR1<-gls(Fueling_slope_females~Yr, data=Stayovers, na.action = na.omit, correlation=corAR1(form= ~Yr))
AIC(M2_gls, M2_gls_AR1)## df AIC
## M2_gls 3 68.96932
## M2_gls_AR1 4 69.51784###########
# lugworms
M3<-lm(lugworms_ad~Yr, data=Stayovers)
summary(M3)##
## Call:
## lm(formula = lugworms_ad ~ Yr, data = Stayovers)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.3713 -2.7367 0.7778 2.6557 6.5469
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -294.3209 298.3152 -0.987 0.337
## Yr 0.1541 0.1487 1.036 0.314
##
## Residual standard error: 3.986 on 18 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.0563, Adjusted R-squared: 0.003869
## F-statistic: 1.074 on 1 and 18 DF, p-value: 0.3138Stayovers$resid_lug<-NA
Stayovers$resid_lug[!is.na(Stayovers$lugworms_ad)]<-residuals(M3)
acf(Stayovers$resid_lug, na.action = na.pass)
M3_gls<-gls(lugworms_ad~Yr, data=Stayovers, na.action = na.omit)
M3_gls_AR1<-gls(lugworms_ad~Yr, data=Stayovers, na.action = na.omit, correlation=corAR1(form = ~Yr))
AIC(M3_gls, M3_gls_AR1)## df AIC
## M3_gls 3 116.4316
## M3_gls_AR1 4 117.0185########################################
# fueling by lugworms
M4<-lm(Fueling_slope_females~lugworms_ad, data=Stayovers)
Stayovers$resid_M4<-NA
Stayovers$resid_M4[!is.na(Stayovers$lugworms_ad) & !is.na(Stayovers$Fueling_slope_females)]<-residuals(M4)
acf(Stayovers$resid_M4, na.action = na.pass)
M4_gls<-gls(Fueling_slope_females~lugworms_ad, data=Stayovers, na.action = na.omit)
M4_gls_AR1<-gls(Fueling_slope_females~lugworms_ad, data=Stayovers,
na.action = na.omit, correlation = corAR1(form = ~ Yr))
AIC(M4_gls, M4_gls_AR1)## df AIC
## M4_gls 3 53.39698
## M4_gls_AR1 4 55.30244###############
# autocorrelation never improved model fit.9.5 Male plot
Supplementary Materials Fig S5.
if (Save_pdf) pdf(file='Males_plot.pdf', useDingbats=FALSE, width=6, height=9)
par(mfrow=c(3,2))
par(mar=c(4, 5,2,2 ))
# set up colors
Late_color<-'#CD6632'
Early_color<-'#6FC067'
axes.labels.cex=1.2
# Panel 1 - stayover change over years
plot(Stayovers$Stayover~Stayovers$Yr, type='p', ylim=c(17, 35), ylab='',
las=1, xlim=c(1992, 2018), xlab='', cex=2, axes=F, yaxs='i', pch=21, col='white', bg='black')
axis(1, cex.axis=axes.labels.cex)
axis(4, labels=FALSE)
par(xpd=TRUE)
mtext(side=1, 'year', line=3, cex=axes.labels.cex)
par(xpd=FALSE)
axis(2, cex.axis=axes.labels.cex, las=1)
mtext(side=2, 'refueling time (d)', line=3,cex=axes.labels.cex,)
abline(lm(Stayover~Yr, data=Stayovers), lwd=2)
# 95%CI
XX<-seq(1992,2018, by=0.1)
Pred<-predict(M1, newdata=data.frame(Yr=XX), se.fit=TRUE)
UCI<-Pred$fit+1.96*Pred$se.fit
LCI<-Pred$fit-1.96*Pred$se.fit
lines(UCI~XX, lty=2)
lines(LCI~XX, lty=2)
box()
Stopover_before_now<-predict(M1, newdata=data.frame(Yr=c(1995, 2016)))
points(Stopover_before_now[1]~I(1995), pch=21, bg=Early_color, col='white', cex=4, lwd=3)
points(Stopover_before_now[2]~I(2016), pch=22, col=grey(0.45), bg=Late_color, cex=4, lwd=3)
#---------------------------------------
# panel 2
Prediction.matrix_cur_males<-subset(Prediction.matrix, fueling_slope_males>=3.5 &
Stopover<=35 & fueling_slope_males<=8.5 & Stopover>17)
Male_Image<-as.image(x=Prediction.matrix_cur_males[,c(7,6)],
Z=Prediction.matrix_cur_males$Prediction_males, nx=200, ny=200)
plot(0,0, type='n', xaxt='n', yaxt='n', xlab='', ylab='', bty='n') # make an empty plot
par(xpd=TRUE)
mtext("refueling rate (g/d)", side=1, line=3, cex=axes.labels.cex)
par(xpd=FALSE)
par(new=T)
par(mar=c(4, 5,2,2 ))
image(Male_Image,
col=Colors(256), xlab='', ylab='', las=1, ylim=c(17,35), xlim=c(3.5,8.5),
cex.axis=axes.labels.cex)
contour(Male_Image, add=T, levels=plogis(Design_males_s%*%as.matrix(summary(Model)$beta[,1, drop=FALSE]))[1],
lwd=3, col=grey(0.5),
labels=round(plogis(Design_males_s%*%as.matrix(summary(Model)$beta[,1, drop=FALSE]))[1],2),
labcex = 1)
contour(Male_Image,add=T, levels=c(0.8, 0.85, 0.9, 0.95), col=grey(0.5), labcex=1)
points(Stayover~Fueling_slope_males, data=Stayovers, pch=21, col='white', bg='black', cex=2)
box()
ScaleX<-seq(Slope_males_before_now[1], 6.46, length.out=11)
points(rep(Stopover_before_now[2], 11)~ScaleX, pch='|')
points(Stopover_before_now[1]~Slope_males_before_now[1], pch=21, bg=Early_color, col='white', cex=4, lwd=3)
points(Stopover_before_now[2]~Slope_males_before_now[1], pch=22, bg='white',col=Late_color, cex=4, lwd=3)
points(Stopover_before_now[2]~c(6.46), pch=22, bg='black', col=Late_color, cex=4, lwd=3)
points(Stopover_before_now[2]~Slope_males_before_now[2], pch=22, col=grey(0.45), bg=Late_color, cex=4, lwd=3)
##-------------------------------------
# Now panel 3
par(mar=c(4, 5,2,2 ))
plot(Stayovers$Fueling_slope_males~Stayovers$Yr, type='p', ylab='',
las=1, xlim=c(1992, 2018), xlab='', axes=F, yaxs='i',
cex=2, pch=21, bg='black', col='white', las=1, ylim=c(3.5,8.5))
axis(1, cex.axis=axes.labels.cex)
axis(3, labels=FALSE)
par(xpd=TRUE)
mtext(side=1, 'year', line=3, cex=axes.labels.cex)
par(xpd=FALSE)
axis(2, las=1, cex.axis=axes.labels.cex)
mtext(side=2, 'refueling rate (g/d)', line=3,cex=axes.labels.cex)
abline(lm(Fueling_slope_males~Yr, data=Stayovers), lwd=2)
# 95%CI
XX<-seq(1992,2018, by=0.1)
Pred<-predict(M3, newdata=data.frame(Yr=XX), se.fit=TRUE)
UCI<-Pred$fit+1.96*Pred$se.fit
LCI<-Pred$fit-1.96*Pred$se.fit
lines(UCI~XX, lty=2)
lines(LCI~XX, lty=2)
box()
points(Slope_males_before_now[1]~I(1995), pch=21, bg=Early_color, col='white', cex=4, lwd=3)
ScaleY<-seq(Slope_males_before_now[1], 6.46, length.out=11)
points(ScaleY~rep(2016, 11), pch=151)
points(Stopover_before_now[1]~Slope_males_before_now[1], pch=21, bg=Early_color, col='white', cex=4, lwd=3)
points(Slope_males_before_now[1]~I(2016), pch=22, bg='white',col=Late_color, cex=4, lwd=3)
points(I(6.46)~I(2016), pch=22, bg='black', col=Late_color, cex=4, lwd=3)
points(Slope_males_before_now[2]~I(2016), pch=22, col=grey(0.45), bg=Late_color, cex=4, lwd=3)
points(Stayovers$Fueling_slope_males~Stayovers$Yr, cex=2, pch=21, bg='black', col='white')
##-------------------------------------
# Now panel 4
par(mar=c(4, 5,2,2 ))
par(xaxs='i')
plot(Stayovers$Stayover~Stayovers$Fueling_slope_males, type='p', ylim=c(17,35), xlim=c(3.5,8.5) ,ylab='',
las=1, xlab='', cex=2, axes=F, yaxs='i', pch=21, col='white', bg='black')
axis(1, cex.axis=axes.labels.cex)
axis(3, labels=FALSE)
par(xpd=TRUE)
mtext(side=1, 'refueling rate (g/d)', line=3, cex=axes.labels.cex)
par(xpd=FALSE)
axis(2,cex.axis=axes.labels.cex, las=1)
mtext(side=2, 'refueling time (d)', line=3,cex=axes.labels.cex)
abline(lm(Stayover~Fueling_slope_males, data=Stayovers), lwd=2)
# 95%CI
XX<-seq(3.5,8.5, by=0.1)
Pred<-predict(M6, newdata=data.frame(Fueling_slope_males=XX), se.fit=TRUE)
UCI<-Pred$fit+1.96*Pred$se.fit
LCI<-Pred$fit-1.96*Pred$se.fit
lines(UCI~XX, lty=2)
lines(LCI~XX, lty=2)
box()
points(Stopover_before_now[1]~Slope_males_before_now[1], pch=21, bg=Early_color, col='white', cex=4, lwd=3)
points(Stopover_before_now[2]~Slope_males_before_now[2], pch=22, col=grey(0.45), bg=Late_color, cex=4, lwd=3)
##-------------------------------------
# Now panel 5
par(mar=c(4, 5,2,2 ))
par(xaxs='r')
plot(lugworms_ad~Yr, data=Lugworms, type='p', ylab='',
las=1, xlim=c(1992, 2018), xlab='', axes=F, yaxs='i',
cex=2, pch=21, bg='black', col='white',
ylim=c(min(Lugworms$lugworms_ad, na.rm=TRUE)-0.5, max(Lugworms$lugworms_ad, na.rm=TRUE)+0.5))
axis(1, cex.axis=axes.labels.cex)
axis(3, labels=FALSE)
par(xpd=TRUE)
mtext(side=1, 'year', line=3, cex=axes.labels.cex)
par(xpd=FALSE)
axis(2, las=1, cex.axis=axes.labels.cex)
mtext(side=2, 'lugworm density (ind./sq.m)', line=3,cex=axes.labels.cex)
abline(lm(lugworms_ad~Yr, data=Lugworms), lwd=2)
# 95%CI
XX<-seq(1992,2018, by=0.1)
Pred<-predict(M7, newdata=data.frame(Yr=XX), se.fit=TRUE)
UCI<-Pred$fit+1.96*Pred$se.fit
LCI<-Pred$fit-1.96*Pred$se.fit
lines(UCI~XX, lty=2)
lines(LCI~XX, lty=2)
box()
##-------------------------------------
# Now panel 6
par(mar=c(4, 5,2,2 ))
plot(fueling_model$BUGSoutput$mean$Alpha_t[,1]~Lugworms$lugworms_ad,
ylim=c(3.5,8.5),ylab='',
xlab='' ,las=1, type='n', axes=FALSE)
axis(2, cex.axis=axes.labels.cex, las=1)
mtext(side=2, 'refueling rate (g/d)', line=3, cex=axes.labels.cex)
axis(1, las=1, cex.axis=axes.labels.cex)
par(xpd=TRUE)
mtext(side=1, 'lugworm density (ind./sq.m)', line=3,cex=axes.labels.cex)
par(xpd=FALSE)
# 95%CI
XX<-seq(min(Lugworms$lugworms_ad), max(Lugworms$lugworms_ad), 0.1)
Y_m<- t(sapply(XX, FUN=function(x)
quantile(x*11 /(x+fueling_model$BUGSoutput$sims.list$K_sat_m),
probs=c(0.025, 0.5,0.975))))
lines(Y_m[,1]~XX, lty=2)
lines(Y_m[,2]~XX, lwd=2)
lines(Y_m[,3]~XX, lty=2)
points(fueling_model$BUGSoutput$mean$Alpha_t[,1]~Lugworms$lugworms_ad,
pch=21, cex=2, bg='black', col='white')
box()
#--------------------------------------------------------------------------
if (Save_pdf) dev.off()10. Plot Figure S1, Phenological data over time
10.1 Load data
# read Snow
Snow<-read.csv('./results/Snow_remote.csv', as.is=TRUE)
# read arrivals to Wadden Sea (predicted by arrival model)
Arr_Wadden<-read.csv('./results/Arrivals_migration_model.csv')
if (Use_local_data) {
# read breeding data
Nests<-read.csv('./data/Nests.csv', as.is=TRUE)
# read crane flies
Tipula<-read.csv('./data/Tipula.csv', as.is=TRUE)
# read arrivals to Khatanga
Arrivals_Khat<-read.csv('./data/Arrivals_Khat.csv', as.is=TRUE)
} else {
# read breeding data
Nests<-read.csv('https://git.io/vXu3y', as.is=TRUE)
# read crane flies
Tipula<-read.csv('https://git.io/vXu3x', as.is=TRUE)
# read arrivals to Khatanga
Arrivals_Khat<-read.csv('https://git.io/vXusf', as.is=TRUE)
}
names(Snow)<-c('Yr', 'yday', 'yday_South')
Nests$Date<-as.Date(Nests$Date, tz='GMT', format='%Y-%m-%d' )
Tipula$Date<-as.Date(Tipula$Date, tz='GMT', format='%Y-%m-%d' )
Arrivals_Khat$Date<-as.Date(Arrivals_Khat$Date, tz='GMT', format='%Y-%m-%d' )
names(Arr_Wadden)[2]<-'yday'10.2 Plot the figure
if (Save_pdf) pdf('Taimyr_and_wadden.pdf', width=5, height=10, useDingbats=FALSE)
par(mar=c(4,4,2,6))
Snow_global_color<-'#bae4bc'
Snow_local_color<-'#8BAB8D'
Tipula_color<-'#8783CE'
Khat_arr_color<-'#7B3575'
Nest_color<-'#CC79A7'
Arr_WZ_col<-'#7B3575'
Color_BA<-'#7B3575'
# snowmelt
plot(yday~Yr, data=Snow, type='n', ylim=c(115, 190), las=1, xlim=c(1993, 2016),axes=FALSE, ylab='')
axis(1)
box()
Snow<-subset(Snow, Yr %in% c(1992:2017))
polygon(y=c(Snow$yday, rep(200, length(Snow$yday))),
x=c(Snow$Yr, rev(Snow$Yr)), cex=4,
bg=grey(0.5), lwd=2, col=Snow_global_color, border='#7bccc4')
points(yday_South~Yr, data=Snow, bg=Snow_local_color, cex=2, col=NULL, pch=21)
# time lines
abline(h=121, col=grey(0.5))
abline(h=152, col=grey(0.5))
abline(h=182, col=grey(0.5))
abline(h=121+15, col=grey(0.5), lty=2)
abline(h=152+15, col=grey(0.5), lty=2)
# and lables for them
par(xpd=TRUE)
text(y=121, x=2019, '1 May')
text(y=152, x=2019, '1 June')
text(y=182, x=2019, '1 July')
par(xpd=FALSE)
# snow
abline(lm(yday_South~Yr, data=subset(Snow, Yr %in% c(1994:2016))), col=Snow_local_color, lwd=2)
# tipula #8783ce
points(yday~Yr, data=subset(Tipula, Region=='South'), pch=21, bg=Tipula_color, lwd=1, cex=2, col=NULL )
lm_tipula<-lm(yday~Yr, data=subset(Tipula, Region=='South'))
lines(predict(lm_tipula, newdata=data.frame(Yr=c(1992,2017)))~
c(1992,2017),lwd=2, col=Tipula_color)
# nests
points(yday~Yr, data=subset(Nests, Region=='South'), pch=21, cex=2, lwd=1, bg=Nest_color, col=NULL)
lm_nests<-lm(yday~Yr, data=subset(Nests, Region=='South'))
lines(predict(lm_nests, newdata=data.frame(Yr=c(1992,2016)))~
c(1992,2017), lwd=2, col=Nest_color)
# Khat arrivals
abline(lm(yday~Yr, data=Arrivals_Khat),col=Khat_arr_color, lwd=2)
points(yday~Yr, data=subset(Arrivals_Khat, Yr>=1993), pch=21, cex=2, bg=Khat_arr_color, lwd=1, col=NULL)
abline(lm(yday~Yr, data=Arr_Wadden), col=Arr_WZ_col, lwd=2)
points(yday~Yr, data=Arr_Wadden, bg=Arr_WZ_col, pch=21, cex=2, col=grey(0.9))
box()
if (Save_pdf) dev.off()11 Plot capture data with slopes (Figures S3 and S4)
11.1 Load data
# a. get arrivals to the Wadden Sea
Arrivals_WS<-read.csv('./results/Arrivals_migration_model.csv')
# b. get arrival mass
Arrival_weights<-read.csv('./results/Arrival_weights.csv')
# c. get annual slopes
Fueling_slopes<-read.csv('./results/Fueling_slopes.csv')
if (Use_local_data) {
# d. get departures from Wadden Sea through arrival to Taimyr
Arrival_Taimyr<-read.csv('./data/Arrivals_Khat.csv')
} else {
# d. get departures from Wadden Sea through arrival to Taimyr
Arrival_Taimyr<-read.csv('https://git.io/vXEbj')
}11.2 Prepare data for plotting
#names(Arrivals_WS)[2]<-'yday'
Arrivals<-merge(Arrivals_WS, Arrival_weights)
Arrivals$Left_95<-Arrivals$arrival-1.96*Arrivals$sigma
Arrivals$Right_95<-Arrivals$arrival+1.96*Arrivals$sigma
Arrivals$Males_low<-Arrivals$Male_arrival_weight-1.96*Arrivals$Male_arrival_weight_sigma
Arrivals$Males_up<-Arrivals$Male_arrival_weight+1.96*Arrivals$Male_arrival_weight_sigma
Arrivals$Females_low<-Arrivals$Female_arrival_weight-1.96*Arrivals$Female_arrival_weight_sigma
Arrivals$Females_up<-Arrivals$Female_arrival_weight+1.96*Arrivals$Female_arrival_weight_sigma
Arrivals<-merge(Arrivals, Fueling_slopes[,1:3])
Arrivals$Male_fueling_intercept<- Arrivals$Male_arrival_weight-
Arrivals$Fueling_slope_males*Arrivals$arrival
Arrivals$Female_fueling_intercept<-Arrivals$Female_arrival_weight-
Arrivals$Fueling_slope_females*Arrivals$arrival
Departure_Wadden<-Arrival_Taimyr[,c(1,3)]
Departure_Wadden$Departure<-Departure_Wadden$yday-5.5
Departure_Wadden<-Departure_Wadden[,-2]
Arrivals<-merge(Arrivals,Departure_Wadden, all.x=TRUE )
Arrivals$Male_Departure_Weight<-Arrivals$Male_fueling_intercept+
Arrivals$Fueling_slope_males*Arrivals$Departure
Arrivals$Female_Departure_Weight<-Arrivals$Female_fueling_intercept+
Arrivals$Fueling_slope_females*Arrivals$Departure
Capture_data_males<-subset(Capture_data, Sex==1 & Yr>1995)
Capture_data_females<-subset(Capture_data, Sex==2 & Yr>1995)
Lwd<-111.3 Plot figure for males
library(ggplot2)
p<-ggplot(data = Capture_data_males, aes(x = yday, y = Mass)) +
ylim(150, 450) +
geom_point(col='#fdb194', data=subset(Capture_data_males, Site!='Castricum')) +
geom_point(col='#F0E442', data=subset(Capture_data_males, Site=='Castricum'))+
facet_wrap( ~ Yr) +
# arrival time sd
geom_segment(data=Arrivals, aes(x = Left_95, y = Male_arrival_weight,
xend=Right_95, yend=Male_arrival_weight), col='#009E73', lwd=Lwd) +
geom_segment(data=Arrivals, aes(x = arrival, y = Males_low,
xend=arrival, yend=Males_up), col='#E69F00', lwd=Lwd) +
geom_abline(data=Arrivals, aes(intercept=Male_fueling_intercept,
slope=Fueling_slope_males), col='#D55E00', lwd=Lwd) +
geom_vline(data=subset(Arrivals), aes(xintercept = Departure), col="#56B4E9", lwd=Lwd) +
geom_hline(col='#0072B2', data=Arrivals, aes(yintercept = Male_Departure_Weight), lwd=Lwd)+
labs(title='Males')
suppressWarnings(print(p))
if (Save_pdf) ggsave(p, file='Males_captures.pdf',
device='pdf', width=10, height=10, useDingbats=FALSE)11.4 Plot figure for females
p1<-ggplot(data = Capture_data_females, aes(x = yday, y = Mass)) +
geom_point(col='#fdb194', data=subset(Capture_data_females, Site!='Castricum')) +
geom_point(col='#F0E442', data=subset(Capture_data_females, Site=='Castricum'))+
facet_wrap( ~ Yr) +
# arrival time sd
geom_segment(data=Arrivals, aes(x = Left_95, y = Female_arrival_weight,
xend=Right_95, yend=Female_arrival_weight), col='#009E73', lwd=Lwd) +
geom_segment(data=Arrivals, aes(x = arrival, y = Females_low,
xend=arrival, yend=Females_up), col='#E69F00', lwd=Lwd) +
geom_abline(data=Arrivals, aes(intercept=Female_fueling_intercept,
slope=Fueling_slope_females), col='#D55E00', lwd=Lwd) +
geom_vline(data=subset(Arrivals), aes(xintercept = Departure), col="#56B4E9", lwd=Lwd) +
geom_hline(col='#0072B2', data=Arrivals, aes(yintercept = Female_Departure_Weight), lwd=Lwd)+
labs(title='Females')
suppressWarnings(print(p1))
if (Save_pdf) ggsave(p1, file='Females_captures.pdf',
device='pdf', width=10, height=10, useDingbats=FALSE)11.5 Plot legend
Leg_frame<-data.frame(x=c(-1, 1), y=c(-1, 1))
X_labels<--0.85
X_text<--0.6
p3<-ggplot(Leg_frame, aes())+
ylim(-1,1)+
xlim(-1,1)+
theme(
text = element_blank(),
line = element_blank(),
title = element_blank())+
geom_point(col='#fdb194', data=data.frame(x=X_labels, y=0.75), aes(x=x, y=y), cex=3)+
annotate("text", x=X_text, y = 0.75, label = "Wadden captures",hjust = 0)+
geom_point(col='#F0E442', data=data.frame(x=X_labels, y=0.5), aes(x=x, y=y), cex=3)+
annotate("text", x=X_text, y = 0.5, label = "Castricum captures",hjust = 0)+
annotate("segment", x=X_labels-0.1, xend=X_labels+0.1, y=0.25, yend=0.25,
col='#009E73', lwd=1.5)+
annotate("text", x=X_text, y = 0.25, label = "arrival period (95%)",hjust = 0)+
annotate("segment", x=X_labels, xend=X_labels, y=0-0.075, yend=0+0.075,
col='#E69F00', lwd=1.5)+
annotate("text", x=X_text, y = 0, label = "arrival mass (95%)",hjust = 0) +
annotate("segment", x=X_labels-0.1, xend=X_labels+0.1, y=-0.25-0.075, yend=-0.25+0.075,
col='#D55E00', lwd=1.5)+
annotate("text", x=X_text, y = -0.25, label = "refueling slope", hjust = 0) +
annotate("segment", x=X_labels, xend=X_labels, y=-0.5-0.075, yend=-0.5+0.075,
col='#56B4E9', lwd=1.5)+
annotate("text", x=X_text, y = -0.5, label = "departure time",hjust = 0) +
annotate("segment", x=X_labels-0.1, xend=X_labels+0.1, y=-0.75, yend=-0.75,
col='#0072B2', lwd=1.5)+
annotate("text", x=X_text, y = -0.75, label = "departure mass",hjust = 0)
print(p3)
if (Save_pdf) ggsave(p3, file='Legend.pdf', device='pdf', width=2,
height=2, useDingbats=FALSE)References
- Bijleveld, A. I. et al. Designing a benthic monitoring programme with multiple conflicting objectives. Methods Ecol. Evol. 3, 526–536 (2012).
- Brown, R. D. & Robinson, D. A. Northern Hemisphere spring snow cover variability and change over 1922-2010 including an assessment of uncertainty. The Cryosphere 5, 219-229 (2011).
- Compton, T. J. et al. Distinctly variable mudscapes: Distribution gradients of intertidal macrofauna across the Dutch Wadden Sea. J. Sea Res. 82, 103–116 (2013).
- van Gils, J. A. et al. Body shrinkage due to Arctic warming reduces red knot fitness in tropical wintering range. Science 352, 819-821 (2016).
- Kéry, M. & Schaub, M. Bayesian population analysis using WinBUGS a hierarchical perspective. (Academic Press, 2012).
- Laake, J. L. RMark: An R interface for analysis of capture-recapture data with MARK. 25 (Alaska Fish. Sci. Cent., NOAA, Natl. Mar. Fish. Serv., 2013).
- Lappo, E., Tomkovich, P. & Syroechkovskiy, E. Atlas of breeding waders in the Russian Arctic. (2012).
- Lebreton, J.-D., Burnham, K. P., Clobert, J. & Anderson, D. R. Modeling survival and testing biological hypotheses using marked animals: a unified approach with case studies. Ecol. Monogr. 62, 67-118 (1992).
- Lindén, A. & Mäntyniemi, S. Using the negative binomial distribution to model overdispersion in ecological count data. Ecology 92, 1414-1421 (2011).
- Plummer, M. JAGS: A program for analysis of Bayesian graphical models using Gibbs sampling. (2003).
- Su, Y.-S. & Masanao Yajima. R2jags: Using R to Run ‘JAGS’. (2015).
- White, G. C. & Burnham, K. P. Program MARK: survival estimation from populations of marked animals. Bird Study 46, S120-S139 (1999).
P.S. The html file from this markdown can be recereated with folloing code
download.file('https://git.io/vXaM6', 'tmp.rmd', cacheOK = FALSE)
rmarkdown::render('tmp.rmd', output_format = 'html_document',
output_options=list(toc=TRUE, toc_float=list(collapsed=FALSE)),
encoding='utf-8')