##install.packages("rjags")
##install.packages("runjags")
#setwd and read table
setwd("~/Documents/Lab/summer18")
cat<-read.csv("caterpillars_cleaned.csv")
pop<-cat$Cat_Pop
plant<-cat$Plant
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="-")
#load rjags
library(rjags)
library(runjags)
##Eclosion for plant and cat
#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)
#Put info into list
dataList<-list(y=y, ytotal=ytotal, pop=as.numeric(pop), plant=as.numeric(plant),ptotal=ptotal, inter=inter, intertotal=intertotal)
dataList
#Make model string adn send it as txt file to JAGS, theta is prob of survival to eclosion
modelString="
model {
for (i in 1:ytotal){
y[i]~dbern(p[i])
logit(p[i])=betapop[pop[i]]+betaplant[plant[i]] + betainter[inter[i]]
}
for (j in 1:ptotal){
betapop[j]~dnorm(u,t)
}
for (k in 1:(ptotal-1)){
betaplant[k]~dnorm(0,tpl)
}
for (b in 1:(intertotal-1)){
betainter[b]~dnorm(0,tb)
}
## sum-to-zero constraint
betaplant[ptotal]= -1 * sum(betaplant[1:(ptotal-1)])
betainter[intertotal]= -1 * sum(betainter[1:(intertotal-1)])
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)
## 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=500)
#Generate MCMC samples that represent posterior dist
codasamples<-coda.samples(jagsmodel, variable.names=c("betapop","betaplant","muPopPlant"), n.iter =3334)
summary(codasamples)
eclosion<-summary(codasamples)[[2]][13:48,]
plot(1:36,eclosion[,3],ylim=c(0,.8),pch=21, cex.axis=0.75, xaxt = "n", main="Eclosion Interactions", xlab="Cat-Plant Populations", ylab="Eclosion proportion")
axis(1, at=1:36, labels=interplot[1:36], las=2, cex.axis=.48)
segments(1:36,eclosion[,1], 1:36, eclosion[,5])
##Pupation for plant and cat
#get pupation data in list form
#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)
#Put info into list
dataList<-list(y=y, ytotal=ytotal, pop=as.numeric(pop), plant=as.numeric(plant),ptotal=ptotal, inter=inter, intertotal=intertotal)
dataList
#Make model string adn send it as txt file to JAGS, theta is prob of survival to eclosion
modelString="
model {
for (i in 1:ytotal){
y[i]~dbern(p[i])
logit(p[i])=betapop[pop[i]]+betaplant[plant[i]] + betainter[inter[i]]
}
for (j in 1:ptotal){
betapop[j]~dnorm(u,t)
}
for (k in 1:(ptotal-1)){
betaplant[k]~dnorm(0,tpl)
}
for (b in 1:(intertotal-1)){
betainter[b]~dnorm(0,tb)
}
## sum-to-zero constraint
betaplant[ptotal]= -1 * sum(betaplant[1:(ptotal-1)])
betainter[intertotal]= -1 * sum(betainter[1:(intertotal-1)])
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)
## 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=500)
#Generate MCMC samples that represent posterior dist
codasamples<-coda.samples(jagsmodel, variable.names=c("betapop","betaplant","muPopPlant"), n.iter =3334)
summary(codasamples)
pupation<-summary(codasamples)[[2]][13:48,]
plot(1:36,pupation[,3],ylim=c(0,.9),pch=21, cex.axis=0.75, xaxt = "n", main="Pupation Interactions", xlab="Cat-Plant Populations", ylab="Pupation proportion")
axis(1, at=1:36, labels=interplot[1:36], las=2, cex.axis=.48)
segments(1:36,pupation[,1],1:36,pupation[,5])
####8 day weight for plant and cat and interaction
#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
ytotal<-length(y)
ptotal<-length(unique(pop))
intertotal<-length(inter)
#Put info into list
dataList<-list(y=y, ytotal=ytotal, pop=as.numeric(pop), plant=as.numeric(plant),ptotal=ptotal, inter=inter, intertotal=intertotal)
dataList
#Make model string adn 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]]
}
for (j in 1:ptotal){
betapop[j]~dnorm(u,t)
}
for (k in 1:(ptotal-1)){
betaplant[k]~dnorm(0,tpl)
}
for (b in 1:(intertotal-1)){
betainter[b]~dnorm(0,tb)
}
## sum-to-zero constraint
betaplant[ptotal]= -1 * sum(betaplant[1:(ptotal-1)])
betainter[intertotal]= -1 * sum(betainter[1:(intertotal-1)])
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)
## 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=500)
#Generate MCMC samples that represent posterior dist
codasamples<-coda.samples(jagsmodel, variable.names=c("betapop","betaplant","muPopPlant"), n.iter =3334)
summary(codasamples)
weight8<-summary(codasamples)[[2]][13:48,]
plot(1:36,weight8[,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=interplot[1:36], las=2, cex.axis=.48)
segments(1:36,weight8[,1],1:36,weight8[,5])
####14 day weight for plant and cat
#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
ytotal<-length(y)
ptotal<-length(unique(pop))
intertotal<-length(inter)
#Put info into list
dataList<-list(y=y, ytotal=ytotal, pop=as.numeric(pop), plant=as.numeric(plant),ptotal=ptotal, inter=inter, intertotal=intertotal)
dataList
#Make model string adn 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]]
}
for (j in 1:ptotal){
betapop[j]~dnorm(u,t)
}
for (k in 1:(ptotal-1)){
betaplant[k]~dnorm(0,tpl)
}
for (b in 1:(intertotal-1)){
betainter[b]~dnorm(0,tb)
}
## sum-to-zero constraint
betaplant[ptotal]= -1 * sum(betaplant[1:(ptotal-1)])
betainter[intertotal]= -1 * sum(betainter[1:(intertotal-1)])
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)
## 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=500)
#Generate MCMC samples that represent posterior dist
codasamples<-coda.samples(jagsmodel, variable.names=c("betapop","betaplant","muPopPlant"), n.iter =3334)
summary(codasamples)
weight16<-summary(codasamples)[[2]][13:48,]
plot(1:36,weight16[,3],ylim=c(0,100),pch=21, cex.axis=0.75, xaxt = "n", main="16 Day Weight Interactions", xlab="Cat-Plant Populations", ylab="16 Day Weight (mg) ")
axis(1, at=1:36, labels=interplot[1:36], las=2, cex.axis=.48)
segments(1:36,weight16[,1],1:36,weight16[,5])