Post date: Sep 06, 2016 4:13:26 PM
See files here: /home/zgompert/Local/timemaFluctMorph/
R code
## read in the phenotypic data
dat<-read.table("Tcris_master_locality_host_year_17.csv",header=TRUE,sep=",")
dat<-dat[-513,]
## format data
y<-dat$melanistic
n<-dat$total
elev<-dat$elevation_Clarissa
host<-as.numeric(dat$host)-1 ## A = 0, C = 1
#pop<-as.numeric(dat$location)
pop<-as.numeric(as.factor(paste(dat$location,dat$host,sep="-"))) ## version with host as part of pop
yr<-dat$year
temp<-dat$mean_MarAprMay_temp
mntn<-dat$refugio_yn
df<-data.frame(y,n,temp,pop,yr,host,elev,mntn)
#miss<-which(is.na(apply(df,1,sum))==TRUE)
#df<-df[-miss,]
df$pop<-as.numeric(as.factor(df$pop)) ## renumber pops
df$temp<-(df$temp-mean(df$temp))/sd(df$temp) ## scale and center temp
df$yr<-df$yr-mean(df$yr) ## center yr
df$elev<-(df$elev-mean(df$elev))/sd(df$elev) ## scale and center elev
D<-c(as.list(df),list(N=dim(df)[1],Npop=length(unique(df$pop))))
D$host<-round(tapply(X=D$host,INDEX=D$pop,mean),0)
D$elev<-tapply(X=D$elev,INDEX=D$pop,mean)
D$mntn<-round(tapply(X=D$mntn,INDEX=D$pop,mean),0)
## basic model with temp, elevation and host
library(rjags)
mod<-jags.model(file="hiermod-simple",data=D,n.chains=3)
update(mod,n.iter=10000)
out<-coda.samples(model=mod,variable.names=c("alpha","beta","theta","epsilon","a","b","c","taue"),n.iter=25000,thin=10)
odic<-dic.samples(model=mod,n.iter=25000,thin=10,type="pD")
save(list=ls(),file="hmod-simple.rdat")
## as above, but includes affect of being on refugio or not
library(rjags)
mod<-jags.model(file="hiermod-refugio",data=D,n.chains=3)
update(mod,n.iter=10000)
out<-coda.samples(model=mod,variable.names=c("alpha","beta","theta","epsilon","a","b","c","d","taue"),n.iter=25000,thin=10)
odic<-dic.samples(model=mod,n.iter=25000,thin=10,type="pD")
save(list=ls(),file="hmod-refugio.rdat")
Results: