Post date: Feb 25, 2014 4:6:11 PM
I used rjags to fit Bayesian regression models for the phenotypic data. I didn't transform any of the data. The models for the survival and weight data had an intercept, plant effect, population effect, and interaction, as well as a family random effects. I ran three chains with 5000 iteration burnins, 100000 iterations of sampling, and thinning intervals of 2. ESS and trace plots looked great.
Here are the key parameter estimates:
Bayesian regression for phenotypic data
survival, pop and interaction significant
median 95\% ETPI
pop -0.5577845 (-1.0005028,-0.1481830)
plant -0.1226343 (-0.5146778,0.2733833)
interaction 3.5886160 (2.9514036,4.2501347)
adult weight, plant significant, 6.4 gram increase on As vs Me
median 95\% ETPI
pop 0.05998483 (-0.4852501,0.6110942)
plant 6.43577523 (5.9512288,6.9168216)
interaction -0.01273623 (-0.6589591,0.6322003)
The R code I used which is from projects/lycaeides_hostplant/phenotypicAnalyses/analyzePheno.R is below.
library(rjags)
## read data
ph<-read.table("performance.txt",header=TRUE)
ph<-ph[ph$family!="?",]
## format survival data
## gla = 0, sla = 1; ms = 0, ac = 1
data<-list(y=as.numeric(ph$survivorship),pop=as.numeric(ph$pop)-1,plant=2-as.numeric(ph$plant),fam=as.numeric(as.character(ph$family)),n=dim(ph)[1],nfam=51)
## specify jags model
jagsModel="
model {
for (i in 1:n){
y[i] ~ dbern(p[i])
logit(p[i])<-beta[1] + beta[2]*pop[i] + beta[3]*plant[i] + beta[4]*pop[i]*plant[i] + alpha0[fam[i]] * (1 - pop[i]) + alpha1[fam[i]] * pop[i]
}
## priors for regression coefficients
for (k in 1:4) {
beta[k] ~ dnorm(0,0.001)
}
## random effect prior
for (j in 1:(nfam-1)){
alpha0[j] ~ dnorm(0,tau0)
alpha1[j] ~ dnorm(0,tau1)
}
alpha0[nfam]<--1*sum(alpha0[1:(nfam-1)])
alpha1[nfam]<--1*sum(alpha1[1:(nfam-1)])
## precision priors
tau0 ~ dgamma(1,0.01)
tau1 ~ dgamma(1,0.01)
## monitor expected values, and alpha sum
ex[1]<-1/(1 + exp(-1 * beta[1])) ## gla on ms
ex[2]<-1/(1 + exp(-1 * (beta[1] + beta[3]))) ## gla on ac
ex[3]<-1/(1 + exp(-1 * (beta[1] + beta[2]))) ## sla on ms
ex[4]<-1/(1 + exp(-1 * (beta[1] + beta[2] + beta[3] + beta[4]))) ## sla on ac
sa0<-sum(alpha0)
sa1<-sum(alpha1)
}
"
## fit model
model<-jags.model(textConnection(jagsModel),data=data,n.chains=3)
update(model,n.iter=5000)
mcmcSamples<-coda.samples(model=model,variable.names=c("beta","sa0","sa1","ex"),n.iter=10000,thin=2)
## summarize and simple plot median and 95% ETPI, population, plant, and pxp significant
out<-summary(mcmcSamples)
##plot(1:4,out$quantiles[5:8,3],ylim=c(0,1))
##segments(1:4,out$quantiles[5:8,1],1:4,out$quantiles[5:8,5])
## format weight data
## gla = 0, sla = 1; ms = 0, ac = 1
x<-which(as.numeric(ph$survivorship)==1) ## retain individuals that lived
data<-list(y=as.numeric(ph$weight[x]),pop=as.numeric(ph$pop)[x]-1,plant=2-as.numeric(ph$plant)[x],fam=as.numeric(as.character(ph$family))[x],n=length(x),nfam=51)
## specify jags model
jagsModel="
model {
for (i in 1:n){
y[i] ~ dnorm(mu[i],tau)
mu[i]<-beta[1] + beta[2]*pop[i] + beta[3]*plant[i] + beta[4]*pop[i]*plant[i] + alpha0[fam[i]] * (1 - pop[i]) + alpha1[fam[i]] * pop[i]
}
## priors for regression coefficients
for (k in 1:4) {
beta[k] ~ dnorm(0,0.001)
}
## random effect prior
for (j in 1:(nfam-1)){
alpha0[j] ~ dnorm(0,tau0)
alpha1[j] ~ dnorm(0,tau1)
}
alpha0[nfam]<--1*sum(alpha0[1:(nfam-1)])
alpha1[nfam]<--1*sum(alpha1[1:(nfam-1)])
## precision priors
tau0 ~ dgamma(1,0.01)
tau1 ~ dgamma(1,0.01)
tau ~ dgamma(1,0.01)
## monitor expected values, and alpha sum
ex[1]<-beta[1] ## gla on ms
ex[2]<-beta[1] + beta[3]## gla on ac
ex[3]<-beta[1] + beta[2] ## sla on ms
ex[4]<-beta[1] + beta[2] + beta[3] + beta[4] ## sla on ac
sa0<-sum(alpha0)
sa1<-sum(alpha1)
}
"
## fit model
model<-jags.model(textConnection(jagsModel),data=data,n.chains=3)
update(model,n.iter=5000)
mcmcSamplesWgt<-coda.samples(model=model,variable.names=c("beta","sa0","sa1","ex"),n.iter=10000,thin=2)
## summarize and simple plot median and 95% ETPI, population, plant, and pxp significant
outWgt<-summary(mcmcSamplesWgt)
plot(1:4,outWgt$quantiles[5:8,3],ylim=c(0,1))
segments(1:4,outWgt$quantiles[5:8,1],1:4,outWgt$quantiles[5:8,5])
postscript("experiment.eps",width=4,height=8)
par(mfrow=c(2,1))
par(pty='s')
par(mar=c(4.5,5,1,0.5))
# survival
plot(c(0.15,0.85),out$statistics[c(5:6),1],ylim=c(0,1),xlim=c(0,1),type='b',lty=2,lwd=1.5,pch=16,axes=FALSE,xlab="",ylab="survival probability",cex.lab=1.5,cex=1.1)
points(c(0.15,0.85),out$statistics[c(7:8),1],type='b',lty=2,lwd=1.5,pch=23,cex=1.1)
lb<-out$statistics[c(5:8),1]-out$statistics[c(5:8),2]
ub<-out$statistics[c(5:8),1]+out$statistics[c(5:8),2]
segments(c(0.15,0.85,0.15,0.85),lb,c(0.15,0.85,0.15,0.85),ub,lwd=2)
axis(1,at=c(0.15,0.85),c("Medicago","Astragalus"),cex.axis=1.5)
axis(2,cex.axis=1.3)
box()
legend(0.05,0.95,c("GLA","SLA"),pch=c(16,23),bty='n',cex=1.1)
# weight
plot(c(0.15,0.85),outWgt$statistics[c(5:6),1],ylim=c(2,9),xlim=c(0,1),type='b',lty=2,lwd=1.5,pch=16,axes=FALSE,xlab="plant treatment",ylab="adult weight (grams)",cex.lab=1.5,cex=1.1)
points(c(0.15,0.85),outWgt$statistics[c(7:8),1],type='b',lty=2,lwd=1.5,pch=23,cex=1.1)
lb<-outWgt$statistics[c(5:8),1]-outWgt$statistics[c(5:8),2]
ub<-outWgt$statistics[c(5:8),1]+outWgt$statistics[c(5:8),2]
segments(c(0.15,0.85,0.15,0.85),lb,c(0.15,0.85,0.15,0.85),ub,lwd=2)
axis(1,at=c(0.15,0.85),c("Medicago","Astragalus"),cex.axis=1.5)
axis(2,cex.axis=1.3)
box()
#legend(0.05,0.95,c("GLA","SLA"),pch=c(16,23),bty='n',cex=1.1)
dev.off()