Code from:

Wood, J.L.A, Yates, M.C., and D.J. Fraser. 2016. Are heritability and selection related to population size
in nature? Meta-analysis and conservation implications. Evolutionary Applications.

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


library(MCMCglmm)

############SELECTION ANALYSIS
############CODE USED TO ESTIMATE THE MAGNITUDE OF SELECTION USING THE FOLDED NORMAL DISTRIBUTION
############POPULATION SIZE BIN ONLY (META-ANALYSIS)


prior1=list(R=list(V=diag(4),n=0.002),G=list(G1=list(V=diag(4),n=0.002),G2=list(V=diag(4),n=0.002)))
mev=grad.lin.SE^2

model1=MCMCglmm(grad.lin.value~pop.bin,mev=mev,data=data,prior=prior1,random=~idh(pop.bin):study.ID+idh(pop.bin):pop.ID,rcov=~idh(pop.bin):units,nitt=600000,burnin=100000,thin=50,family="gaussian")
posterior.mode(model1$VCV)
posterior.mode(model1$Sol)
HPDinterval(model1$Sol)
summary(model1)
plot(model1$Sol)
plot(model1$VCV)
autocorr(model1$VCV)

effectiveSize(model1$Sol)
effectiveSize(model1$VCV)
heidel.diag(model1$VCV)

model1$Sol[1:10,]
model1$VCV[1:10,]

m1 = c(1, 5, 10)

mu.fnorm<-function(mu,sigma){dnorm(mu,0,sigma)*2*sigma^2+mu*(2*pnorm(mu,0,sigma)-1)}
newModeModel1 = mu.fnorm(model1$Sol[,1], sqrt(rowSums(model1$VCV[,m1])))
HPDinterval(newModeModel1)
posterior.mode(newModeModel1)

var.fnorm<-function(mu,sigma){mu^2+sigma^2-(sigma*sqrt(2/pi)*exp((-1*mu^2)/(2*sigma^2))+mu*(1-2*pnorm(-1*mu/sigma,0,1)))^2}
newVarModel1 = var.fnorm(model1$Sol[,1], sqrt(rowSums(model1$VCV[,m1])))
HPDinterval(newVarModel1)
posterior.mode(newVarModel1)





###########POPULATION SIZE BIN ONLY (UNWEIGHTED ANALYSIS)


prior1=list(R=list(V=diag(4),n=0.002),G=list(G1=list(V=diag(4),n=0.002),G2=list(V=diag(4),n=0.002)))

model1=MCMCglmm(grad.lin.value~pop.bin,data=data,prior=prior1,random=~idh(pop.bin):study.ID+idh(pop.bin):pop.ID,rcov=~idh(pop.bin):units,nitt=600000,burnin=100000,thin=50,family="gaussian")
posterior.mode(model1$VCV)
posterior.mode(model1$Sol)
HPDinterval(model1$Sol)
summary(model1)
plot(model1$Sol)
plot(model1$VCV)
autocorr(model1$VCV)

effectiveSize(model1$Sol)
effectiveSize(model1$VCV)
heidel.diag(model1$VCV)

model1$Sol[1:10,]
model1$VCV[1:10,]

m1 = c(1, 5, 10)

mu.fnorm<-function(mu,sigma){dnorm(mu,0,sigma)*2*sigma^2+mu*(2*pnorm(mu,0,sigma)-1)}
newModeModel1 = mu.fnorm(model1$Sol[,1], sqrt(rowSums(model1$VCV[,m1])))
HPDinterval(newModeModel1)
posterior.mode(newModeModel1)

var.fnorm<-function(mu,sigma){mu^2+sigma^2-(sigma*sqrt(2/pi)*exp((-1*mu^2)/(2*sigma^2))+mu*(1-2*pnorm(-1*mu/sigma,0,1)))^2}
newVarModel1 = var.fnorm(model1$Sol[,1], sqrt(rowSums(model1$VCV[,m1])))
HPDinterval(newVarModel1)
posterior.mode(newVarModel1)




############POPULATION SIZE BIN X TRAIT CLASS INTERACTION (META-ANALYSIS)

prior4=list(R=list(V=diag(8),n=0.002),G=list(G1=list(V=diag(8),n=0.002),G2=list(V=diag(8),n=0.002)))
mev=grad.lin.SE^2

model4=MCMCglmm(grad.lin.value~pop.bin*trait.class,mev=mev,data=data,prior=prior4,random=~idh(pop.bin:trait.class):study.ID+idh(pop.bin:trait.class):pop.ID,rcov=~idh(pop.bin:trait.class):units,nitt=600000,burnin=100000,thin=50,family="gaussian")
posterior.mode(model4$VCV)
posterior.mode(model4$Sol)
HPDinterval(model4$Sol)
summary(model4)
plot(model4$Sol)
plot(model4$VCV)
autocorr(model4$VCV)

effectiveSize(model1$Sol)
effectiveSize(model1$VCV)
heidel.diag(model1$VCV)

model4$Sol[1:10,]
model4$VCV[1:10,]

m4 = c(1, 9, 18)

mu.fnorm<-function(mu,sigma){dnorm(mu,0,sigma)*2*sigma^2+mu*(2*pnorm(mu,0,sigma)-1)}
newModeModel4 = mu.fnorm(model4$Sol[,1], sqrt(rowSums(model4$VCV[,m4])))
HPDinterval(newModeModel4)
posterior.mode(newModeModel4)

var.fnorm<-function(mu,sigma){mu^2+sigma^2-(sigma*sqrt(2/pi)*exp((-1*mu^2)/(2*sigma^2))+mu*(1-2*pnorm(-1*mu/sigma,0,1)))^2}
newVarModel4 = var.fnorm(model4$Sol[,1], sqrt(rowSums(model4$VCV[,m4])))
HPDinterval(newVarModel4)
posterior.mode(newVarModel4)




