#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