The R code provided below performs the tests of the slope of phenology (for Anthocharis cardamines, Cardamine pratensis or Alliaria petiolata) on temprature over space and time (PART 1) and the multiple regression of A. cardamines phenology on the two host plants and temperatur eover space and time (PART 2).

The tab delimited file of phenology and average temperature data can be downloaded from Dryad (doi:10.5061/dryad.733d9).This file should be saved on your computer and read in.

phenol<-read.table("...../butterfly.dat",sep="\t",header=T)
#Note that the dataset provided here only contains the average temperature data for the (latitudinally-invarant) periods that were found to return the highest R2 for each of the focal species. If you require temperature data for alternative time periods, including from daylength inititiated time wiondows please contact Albert Phillimore (albert.phillimore@ed.ac.uk).

library(MCMCglmm)
#Please refer to the vignette written by Jarrod Hadfield and available on CRAN (http://cran.r-project.org/web/packages/MCMCglmm/vignettes/CourseNotes.pdf) for a detailed introduction to using MCMCglmm. Many thanks to Jarrod for substantial advice on fitting the models described below.

#######################################################################################
#We will use this function later to calculate a psuedo R2. See Phillimore A.B., Hadfield J.D., Jones O.R. & Smithers R.J. (2010). Differences in spawning date between populations of common frog reveal local adaptation. Proceedings of the National Academy of Sciences, 107, 8292-8297.

get.R2<-function(VCVmatrix,rand.effects){
	numerator<-rep(0,dim(VCVmatrix)[1])
	denominator<-rep(0,dim(VCVmatrix)[1])
	for(reps in 1:rand.effects){
		numerator<-numerator+(VCVmatrix[,(reps*4)-2]^2)/VCVmatrix[,(reps*4)]
		denominator<-denominator+VCVmatrix[,(reps*4)-3]
		}

		R2vals<-mcmc(numerator/denominator)
	return(list("R2"=median(R2vals),"HPDinterval"=HPDinterval(R2vals,0.95)))
}


#####################################################################################
#1. analyses with species phenology and average temperature as a bivariate response.
#####################################################################################

###############Anthocharis_cardamines################################################

Ac_data<-phenol[which(is.na(phenol$Ac_phenol)==FALSE),]

#MCMCglmm with default priors. 
Ac_model<-MCMCglmm(cbind(Ac_phenol,Avgtemp_61to125)~trait-1, random=~us(trait):Grid_150km+us(trait):Year,rcov=~us(trait):units,family=c("gaussian","gaussian"),nitt=23000,burnin=3000,data=Ac_data)

####model diagnostics####

#What degree of autocorrelation is there between succesive samples of the posterior distribution?
autocorr(Ac_model$Sol)
autocorr(Ac_model$VCV)

#Also check the trace files for the estimated means and variances.
plot(Ac_model$Sol)
plot(Ac_model$VCV)

#A summary table for the model.
summary(Ac_model)
#Note in all cases we want the effective sample size to be large (approaching 2000).


####Model fit and slope estimates####

get.R2(Ac_model$VCV,3)
#Note that this pseudo R2 takes into account the proportion of the variance in phenology that temperature explains (i) across grid cells, (ii) across years and (iii) within grid cells and years.

#Median spatial slope and the 95% highest posterior density
median(Ac_model$VCV[,2]/Ac_model$VCV[,4])
HPDinterval(Ac_model$VCV[,2]/Ac_model$VCV[,4])

#Median temporal slope and the 95% highest posterior density
median(Ac_model$VCV[,6]/Ac_model$VCV[,8])
HPDinterval(Ac_model$VCV[,6]/Ac_model$VCV[,8])

#Median residual slope and the 95% highest posterior density
median(Ac_model$VCV[,10]/Ac_model$VCV[,12])
HPDinterval(Ac_model$VCV[,10]/Ac_model$VCV[,12])

#Test for slope difference (delta b = spatial slope - temporal slope)
median(Ac_model$VCV[,2]/Ac_model$VCV[,4]-Ac_model$VCV[,6]/Ac_model$VCV[,8])
HPDinterval(Ac_model$VCV[,2]/Ac_model$VCV[,4]-Ac_model$VCV[,6]/Ac_model$VCV[,8])
#an HPD interval that does not include zero corresponds to significant evidence for local adaptation.

