Electronic Supplementary Material (ESM 1) For

6

Wells et al., electronic supplementary material: BUGS model code

Electronic supplementary material (ESM 1) for

Inferring host specificity and network formation through agent-based models: tick – mammal interactions in Borneo

Konstans Wells · Robert B. O’Hara · Martin Pfeiffer · Maklarin B. Lakim · Trevor N. Petney · Lance A. Durden

Contact:

Description: Model code BUGS software (e.g. OpenBUGS as freely available at http://openbugs.info/w/) for the multi-species hierarchical model, estimating infestation probabilities for species-pair of ticks with different hosts and calculating indices of species specificity and random matrices of species interactions.

model

{

for(h in 1: NHostInd){

for(j in 1: NTickSpec){

# logistic regression for modeling individual-level interaction probability between host and ticks under consideration of covariates:

logit(psi[h,j]) <- vTickSp[j]

+ b.hsp[hostspec[h]]

+ b.htint[j, hostspec[h]]

+ b.hfa[j, hostfam[h]]

+ b.hor[j, hostord[h]]

+ b.hstrat[j, hoststrata[h]]

+ b.ar[j, area[h]]

+ b.ft[j, foresttype[h]]

+ b.y[j, year[h]]

+ b.m[j, month[h]]

+ b.hw[j] * hostweight[h]

+ b.locabund[j] * localabundance[h]

}

}

for(h in 1: NHostInd){

for(j in 1: NTickSpec){

# link of model to oberserved data with obs[i,j]being the observed interaction between host individual h and tick species j

obs[h,j] ~ dbern(psi[h,j])

}

}

# Priors

# Imputing missing data of host weight

for(h in 1: NHostInd){

hostweight[h] ~ dnorm(0,1)

}

for(i in 1: NHostSpec){

b.hsp.star[i] <- b.hsp[i] - mean(b.hsp[])

# note: mu for mu.hsp[i] allowed to vary across host species

b.hsp[i] ~ dnorm(mu.hsp[i], tau.hsp)

mu.hsp[i] ~ dnorm(0,1)

}

for(j in 1: NTickSpec){

vTickSp.Star[j] <- vTickSp[j] - mean(vTickSp[])

# note: mu for mu.TickSp[j] allowed to vary across tick species

vTickSp[j] ~ dnorm(mu.TickSp[j], tau.TickSp)

mu.TickSp[j] ~ dnorm(0,1)

for(i in 1: NHostSpec){

b.htint.star[j,i] <- b.htint[j,i] - mean(b.htint[,])

b.htint[j,i] ~ dnorm(mu.htint, tau.htint)

}

for(f in 1: NHostFam){

b.hfa.star[j,f] <- b.hfa[j,f] - mean(b.hfa[j,])

b.hfa[j,f] ~ dnorm(mu.hfa, tau.hfa)

}

for(o in 1: NHostOrd){

b.hor.star[j,o] <- b.hor[j,o] - mean(b.hor[j,])

b.hor[j,o] ~ dnorm(mu.hor, tau.hor)

}

for(s in 1: NHostStrata){

b.hstrat.star[j,s] <- b.hstrat[j,s] - mean(b.hstrat[j,])

b.hstrat[j,s] ~ dnorm(mu.hstrat, tau.hstrat)

}

for(a in 1: NArea){

b.ar.star[j,a] <- b.ar[j,a] - mean(b.ar[j,])

b.ar[j,a] ~ dnorm(mu.ar, tau.ar)

}

for(f in 1: NForTyp){

b.ft.star[j,f] <- b.ft[j,f] - mean(b.ft[j,])

b.ft[j,f] ~ dnorm(mu.ft, tau.ft)

}

for(y in 1: NYear){

b.y.star[j,y] <- b.y[j,y] - mean(b.y[j,])

b.y[j,y] ~ dnorm(mu.y, tau.y)

}

for(m in 1: NMonth){

b.m.star[j,m] <- b.m[j,m] - mean(b.m[j,])

b.m[j,m] ~ dnorm(mu.m, tau.m)

}

b.hw[j] ~ dnorm(mu.hw, tau.hw)

b.locabund[j] ~ dnorm(mu.locabund, tau.locabund)

}

tau.TickSp <- pow(sd.TickSp, -2); sd.TickSp ~ dunif(0,10)

tau.hsp <- pow(sd.hsp, -2); sd.hsp ~ dunif(0,10)

tau.htint <- pow(sd.htint, -2); sd.htint ~ dunif(0,10)

tau.hfa <- pow(sd.hfa, -2); sd.hfa ~ dunif(0,10)

tau.hor <- pow(sd.hor, -2); sd.hor ~ dunif(0,10)

tau.hstrat <- pow(sd.hstrat, -2); sd.hstrat ~ dunif(0,10)

tau.ar <- pow(sd.ar, -2); sd.ar ~ dunif(0,10)

tau.ft <- pow(sd.ft, -2); sd.ft ~ dunif(0,10)

tau.y <- pow(sd.y, -2); sd.y ~ dunif(0,10)

tau.m <- pow(sd.m, -2); sd.m ~ dunif(0,10)

tau.hw <- pow(sd.hw, -2); sd.hw ~ dunif(0,10)

tau.locabund <- pow(sd.locabund, -2); sd.locabund ~ dunif(0,10)

mu.htint ~ dnorm(0,1)

mu.hfa ~ dnorm(0,1)

mu.hor ~ dnorm(0,1)

mu.hstrat ~ dnorm(0,1)

mu.ar ~ dnorm(0,1)

mu.ft ~ dnorm(0,1)

mu.y ~ dnorm(0,1)

mu.m ~ dnorm(0,1)

mu.hw ~ dnorm(0,1)

mu.locabund ~ dnorm(0,1)

################

# Calculate variance components

################

Var.vTickSp <- pow(sd(vTickSp.Star[]), 2)

Var.b.hsp <- pow(sd(b.hsp.star[]), 2)

propVar.vTickSp <- Var.vTickSp / mean(Var.psi.total[])

propVar.b.hsp <- Var.b.hsp / mean(Var.psi.total[])

for(j in 1: NTickSpec){

Var.b.htint[j] <- pow(sd(b.htint.star[j,]), 2)

Var.b.hfa[j] <- pow(sd(b.hfa.star[j,]), 2)

Var.b.hor[j] <- pow(sd(b.hor.star[j,]), 2)

Var.b.hstrat[j] <- pow(sd(b.hstrat.star[j,]), 2)

Var.b.ar[j] <- pow(sd(b.ar.star[j,]), 2)

Var.b.ft[j] <- pow(sd(b.ft.star[j,]), 2)

Var.b.y[j] <- pow(sd(b.y.star[j,]), 2)

Var.b.m[j] <- pow(sd(b.m.star[j,]), 2)

Var.b.hw[j] <- b.hw[j] * b.hw[j]

Var.b.locabund[j] <- b.locabund[j] * b.locabund[j]

Var.psi.total[j] <- Var.vTickSp + Var.b.hsp + Var.b.htint[j] + Var.b.hfa[j] + Var.b.hor[j] + Var.b.hstrat[j] + Var.b.ar[j] + Var.b.ft[j] + Var.b.y[j] + Var.b.m[j] + Var.b.hw[j] + Var.b.locabund[j]

propVar.b.htint[j] <- Var.b.htint[j] / Var.psi.total[j]

propVar.b.hfa[j] <- Var.b.hfa[j] / Var.psi.total[j]

propVar.b.hor[j] <- Var.b.hor[j] / Var.psi.total[j]

propVar.b.hstrat[j] <- Var.b.hstrat[j] / Var.psi.total[j]

propVar.b.ar[j] <- Var.b.ar[j] / Var.psi.total[j]

propVar.b.ft[j] <- Var.b.ft[j] / Var.psi.total[j]

propVar.b.y[j] <- Var.b.y[j] / Var.psi.total[j]

propVar.b.m[j] <- Var.b.m[j] / Var.psi.total[j]

propVar.b.hw[j] <- Var.b.hw[j] / Var.psi.total[j]

propVar.b.locabund[j] <- Var.b.locabund[j] / Var.psi.total[j]

}

################

# Calculate mean infestation probability / prevalence for each species pair

################

for (i in 1: NHostSpec){

for(j in 1: NTickSpec){

logit(psiMean.sp[i,j]) <- vTickSp[j]

+ b.hsp[hostspec.sp[i]]

+ b.htint[j,hostspec.sp[i]]

+ b.hfa[j, hostfam.sp[i]]

+ b.hor[j, hostord.sp[i]]

+ b.hstrat[j, hoststrata.sp[i]]

+ mean(b.ar[j,])

+ mean(b.ft[j,])

+ mean(b.y[j,])

+ mean(b.m[j,])

+ b.hw[j] * hostweight.sp[i]

+ b.locabund[j] * localabundance.sp[i]

}

}

################

# Calculate indices of host/species specificity and interactions based on Rao’s quadratic entropy

################

NHostTotal <- sum(NAbundSpec[])

for(j in 1: NTickSpec){

Q.tick[j] <- sum(Q.ind[j,,])

Q.tax.tick[j] <- sum(Q.tax.ind[j,,])

Q.abund.tick[j] <- sum(Q.abund.ind[j,,])

Q.tax.abund.tick[j] <- sum(Q.tax.abund.ind[j,,])

# Ratio of sum of weighted taxonomic distinctness to weighting factor

S.TD[j] <- Q.tax.tick[j] / Q.tick[j]

# Ratio of sum of weighted host availability to weighting factor

S.AB[j] <- Q.abund.tick[j] / Q.tick[j]

S.TD.2[j] <- Q.tax.abund.tick[j] / Q.abund.tick[j]

S.AB.2[j] <- Q.tax.abund.tick[j] / Q.tax.tick[j]

S.TDAB[j] <- Q.tax.abund.tick[j] / Q.tick[j]

# Shannon index of host diversity

Shannon.tick[j] <- -1 * sum(Shannon.ind[j,])

for (i in 1: NHostSpec){

Shannon.ind[i,j] <- psiMean.sp[i,j] * log(psiMean.sp[i,j])

for (s in 1: NHostSpec){

# Calculate quadradic entropy term fpr pairs of host species

Q.ind[j,i,s] <- psiMean.sp[i,j] * psiMean.sp[s,j]

# quadradic entropy term weighted by taxonomic distance for pairs of host species

Q.tax.ind[j,i,s] <- psiMean.sp[i,j] * psiMean.sp[s,j] * host.taxdis[i,s]

# quadradic entropy term weighted by host availability for pairs of host species

Q.abund.ind[j,i,s] <- psiMean.sp[i,j] * psiMean.sp[s,j] * NAbundSpec[i] * NAbundSpec[s] / (NHostTotal * NHostTotal)

# quadradic entropy term weighted by taxonomic distance AND host availability for pairs of host species

Q.tax.abund.ind[j,i,s] <- psiMean.sp[i,j] * psiMean.sp[s,j] * host.taxdis[i,s] * NAbundSpec[i] * NAbundSpec[s] / (NHostTotal * NHostTotal)

}

}

}

################

## Calculate random interactions matrices based on binomial distribution and two sampling scenarios (“FieldScenario” and “ExhaustiveSample”)

################

for(j in 1: NTickSpec){

# Number of expected host associations under field-based scenario (“FieldScenario”)

NoHostAssoc.NAbundSpec[j] <- sum(assoc.NAbundSpec[,j])

# Number of expected host associations under exhaustive sampling scenario (“ExhaustiveSample”)

NoHostAssoc.Equal[j] <- sum(assoc.Equal[,j])

NoHostAssoc.pa[j] <- sum(a.pa[,j])

for (i in 1: NHostSpec){

# Random draw of host-tick associations based on average interaction frequencies and relative observed abundances (“FieldScenario”)

F.NAbundSpec[i,j] ~ dbin(psiMean.sp[i,j], NAbundSpec[i])

assoc.NAbundSpec[i,j] <- step(F.NAbundSpec[i,j]-1)

# Random draw of host-tick associations based on average interaction frequencies and equally large abundances of all host species “ExhaustiveSample”)

F.Equal[i,j] ~ dbin(psiMean.sp[i,j], 1000)

assoc.Equal[i,j] <- step(F.Equal[i,j]-1)

a.pa[i,j] ~ dbern(psiMean.sp[i,j])

}

}

for (i in 1: NHostSpec){

NoTickAssoc.pa[i] <- sum(a.pa[i,])

}

}