This is for plant and caterpillar affecting caterpillar trait of interest, but no interaction between plant and cat population.
##install.packages("rjags")
##install.packages("runjags")
#setwd and read table
setwd("~/Documents/Lab/summer18")
cat<-read.csv("caterpillars_cleaned.csv")
#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
pop<-cat$Cat_Pop
plant<-cat$Plant
#get total number of y and p
ytotal<-length(y)
ptotal<-length(unique(pop))
#Put info into list
dataList<-list(y=y, ytotal=ytotal, pop=as.numeric(pop), plant=as.numeric(plant),ptotal=ptotal)
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]]
}
for (j in 1:ptotal){
betapop[j]~dnorm(u,t)
}
for (j in 1:(ptotal-1)){
betaplant[j]~dnorm(0,tpl)
}
## sum-to-zero constraint
betaplant[ptotal]= -1 * sum(betaplant[1:(ptotal-1)])
u~dnorm(0,.00001)
t~dgamma(0.01,0.01)
tpl~dgamma(0.01,0.01)
## variables to monitor only
for (j in 1:ptotal){
ibetapop[j]=1/(1 + exp(-betapop[j]))
for(k in 1:ptotal){
ibetapoppl[k+(ptotal * (j-1))]=1/(1+exp(-(betapop[j] + betaplant[k])))
}
}
}
"
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("ibetapop","ibetapoppl"), n.iter =3334)
summary(codasamples)
##Pupation for plant and cat
#get pupation data in list form
#set y and p (plant) values
y<- cat$Binary_Pupation
pop<-cat$Cat_Pop
plant<-cat$Plant
#get total number of y and p
ytotal<-length(y)
ptotal<-length(unique(pop))
#Put info into list
dataList<-list(y=y, ytotal=ytotal, pop=as.numeric(pop), plant=as.numeric(plant),ptotal=ptotal)
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]]
}
for (j in 1:ptotal){
betapop[j]~dnorm(u,t)
}
for (j in 1:(ptotal-1)){
betaplant[j]~dnorm(0,tpl)
}
## sum-to-zero constraint
betaplant[ptotal]= -1 * sum(betaplant[1:(ptotal-1)])
u~dnorm(0,.00001)
t~dgamma(0.01,0.01)
tpl~dgamma(0.01,0.01)
## variables to monitor only
for (j in 1:ptotal){
ibetapop[j]=1/(1 + exp(-betapop[j]))
for(k in 1:ptotal){
ibetapoppl[k+(ptotal * (j-1))]=1/(1+exp(-(betapop[j] + betaplant[k])))
}
}
}
"
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("ibetapop","ibetapoppl"), n.iter =3334)
summary(codasamples)
####8 day weight for plant and cat
#get 8 day weight data in list form
#set y and p (plant) values
y<- cat$day_8_weight_mg
pop<-cat$Cat_Pop
plant<-cat$Plant
#get total number of y and p
ytotal<-length(y)
ptotal<-length(unique(pop))
#Put info into list
dataList<-list(y=y, ytotal=ytotal, pop=as.numeric(pop), plant=as.numeric(plant),ptotal=ptotal)
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]~dnorm(mn[i],tau)
mn[i]=betapop[pop[i]]+betaplant[plant[i]]
}
for (j in 1:ptotal){
betapop[j]~dnorm(u,t)
}
for (j in 1:(ptotal-1)){
betaplant[j]~dnorm(0,tpl)
}
## sum-to-zero constraint
betaplant[ptotal]= -1 * sum(betaplant[1:(ptotal-1)])
u~dnorm(0,.00001)
t~dgamma(0.01,0.01)
tau~dgamma(0.01,0.01)
tpl~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]
}
}
}
"
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)
####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
pop<-cat$Cat_Pop
plant<-cat$Plant
#get total number of y and p
ytotal<-length(y)
ptotal<-length(unique(pop))
#Put info into list
dataList<-list(y=y, ytotal=ytotal, pop=as.numeric(pop), plant=as.numeric(plant),ptotal=ptotal)
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]~dnorm(mn[i],tau)
mn[i]=betapop[pop[i]]+betaplant[plant[i]]
}
for (j in 1:ptotal){
betapop[j]~dnorm(u,t)
}
for (j in 1:(ptotal-1)){
betaplant[j]~dnorm(0,tpl)
}
## sum-to-zero constraint
betaplant[ptotal]= -1 * sum(betaplant[1:(ptotal-1)])
u~dnorm(0,.00001)
t~dgamma(0.01,0.01)
tau~dgamma(0.01,0.01)
tpl~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]
}
}
}
"
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("ibetapop","ibetapoppl"), n.iter =3334)
summary(codasamples)