############POPULATION SIZE BIN X TRAIT CLASS INTERACTION (UNWEIGHTED ANALYSIS)

prior4=list(R=list(V=diag(8),n=0.002),G=list(G1=list(V=diag(8),n=0.002),G2=list(V=diag(8),n=0.002)))

model4=MCMCglmm(grad.lin.value~pop.bin*trait.class,data=data,prior=prior4,random=~idh(pop.bin:trait.class):study.ID+idh(pop.bin:trait.class):pop.ID,rcov=~idh(pop.bin:trait.class):units,nitt=600000,burnin=100000,thin=50,family="gaussian")
posterior.mode(model4$VCV)
posterior.mode(model4$Sol)
HPDinterval(model4$Sol)
summary(model4)

plot(model4$Sol)
plot(model4$VCV)

model4$Sol[1:10,]
model4$VCV[1:10,]

m4 = c(1, 9, 18)

mu.fnorm<-function(mu,sigma){dnorm(mu,0,sigma)*2*sigma^2+mu*(2*pnorm(mu,0,sigma)-1)}
newModeModel4 = mu.fnorm(model4$Sol[,1], sqrt(rowSums(model4$VCV[,m4])))
HPDinterval(newModeModel4)
posterior.mode(newModeModel4)

var.fnorm<-function(mu,sigma){mu^2+sigma^2-(sigma*sqrt(2/pi)*exp((-1*mu^2)/(2*sigma^2))+mu*(1-2*pnorm(-1*mu/sigma,0,1)))^2}
newVarModel4 = var.fnorm(model4$Sol[,1], sqrt(rowSums(model4$VCV[,m4])))
HPDinterval(newVarModel4)
posterior.mode(newVarModel4)




############POPULATION SIZE BIN X TAXA INTERACTION (META-ANALYSIS)

prior1=list(R=list(V=diag(8),n=0.002),G=list(G1=list(V=diag(8),n=0.002),G2=list(V=diag(8),n=0.002)))
mev=grad.lin.SE^2

model1=MCMCglmm(grad.lin.value~pop.bin*taxon,mev=mev,data=data,prior=prior1,random=~idh(pop.bin:taxon):study.ID+idh(pop.bin:taxon):pop.ID,rcov=~idh(pop.bin:taxon):units,nitt=600000,burnin=100000,thin=50,family="gaussian")
posterior.mode(model1$VCV)
posterior.mode(model1$Sol)
HPDinterval(model1$Sol)
summary(model1)
plot(model1$Sol)
plot(model1$VCV)
autocorr(model1$Sol)

effectiveSize(model1$Sol)
effectiveSize(model1$VCV)
heidel.diag(model1$VCV)
heidel.diag(model1$Sol)

model1$Sol[1:10,]
model1$VCV[1:10,]

m1 = c(1, 9, 18)

mu.fnorm<-function(mu,sigma){dnorm(mu,0,sigma)*2*sigma^2+mu*(2*pnorm(mu,0,sigma)-1)}
newModeModel1 = mu.fnorm(model1$Sol[,1], sqrt(rowSums(model1$VCV[,m1])))
HPDinterval(newModeModel1)
posterior.mode(newModeModel1)

var.fnorm<-function(mu,sigma){mu^2+sigma^2-(sigma*sqrt(2/pi)*exp((-1*mu^2)/(2*sigma^2))+mu*(1-2*pnorm(-1*mu/sigma,0,1)))^2}
newVarModel1 = var.fnorm(model1$Sol[,1], sqrt(rowSums(model1$VCV[,m1])))
HPDinterval(newVarModel1)
posterior.mode(newVarModel1)



############POPULATION SIZE BIN X TAXA INTERACTION (UNWEIGHTED ANALYSIS)

prior1=list(R=list(V=diag(8),n=0.002),G=list(G1=list(V=diag(8),n=0.002),G2=list(V=diag(8),n=0.002)))

model1=MCMCglmm(grad.lin.value~pop.bin*taxon,data=data,prior=prior1,random=~idh(pop.bin:taxon):study.ID+idh(pop.bin:taxon):pop.ID,rcov=~idh(pop.bin:taxon):units,nitt=600000,burnin=100000,thin=50,family="gaussian")
posterior.mode(model1$VCV)
posterior.mode(model1$Sol)
HPDinterval(model1$Sol)
summary(model1)
plot(model1$Sol)
plot(model1$VCV)
autocorr(model1$Sol)

effectiveSize(model1$Sol)
effectiveSize(model1$VCV)
heidel.diag(model1$VCV)
heidel.diag(model1$Sol)

model1$Sol[1:10,]
model1$VCV[1:10,]

m1 = c(1, 9, 18)

mu.fnorm<-function(mu,sigma){dnorm(mu,0,sigma)*2*sigma^2+mu*(2*pnorm(mu,0,sigma)-1)}
newModeModel1 = mu.fnorm(model1$Sol[,1], sqrt(rowSums(model1$VCV[,m1])))
HPDinterval(newModeModel1)
posterior.mode(newModeModel1)

var.fnorm<-function(mu,sigma){mu^2+sigma^2-(sigma*sqrt(2/pi)*exp((-1*mu^2)/(2*sigma^2))+mu*(1-2*pnorm(-1*mu/sigma,0,1)))^2}
newVarModel1 = var.fnorm(model1$Sol[,1], sqrt(rowSums(model1$VCV[,m1])))
HPDinterval(newVarModel1)
posterior.mode(newVarModel1)