###############Cardamine pratensis################################################

Cp_data<-phenol[which(is.na(phenol$Cp_phenol)==FALSE),]

Cp_model<-MCMCglmm(cbind(Cp_phenol,Avgtemp_61to95)~trait-1, random=~us(trait):Grid_150km+us(trait):Year,rcov=~us(trait):units,family=c("gaussian","gaussian"),nitt=23000,burnin=3000,data=Cp_data)


####model diagnostics####

#What degree of autocorrelation is there between succesive samples of the posterior distribution?
autocorr(Cp_model$Sol)
autocorr(Cp_model$VCV)

#Also check the trace files for the estimated means and variances.
plot(Cp_model$Sol)
plot(Cp_model$VCV)

#A summary table for the model.
summary(Cp_model)
#Note in all cases we want the effective sample size to be large (approaching 2000).


####Model fit and slope estimates####

get.R2(Cp_model$VCV,3)
#Note that this pseudo R2 takes into account the proportion of the variance in phenology that temperature explains (i) across grid cells, (ii) across years and (iii) within grid cells and years.

#Median spatial slope and the 95% highest posterior density
median(Cp_model$VCV[,2]/Cp_model$VCV[,4])
HPDinterval(Cp_model$VCV[,2]/Cp_model$VCV[,4])

#Median temporal slope and the 95% highest posterior density
median(Cp_model$VCV[,6]/Cp_model$VCV[,8])
HPDinterval(Cp_model$VCV[,6]/Cp_model$VCV[,8])

#Median residual slope and the 95% highest posterior density
median(Cp_model$VCV[,10]/Cp_model$VCV[,12])
HPDinterval(Cp_model$VCV[,10]/Cp_model$VCV[,12])

#Test for slope difference (delta b = spatial slope - temporal slope)
median(Cp_model$VCV[,2]/Cp_model$VCV[,4]-Cp_model$VCV[,6]/Cp_model$VCV[,8])
HPDinterval(Cp_model$VCV[,2]/Cp_model$VCV[,4]-Cp_model$VCV[,6]/Cp_model$VCV[,8])

###############Alliara petiolata################################################

Ap_data<-phenol[which(is.na(phenol$Ap_phenol)==FALSE),]

Ap_model<-MCMCglmm(cbind(Ap_phenol,Avgtemp_46to95)~trait-1, random=~us(trait):Grid_150km+us(trait):Year,rcov=~us(trait):units,family=c("gaussian","gaussian"),nitt=23000,burnin=3000,data=Ap_data)


####model diagnostics####

#What degree of autocorrelation is there between succesive samples of the posterior distribution?
autocorr(Ap_model$Sol)
autocorr(Ap_model$VCV)

#Also check the trace files for the estimated means and variances.
plot(Ap_model$Sol)
plot(Ap_model$VCV)

#A summary table for the model.
summary(Ap_model)
#Note in all cases we want the effective sample size to be large (approaching 2000).


####Model fit and slope estimates####

get.R2(Ap_model$VCV,3)
#Note that this pseudo R2 takes into account the proportion of the variance in phenology that temperature explains (i) across grid cells, (ii) across years and (iii) within grid cells and years.

#Median spatial slope and the 95% highest posterior density
median(Ap_model$VCV[,2]/Ap_model$VCV[,4])
HPDinterval(Ap_model$VCV[,2]/Ap_model$VCV[,4])

#Median temporal slope and the 95% highest posterior density
median(Ap_model$VCV[,6]/Ap_model$VCV[,8])
HPDinterval(Ap_model$VCV[,6]/Ap_model$VCV[,8])

#Median residual slope and the 95% highest posterior density
median(Ap_model$VCV[,10]/Ap_model$VCV[,12])
HPDinterval(Ap_model$VCV[,10]/Ap_model$VCV[,12])

