---
title: "Summer 18"
author: "Saley"
date: "3/29/2019"
output: pdf_document
---
**Running bayesian model for 8 day weight**
```{r setup, include=FALSE, results="hide"}
knitr::opts_chunk$set(echo = TRUE)
#Load all necessary packages
library(rjags)
library(coda)
library(dplyr)
library(lme4)
library(blme)
library(RLRsim)
library(randomForest)
#For model with populations and family#
setwd("~/Documents/Lab/summer18")
cat<-read.csv("caterpillars_cleaned.csv")
pop<-cat$Cat_Pop
plant<-cat$Plant
familynames<-cat$family
family<-as.numeric(familynames)
familytotal<-length(unique(family))
inter<-as.numeric(as.factor(paste(cat$Cat_Pop,cat$Plant,sep="_")))
a<-as.character(unique(pop))
a<-sort(a)
plant<-cat$Plant
b<-as.character(unique(plant))
b<-sort(b)
interplot<-paste(rep(a, each=6), b, sep="-")
interplotmel<-interplot[c(1:12,19:30)]
interplotcol<-interplot[c(13:18)]
interplotply<-interplot[c(31:36)]
interplotall<-c(interplotmel,interplotcol,interplotply)
interplotall
#For model with species and populations
spec_cat<-read.csv("caterpillars_cleaned_species.csv")
spec_pop<-spec_cat$Cat_Pop
spec_plant<-spec_cat$Plant
spec_species<-spec_cat$species
spec_population<-as.numeric(spec_pop)
spec_poptotal<-length(unique(spec_population))
spec_inter<-as.numeric(as.factor(paste(spec_cat$species,spec_cat$Plant,sep="_")))
a<-as.character(unique(spec_species))
spec_plant<-spec_cat$Plant
b<-as.character(unique(spec_plant))
b<-sort(b)
spec_interplot<-paste(rep(a, each=6), b, sep="-")
spec_interplot
```
*This hierarchical model accounts for interactions between plant and caterpillar populations. In the figures, the caterpillar population comes before the hyphen and the plant population comes after the hyphen. In the figure with all species, PLY has the highest weight, followed by colias, and then melissa. You can also see that there is a similar pattern for each population of caterpillars in that APLL has the highest weight.For the melissa only figure it is clear that there is nothing significantly different as all the 95% CI overlap. The Gelman–Rubin Diagnostic shows that convergence was good for the model by analyzing the difference between the chains. *
```{r, eval=TRUE}
#get 8 day weight data in list form
#set y and p (plant) values
y<- cat$day_8_weight_mg
#get total number of y and p and family
ytotal<-length(y)
ptotal<-length(unique(pop))
intertotal<-length(inter)
#Put info into list
dataList<-list(y=y, family=family, ytotal=ytotal, pop=as.numeric(pop), plant=as.numeric(plant),ptotal=ptotal, inter=inter, intertotal=intertotal)
#Make model string and send it as txt file to JAGS
modelString="
model {
for (i in 1:ytotal){
y[i]~dnorm(mn[i],tau)
mn[i]=betapop[pop[i]]+betaplant[plant[i]] + betainter[inter[i]] + betafamily[family[i]]
}
for (j in 1:ptotal){
betapop[j]~dnorm(u,t)
}
for (k in 1:(ptotal-1)){
betaplant[k]~dnorm(0,tpl)
}
betaplant[ptotal]= -1 * sum(betaplant[1:(ptotal-1)])
for (b in 1:(intertotal-1)){
betainter[b]~dnorm(0,tb)
}
betainter[intertotal]= -1 * sum(betainter[1:(intertotal-1)])
## BST has 7 families
for (c in 1:(7-1)){
betafamily[c]~dnorm(0,t1)
}
betafamily[7]<--1 * sum(betafamily[1:(7-1)])
## BWP only has 1
betafamily[8]<-0
## COL is 9 - 17
for(c in 9:(17-1)){
betafamily[c]~dnorm(0,t2)
}
betafamily[17]<--1 * sum(betafamily[9:(17-1)])
##HWR has 7 families
for(c in 18:(24-1)){
betafamily[c]~dnorm(0,t3)
}
betafamily[24]<--1 * sum(betafamily[18:(24-1)])
##JJT has 7 families
for(c in 25:(31-1)){
betafamily[c]~dnorm(0,t4)
}
betafamily[31]<--1 * sum(betafamily[25:(31-1)])
##PLY has 1 family
betafamily[32]<--0
u~dnorm(0,.00001)
t~dgamma(0.01,0.01)
tau~dgamma(0.01,0.01)
tpl~dgamma(0.01,0.01)
tb~dgamma(0.01,0.01)
t1~dgamma(0.01,0.01)
t2~dgamma(0.01,0.01)
t3~dgamma(0.01,0.01)
t4~dgamma(0.01,0.01)
## variables to monitor only
for (j in 1:ptotal){
for(k in 1:ptotal){
muPopPlant[k+(ptotal * (j-1))]=betapop[j] + betaplant[k]+ betainter[k+(ptotal * (j-1))]
}
}
}
"
writeLines(modelString, con="TEMPmodel.txt")
##Generate Chains
#Get info into JAGs and let it generate MCMC sample from posterior
jagsmodel<- jags.model(file="TEMPmodel.txt", data=dataList, n.chains=3, n.adapt=500)
#Run chains to acomplish burn in
update(jagsmodel, n.iter=3000)
#Generate MCMC samples that represent posterior dist
codasamples<-coda.samples(jagsmodel, variable.names=c("betapop","betaplant","muPopPlant"), n.iter =3334)
summary(codasamples)
weight8col<-summary(codasamples)[[2]][25:30,]
weight8ply<-summary(codasamples)[[2]][43:48,]
weight8mel<-summary(codasamples)[[2]][c(13:24,31:42),]
weight8all<-rbind(weight8mel,weight8col,weight8ply)
#Plot all
plot(1:36,weight8all[,3],ylim=c(0,15),pch=21, cex.axis=0.75, xaxt = "n", main="8 Day Weight Interactions", xlab="Cat-Plant Populations", ylab="8 Day Weight (mg)")
axis(1, at=1:36, labels=interplotall[1:36], las=2, cex.axis=.48)
segments(1:36,weight8all[,1],1:36,weight8all[,5])
#Plot Mel
plot(1:24,weight8mel[,3],ylim=c(0,5),pch=21, cex.axis=0.75, xaxt = "n", main="8 Day Weight Interactions Melissa", xlab="Cat-Plant Populations", ylab="8 Day Weight (mg)")
axis(1, at=1:25, labels=interplotmel[1:25], las=2, cex.axis=.48)
segments(1:24,weight8mel[,1],1:24,weight8mel[,5])
gelman.diag(codasamples, confidence=0.95, transform=FALSE,multivariate = FALSE)
#Heatmap
matrixsamples8<-matrix(NA,nrow=3334*3,ncol=36)
for(i in 1:36){
matrixsamples8[,i]<-c(codasamples[[1]][,i+12],codasamples[[2]][,i+12],codasamples[[3]][,i+12])
}
dim(matrixsamples8)
heatmap_8<-matrix(NA,nrow=36,ncol=36)
for(i in 1:36){for(j in 1:36){
if(i != j){
heatmap_8[i,j]<-mean(matrixsamples8[,i] > matrixsamples8[,j])
}
}}
rownames(heatmap_8)<- interplotall
colnames(heatmap_8)<- interplotall
heatmap.2(heatmap_8, Rowv=NA, Colv=NA, main = "Weight 8", trace="none", cexRow=.45, dendrogram="none", density.info="none",cexCol =.45, legend())
```
**8 day weight- 3 species and Mel pop nested**
```{r, eval=TRUE}
#setup
#For model with species and populations
spec_cat<-read.csv("caterpillars_cleaned_species.csv")
spec_pop<-spec_cat$Cat_Pop
spec_plant<-spec_cat$Plant
spec_species<-spec_cat$species
spec_population<-as.numeric(spec_pop)
spec_poptotal<-length(unique(spec_population))
spec_inter<-as.numeric(as.factor(paste(spec_cat$species,spec_cat$Plant,sep="_")))
a<-as.character(unique(spec_species))
spec_plant<-spec_cat$Plant
b<-as.character(unique(spec_plant))
b<-sort(b)
spec_interplot<-paste(rep(a, each=6), b, sep="-")
spec_interplot
#get 8 day weight data in list form
#set y and p (plant) values
y<- spec_cat$day_8_weight_mg
#get total number of y and p and family
ytotal<-length(y)
spectotal<-length(unique(spec_species))
planttotal<-length(unique(spec_plant))
poptotal<-length(unique(spec_pop))
spec_intertotal<-length(spec_inter)
#Put info into list
dataList<-list(y=y, ytotal=ytotal, species=as.numeric(spec_species),spectotal=spectotal, plant=as.numeric(spec_plant),planttotal=planttotal, poptotal=poptotal, spec_pop=spec_pop, spec_inter=spec_inter, spec_intertotal=spec_intertotal)
dataList
pop
#Make model string and send it as txt file to JAGS
modelString="
#plant total
model {
for (i in 1:ytotal){
y[i]~dnorm(mn[i],tau)
mn[i]=betaspec[species[i]]+betaplant[plant[i]] + betainter[spec_inter[i]]+ betapop[spec_pop[i]]
}
for (j in 1:spectotal){
betaspec[j]~dnorm(u,t)
}
for (k in 1:(planttotal-1)){
betaplant[k]~dnorm(0,tpl)
}
betaplant[planttotal]= -1 * sum(betaplant[1:(planttotal-1)])
for (b in 1:(spec_intertotal-1)){
betainter[b]~dnorm(0,tb)
}
betainter[spec_intertotal]= -1 * sum(betainter[1:(spec_intertotal-1)])
## BST has 4 populations
for (c in 1:(4-1)){
betapop[c]~dnorm(0,t1)
}
betapop[4]<--1 * sum(betapop[1:(4-1)])
## COL only has 1
betapop[5]<-0
## PLY only has 1
betapop[6]<-0
u~dnorm(0,.00001)
t~dgamma(0.01,0.01)
tau~dgamma(0.01,0.01)
tpl~dgamma(0.01,0.01)
tb~dgamma(0.01,0.01)
t1~dgamma(0.01,0.01)
## variables to monitor only
for (j in 1:spectotal){
for(k in 1:planttotal){
muSpecPlant[k+(planttotal * (j-1))]=betaspec[j] + betaplant[k]+ betainter[k+(poptotal * (j-1))]
}
}
}
"
writeLines(modelString, con="TEMPmodel.txt")
##Generate Chains
#Get info into JAGs and let it generate MCMC sample from posterior
jagsmodel<- jags.model(file="TEMPmodel.txt", data=dataList, n.chains=3, n.adapt=500)
#Run chains to acomplish burn in
update(jagsmodel, n.iter=3000)
#Generate MCMC samples that represent posterior dist
#codasamples<-coda.samples(jagsmodel, variable.names=c("betaspec","betaplant","muSpecPlant","betapop"), n.iter =3334)
codasamples_spec_plant<-coda.samples(jagsmodel, variable.names=c("muSpecPlant"), n.iter =3334)
weight8col<-summary(codasamples_spec_plant)[[2]][1:6,]
weight8col
weight8ply<-summary(codasamples_spec_plant)[[2]][13:18,]
weight8mel<-summary(codasamples_spec_plant)[[2]][7:12,]
weight8all<-rbind(weight8mel,weight8col,weight8ply)
plot(1:18,weight8all[,3],ylim=c(0,25),pch=21, cex.axis=0.75, xaxt = "n", main="8 Day Weight Interactions", xlab="Species-Plant Populations", ylab="8 Day Weight (mg)")
axis(1, at=1:18, labels=spec_interplot[1:18], las=2, cex.axis=.48)
segments(1:18,weight8all[,1],1:18,weight8all[,5])
gelman.diag(codasamples_spec_plant, confidence=0.95, transform=FALSE,multivariate = FALSE)
#Heatmap
matrixsamples8<-matrix(NA,nrow=3334*3,ncol=18)
for(i in 1:18){
matrixsamples8[,i]<-c(codasamples_spec_plant[[1]][,i+12],codasamples[[2]][,i+12],codasamples[[3]][,i+12])
}
dim(matrixsamples8)
heatmap_8<-matrix(NA,nrow=36,ncol=36)
for(i in 1:36){for(j in 1:36){
if(i != j){
heatmap_8[i,j]<-mean(matrixsamples8[,i] > matrixsamples8[,j])
}
}}
rownames(heatmap_8)<- interplotall
colnames(heatmap_8)<- interplotall
heatmap.2(heatmap_8, Rowv=NA, Colv=NA, main = "Weight 8", trace="none", cexRow=.45, dendrogram="none", density.info="none",cexCol =.45, legend())
```
**14 Day Weight**
*This hierarchical model accounts for interactions between plant and caterpillar populations. In the figures, the caterpillar population comes before the hyphen and the plant population comes after the hyphen. In the figure with all species, PLY has the highest weight, followed by colias, and then melissa. You can also see that there is a similar pattern for each population of caterpillars in that APLL has the highest weight.For the melissa only figure it is clear that there is nothing significantly different as all the 95% CI overlap. The Gelman–Rubin Diagnostic shows that convergence was good for the model by analyzing the difference between the chains. *
```{r, eval=TRUE}
##14 day weight**
#get 14 day weight data in list form
#set y and p (plant) values
y<- cat$day_14_weight_mg
#get total number of y and p and family
ytotal<-length(y)
ptotal<-length(unique(pop))
intertotal<-length(inter)
familytotal<-length(unique(family))
#Put info into list
dataList<-list(y=y, family=family, ytotal=ytotal, pop=as.numeric(pop), plant=as.numeric(plant),ptotal=ptotal, inter=inter, intertotal=intertotal)
#Make model string and send it as txt file to JAGS
modelString="
model {
for (i in 1:ytotal){
y[i]~dnorm(mn[i],tau)
mn[i]=betapop[pop[i]]+betaplant[plant[i]] + betainter[inter[i]] + betafamily[family[i]]
}
for (j in 1:ptotal){
betapop[j]~dnorm(u,t)
}
for (k in 1:(ptotal-1)){
betaplant[k]~dnorm(0,tpl)
}
betaplant[ptotal]= -1 * sum(betaplant[1:(ptotal-1)])
for (b in 1:(intertotal-1)){
betainter[b]~dnorm(0,tb)
}
betainter[intertotal]= -1 * sum(betainter[1:(intertotal-1)])
## BST has 7 families
for (c in 1:(7-1)){
betafamily[c]~dnorm(0,t1)
}
betafamily[7]<--1 * sum(betafamily[1:(7-1)])
## BWP only has 1
betafamily[8]<-0
## COL is 9 - 17
for(c in 9:(17-1)){
betafamily[c]~dnorm(0,t2)
}
betafamily[17]<--1 * sum(betafamily[9:(17-1)])
##HWR has 7 families
for(c in 18:(24-1)){
betafamily[c]~dnorm(0,t3)
}
betafamily[24]<--1 * sum(betafamily[18:(24-1)])
##JJT has 7 families
for(c in 25:(31-1)){
betafamily[c]~dnorm(0,t4)
}
betafamily[31]<--1 * sum(betafamily[25:(31-1)])
##PLY has 1 family
betafamily[32]<--0
u~dnorm(0,.00001)
t~dgamma(0.01,0.01)
tau~dgamma(0.01,0.01)
tpl~dgamma(0.01,0.01)
tb~dgamma(0.01,0.01)
t1~dgamma(0.01,0.01)
t2~dgamma(0.01,0.01)
t3~dgamma(0.01,0.01)
t4~dgamma(0.01,0.01)
## variables to monitor only
for (j in 1:ptotal){
for(k in 1:ptotal){
muPopPlant[k+(ptotal * (j-1))]=betapop[j] + betaplant[k]+ betainter[k+(ptotal * (j-1))]
}
}
}
"
writeLines(modelString, con="TEMPmodel.txt")
##Generate Chains
#Get info into JAGs and let it generate MCMC sample from posterior
jagsmodel<- jags.model(file="TEMPmodel.txt", data=dataList, n.chains=3, n.adapt=500)
#Run chains to acomplish burn in
update(jagsmodel, n.iter=3000)
#Generate MCMC samples that represent posterior dist
codasamples<-coda.samples(jagsmodel, variable.names=c("betapop","betaplant","muPopPlant"), n.iter =3334)
summary(codasamples)
weight14col<-summary(codasamples)[[2]][25:30,]
weight14ply<-summary(codasamples)[[2]][43:48,]
weight14mel<-summary(codasamples)[[2]][c(13:24,31:42),]
weight14all<-rbind(weight16mel,weight16col,weight16ply)
#Plot all
plot(1:36,weight14all[,3],ylim=c(0,110),pch=21, cex.axis=0.75, xaxt = "n", main="14 Day Weight Interactions", xlab="Cat-Plant Populations", ylab="16 Day Weight (mg)")
axis(1, at=1:36, labels=interplotall[1:36], las=2, cex.axis=.48)
segments(1:36,weight14all[,1],1:36,weight14all[,5])
#Plot Mel
plot(1:24,weight14mel[,3],ylim=c(0,30),pch=21, cex.axis=0.75, xaxt = "n", main="14 Day Weight Interactions Melissa", xlab="Cat-Plant Populations", ylab="16 Day Weight (mg)")
axis(1, at=1:25, labels=interplotmel[1:25], las=2, cex.axis=.48)
segments(1:24,weight14mel[,1],1:24,weight14mel[,5])
gelman.diag(codasamples, confidence=0.95, transform=FALSE,multivariate = FALSE)
#Heatmap
matrixsamples14<-matrix(NA,nrow=3334*3,ncol=36)
for(i in 1:36){
matrixsamples14[,i]<-c(codasamples[[1]][,i+12],codasamples[[2]][,i+12],codasamples[[3]][,i+12])
}
dim(matrixsamples14)
heatmap_14<-matrix(NA,nrow=36,ncol=36)
for(i in 1:36){for(j in 1:36){
if(i != j){
heatmap_14[i,j]<-mean(matrixsamples14[,i] > matrixsamples14[,j])
}
}}
rownames(heatmap_14)<- interplotall
colnames(heatmap_14)<- interplotall
heatmap.2(heatmap_14, Rowv=NA, Colv=NA, main = "Weight 14", trace="none", cexRow=.45, dendrogram="none", density.info="none",cexCol =.45, legend())
```
**Pupation**
*COL had the highest pupation followed by BST with everything else being extremely low. Again it looks like APLL had the highest pupation proportion*
```{r, eval=TRUE}
#set y and p (plant) values
y<- cat$Binary_Pupation
#get total number of y and p
ytotal<-length(y)
ptotal<-length(unique(pop))
intertotal<-length(inter)
familytotal<-length(unique(family))
#Put info into list
dataList<-list(y=y, family=family, ytotal=ytotal, pop=as.numeric(pop), plant=as.numeric(plant),ptotal=ptotal, inter=inter, intertotal=intertotal)
#Make model string and send it as txt file to JAGS
modelString="
model {
for (i in 1:ytotal){
y[i]~dbern(p[i])
logit(p[i])=betapop[pop[i]]+betaplant[plant[i]] + betainter[inter[i]]+betafamily[family[i]]
}
for (j in 1:ptotal){
betapop[j]~dnorm(u,t)
}
for (k in 1:(ptotal-1)){
betaplant[k]~dnorm(0,tpl)
}
betaplant[ptotal]= -1 * sum(betaplant[1:(ptotal-1)])
for (b in 1:(intertotal-1)){
betainter[b]~dnorm(0,tb)
}
betainter[intertotal]= -1 * sum(betainter[1:(intertotal-1)])
## BST has 7 families
for (c in 1:(7-1)){
betafamily[c]~dnorm(0,t1)
}
betafamily[7]<--1 * sum(betafamily[1:(7-1)])
## BWP only has 1
betafamily[8]<-0
## COL is 9 - 17
for(c in 9:(17-1)){
betafamily[c]~dnorm(0,t2)
}
betafamily[17]<--1 * sum(betafamily[9:(17-1)])
##HWR has 7 families
for(c in 18:(24-1)){
betafamily[c]~dnorm(0,t3)
}
betafamily[24]<--1 * sum(betafamily[18:(24-1)])
##JJT has 7 families
for(c in 25:(31-1)){
betafamily[c]~dnorm(0,t4)
}
betafamily[31]<--1 * sum(betafamily[25:(31-1)])
##PLY has 1 family
betafamily[32]<--0
u~dnorm(0,.00001)
t~dgamma(0.01,0.01)
tau~dgamma(0.01,0.01)
tpl~dgamma(0.01,0.01)
tb~dgamma(0.01,0.01)
t1~dgamma(0.01,0.01)
t2~dgamma(0.01,0.01)
t3~dgamma(0.01,0.01)
t4~dgamma(0.01,0.01)
## variables to monitor only
for (j in 1:ptotal){
for(k in 1:ptotal){
muPopPlant[k+(ptotal * (j-1))]=1/(1+exp(-(betapop[j] + betaplant[k]+ betainter[k+(ptotal * (j-1))])))
}
}
}
"
writeLines(modelString, con="TEMPmodel.txt")
##Generate Chains
#Get info into JAGs and let it generate MCMC sample from posterior
jagsmodel<- jags.model(file="TEMPmodel.txt", data=dataList, n.chains=3, n.adapt=500)
#Run chains to acomplish burn in
update(jagsmodel, n.iter=3000)
#Generate MCMC samples that represent posterior dist
codasamples<-coda.samples(jagsmodel, variable.names=c("betapop","betaplant","muPopPlant"), n.iter =3334)
codasamples_family<-coda.samples(jagsmodel, variable.names=c("betapop","betaplant","muPopPlant","betafamily"), n.iter =3334)
summary(codasamples)
pupcol<-summary(codasamples)[[2]][25:30,]
pupply<-summary(codasamples)[[2]][43:48,]
pupmel<-summary(codasamples)[[2]][c(13:24,31:42),]
pupall<-rbind(pupmel,pupcol,pupply)
plot(1:36,pupall[,3],ylim=c(0,1),pch=21, cex.axis=0.75, xaxt = "n", main="Pupation Interactions", xlab="Cat-Plant Populations", ylab="Pupation proportion")
axis(1, at=1:36, labels=interplotall[1:36], las=2, cex.axis=.48)
segments(1:36,pupall[,1],1:36,pupall[,5])
#Plot Mel
plot(1:24,pupmel[,3],ylim=c(0,0.4),pch=21, cex.axis=0.75, xaxt = "n", main="Pupation Interactions Melissa", xlab="Cat-Plant Populations", ylab="Pupation proportion")
axis(1, at=1:25, labels=interplotmel[1:25], las=2, cex.axis=.48)
segments(1:24,pupmel[,1],1:24,pupmel[,5])
gelman.diag(codasamples, confidence=0.95, transform=FALSE,multivariate = FALSE)
#Heatmap
matrixsamplespup<-matrix(NA,nrow=3334*3,ncol=36)
for(i in 1:36){
matrixsamplespup[,i]<-c(codasamples[[1]][,i+12],codasamples[[2]][,i+12],codasamples[[3]][,i+12])
}
dim(matrixsamplespup)
heatmap_pup<-matrix(NA,nrow=36,ncol=36)
for(i in 1:36){for(j in 1:36){
if(i != j){
heatmap_pup[i,j]<-mean(matrixsamplespup[,i] > matrixsamplespup[,j])
}
}}
rownames(heatmap_pup)<- interplotall
colnames(heatmap_pup)<- interplotall
heatmap.2(heatmap_pup, Rowv=NA, Colv=NA, main = "Pupation", trace="none", cexRow=.45, dendrogram="none", density.info="none",cexCol =.45, legend())
```
**Eclosion**
*Not all that many caterpillar made it to eclosion. Colias had the highest, with BST following. This time caterpillars on APLL were not clearly more sucessful*
```{r, eval=TRUE}
#get eclosion data in list form
#set y and p (plant) values
y<- cat$Binary_Eclosion
interplot<-paste(rep(a, each=6), b)
#get total number of y and p
ytotal<-length(y)
ptotal<-length(unique(pop))
intertotal<-length(inter)
familytotal<-length(unique(family))
#Put info into list
dataList<-list(y=y, family=family, ytotal=ytotal, pop=as.numeric(pop), plant=as.numeric(plant),ptotal=ptotal, inter=inter, intertotal=intertotal)
#Make model string and send it as txt file to JAGS
modelString="
model {
for (i in 1:ytotal){
y[i]~dbern(p[i])
logit(p[i])=betapop[pop[i]]+betaplant[plant[i]] + betainter[inter[i]]+betafamily[family[i]]
}
for (j in 1:ptotal){
betapop[j]~dnorm(u,t)
}
for (k in 1:(ptotal-1)){
betaplant[k]~dnorm(0,tpl)
}
betaplant[ptotal]= -1 * sum(betaplant[1:(ptotal-1)])
for (b in 1:(intertotal-1)){
betainter[b]~dnorm(0,tb)
}
betainter[intertotal]= -1 * sum(betainter[1:(intertotal-1)])
## BST has 7 families
for (c in 1:(7-1)){
betafamily[c]~dnorm(0,t1)
}
betafamily[7]<--1 * sum(betafamily[1:(7-1)])
## BWP only has 1
betafamily[8]<-0
## COL is 9 - 17
for(c in 9:(17-1)){
betafamily[c]~dnorm(0,t2)
}
betafamily[17]<--1 * sum(betafamily[9:(17-1)])
##HWR has 7 families
for(c in 18:(24-1)){
betafamily[c]~dnorm(0,t3)
}
betafamily[24]<--1 * sum(betafamily[18:(24-1)])
##JJT has 7 families
for(c in 25:(31-1)){
betafamily[c]~dnorm(0,t4)
}
betafamily[31]<--1 * sum(betafamily[25:(31-1)])
##PLY has 1 family
betafamily[32]<--0
u~dnorm(0,.00001)
t~dgamma(0.01,0.01)
tau~dgamma(0.01,0.01)
tpl~dgamma(0.01,0.01)
tb~dgamma(0.01,0.01)
t1~dgamma(0.01,0.01)
t2~dgamma(0.01,0.01)
t3~dgamma(0.01,0.01)
t4~dgamma(0.01,0.01)
## variables to monitor only
for (j in 1:ptotal){
for(k in 1:ptotal){
muPopPlant[k+(ptotal * (j-1))]=1/(1+exp(-(betapop[j] + betaplant[k]+ betainter[k+(ptotal * (j-1))])))
}
}
}
"
writeLines(modelString, con="TEMPmodel.txt")
##Generate Chains
#Get info into JAGs and let it generate MCMC sample from posterior
jagsmodel<- jags.model(file="TEMPmodel.txt", data=dataList, n.chains=3, n.adapt=500)
#Run chains to acomplish burn in
update(jagsmodel, n.iter=3000)
#Generate MCMC samples that represent posterior dist
codasamples<-coda.samples(jagsmodel, variable.names=c("betapop","betaplant","muPopPlant"), n.iter =3334)
codasamples_family<-coda.samples(jagsmodel, variable.names=c("betapop","betaplant","muPopPlant","betafamily"), n.iter =3334)
summary(codasamples)
eclocol<-summary(codasamples)[[2]][25:30,]
ecloply<-summary(codasamples)[[2]][43:48,]
eclomel<-summary(codasamples)[[2]][c(13:24,31:42),]
ecloall<-rbind(eclomel,eclocol,ecloply)
plot(1:36,ecloall[,3],ylim=c(0,.85),pch=21, cex.axis=0.75, xaxt = "n", main="Eclosion Interactions", xlab="Cat-Plant Populations", ylab="Eclosion proportion")
axis(1, at=1:36, labels=interplotall[1:36], las=2, cex.axis=.48)
segments(1:36,ecloall[,1], 1:36, ecloall[,5])
#Plot Mel
plot(1:24,eclomel[,3],ylim=c(0,.1),pch=21, cex.axis=0.75, xaxt = "n", main="Eclosion Interactions Melissa", xlab="Cat-Plant Populations", ylab="Eclosion proportion")
axis(1, at=1:25, labels=interplotmel[1:25], las=2, cex.axis=.48)
segments(1:24,eclomel[,1],1:24,eclomel[,5])
gelman.diag(codasamples, confidence=0.95, transform=FALSE,multivariate = FALSE)
#Heatmap mel
matrixsampleseclo<-matrix(NA,nrow=3334*3,ncol=36)
for(i in 1:36){
matrixsampleseclo[,i]<-c(codasamples[[1]][,i+12],codasamples[[2]][,i+12],codasamples[[3]][,i+12])
}
dim(matrixsampleseclo)
heatmap_eclo<-matrix(NA,nrow=36,ncol=36)
for(i in 1:36){for(j in 1:36){
if(i != j){
heatmap_eclo[i,j]<-mean(matrixsampleseclo[,i] > matrixsampleseclo[,j])
}
}}
rownames(heatmap_eclo)<- interplotall
colnames(heatmap_eclo)<- interplotall
heatmap.2(heatmap_eclo, Rowv=NA, Colv=NA, main = "Eclosion", trace="none", cexRow=.45, dendrogram="none", density.info="none",cexCol =.45, legend())
```
```{r, eval=TRUE, include=FALSE}
##Chem
setwd("~/Documents/Lab/summer18")
chem<- read.csv("plant_chemistry.csv", header=T)
cat<- read.csv("caterpillars_cleaned.csv", header=T)
trait<-read.csv("all_plant_traits.csv", header=T)
##Means of chemicals and plant traits by pop
avgchem<-aggregate(chem[, 8:169], list(chem$population), mean)
colnames(avgchem)[colnames(avgchem)=="Group.1"] <- "Plant"
plant_trait<-trait[complete.cases(plant_trait), ]
avgplant<-aggregate(plant_trait[,-1], list(plant_trait$pop), mean)
colnames(avgplant)[colnames(avgplant)=="Group.1"] <- "Plant"
avgplant<-avgplant[,-2]
#install.packages("dplyr")
##Rename to match column
colnames(avgchem)[colnames(avgchem)=="Group.1"] <- "Plant"
##Combine avg chem with caterpillar data
cat_chem_all<-left_join(cat,avgchem)
cat_chem_plant_all<-left_join(cat_chem_all,avgplant)
mel_chem<-subset(cat_chem_plant_all, Cat_Pop=="JJT"|Cat_Pop=="BST"|Cat_Pop=="BWP"|Cat_Pop=="HWR")
mel_chem_plant<-mel_chem[complete.cases(mel_chem), ]
reml<-cat_chem_plant_all[,-c(2,4,6:9, 21:25)]
reml_cat_pop_fam<-reml[,c(2:3,6,8,11,13)]
##For random forest plant traits
plantchem<-chem[,-c(1:3,5:7)]
dput(names(plantchem))
colnames(plantchem)
```
**CHEM PCA**
```{r, eval=TRUE, include=FALSE}
setwd("~/Documents/Lab/summer18")
chem<- read.csv("plant_chemistry.csv", header=T)
plant_trait<-read.csv("all_plant_traits.csv", header=T)
library(RColorBrewer)
#scaled
cs<-brewer.pal(n=6,"Set1")
chem.pca <- prcomp(chem[,c(4:165)],scale. = TRUE)
plot(chem.pca$x[,1],chem.pca$x[,2],pch=20,col=cs[as.numeric(chem$population)], xlab="PC1", ylab="PC2")
legend("topleft",legend = c("Population", "Family","c","c","s","f"), pch = 15, col=cs[as.numeric(chem$population)])
#Not Scaled
chem.pca <- prcomp(chem[,c(8:169)],scale. = FALSE)
summary(chem.pca)
```
**REML for Plant Chemistry**
*Finished with plot for %pop*
```{r, eval=TRUE, results="hide", warning=FALSE}
setwd("~/Documents/Lab/summer18")
chem<- read.csv("plant_chemistry.csv", header=T)
x<-colnames(chem)
chemnames<-x[4:165]
#install.packages("lme4")
library(lme4)
#install.packages("RLRsim")
library("RLRsim")
#install.packages("blme")
library(blme)
chem_reml<-chem
chem_reml[,1:165] <- sapply(chem_reml[,1:165],as.numeric)
dim(chem_reml)
head(chem_reml)
Nx<-dim(chem_reml)[2]
Nx
#number of populations
L<-length(unique(chem_reml[,1])) ## no. lines
L
#names of populations
ul<-unique(chem_reml[,1])
ul
#Make blank matrix to fill in
vars_chem<-matrix(NA,nrow=Nx-3,ncol=4)
#Chem nested
for(i in 4:Nx){
nas<-which(is.na(chem_reml[,i])==FALSE)
#go through columns name y becomes column
dat<-data.frame(y=chem_reml[nas,i],pop=as.numeric(as.factor(chem_reml[nas,1])))
#Go through columns in lmer
est<-lmer(y ~ 1 + (1 | pop),data=dat,verbose=1,REML=TRUE)
vc<-as.data.frame(VarCorr(est))
vars_chem[i-3,1]<-vc$vcov[1]/sum(vc$vcov)
out<-exactRLRT(est)
#Output p value
vars_chem[i-3,3]<-out$p.value
#output statistic
vars_chem[i-3,2]<-out$statistic
est<-blmer(y ~ 1 + (1 | pop),data=dat,REML=TRUE)
vc<-as.data.frame(VarCorr(est))
vars_chem[i-3,4]<-vc$vcov[1]/sum(vc$vcov)
}
chemrem<-round(vars_chem,3)
colnames(chemrem) <- c("%pop","STAT","p-value","Bayes%fam")
chemrem
plot(1:162,chemrem[,1],ylim=c(0,.15),pch=21, cex.axis=0.75, xaxt = "n", main="REML for Populations of Plant", xlab="Chemical", ylab="%Pop")
axis(1, at=1:162, labels=chemnames, las=2, cex.axis=.4)
```
*Plant Trait REML*
```{r, eval=TRUE, results="hide", warning=FALSE}
trait<-read.csv("all_plant_traits.csv", header=T)
library(lme4)
library(blme)
library(RLRsim)
trait_reml<-trait
trait_names<-colnames(trait_reml)
trait_names<-trait_names[4:14]
trait_reml[,1:14] <- sapply(trait_reml[,1:14],as.numeric)
dim(trait_reml)
head(trait_reml)
Nx<-dim(trait_reml)[2]
Nx
#number of populations
L<-length(unique(trait_reml[,3])) ## no. lines
L
#names of populations
ul<-unique(trait_reml[,3])
ul
#Make blank matrix to fill in
vars_trait<-matrix(NA,nrow=Nx-3,ncol=8)
for(i in 4:Nx){
nas<-which(is.na(trait_reml[,i])==FALSE)
dat<-data.frame(y=trait_reml[nas,i],pop=as.numeric(as.factor(trait_reml[nas,3])),fam=as.numeric(as.factor(trait_reml[nas,2])))
est1<-lmer(y ~ 1 + (1|pop)+(1|pop:fam) ,data=dat,verbose=1,REML=TRUE)
estX<-lmer(y ~ 1 + (1|pop:fam) ,data=dat,verbose=1,REML=TRUE)
est0<-lmer(y ~ 1 + (1 | pop) ,data=dat,verbose=1,REML=TRUE)
vc<-as.data.frame(VarCorr(est1))
vars_trait[i-3,1]<-vc$vcov[1]/sum(vc$vcov)
vars_trait[i-3,2]<-vc$vcov[2]/sum(vc$vcov)
out<-exactRLRT(est0)
vars_trait[i-3,4]<-out$p.value
vars_trait[i-3,3]<-out$statistic
out<-exactRLRT(m=estX,mA=est1,m0=est0)
vars_trait[i-3,6]<-out$p.value
vars_trait[i-3,5]<-out$statistic
}
par(xpd=FALSE)
barplot(t(vars_trait[,1:2]),horiz=FALSE,beside =FALSE, main="REML Plant Traits", ylab="%Effect",col=c("black","grey"))
axis(1, at=1:11, labels=trait_names,las=2, cex.axis=.5, xpd = NA)
legend("topright",legend = c("Population", "Family"), pch = 15, col=c("black","grey"))
```
**Plant Population from Plant Traits**
*Only done for chemicals- would have to use the average for plant physically measured traits. It looks like plant chemical concentrations are not good predictor of plant populations (really, really high OOB)*
```{r, eval=TRUE}
library(randomForest)
#Run RF
plant_chem_forrf<-chem[,-c(2:3)]
plantpop_rf<-randomForest(population~., data=plant_chem_forrf, importance=TRUE, proximity=TRUE, ntree=300,mtry=7)
plantpop_rf
importance(plantpop_rf)
varImpPlot(plantpop_rf)
```