#Appendix S2: JAGS (BUGS) code used to implement the model in Appendix S1 from the paper
#Gene flow relates to evolutionary divergence among populations at the range margin

## BUGS code (implemented in JAGS; Plummer 2003)
# the code is divided into two models
# the first model (Appendix S2a) is the calibration model that used morphological data from the same individuals 
# to convert measurements made by callipers to those made on digital images
# the outputs from this model are then fed into the GLMM model (Appendix S2b) as a measure of observational uncertainty 

# models were written and implemented by ML 

###### APPENDIX S2A #########
 ### calibration model ###

#data supplied to model
#nobs.A = number of observations
#digital.A = measures of trait A from digital images
#caliper.A = measures of trait A from callipers from same individual  
#convert.A = the conversion factor (as a simple additive measure) we are estimating

model{

#likelihood and process model for morphological trait A	
for(i in 1:nobs.A){

digital.A[i]~dnorm(mu.A[i],tau.A)
mu.A[i]<-caliper.A[i]+convert.A
	
}
	 
#priors
tau.A<-1/sigma.A^2 
sigma.A~dunif(0,10)
convert.A~dnorm(0,0.001)

}  #close model 1a

###### APPENDIX 2B #########
 ### GLMM model relating trait size to isolation, origin, sex & year ###

#nobs = number of observations
#y = observations of the trait of interest
#z = 'true' unobservable value of the trait (to this the observational error of 'convert' is added from model S1A) 
#year = data on year to account for between year environmental effects on trait sizes
#iso = if from continuous or isolated populations
#field = if collected from the wild or lab-raised
#sex = male or female (NOTE this term and its interactions removed for sex-specific traits)
#site = a 4 category variable defining the 4 sites where individuals were collected
#n.sites = the number of sites (4)

model{
	
for(i in 1:nobs){

#likelihood
	
	y[i]~dnorm(z[i]+convert[i],tau.convert[i])

	z[i]~dnorm(mu[i],tau)

#process model
	
	mu[i]<-alpha[site[i]] 
	+ b.year*year[i] 
	+ b.field*field[i] 
	+ b.iso*iso[i] 
	+ b.sex*sex[i] 
	+ int.field.iso*field[i]*iso[i] 
	+ int.field.sex*field[i]*sex[i] 
	+ int.iso.sex*iso[i]*sex[i] 
	+ int.3way*field[i]*sex[i]*iso[i] 
	
	tau.convert[i]<-1/sigma.convert[i]^2
	
}

for(j in 1:n.sites){

alpha[j] ~ dnorm(5, 0.001)

}
	
#priors	
b.year~dnorm(0, 0.001)
b.iso~dnorm(0, 0.001)	
b.field~dnorm(0, 0.001)
b.sex~dnorm(0, 0.001)
int.field.iso~dnorm(0, 0.001)
int.field.sex~dnorm(0, 0.001)	
int.iso.sex~dnorm(0, 0.001)
int.3way~dnorm(0, 0.001)
tau<-1/sigma^2
sigma~dunif(0,20)	

} #close model 1b