#Test for slope difference (delta b = spatial slope - temporal slope)
median(Ap_model$VCV[,2]/Ap_model$VCV[,4]-Ap_model$VCV[,6]/Ap_model$VCV[,8])
HPDinterval(Ap_model$VCV[,2]/Ap_model$VCV[,4]-Ap_model$VCV[,6]/Ap_model$VCV[,8])




#####################################################################################
#2. Bivariate and multivariate responses with butterfly and host both as responses
#####################################################################################
#first we scale the data to have mean = 0 and variance = 1.
phenol$Ac_phenol_scale<-scale(phenol$Ac_phenol)
phenol$Cp_phenol_scale<-scale(phenol$Cp_phenol)
phenol$Ap_phenol_scale<-scale(phenol$Ap_phenol)
phenol$temp_scale<-scale(phenol$Avgtemp_61to125)

#We will use D later to place estimated (co)variances back onto the original scale.
D<-diag(c(attr(phenol$Ac_phenol_scale, "scaled:scale"),attr(phenol$Cp_phenol_scale, "scaled:scale"),attr(phenol$Ap_phenol_scale, "scaled:scale"),attr(phenol$temp_scale, "scaled:scale")))



#As you will see if you look at the data for phenology, many observations are missing. In this framework missing data are treated as missing at random conditional on the model. 
#model 1 - using default priors. This model sometimes takes a couple of attempts to run.
model<-1
while(length(model)==1){
try(model<-MCMCglmm(cbind(Ac_phenol_scale,Cp_phenol_scale,Ap_phenol_scale,temp_scale)~trait-1, random=~us(trait):Grid_150km+us(trait):Year,rcov=~us(trait):units,family=c("gaussian","gaussian","gaussian","gaussian"),nitt=23000,data=phenol))}

#Selecting appropriate priors for this model that were uninformative with respect to covariances was challenging. In the end the default priors were preferred. Reassuringly, very similar results twere obtained running the equivalent model in a frequentist framework using AsReml3 (Gilmore, A. R., B. J. Gogel, B. R. Cullis, and R. Thompson. 2009. Asreml user guide release 3.0. VSN International Ltd, Hemel Hempsted).


####################
#Here we project estimated (co)variance components back onto the original scale.
selectedmodel<-model
originalscaleVCV<-selectedmodel$VCV
for(x in 1:length(originalscaleVCV[,1])){
	C1<-matrix(originalscaleVCV[x,1:16],nrow=4)
	C2<-matrix(originalscaleVCV[x,17:32],nrow=4)
	C3<-matrix(originalscaleVCV[x,33:48],nrow=4)
	
	originalscaleVCV[x,]<-c(as.numeric(D%*%C1%*%D),as.numeric(D%*%C2%*%D),as.numeric(D%*%C3%*%D))
}


####Model checking######
autocorr(model$Sol)
autocorr(model$VCV)

plot(model$Sol)
plot(model$VCV)
#Most of these look ok. However, at the residual level the covariance between the phenology of different species shows a high level of autocorrelation in successive MCMC iterations. This is becuase there is no information in the data regarding the covariance between the two and all information comes from the prior.

summary(model)
#Note the small effective sample size for covariances among the residuals between species. 

########################
chosenmodelVCV<-originalscaleVCV
#chosenmodelVCV<-model1$VCV
#Specifiy whether you want to conduct remaining analyses on the scaled data used in the mdoel, or on the data that have been converted back to the original scale.

#obtain the posterior median covariance matrices over space and time.

#Space matrix - posterior median of VCV matrix among species and temperature.

spatial<-matrix(c(median(chosenmodelVCV[,1]),median(chosenmodelVCV[,2]),median(chosenmodelVCV[,3]),median(chosenmodelVCV[,4]),median(chosenmodelVCV[,5]),median(chosenmodelVCV[,6]),median(chosenmodelVCV[,7]),median(chosenmodelVCV[,8]),median(chosenmodelVCV[,9]),median(chosenmodelVCV[,10]),median(chosenmodelVCV[,11]),median(chosenmodelVCV[,12]),median(chosenmodelVCV[,13]),median(chosenmodelVCV[,14]),median(chosenmodelVCV[,15]),median(chosenmodelVCV[,16])),4,4)


