Post date: Apr 27, 2017 4:19:22 PM
We wanted to examine stripe frequency over time and the stability of differences in stripe frequency between A and C for the full 25 year time series. This is complicated by the fact that the sampling of sites over time is quite sparse; many sites were sampled once or a few times and several were sampled much more (but still not every year). Thus, simple annual averages confound spatial and temporal variation. To get around this I used a hierarchical Bayesian model that included site means and annual offsets for each host. The JAGS model and key R code from plotFreqs.R is below. The results, code, etc. are in king:/uufs/chpc.utah.edu/common/home/u6000989/projects/timema_fluct/timemaFluctMorph/.
JAGS model:
model{
## binomial sampling dist
for (i in 1:NoA){
yA[i] ~ dbinom(pA[i],nA[i])
logit(pA[i])<-beta0A[sitesA[i]] + beta1A[yrA[i]]
}
for (i in 1:NoC){
yC[i] ~ dbinom(pC[i],nC[i])
logit(pC[i])<-beta0C[sitesC[i]] + beta1C[yrC[i]]
}
## hierarchical priors on site and year and error effects
for(j in 1:NsiteA){
beta0A[j] ~ dnorm(muA,tA[1])
}
for(j in 1:NsiteC){
beta0C[j] ~ dnorm(muC,tC[1])
}
for(j in 1:(Nyr-1)){ ## sum2zero
beta1A[j] ~ dnorm(0,tA[2])
beta1C[j] ~ dnorm(0,tC[2])
}
beta1A[Nyr]<--1 * sum(beta1A[1:(Nyr-1)])
beta1C[Nyr]<--1 * sum(beta1C[1:(Nyr-1)])
#for(i in 1:(N-1)){
# e[i] ~ dnorm(0,t[3])
#}
#e[N]<- -1 * sum(e[1:(N-1)])
## hyperpriors
muA ~ dnorm(0,1e-6)
muC ~ dnorm(0,1e-6)
for(k in 1:2){
tA[k] ~ dgamma(0.01,0.001)
tC[k] ~ dgamma(0.01,0.001)
}
## derived quantities
for(j in 1:Nyr){
dif[j]<-beta1A[j]-beta1C[j]
yrMnA[j]<-1/(1 + exp(-1 * (muA+beta1A[j])))
yrMnC[j]<-1/(1 + exp(-1 * (muC+beta1C[j])))
yrDif[j]<-yrMnA[j]-yrMnC[j]
}
}
Key R code:
dat<-read.csv("Tcris_master_26.csv")
mt0<-which(dat$refugio_yn==0)
dat0<-dat[mt0,]
cs<-c("orange","blue")
A<-which(dat0$host=='A')
C<-which(dat0$host=='C')
Amns<-tapply(X=dat0$percent_striped_no_mel[A],INDEX=round(dat0$elevation_Clarissa_3,-2)[A],mean,na.rm=TRUE)
Cmns<-tapply(X=dat0$percent_striped_no_mel[C],INDEX=round(dat0$elevation_Clarissa_3,-2)[C],mean,na.rm=TRUE)
ulocs<-as.character(unique(dat0$location))
keep<-which(table(as.character(dat0$location)) > 5)
x<-which(dat0$location %in% names(keep))
A<-x[which(dat0$host[x]=='A')]
C<-x[which(dat0$host[x]=='C')]
Ayr<-tapply(X=dat0$percent_striped_no_mel[A],INDEX=dat0$year[A],mean,na.rm=TRUE)
Cyr<-tapply(X=dat0$percent_striped_no_mel[C],INDEX=dat0$year[C],mean,na.rm=TRUE)
x<-which(dat0$location %in% names(keep))
y<-dat0$striped[x]
n<-dat0$striped[x]+dat0$unstriped[x]
yr<-as.numeric(as.factor(dat0$year[x]))
A<-which(dat0$host[x]=='A')
C<-which(dat0$host[x]=='C')
yA<-y[A]
yC<-y[C]
nA<-n[A]
nC<-n[C]
yrA<-yr[A]
yrC<-yr[C]
sitesA<-as.numeric(as.factor(as.character(dat0$location[x][A])))
sitesC<-as.numeric(as.factor(as.character(dat0$location[x][C])))
library(rjags)
D<-list(yA=yA,yC=yC,nA=nA,nC=nC,yrA=yrA,yrC=yrC,sitesA=sitesA,sitesC=sitesC,NoA=length(A),NoC=length(C),Nyr=length(unique(yr)),NsiteA=length(unique(sitesA)),NsiteC=length(unique(sitesC)))
mod<-jags.model(file="stripemod-host",data=D,n.chains=3)
update(mod,n.iter=10000)
sam<-coda.samples(mod,n.iter=20000,thin=5,variable.names=c("yrMnA","yrMnC","yrDif","dif"))
pdf("yearMnDifPlotsv2.pdf",width=8,height=11)
par(mar=c(4,5,2,1))
par(mfrow=c(2,1))
plot(o[61:80,3],type='l',col="blue",ylim=c(0,1),cex.lab=1.5,axes=FALSE,xlab="year",ylab="stripe frequency")
axis(1,at=1:20,sort(unique(dat0$year)))
axis(2)
box()
polygon(x=c(1:20,rev(1:20)),y=c(o[61:80,1],rev(o[61:80,5])),col="lightblue",density=80)
lines(o[61:80,3],col="blue",lwd=2)
polygon(x=c(1:20,rev(1:20)),y=c(o[41:60,1],rev(o[41:60,5])),col="orange",density=80)
lines(o[41:60,3],col="orangered",lwd=2)
mn<-mean(o[21:40,3])
plot(o[21:40,3]-mn,type='l',col="green",ylim=c(-0.5,0.5),cex.lab=1.5,axes=FALSE,xlab="year",ylab="difference in stripe frequency")
axis(1,at=1:20,sort(unique(dat0$year)))
axis(2)
box()
polygon(x=c(1:20,rev(1:20)),y=c(o[21:40,1]-mn,rev(o[21:40,5]-mn)),col="green",density=80)
lines(o[21:40,3]-mn,col="forestgreen",lwd=2.5)
abline(h=0,lty=2,lwd=1.5)
dev.off()
And here is the main summary plot. The top panel is the annual mean stripe freq. for A (orange) and C (blue) on the natural scale (lines and 95% CIs). This is again in a hierarchical model that accounts for the uneven sampling. The second panel is then the difference in frequency over time with the mean difference removed. Note that most of the uncertainty comes from the sparse sampling.