### Susu Rytteri, 2021-02-04 ### This file contains the script needed to conduct the analyses of the paper: ### Rytteri S, Kuussaari M & Saastamoinen M 2021: Microclimatic variability buffers butterfly populations ### against increased mortality caused by phenological asynchrony between larvae and their host plants. ### Oikos. https://doi.org/10.1111/oik.07653 ### Spatiotemporal regression models for larval temperature and activity. ### Spatial regression models for larval survival and weight gain, local population growth, ### host plant growth, and larval food availability. # Set the working directory: setwd("") # Install and load R-INLA install.packages("INLA", repos=c(getOption("repos"), INLA="https://inla.r-inla-download.org/R/stable"), dep=TRUE) library(INLA) ### MODELS FOR LARVAL ACTIVITY AND TEMPERATURE ### # Read in the temperature and activity data set: tdata=data.frame(read.table("Rytteri_et_al_2021_temperature_data.txt",header=T)) # Standardize the explanatory variables: MyStd <- function(x) { (x - mean(x)) / sd(x)} tdata$Groupsize.std <- MyStd(tdata$Groupsize) tdata$Air_Temp.std <- MyStd(tdata$Air_Temp) tdata$Larval_Temp.std <- MyStd(tdata$Larval_Temp) tdata$Sun_nest.std <- MyStd(tdata$Sun_nest) # Make a mesh # Define sampling locations in kilometres: tdata$Xkm <- tdata$X / 1000 tdata$Ykm <- tdata$Y / 1000 tLoc <- cbind(tdata$Xkm, tdata$Ykm) tmesh <- inla.mesh.2d(tLoc, max.edge=c(1, 5), cutoff = 0.01) # Maximum edge length of triangles 1km (5km in the buffer zone) and minimum edge length 10 m. # Define the SPDE spde <- inla.spde2.pcmatern(tmesh, alpha = 2, prior.range = c(1,0.5), prior.sigma = c(1,0.25)) # Prior spatial range of random field (rw) is P(rw < 1km) = 0.5 # and marginal sd of random field (sdw) is P(sdw > 1km) = 0.25. # Specify the grouping structure for the temporal autocorrelation: Group <- tdata$Group # We have 17 groups NGroup <- length(unique(Group)) NGroup # Define the spatial random field w.index <- inla.spde.make.index(name = 'w', n.spde = tmesh$n, n.group = NGroup) # Define the projector matrix A2 <- inla.spde.make.A(tmesh, group = Group, loc = tLoc) # Make a stack # Create a data frame with an intercept and covariates for the model of larval activity: N <- nrow(tdata) Xa <- data.frame(Intercept = rep(1, N), Group.std = tdata$Groupsize.std, AmbT.std = tdata$Air_Temp.std, LarvT.std = tdata$Larval_Temp.std, Sun.std = tdata$Sun_nest.std ) Stacka <- inla.stack( tag = "Fit", data = list(y = tdata$Larval_act), A = list(A2, 1), effects = list( w = w.index, #Spatiotemporal correlations Xa = as.data.frame(Xa))) #Covariates # Define the formula of autoregressive correlation in terms of the response variable, covariates and the # spatially and temporally correlated terms fa <- y ~ -1 + Intercept + Group.std + AmbT.std + Sun.std + LarvT.std + f(w, model = spde, group = w.group, control.group = list(model='ar1')) # Run the spatiotemporal model in INLA Ia <- inla(fa, family = "Gaussian", data = inla.stack.data(Stacka), control.compute = list(dic = TRUE, waic = TRUE), control.predictor = list(A = inla.stack.A(Stacka))) summary(Ia) # Summary of the model for larval activity # Create a data frame with an intercept and covariates for the model of larval temperature: N <- nrow(tdata) Xt <- data.frame(Intercept = rep(1, N), Group.std = tdata$Groupsize.std, AmbT.std = tdata$Air_Temp.std, Sun.std = tdata$Sun_nest.std ) Stackt <- inla.stack( tag = "Fit", data = list(y = tdata$Larval_Temp), A = list(A2, 1), effects = list( w = w.index, Xt = as.data.frame(Xt))) ft <- y ~ -1 + Intercept + Group.std + AmbT.std + Sun.std + f(w, model = spde, group = w.group, control.group = list(model='ar1')) It <- inla(ft, family = "Gaussian", data = inla.stack.data(Stackt), control.compute = list(dic = TRUE, waic = TRUE), control.predictor = list(A = inla.stack.A(Stackt))) summary(It) # Summary of the model for larval temperature ### MODELS FOR LARVAL SURVIVAL AND WEIGHT GAIN ### # Read in the survival and weight gain data set: mdata=data.frame(read.table("Rytteri_et_al_2021_monitoring_data.txt",header=T)) mdata$Group5.std <- MyStd(mdata$Group5) mdata$DD.std <- MyStd(mdata$Degree_DayMean) mdata$Prec.std <- MyStd(mdata$Daily_Precip) mdata$Shade.std <- MyStd(mdata$Nest_Shade) mdata$Pred.std <- MyStd(mdata$Predation) mdata$Food.std <- MyStd(mdata$Food_6) mdata$Dir.std <- MyStd(mdata$N_Dir) mdata$Veg.std <- MyStd(mdata$Nest_Veg) mdata$Xkm <- mdata$Lon / 1000 mdata$Ykm <- mdata$Lat / 1000 mLoc <- cbind(mdata$Xkm, mdata$Ykm) mmesh <- inla.mesh.2d(mLoc, max.edge=c(1, 5), cutoff = 0.01) A2 <- inla.spde.make.A(mmesh, loc = mLoc) spde <- inla.spde2.pcmatern(mmesh, alpha = 2, prior.range = c(1,0.5), prior.sigma = c(1,0.25)) w.index <- inla.spde.make.index( name = 'w', n.spde = spde$n.spde, n.group = 1, n.repl = 1) N <- nrow(mdata) Xm <- data.frame(Intercept = rep(1, N), Group5.std = mdata$Group5.std, DD.std = mdata$DD.std, Prec.std = mdata$Prec.std, Shade.std = mdata$Shade.std, Pred.std = mdata$Pred.std, Food.std = mdata$Food.std, Dir.std = mdata$Dir.std, Veg.std = mdata$Veg.std ) ### MODEL FOR SURVIVAL OVER THE SPRING Stacks <- inla.stack( tag = "Fit", data = list(y = mdata$Group_end), A = list(A2, 1), effects = list( w = w.index, Xm = as.data.frame(Xm))) fs <- y ~ -1 + Intercept + Group5.std + DD.std + Prec.std + Shade.std + Pred.std + Food.std + Dir.std + Veg.std + f(w, model = spde) Is <- inla(fs, family = "binomial", Ntrials = mdata$Group5, data = inla.stack.data(Stacks), control.compute = list(dic = TRUE, waic = TRUE), control.predictor = list(A = inla.stack.A(Stacks))) summary(Is) # Summary of the model for larval survival ### MODEL FOR WEIGHT GAIN OVER THE SPRING Stackw <- inla.stack( tag = "Fit", data = list(y = mdata$Growth_total), A = list(A2, 1), effects = list( w = w.index, Xm = as.data.frame(Xm))) fw <- y ~ -1 + Intercept + Group5.std + DD.std + Prec.std +Shade.std + Pred.std + Food.std + Dir.std + Veg.std + f(w, model = spde) Iw <- inla(fw, family = "gaussian", data=inla.stack.data(Stackw), control.compute = list(dic = TRUE), control.predictor = list(A = inla.stack.A(Stackw))) summary(Iw) # Summary of the model for larval weight gain ### MODEL FOR POPULATION GROWTH 2016 ########################################################################### data16 <- read.table(file = "Rytteri_et_al_2021_survey15-16_data.txt", header = TRUE, dec = ".") data16$LgNest15.std <- MyStd(data16$LgNest15) data16$DegreeD_5.std <- MyStd(data16$DegreeD_5) data16$Precipitation_5.std <- MyStd(data16$Precipitation_5) data16$Host_mean.std <- MyStd(data16$Host_mean) data16$Shade_mean.std <- MyStd(data16$Shade_mean) data16$Shade_Stddev.std <- MyStd(data16$Shade_Stddev) data16$CompassP1_mean.std <- MyStd(data16$CompassP1_mean) data16$Vegetation_mean.std <- MyStd(data16$Vegetation_mean) data16$Xkm <- data16$X / 1000 data16$Ykm <- data16$Y / 1000 Loc16 <- cbind(data16$Xkm, data16$Ykm) mesh16 <- inla.mesh.2d(Loc16, max.edge=c(5, 10), cutoff = 0.2) A2 <- inla.spde.make.A(mesh16, loc = Loc16) spde2 <- inla.spde2.pcmatern(mesh16, alpha = 2, prior.range = c(5,0.5), prior.sigma = c(1,0.25)) w2.index <- inla.spde.make.index('w', n.spde = spde2$n.spde) Xm16 <- model.matrix(~ LgNest15.std + DegreeD_5.std + Precipitation_5.std + Host_mean.std + Shade_mean.std + Shade_Stddev.std + CompassP1_mean.std + Vegetation_mean.std, data = data16) N <- nrow(data16) X16 <- data.frame(Nests = Xm16[,2], DD = Xm16[,3], Prec = Xm16[,4], Host = Xm16[,5], Shade = Xm16[,6], Shadedev = Xm16[,7], Compass = Xm16[,8], Veg = Xm16[,9] ) Stack16 <- inla.stack( tag = "Fit", data = list(y = data16$LgPopChange), A = list(1, 1, A2), effects = list( Intercept = rep(1, N), X = X16, w = w2.index)) f16 <- y ~ -1 + Intercept + Nests + DD + Prec+ Host + Shade + Shadedev + Compass + Veg + f(w, model = spde2) I16 <- inla(f16, family = "gaussian", data=inla.stack.data(Stack16), control.compute = list(dic = TRUE), control.predictor = list(A = inla.stack.A(Stack16))) summary(I16) # Summary of the model for population growth rate 2015-2016 ### MODEL FOR POPULATION GROWTH 2017 ########################################################################## data17 <- read.table(file = "Rytteri_et_al_2021_survey16-17_data.txt", header = TRUE, dec = ".") data17$LgNest16.std <- MyStd(data17$LgNest16) data17$DegreeD_5.std <- MyStd(data17$DegreeD_5) data17$Precipitation_5.std <- MyStd(data17$Precipitation_5) data17$Host_mean.std <- MyStd(data17$Host_mean) data17$Shade_mean.std <- MyStd(data17$Shade_mean) data17$Shade_Stddev.std <- MyStd(data17$Shade_Stddev) data17$CompassP1_mean.std <- MyStd(data17$CompassP1_mean) data17$Vegetation_mean.std <- MyStd(data17$Vegetation_mean) data17$Xkm <- data17$X / 1000 data17$Ykm <- data17$Y / 1000 Loc17 <- cbind(data17$Xkm, data17$Ykm) mesh17 <- inla.mesh.2d(Loc17, max.edge=c(5, 10), cutoff = 0.2) A2 <- inla.spde.make.A(mesh17, loc = Loc17) spde2 <- inla.spde2.pcmatern(mesh17, alpha = 2, prior.range = c(5,0.5), prior.sigma = c(1,0.25)) w2.index <- inla.spde.make.index('w', n.spde = spde2$n.spde) Xm17 <- model.matrix(~ LgNest16.std + DegreeD_5.std + Precipitation_5.std + Host_mean.std + Shade_mean.std + Shade_Stddev.std + CompassP1_mean.std + Vegetation_mean.std, data = data17) N <- nrow(data17) X17 <- data.frame(Nests = Xm17[,2], DD = Xm17[,3], Prec = Xm17[,4], Host = Xm17[,5], Shade = Xm17[,6], Shadedev = Xm17[,7], Compass = Xm17[,8], Veg = Xm17[,9] ) Stack17 <- inla.stack( tag = "Fit", data = list(y = data17$LgPopChange), A = list(1, 1, A2), effects = list( Intercept = rep(1, N), X = X17, w = w2.index)) f17 <- y ~ -1 + Intercept + Nests + DD + Prec+ Host + Shade + Shadedev + Compass + Veg + f(w, model = spde2) I17 <- inla(f17, family = "gaussian", data=inla.stack.data(Stack17), control.compute = list(dic = TRUE), control.predictor = list(A = inla.stack.A(Stack17))) summary(I17) # Summary of the model for population growth rate 2016-2017 ### MODEL FOR HOST PLANT GROWTH OVER SPRING ###################################################################### datah <- read.table(file = "Rytteri_et_al_2021_host_plant_data.txt", header = TRUE, dec = ".") datah$DD.std <- MyStd(datah$Degree_DayMean) datah$Prec.std <- MyStd(datah$Daily_Precip) datah$Dir.std <- MyStd(datah$N_Dir) datah$Veg.std <- MyStd(datah$Nest_Veg) datah$Shade.std <- MyStd(datah$Nest_Shade) datah$Xkm <- datah$Lon / 1000 datah$Ykm <- datah$Lat / 1000 Loch <- cbind(datah$Xkm, datah$Ykm) meshh <- inla.mesh.2d(Loch, max.edge=c(1, 5), cutoff = 0.01) A2 <- inla.spde.make.A(meshh, loc = Loch) spde2 <- inla.spde2.pcmatern(meshh, alpha = 2, prior.range = c(1,0.5), prior.sigma = c(1,0.25)) w2.index <- inla.spde.make.index('w', n.spde = spde2$n.spde) Xmh <- model.matrix(~ DD.std + Prec.std + Dir.std + Veg.std + Shade.std, data = datah) N <- nrow(datah) Xh <- data.frame(DD = Xmh[,2], Prec = Xmh[,3], Dir = Xmh[,4], Veg = Xmh[,5], Shade = Xmh[,6] ) fh <- y ~ -1 + Intercept + DD + Prec + Dir + Shade + Veg + f(w, model = spde2) # Host plant growth in the first half of the study season Stackhle <- inla.stack( tag = "Fit", data = list(y = datah$Ch_H_l_e), A = list(1, 1, A2), effects = list( Intercept = rep(1, N), X = Xh, w = w2.index)) Ihle <- inla(fh, family = "gaussian", data=inla.stack.data(Stackhle), control.compute = list(dic = TRUE), control.predictor = list(A = inla.stack.A(Stackhle))) summary(Ihle) # Summary of the model for host plant growth in the first half of the study season # Host plant growth in the second half of the study season Stackhll <- inla.stack( tag = "Fit", data = list(y = datah$Ch_H_l_l), A = list(1, 1, A2), effects = list( Intercept = rep(1, N), X = Xh, w = w2.index)) Ihll <- inla(fh, family = "gaussian", data=inla.stack.data(Stackhll), control.compute = list(dic = TRUE), control.predictor = list(A = inla.stack.A(Stackhll))) summary(Ihll) # Summary of the model for host plant growth in the second half of the study season # Host plant growth of the whole study season Stackhlt <- inla.stack( tag = "Fit", data = list(y = datah$Ch_H_l_tot), A = list(1, 1, A2), effects = list( Intercept = rep(1, N), X = Xh, w = w2.index)) Ihlt <- inla(fh, family = "gaussian", data=inla.stack.data(Stackhlt), control.compute = list(dic = TRUE), control.predictor = list(A = inla.stack.A(Stackhlt))) summary(Ihlt) # Summary of the model for host plant growth in the whole study season ### MODEL FOR LARVAL FOOD AVAILABILITY OVER SPRING ############################################################### dataf <- read.table(file = "Rytteri_et_al_2021_food_availability_data.txt", header = TRUE, dec = ".") dataf$DD.std <- MyStd(dataf$Degree_DayMean) dataf$Prec.std <- MyStd(dataf$Daily_Precip) dataf$Dir.std <- MyStd(dataf$N_Dir) dataf$Veg.std <- MyStd(dataf$Nest_Veg) dataf$Shade.std <- MyStd(dataf$Nest_Shade) dataf$Group.std <- MyStd(dataf$Group5) dataf$Xkm <- dataf$Lon / 1000 dataf$Ykm <- dataf$Lat / 1000 Locf <- cbind(dataf$Xkm, dataf$Ykm) meshf <- inla.mesh.2d(Locf, max.edge=c(1, 5), cutoff = 0.01) A2 <- inla.spde.make.A(meshf, loc = Locf) spde2 <- inla.spde2.pcmatern(meshf, alpha = 2, prior.range = c(1,0.5), prior.sigma = c(1,0.25)) w2.index <- inla.spde.make.index('w', n.spde = spde2$n.spde) Xmf <- model.matrix(~ DD.std + Prec.std + Dir.std + Veg.std + Shade.std + Group.std, data = dataf) N <- nrow(dataf) Xf <- data.frame(DD = Xmf[,2], Prec = Xmf[,3], Dir = Xmf[,4], Veg = Xmf[,5], Shade = Xmf[,6], Group = Xmf[,7] ) ff <- y ~ -1 + Intercept + DD + Prec + Dir + Shade + Veg + Group + f(w, model = spde2) # Change of food availability in the first half of the study season Stackfe <- inla.stack( tag = "Fit", data = list(y = dataf$Ch_F_e), A = list(1, 1, A2), effects = list( Intercept = rep(1, N), X = Xf, w = w2.index)) Ife <- inla(ff, family = "gaussian", data=inla.stack.data(Stackfe), control.compute = list(dic = TRUE), control.predictor = list(A = inla.stack.A(Stackfe))) summary(Ife) # Summary of the model for change of food availability in the first half of the study season # Change of food availability in the second half of the study season Stackfl <- inla.stack( tag = "Fit", data = list(y = dataf$Ch_F_l), A = list(1, 1, A2), effects = list( Intercept = rep(1, N), X = Xf, w = w2.index)) Ifl <- inla(ff, family = "gaussian", data=inla.stack.data(Stackfl), control.compute = list(dic = TRUE), control.predictor = list(A = inla.stack.A(Stackfl))) summary(Ifl) # Summary of the model for change of food availability in the second half of the study season # Change of food availability of the whole study season Stackft <- inla.stack( tag = "Fit", data = list(y = dataf$Ch_F_tot), A = list(1, 1, A2), effects = list( Intercept = rep(1, N), X = Xf, w = w2.index)) Ift <- inla(ff, family = "gaussian", data=inla.stack.data(Stackft), control.compute = list(dic = TRUE), control.predictor = list(A = inla.stack.A(Stackft))) summary(Ift) # Summary of the model for change of food availability in the whole study season