#Time matrix - posterior median of VCV matrix among species and temperature.
temporal<-matrix(c(median(chosenmodelVCV[,17]),median(chosenmodelVCV[,18]),median(chosenmodelVCV[,19]),median(chosenmodelVCV[,20]),median(chosenmodelVCV[,21]),median(chosenmodelVCV[,22]),median(chosenmodelVCV[,23]),median(chosenmodelVCV[,24]),median(chosenmodelVCV[,25]),median(chosenmodelVCV[,26]),median(chosenmodelVCV[,27]),median(chosenmodelVCV[,28]),median(chosenmodelVCV[,29]),median(chosenmodelVCV[,30]),median(chosenmodelVCV[,31]),median(chosenmodelVCV[,32])),4,4)


#####################################################################################
#2a. Analyses of how plant phenology predicts butterfly phenology over space and time. See Fig. 4.
#####################################################################################

#butterfly and Cardamine pratensis

#slope of buterfly phenology on C. pratensis phenology over space
median(chosenmodelVCV[,2]/chosenmodelVCV[,6])
HPDinterval(chosenmodelVCV[,2]/chosenmodelVCV[,6])

#slope of buterfly phenology on C. pratensis phenology over time
median(chosenmodelVCV[,18]/chosenmodelVCV[,22])
HPDinterval(chosenmodelVCV[,18]/chosenmodelVCV[,22])

#slope difference - if the HPD interval does not overlap zero this is consistent with local adaptation of the butterfly phenology to C. pratensis phenology.
median(chosenmodelVCV[,2]/chosenmodelVCV[,6]-chosenmodelVCV[,18]/chosenmodelVCV[,22])
HPDinterval(chosenmodelVCV[,2]/chosenmodelVCV[,6]-chosenmodelVCV[,18]/chosenmodelVCV[,22])


#butterfly and Alliaria petiolata

#slope of buterfly phenology on A. petiolata phenology over space
median(chosenmodelVCV[,3]/chosenmodelVCV[,11])
HPDinterval(chosenmodelVCV[,3]/chosenmodelVCV[,11])

#slope of buterfly phenology on A. petiolata phenology over time
median(chosenmodelVCV[,19]/chosenmodelVCV[,27])
HPDinterval(chosenmodelVCV[,19]/chosenmodelVCV[,27])

#slope difference - if the HPD interval does not overlap zero this is consistent with local adaptation of the butterfly phenology to A. petiolata phenology.
median(chosenmodelVCV[,3]/chosenmodelVCV[,11]-chosenmodelVCV[,19]/chosenmodelVCV[,27])
HPDinterval(chosenmodelVCV[,3]/chosenmodelVCV[,11]-chosenmodelVCV[,19]/chosenmodelVCV[,27])


#####################################################################################
#2b. Multiple regression exploring the effects of temperature and plant phenology on butterfly phenology over space and time.
#####################################################################################

space.regress.coefs<-matrix(nrow=length(chosenmodelVCV[,1]),ncol=3)
time.regress.coefs<-matrix(nrow=length(chosenmodelVCV[,1]),ncol=3)

space.regress.coefs_noAp<-matrix(nrow=length(chosenmodelVCV[,1]),ncol=2)
time.regress.coefs_noAp<-matrix(nrow=length(chosenmodelVCV[,1]),ncol=2)
#chose to remove Alliara petiolata as high spatial correlation with Cardamine pratensis introduces collinearity.

full.cor.mat.space<-matrix(nrow=length(chosenmodelVCV[,1]),ncol=16)
full.cor.mat.time<-matrix(nrow=length(chosenmodelVCV[,1]),ncol=16)

