Post date: Oct 04, 2016 10:27:39 PM
I started a final set of runs with elevation, climPC1, climPC2 or climPC1 and climPC2 as covariates. I am estimating model parameters, calculating DIC and conducting cross-validation (only including samples with n = 25+, 4 subsets for CV). This is running on my desktop computer and results will be in ~/Local/timemaFluctMorph/.
Here is the R code (fluctMorphAnalysis.R):
## read in the phenotypic data
dat<-read.table("Tcris_master_26.csv",header=TRUE,sep=",")
dat<-dat[-452,] ## pooled sample, no elevation data
###### format data #########
y<-dat$melanistic
n<-dat$total
elev<-dat$elevation_Clarissa_3
host<-as.numeric(dat$host)-1 ## A = 0, C = 1
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
## combine and clean-up data
df<-data.frame(y,n,temp,pop,yr,host,elev,mntn,pc1=dat$clim_PC1,pc2=dat$clim_PC2)
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)
#### fit models ####
### elevation ####
library(rjags)
mod<-jags.model(file="hiermod-refugio",data=D,n.chains=3)
update(mod,n.iter=10000)
outElev<-coda.samples(model=mod,variable.names=c("alpha","beta","theta","epsilon","a","b","c","d","taue"),n.iter=25000,thin=10)
odicElev<-dic.samples(model=mod,n.iter=25000,thin=10,type="pD")
save(list=ls(),file="postElev.rdat")
### PC1 ####
library(rjags)
mod<-jags.model(file="hiermod-pc1",data=D,n.chains=3)
update(mod,n.iter=10000)
outElev<-coda.samples(model=mod,variable.names=c("alpha","beta","theta","epsilon","a","b","c","d","taue"),n.iter=25000,thin=10)
odicElev<-dic.samples(model=mod,n.iter=25000,thin=10,type="pD")
save(list=ls(),file="postPc1.rdat")
### PC2 ###
library(rjags)
mod<-jags.model(file="hiermod-pc2",data=D,n.chains=3)
update(mod,n.iter=10000)
outElev<-coda.samples(model=mod,variable.names=c("alpha","beta","theta","epsilon","a","b","c","d","taue"),n.iter=25000,thin=10)
odicElev<-dic.samples(model=mod,n.iter=25000,thin=10,type="pD")
save(list=ls(),file="postPc2.rdat")
### PC1 & PC2 ###
library(rjags)
mod<-jags.model(file="hiermod-pc1p2",data=D,n.chains=3)
update(mod,n.iter=10000)
outElev<-coda.samples(model=mod,variable.names=c("alpha","beta","theta","epsilon","a","b","c","d","e","taue"),n.iter=25000,thin=10)
odicElev<-dic.samples(model=mod,n.iter=25000,thin=10,type="pD")
save(list=ls(),file="postPc1p2.rdat")
####################### cross validation #############################
library(rjags)
a<-which(D$n > 25) ## minimum sample of 25 for cv
v<-rep(1:4,each=54)
vv<-sample(v,length(v),replace=FALSE)
Y<-D$y
## elevation
out<-vector("list",4)
for(i in 1:4){
D$y<-Y
D$y[a[vv==i]]<-NA
mod<-jags.model(file="hiermod-refugio",data=D,n.chains=3)
update(mod,n.iter=10000)
out[[i]]<-coda.samples(model=mod,variable.names=c("y"),n.iter=25000,thin=10)
}
D$y<-Y
save(list=ls(),file="cvElev.rdat")
## pc1
out<-vector("list",4)
for(i in 1:4){
D$y<-Y
D$y[a[vv==i]]<-NA
mod<-jags.model(file="hiermod-pc1",data=D,n.chains=3)
update(mod,n.iter=10000)
out[[i]]<-coda.samples(model=mod,variable.names=c("y"),n.iter=25000,thin=10)
}
D$y<-Y
save(list=ls(),file="cvPc1.rdat")
## pc2
out<-vector("list",4)
for(i in 1:4){
D$y<-Y
D$y[a[vv==i]]<-NA
mod<-jags.model(file="hiermod-pc2",data=D,n.chains=3)
update(mod,n.iter=10000)
out[[i]]<-coda.samples(model=mod,variable.names=c("y"),n.iter=25000,thin=10)
}
D$y<-Y
save(list=ls(),file="cvPc2.rdat")
## pc2
out<-vector("list",4)
for(i in 1:4){
D$y<-Y
D$y[a[vv==i]]<-NA
mod<-jags.model(file="hiermod-pc1p2",data=D,n.chains=3)
update(mod,n.iter=10000)
out[[i]]<-coda.samples(model=mod,variable.names=c("y"),n.iter=25000,thin=10)
}
D$y<-Y
save(list=ls(),file="cvPc1p2.rdat")