# If applying part or all of this code please cite both the Dryad package and the original article. # If you have any questions or concerns, please email me (erik.svensson@biol.lu.se) # Signed: Erik Svensson #Load packages library(Matrix) library(MASS) library(ggplot2) library(MCMCglmm) #LOAD "FieldSelection"-data (edit to match own directory address) dat<-read.csv(file = "~/YOUR DIRECTORY/FieldSelection.csv", sep = ",", header = T) #Bayesian analysis of relative fitnesses and selection coefficients from three years of field data (2015-2017) #The code below will re-generate Fig. 4A and 4B in the manuscript #We also provide code for posthoc-tests for morph differences in relative fitnesses and selection coefficients IJ <- (1/3) * (diag(2) + matrix(1, 2, 2)) #residual covariance matrix 1/J*(I+J), I and J are identity and unit matrices respectively prior = list(R = list(V = IJ, fix = 1)) prior$B = list(mu = rep(0, 4), V = kronecker(IJ, diag(2)) * (1.7 + pi^2/3)) # kronecker prior close to being flat for the 2 way (A vs. I, I vs. O and A vs. O) marginal probabilities within each age prior$G = list(G1=list(V=1, nu=0.002), G2=list(V=1, nu=0.002)) #double check adequacy of univariate inverse Wishart prior with variance set to 1 and low degree of belief m1<- MCMCglmm(Morph ~ at.level(Age, 1):trait + at.level(Age, 2):trait - 1, random = ~us(at.level(trait,1)):Locale + us(at.level(trait,2)):Locale, rcov = ~us(trait):units, data = dat, family = "categorical", prior = prior, nitt = 50000, burnin = 10000, thin = 1000, #increase number of iterations to get the fuzzy worm! pr=TRUE, saveX = TRUE, verbose=TRUE) #nitt = 2100000, burnin = 100000, thin = 1000, #increase number of iterations to get the fuzzy worm! plot(m1) #rescale the intercepts as if the residual covariance matrix was zero Delta <- cbind(c(-1, 1, 0), c(-1, 0, 1)) c2 <- (16 * sqrt(3)/(15 * pi))^2 D <- ginv(Delta %*% t(Delta)) %*% Delta Int1 <- t(apply(m1$Sol[, 1:2], 1, function(x) { D %*% (x/sqrt(1 + c2 * diag(IJ)) ) })) summary(mcmc(exp(Int1)/rowSums(exp(Int1)))) Afreq<- (mcmc(exp(Int1)/rowSums(exp(Int1)))) #morph frequencies in adults Int2 <- t(apply(cbind(m1$Sol[, 3:4]), 1, function(x) { D %*% (x/sqrt(1 + c2 * diag(IJ)) ) })) summary(mcmc(exp(Int2)/rowSums(exp(Int2)))) Jfreq<-(mcmc(exp(Int2)/rowSums(exp(Int2)))) #morph frequencies in juveniles MF<-cbind(Jfreq,Afreq) #matrix of posterior distributions from m1 head(MF) relfit<-Afreq/Jfreq # matrix of relative fitness by morph fit<-c(relfit[,1],relfit[,2],relfit[,3]) Morph<-c(rep("A",nrow(Afreq)), rep("I",nrow(Afreq)),rep("O",nrow(Afreq))) freqpost<-data.frame(Morph,fit) dev.off() P<-ggplot (freqpost, aes( x=fit ,fill=Morph, color=Morph)) + geom_density(alpha=0.5)+ scale_color_manual(values=c("Dodgerblue","darkgoldenrod2", "tomato3"),labels=c("A","I","O"))+ scale_fill_manual(values=c("Dodgerblue","darkgoldenrod2", "tomato3"))+ xlab(label = "Cross-sectional relative morph fitness (w)")+ ylab(label = "Posterior density")+ theme_bw()+ #theme(legend.position="none")+ theme(panel.grid.minor.x=element_blank(), panel.grid.major.x=element_blank())+ theme(panel.grid.minor.y=element_blank(), panel.grid.major.y=element_blank())+ theme(text = element_text(size=16)) P #This will generate Fig. 4A in hte manuscript #Calculate selection coefficient as delta p/(p'*(1-p)-p*r*q)) using the entire posterior. Maybe sampling from the posterior is more appropriate? SelA<-(MF[,4]-MF[,1])/(MF[,4]*(1-MF[,1])-MF[,1]*MF[,2]*MF[,3]) SelI<-(MF[,5]-MF[,2])/(MF[,5]*(1-MF[,2])-MF[,1]*MF[,2]*MF[,3]) SelO<-(MF[,6]-MF[,3])/(MF[,6]*(1-MF[,3])-MF[,1]*MF[,2]*MF[,3]) s<-c(SelA,SelI,SelO) freqpost2<-data.frame(Morph,s) P2<-ggplot (freqpost2, aes( x=s ,fill=Morph, color=Morph)) + geom_density(alpha=0.5)+ scale_color_manual(values=c("Dodgerblue","darkgoldenrod2", "tomato3"))+ scale_fill_manual(values=c("Dodgerblue","darkgoldenrod2", "tomato3"))+ xlab(label = "Cross-sectional selection coefficient (s)")+ ylab(label = "Posterior density")+ theme_bw()+ #theme(legend.position="none")+ theme(panel.grid.minor.x=element_blank(), panel.grid.major.x=element_blank())+ theme(panel.grid.minor.y=element_blank(), panel.grid.major.y=element_blank())+ theme(text = element_text(size=16)) P2 #This will generate Fig. 4B in hte manuscript #Post-hoc tests of difference in selection coefficients from subtracting the posterior distributions #Need to do this for the relative fitness-distributions as well sum((SelA-SelI)<0*2) #How many times out of all runs are SelI greater than SelA? sum((SelO-SelI)<0*2) #How many times out of all runs are SelI greater than SelO? sum((SelA-SelO)<0*2) #How many times out of all runs are SelO greater than SelA? length(SelA) #Length of the vector SelA (the number above should be divided by this number) length(SelO) #Length of the vector SelO (the number above should be divided by this number) length(SelI) #Length of the vector SelI (the number above should be divided by this number) sum((SelI-SelO)<0*2)/40 #Formal posthoc-tests sum((SelO-SelI)<0*2)/40 sum((SelI-SelA)<0*2)/40 sum((SelA-SelI)<0*2)/40 sum((SelO-SelA)<0*2)/40 sum((SelA-SelO)<0*2)/40