for (post.samp in 1:length(chosenmodelVCV[,1])){
	
	space.mat<-matrix(chosenmodelVCV[post.samp,1:16],4,4)
	time.mat<-matrix(chosenmodelVCV[post.samp,17:32],4,4)
	
	###############
	space.regress.coefs[post.samp,]<-solve(space.mat[2:4,2:4], space.mat[1, 2:4])
	time.regress.coefs[post.samp,]<-solve(time.mat[2:4,2:4], time.mat[1, 2:4])
	space.regress.coefs_noAp[post.samp,]<-solve(space.mat[c(2,4),c(2,4)], space.mat[1, c(2,4)])
	time.regress.coefs_noAp[post.samp,]<-solve(time.mat[c(2,4),c(2,4)], time.mat[1, c(2,4)])
    
	###############
	
	full.cor.mat.space[post.samp,]<-as.numeric(cov2cor(space.mat))
	full.cor.mat.time[post.samp,]<-as.numeric(cov2cor(time.mat))
	
}

matrix(nrow=4,ncol=4,apply(full.cor.mat.space,2,median))
matrix(nrow=4,ncol=4,apply(full.cor.mat.time,2,median))
#The median of the posterior sample of correlation matrices

#############################################################
#Below we obtain the coefficients reported in Fig. 4 (includes both plant species and temperature).

#####Spatial######

c(median(mcmc(space.regress.coefs[,1])),median(mcmc(space.regress.coefs[,2])),median(mcmc(space.regress.coefs[,3])))
#the median correlation coefficients for C. pratensis, A. petiolata and temperature over space.

#the 95% HPD for each spatial correlation coefficient.
#C. pratensis as predictor over space
HPDinterval(mcmc(space.regress.coefs[,1]))

#A. petiolata as predictor over space
HPDinterval(mcmc(space.regress.coefs[,2]))

#temperature as predictor over space
HPDinterval(mcmc(space.regress.coefs[,3]))

#####Temporal######

c(median(mcmc(time.regress.coefs[,1])),median(mcmc(time.regress.coefs[,2])),median(mcmc(time.regress.coefs[,3])))
#the median correlation coefficients for C. pratensis, A. petiolata and temperature over time.

#the 95% HPD for each temporal correlation coefficient.
#C. pratensis as predictor over space
HPDinterval(mcmc(time.regress.coefs[,1]))	

#A. petiolata as predictor over time
HPDinterval(mcmc(time.regress.coefs[,2]))

#temperature as predictor over time
HPDinterval(mcmc(time.regress.coefs[,3]))


#####Delta b######

c(median(mcmc(space.regress.coefs[,1]-time.regress.coefs[,1])),median(mcmc(space.regress.coefs[,2]-time.regress.coefs[,2])),median(mcmc(space.regress.coefs[,3]-time.regress.coefs[,3])))

#95% delta b for C. pratensis
HPDinterval(mcmc(space.regress.coefs[,1]-time.regress.coefs[,1]))

#95% delta b for A. petiolata
HPDinterval(mcmc(space.regress.coefs[,2]-time.regress.coefs[,2]))

#95% delta b for temperature
HPDinterval(mcmc(space.regress.coefs[,3]-time.regress.coefs[,3]))

#############################################################
#Below we obtain the coefficients reported in Fig. 4 (includes only C. pratensis and temperature).

#####Spatial######
c(median(mcmc(space.regress.coefs_noAp[,1])),median(mcmc(space.regress.coefs_noAp[,2])))

#C. pratensis as predictor over space
HPDinterval(mcmc(space.regress.coefs_noAp[,1]))

#temperature as predictor over space
HPDinterval(mcmc(space.regress.coefs_noAp[,2]))


#####Temporal######

c(median(mcmc(time.regress.coefs_noAp[,1])),median(mcmc(time.regress.coefs_noAp[,2])))

#C. pratensis as predictor over space
HPDinterval(mcmc(time.regress.coefs_noAp[,1]))

#temperature as predictor over time
HPDinterval(mcmc(time.regress.coefs_noAp[,2]))


#####Delta b######

median(mcmc(space.regress.coefs_noAp[,1]-time.regress.coefs_noAp[,1]))
#95% delta b for C. pratensis
HPDinterval(mcmc(space.regress.coefs_noAp[,1]-time.regress.coefs_noAp[,1]))

median(mcmc(space.regress.coefs_noAp[,2]-time.regress.coefs_noAp[,2]))
#95% delta b for temperature
HPDinterval(mcmc(space.regress.coefs_noAp[,2]-time.regress.coefs_noAp[,2]))