Post date: May 16, 2018 10:45:26 PM
Patrik and Clarissa conducted an experiment that involved transplanting T. chumash to mountain mahogany (MM) or adenostoma plus ceanothus (AC) to test the hypothesis that the more continuous color variation on the former did not select against intermediate Timema whereas the more discrete color variation on the latter did. There were two blocks and for each 20 green, intermediate and melanistic T. chumash were released on MM and AC (note that for the later plants growing together were used).
I fit a multinomial-Dirichlet model for the recapture data using JAGS. This was done with either a constraint forcing the same parameter values for the two blocks and with a model allowing each host and block to have unique relative fitness values. The former model was preferred by DIC and as predicted we found higher relative fitness for intermediates on MM.
The R code (hetExperiment.R) below and all of the other files are at king:/uufs/chpc.utah.edu/common/home/u6000989/projects/timema_genarch/color_exp/. Note that we focused on the non-standardized relative fitness values.
## data
MM1<-c(4,5,2)
MM2<-c(4,8,3)
CA1<-c(2,0,2)
CA2<-c(6,2,4)
y<-rbind(MM1,MM2,CA1,CA2)
n<-apply(y,1,sum)
nc<-3
burn<-1000
mcmc<-9000
thin<-3
mod<-jags.model("model_uncon",data=list(y=y,n=n,N=4,a=c(0.5,0.5,0.5)),n.chains=nc)
update(mod,n.iter=burn)
out<-coda.samples(model=mod,variable.names=c("p"),n.iter=mcmc,thin=thin)
diccon<-dic.samples(model=mod,n.iter=mcmc,thin=thin,type="pD")
## treatment
mod<-jags.model("model_trt",data=list(y=y,n=n,a=c(0.5,0.5,0.5)),n.chains=nc)
update(mod,n.iter=burn)
out<-coda.samples(model=mod,variable.names=c("p1","p2"),n.iter=mcmc,thin=thin)
dictrt<-dic.samples(model=mod,n.iter=mcmc,thin=thin,type="pD")
## unique
mod<-jags.model("model_uni",data=list(y=y,n=n,a=c(0.5,0.5,0.5)),n.chains=nc)
update(mod,n.iter=burn)
out<-coda.samples(model=mod,variable.names=c("p1","p2","p3","p4"),n.iter=mcmc,thin=thin)
dicuni<-dic.samples(model=mod,n.iter=mcmc,thin=thin,type="pD")
## treatment downstream
mod<-jags.model("model_trt",data=list(y=y,n=n,a=c(0.5,0.5,0.5)),n.chains=nc)
update(mod,n.iter=burn)
out<-coda.samples(model=mod,variable.names=c("p1","p2"),n.iter=mcmc,thin=thin)
outM<-matrix(NA,nrow=3000*3,ncol=6)
for(i in 1:6){
outM[,i]<-c(out[[1]][,i],out[[2]][,i],out[[3]][,i])
}
apply(outM,2,quantile,probs=c(.5,.025,.975))
mn<-apply(outM,2,mean)
s<-apply(outM,2,sd)
seup<-mn+s
sedwn<-mn-s
mean(outM[,2] > outM[,5])
## 0.995333
#####################################3
## plot for ms ###
MM<-MM1+MM2
CA<-CA1+CA2
UB<-22
cs<-c("darkgreen","lightyellow3","brown")
nms<-c("Green","Int.","Melan.")
pdf("fig_experiment.pdf",width=6,height=9)
par(mfrow=c(3,2))
par(mar=c(4.5,5.5,2.5,1.5))
barplot(c(20,20,20),col=cs,ylim=c(0,UB),ylab="Number",cex.lab=1.75,names.arg=nms,cex.names=1.3)
title(main="(a) Release MM",cex.main=1.5)
barplot(c(20,20,20),col=cs,ylim=c(0,UB),ylab="Number",cex.lab=1.75,names.arg=nms,cex.names=1.3)
title(main="(b) Release A/C",cex.main=1.5)
barplot(MM,col=cs,ylim=c(0,UB),ylab="Number",cex.lab=1.75,names.arg=nms,cex.names=1.3)
title(main="(c) Recapture MM",cex.main=1.5)
barplot(CA,col=cs,ylim=c(0,UB),ylab="Number",cex.lab=1.75,names.arg=nms,cex.names=1.3)
title(main="(d) Recapture A/C",cex.main=1.5)
x<-barplot(mn[1:3],col=cs,ylim=c(0,1),ylab="Relative fitness",cex.lab=1.75,names.arg=nms,cex.names=1.3)
segments(x,seup[1:3],x,sedwn[1:3])
title(main="(e) Fitness MM",cex.main=1.5)
x<-barplot(mn[4:6],col=cs,ylim=c(0,1),ylab="Relative fitness",cex.lab=1.75,names.arg=nms,cex.names=1.3)
segments(x,seup[4:6],x,sedwn[4:6])
title(main="(f) Fitness A/C",cex.main=1.5)
dev.off()
###############################333
outMstan<-outM
for(i in 1:3){
outMstan[,i]<-outM[,i]/outM[,1]
}
for(i in 4:6){
outMstan[,i]<-outM[,i]/outM[,4]
}
mean(outMstan[,2] > outMstan[,5])
## 0.988
apply(outMstan,2,quantile,probs=c(.5,.025,.975))
# [,1] [,2] [,3] [,4] [,5] [,6]
#50% 1 0.4892992 0.1938114 1 0.12984869 0.3658709
#2.5% 1 0.3101914 0.0756813 1 0.02591303 0.1675710
#97.5% 1 0.6740828 0.3628521 1 0.33833922 0.6060812