Post date: Aug 05, 2017 7:20:18 PM
The data, code and figures are all in /uufs/chpc.utah.edu/common/home/u6000989/projects/timema_fluct/predictingEvolution/. Separate R files exist for some figures, but the main analysis is in predictingEvoAnalyses.R and uses functions from predictingEvoFuncs.R. Here is the full R code:
## forecasting and predicting evolutionary change across data sets
## load functions and libraries
library(rjags)
source("predictingEvoFuncs.R")
########################## timema stripe #############################3
dat<-read.csv("Tcris_master_26.csv")
hv<-which(dat$location=="HV")
y<-tapply(X=dat$striped[hv],INDEX=dat$year[hv],sum)
n<-tapply(X=dat$striped[hv]+dat$unstriped[hv],INDEX=dat$year[hv],sum)
## add 2017 data
y<-c(y,90+280)
n<-c(n,164+90+106+280)
## freq
p<-y[-c(1:2)]/n[-c(1:2)]
## dp
dp<-p[2:18]-p[1:17]
dat<-list(y=dp,N=17)
## determine best arma model based on DIC
dicTcrStripe1<-fitArmaModels(dat=dat,niter=100000,thin=50,nchains=10)
dicTcrStripe2<-fitArmaModels(dat=dat,niter=100000,thin=50,nchains=10)
## best model was 6 = ARMA(1,2)
cvTcrStripe1<-crossValidation(dat=dat,niter=100000,thin=50,burn=50000,nchains=10,mymod="dp-arma12model",init.theta=1,init.phi=2)
cvTcrStripe2<-crossValidation(dat=dat,niter=100000,thin=50,burn=50000,nchains=10,mymod="dp-arma12model",init.theta=1,init.phi=2)
## forecasting
drop<-3:10
forecastTcrStripe1<-forecastChange(dat=dat,niter=100000,thin=50,burn=50000,nchains=10,mymod="dp-arma12model",init.theta=1,init.phi=2,drop=drop)
forecastTcrStripe2<-forecastChange(dat=dat,niter=100000,thin=50,burn=50000,nchains=10,mymod="dp-arma12model",init.theta=1,init.phi=2,drop=drop)
## make plots
yr<-2000:2017
pdf("p-stripeTimeSeriesRaw.pdf",width=6,height=5)
par(mar=c(5,5,0.5,0.5))
plot(yr,p,type='b',ylim=c(0,1),pch=19,cex.lab="1.5",xlab="year",ylab="stripe frequency")
dev.off()
yr<-2000:2016
pdf("dp-stripeTimeSeriesRaw.pdf",width=6,height=5)
par(mar=c(5,5,0.5,0.5))
plot(yr,dp,type='b',ylim=c(-0.33,0.27),pch=19,cex.lab="1.5",xlab="year",ylab="change in stripe frequency")
abline(h=0,lty=2,col="darkgray")
dev.off()
pred<-(cvTcrStripe1[[1]][,3]+cvTcrStripe2[[1]][,3])/2
o<-lm(dat$y[-c(1:2)]~pred[-c(1:2)])
pdf("dp-stripeObsVpredict.pdf",width=5,height=5)
par(mar=c(5,5,0.5,0.5))
plot(pred,dp,pch=19,xlab="predicted change",ylab="observed change",cex.lab=1.4)
abline(o$coefficients,col="darkgray",lwd=1.5)
dev.off()
# Estimate Std. Error t value Pr(>|t|)
#(Intercept) -0.005463 0.027032 -0.202 0.842976
#pred[-c(1:2)] 0.938310 0.171426 5.474 0.000107 ***
#Residual standard error: 0.1047 on 13 degrees of freedom
#Multiple R-squared: 0.6974, Adjusted R-squared: 0.6741
#F-statistic: 29.96 on 1 and 13 DF, p-value: 0.0001068
forecastTcrStripeComb<-(forecastTcrStripe1[[1]]+forecastTcrStripe2[[1]])/2
apply(forecastTcrStripeComb,2,median)
#[1] 0.8618326 0.9282905 0.5872508 0.9863070
# median r^2, r, lb r, ub r
pdf("dp-stripeNyrPredict.pdf",width=5,height=5)
par(mar=c(5,5,.5,.5))
plot(drop,forecastTcrStripeComb[,1],type='b',xlab="years predicted",ylab=expression(paste("predictive ",r^2)),cex.lab=1.4,ylim=c(0,1),pch=19)
dev.off()
pdf("dp-stripeNyrPredictCorCI.pdf",width=5,height=5)
par(mar=c(5,5,.5,.5))
plot(drop,forecastTcrStripeComb[,2],type='n',xlab="years predicted",ylab="predictive correlation",cex.lab=1.4,ylim=c(-1,1),pch=19)
polygon(c(drop[-1],rev(drop[-1])),c(forecastTcrStripeComb[-1,3],rev(forecastTcrStripeComb[-1,4])),col="gray",density=50)
points(drop,forecastTcrStripeComb[,2],type='b',pch=19)
abline(h=0,lty=2,col="darkblue")
dev.off()
save(list=ls(),file="predictTcrStr.rdat")
########################## timema color #############################3
dat<-read.csv("Tcris_master_26.csv")
hv<-which(dat$location=="HV")
y<-tapply(X=dat$melanistic[hv],INDEX=dat$year[hv],sum)
x<-tapply(X=dat$mean_FebMarApr_temp[hv],INDEX=dat$year[hv],mean)
n<-tapply(X=dat$total[hv],INDEX=dat$year[hv],sum)
## add 2017 data
y<-c(y,233+90)
n<-c(n,164+90+106+280+233+90)
## freq
p<-y[-c(1:2)]/n[-c(1:2)]
## dp
dp<-p[2:18]-p[1:17]
dat<-list(y=dp,N=17)
## determine best arma model based on DIC
dicTcrColor1<-fitArmaModels(dat=dat,niter=100000,thin=50,nchains=10)
dicTcrColor2<-fitArmaModels(dat=dat,niter=100000,thin=50,nchains=10)
## best model was 6 = ARMA(1,2)
cvTcrColor1<-crossValidation(dat=dat,niter=100000,thin=50,burn=50000,nchains=10,mymod="dp-arma12model",init.theta=1,init.phi=2)
cvTcrColor2<-crossValidation(dat=dat,niter=100000,thin=50,burn=50000,nchains=10,mymod="dp-arma12model",init.theta=1,init.phi=2)
## forecasting
drop<-3:10
forecastTcrColor1<-forecastChange(dat=dat,niter=100000,thin=50,burn=50000,nchains=10,mymod="dp-arma12model",init.theta=1,init.phi=2,drop=drop)
forecastTcrColor2<-forecastChange(dat=dat,niter=100000,thin=50,burn=50000,nchains=10,mymod="dp-arma12model",init.theta=1,init.phi=2,drop=drop)
## make plots
yr<-2000:2017
pdf("p-colorTimeSeriesRaw.pdf",width=6,height=5)
par(mar=c(5,5,0.5,0.5))
plot(yr,p,type='b',ylim=c(0,1),pch=19,cex.lab="1.5",xlab="year",ylab="color frequency")
dev.off()
yr<-2000:2016
pdf("dp-colorTimeSeriesRaw.pdf",width=6,height=5)
par(mar=c(5,5,0.5,0.5))
plot(yr,dp,type='b',ylim=c(-0.33,0.27),pch=19,cex.lab="1.5",xlab="year",ylab="change in color frequency")
abline(h=0,lty=2,col="darkgray")
dev.off()
pred<-(cvTcrColor1[[1]][,3]+cvTcrColor2[[1]][,3])/2
o<-lm(dat$y[-c(1:2)]~pred[-c(1:2)])
pdf("dp-colorObsVpredict.pdf",width=5,height=5)
par(mar=c(5,5,0.5,0.5))
plot(pred,dp,pch=19,xlab="predicted change",ylab="observed change",cex.lab=1.4)
abline(o$coefficients,col="darkgray",lwd=1.5)
dev.off()
# Estimate Std. Error t value Pr(>|t|)
#(Intercept) 0.03373 0.02652 1.272 0.226
#pred[-c(1:2)] -1.24295 1.02359 -1.214 0.246
#Residual standard error: 0.06903 on 13 degrees of freedom
#Multiple R-squared: 0.1019, Adjusted R-squared: 0.03279
#F-statistic: 1.475 on 1 and 13 DF, p-value: 0.2462
forecastTcrColorComb<-(forecastTcrColor1[[1]]+forecastTcrColor2[[1]])/2
apply(forecastTcrColorComb,2,median)
#[1] 0.1388707 -0.2959806 -0.7044798 0.5042099
# median r^2, r, lb r, ub r
pdf("dp-colorNyrPredict.pdf",width=5,height=5)
par(mar=c(5,5,.5,.5))
plot(drop,forecastTcrColorComb[,1],type='b',xlab="years predicted",ylab=expression(paste("predictive ",r^2)),cex.lab=1.4,ylim=c(0,1),pch=19)
dev.off()
pdf("dp-colorNyrPredictCorCI.pdf",width=5,height=5)
par(mar=c(5,5,.5,.5))
plot(drop,forecastTcrColorComb[,2],type='n',xlab="years predicted",ylab="predictive correlation",cex.lab=1.4,ylim=c(-1,1),pch=19)
polygon(c(drop[-1],rev(drop[-1])),c(forecastTcrColorComb[-1,3],rev(forecastTcrColorComb[-1,4])),col="gray",density=50)
points(drop,forecastTcrColorComb[,2],type='b',pch=19)
abline(h=0,lty=2,col="darkblue")
dev.off()
save(list=ls(),file="predictTcrCol.rdat")
######################### G. fortis, body size #########################
dat<-read.csv("fortis.csv")
y1<-dat$PC1.body[2:41] ## pc1 body size
y2<-dat$PC1.beak[2:41] ## pc1 beak size
y3<-dat$PC2.beak[2:41] ## pc2 beak size
p<-y1
## dp
dp<-p[2:40]-p[1:39]
dat<-list(y=dp,N=39)
## determine best arma model based on DIC
dicGfortBody1<-fitArmaModels(dat=dat,niter=100000,thin=50,nchains=10)
dicGfortBody2<-fitArmaModels(dat=dat,niter=100000,thin=50,nchains=10)
## best model was 1, ARMA(0,0); then 7, ARMA(2,2)
bestMod<-"dp-arma22model"
theta<-2
phi<-2
cvGfortBody1<-crossValidation(dat=dat,niter=100000,thin=50,burn=50000,nchains=10,mymod=bestMod,init.theta=theta,init.phi=phi)
cvGfortBody2<-crossValidation(dat=dat,niter=100000,thin=50,burn=50000,nchains=10,mymod=bestMod,init.theta=theta,init.phi=phi)
## forecasting
drop<-3:10
forecastGfortBody1<-forecastChange(dat=dat,niter=100000,thin=50,burn=50000,nchains=10,mymod=bestMod,init.theta=theta,init.phi=phi,drop=drop)
forecastGfortBody2<-forecastChange(dat=dat,niter=100000,thin=50,burn=50000,nchains=10,mymod=bestMod,init.theta=theta,init.phi=phi,drop=drop)
## make plots
yr<-1973:2012
pdf("p-GfortBodyTimeSeriesRaw.pdf",width=6,height=5)
par(mar=c(5,5,0.5,0.5))
plot(yr,p,type='b',pch=19,cex.lab="1.5",xlab="year",ylab="body size")
dev.off()
yr<-1973:2011
pdf("dp-GfortBodyTimeSeriesRaw.pdf",width=6,height=5)
par(mar=c(5,5,0.5,0.5))
plot(yr,dp,type='b',pch=19,cex.lab="1.5",xlab="year",ylab="change in body size")
abline(h=0,lty=2,col="darkgray")
dev.off()
pred<-(cvGfortBody1[[1]][,3]+cvGfortBody2[[1]][,3])/2
o<-lm(dat$y[-c(1:2)]~pred[-c(1:2)])
pdf("dp-GfortBodyObsVpredict.pdf",width=5,height=5)
par(mar=c(5,5,0.5,0.5))
plot(pred,dp,pch=19,xlab="predicted change",ylab="observed change",cex.lab=1.4)
abline(o$coefficients,col="darkgray",lwd=1.5)
dev.off()
# Estimate Std. Error t value Pr(>|t|)
#(Intercept) -0.07142 0.03189 -2.240 0.03156 *
#pred[-c(1:2)] -1.89589 0.54332 -3.489 0.00133 **
#Residual standard error: 0.1682 on 35 degrees of freedom
#Multiple R-squared: 0.2581, Adjusted R-squared: 0.2369
#F-statistic: 12.18 on 1 and 35 DF, p-value: 0.001327
forecastGfortBodyComb<-(forecastGfortBody1[[1]]+forecastGfortBody2[[1]])/2
apply(forecastGfortBodyComb,2,median)
# 0.2565920 0.2593157 -0.6401642 0.7617514
# median r^2, r, lb r, ub r
pdf("dp-GfortBodyNyrPredict.pdf",width=5,height=5)
par(mar=c(5,5,.5,.5))
plot(drop,forecastGfortBodyComb[,1],type='b',xlab="years predicted",ylab=expression(paste("predictive ",r^2)),cex.lab=1.4,ylim=c(0,1),pch=19)
dev.off()
pdf("dp-GfortBodyNyrPredictCorCI.pdf",width=5,height=5)
par(mar=c(5,5,.5,.5))
plot(drop,forecastGfortBodyComb[,2],type='n',xlab="years predicted",ylab="predictive correlation",cex.lab=1.4,ylim=c(-1,1),pch=19)
polygon(c(drop[-1],rev(drop[-1])),c(forecastGfortBodyComb[-1,3],rev(forecastGfortBodyComb[-1,4])),col="gray",density=50)
points(drop,forecastGfortBodyComb[,2],type='b',pch=19)
abline(h=0,lty=2,col="darkblue")
dev.off()
save(list=ls(),file="predictGfortBodySize.rdat")
################ G. fortis, beak size 1 #############################3
dat<-read.csv("fortis.csv")
y1<-dat$PC1.body[2:41] ## pc1 body size
y2<-dat$PC1.beak[2:41] ## pc1 beak size
y3<-dat$PC2.beak[2:41] ## pc2 beak size
p<-y2
## dp
dp<-p[2:40]-p[1:39]
dat<-list(y=dp,N=39)
## determine best arma model based on DIC
dicGfortBeakA1<-fitArmaModels(dat=dat,niter=100000,thin=50,nchains=10)
dicGfortBeakA2<-fitArmaModels(dat=dat,niter=100000,thin=50,nchains=10)
## best model was 1, ARMA(0,0); then 3, ARMA(0,1)
bestMod<-"dp-arma01model"
theta<-NA
phi<-1
cvGfortBeakA1<-crossValidation(dat=dat,niter=100000,thin=50,burn=50000,nchains=10,mymod=bestMod,init.theta=theta,init.phi=phi)
cvGfortBeakA2<-crossValidation(dat=dat,niter=100000,thin=50,burn=50000,nchains=10,mymod=bestMod,init.theta=theta,init.phi=phi)
## forecasting
drop<-3:10
forecastGfortBeakA1<-forecastChange(dat=dat,niter=100000,thin=50,burn=50000,nchains=10,mymod=bestMod,init.theta=theta,init.phi=phi,drop=drop)
forecastGfortBeakA2<-forecastChange(dat=dat,niter=100000,thin=50,burn=50000,nchains=10,mymod=bestMod,init.theta=theta,init.phi=phi,drop=drop)
## make plots
yr<-1973:2012
pdf("p-GfortBeak1TimeSeriesRaw.pdf",width=6,height=5)
par(mar=c(5,5,0.5,0.5))
plot(yr,p,type='b',pch=19,cex.lab="1.5",xlab="year",ylab="beak PC1")
dev.off()
yr<-1973:2011
pdf("dp-GfortBeak1TimeSeriesRaw.pdf",width=6,height=5)
par(mar=c(5,5,0.5,0.5))
plot(yr,dp,type='b',pch=19,cex.lab="1.5",xlab="year",ylab="change in beak PC1")
abline(h=0,lty=2,col="darkgray")
dev.off()
pred<-(cvGfortBeakA1[[1]][,3]+cvGfortBeakA2[[1]][,3])/2
o<-lm(dat$y[-c(1:2)]~pred[-c(1:2)])
pdf("dp-GfortBeak1ObsVpredict.pdf",width=5,height=5)
par(mar=c(5,5,0.5,0.5))
plot(pred,dp,pch=19,xlab="predicted change",ylab="observed change",cex.lab=1.4)
abline(o$coefficients,col="darkgray",lwd=1.5)
dev.off()
# Estimate Std. Error t value Pr(>|t|)
#(Intercept) -0.1872 0.0560 -3.344 0.001980 **
#pred[-c(1:2)] -6.0475 1.6519 -3.661 0.000822 ***
#Residual standard error: 0.1908 on 35 degrees of freedom
#Multiple R-squared: 0.2769, Adjusted R-squared: 0.2562
#F-statistic: 13.4 on 1 and 35 DF, p-value: 0.0008223
forecastGfortBeakAComb<-(forecastGfortBeakA1[[1]]+forecastGfortBeakA2[[1]])/2
apply(forecastGfortBeakAComb,2,median)
# 0.1405218 -0.0460291 -0.6096585 0.7513293
# median r^2, r, lb r, ub r
pdf("dp-GfortBeak1NyrPredict.pdf",width=5,height=5)
par(mar=c(5,5,.5,.5))
plot(drop,forecastGfortBeakAComb[,1],type='b',xlab="years predicted",ylab=expression(paste("predictive ",r^2)),cex.lab=1.4,ylim=c(0,1),pch=19)
dev.off()
pdf("dp-GfortBeak1NyrPredictCorCI.pdf",width=5,height=5)
par(mar=c(5,5,.5,.5))
plot(drop,forecastGfortBeakAComb[,2],type='n',xlab="years predicted",ylab="predictive correlation",cex.lab=1.4,ylim=c(-1,1),pch=19)
polygon(c(drop[-1],rev(drop[-1])),c(forecastGfortBeakAComb[-1,3],rev(forecastGfortBeakAComb[-1,4])),col="gray",density=50)
points(drop,forecastGfortBeakAComb[,2],type='b',pch=19)
abline(h=0,lty=2,col="darkblue")
dev.off()
save(list=ls(),file="predictGfortBeakSizeA.rdat")
######################## G. fortis, beak size 2 ##############################
dat<-read.csv("fortis.csv")
y1<-dat$PC1.body[2:41] ## pc1 body size
y2<-dat$PC1.beak[2:41] ## pc1 beak size
y3<-dat$PC2.beak[2:41] ## pc2 beak size
p<-y3
## dp
dp<-p[2:40]-p[1:39]
dat<-list(y=dp,N=39)
## determine best arma model based on DIC
dicGfortBeakB1<-fitArmaModels(dat=dat,niter=100000,thin=50,nchains=10)
dicGfortBeakB2<-fitArmaModels(dat=dat,niter=100000,thin=50,nchains=10)
## best model was 3, ARMA(0,1)
bestMod<-"dp-arma01model"
theta<-NA
phi<-1
cvGfortBeakB1<-crossValidation(dat=dat,niter=100000,thin=50,burn=50000,nchains=10,mymod=bestMod,init.theta=theta,init.phi=phi)
cvGfortBeakB2<-crossValidation(dat=dat,niter=100000,thin=50,burn=50000,nchains=10,mymod=bestMod,init.theta=theta,init.phi=phi)
## forecasting
drop<-3:10
forecastGfortBeakB1<-forecastChange(dat=dat,niter=100000,thin=50,burn=50000,nchains=10,mymod=bestMod,init.theta=theta,init.phi=phi,drop=drop)
forecastGfortBeakB2<-forecastChange(dat=dat,niter=100000,thin=50,burn=50000,nchains=10,mymod=bestMod,init.theta=theta,init.phi=phi,drop=drop)
## make plots
yr<-1973:2012
pdf("p-GfortBeak2TimeSeriesRaw.pdf",width=6,height=5)
par(mar=c(5,5,0.5,0.5))
plot(yr,p,type='b',pch=19,cex.lab="1.5",xlab="year",ylab="beak PC2")
dev.off()
yr<-1973:2011
pdf("dp-GfortBeak2TimeSeriesRaw.pdf",width=6,height=5)
par(mar=c(5,5,0.5,0.5))
plot(yr,dp,type='b',pch=19,cex.lab="1.5",xlab="year",ylab="change in beak PC2")
abline(h=0,lty=2,col="darkgray")
dev.off()
pred<-(cvGfortBeakB1[[1]][,3]+cvGfortBeakB2[[1]][,3])/2
o<-lm(dat$y[-c(1:2)]~pred[-c(1:2)])
pdf("dp-GfortBeak2ObsVpredict.pdf",width=5,height=5)
par(mar=c(5,5,0.5,0.5))
plot(pred,dp,pch=19,xlab="predicted change",ylab="observed change",cex.lab=1.4)
abline(o$coefficients,col="darkgray",lwd=1.5)
dev.off()
# Estimate Std. Error t value Pr(>|t|)
#(Intercept) 0.002860 0.005234 0.546 0.588
#pred[-c(1:2)] 0.462132 0.423803 1.090 0.283
#Residual standard error: 0.02893 on 35 degrees of freedom
#Multiple R-squared: 0.03286, Adjusted R-squared: 0.005224
#F-statistic: 1.189 on 1 and 35 DF, p-value: 0.283
forecastGfortBeakBComb<-(forecastGfortBeakB1[[1]]+forecastGfortBeakB2[[1]])/2
apply(forecastGfortBeakBComb,2,median)
# 0.05675793 0.06829066 -0.70978029 0.82428742
# median r^2, r, lb r, ub r
pdf("dp-GfortBeak2NyrPredict.pdf",width=5,height=5)
par(mar=c(5,5,.5,.5))
plot(drop,forecastGfortBeakBComb[,1],type='b',xlab="years predicted",ylab=expression(paste("predictive ",r^2)),cex.lab=1.4,ylim=c(0,1),pch=19)
dev.off()
pdf("dp-GfortBeak2NyrPredictCorCI.pdf",width=5,height=5)
par(mar=c(5,5,.5,.5))
plot(drop,forecastGfortBeakBComb[,2],type='n',xlab="years predicted",ylab="predictive correlation",cex.lab=1.4,ylim=c(-1,1),pch=19)
polygon(c(drop[-1],rev(drop[-1])),c(forecastGfortBeakBComb[-1,3],rev(forecastGfortBeakBComb[-1,4])),col="gray",density=50)
points(drop,forecastGfortBeakBComb[,2],type='b',pch=19)
abline(h=0,lty=2,col="darkblue")
dev.off()
save(list=ls(),file="predictGfortBeakSizeB.rdat")
######################### G. scandens, body size #########################
dat<-read.csv("scandens.csv")
y1<-dat$PC1.body[1:40] ## pc1 body size
y2<-dat$PC1.beak[1:40] ## pc1 beak size
y3<-dat$PC2.beak[1:40] ## pc2 beak size
p<-y1
## dp
dp<-p[2:40]-p[1:39]
dat<-list(y=dp,N=39)
## determine best arma model based on DIC
dicGscanBody1<-fitArmaModels(dat=dat,niter=100000,thin=50,nchains=10)
dicGscanBody2<-fitArmaModels(dat=dat,niter=100000,thin=50,nchains=10)
## best model was 6, ARMA(1,2)
bestMod<-"dp-arma12model"
theta<-1
phi<-2
cvGscanBody1<-crossValidation(dat=dat,niter=100000,thin=50,burn=50000,nchains=10,mymod=bestMod,init.theta=theta,init.phi=phi)
cvGscanBody2<-crossValidation(dat=dat,niter=100000,thin=50,burn=50000,nchains=10,mymod=bestMod,init.theta=theta,init.phi=phi)
## forecasting
drop<-3:10
forecastGscanBody1<-forecastChange(dat=dat,niter=100000,thin=50,burn=50000,nchains=10,mymod=bestMod,init.theta=theta,init.phi=phi,drop=drop)
forecastGscanBody2<-forecastChange(dat=dat,niter=100000,thin=50,burn=50000,nchains=10,mymod=bestMod,init.theta=theta,init.phi=phi,drop=drop)
## make plots
yr<-1973:2012
pdf("p-GscanBodyTimeSeriesRaw.pdf",width=6,height=5)
par(mar=c(5,5,0.5,0.5))
plot(yr,p,type='b',pch=19,cex.lab="1.5",xlab="year",ylab="body size")
dev.off()
yr<-1973:2011
pdf("dp-GscanBodyTimeSeriesRaw.pdf",width=6,height=5)
par(mar=c(5,5,0.5,0.5))
plot(yr,dp,type='b',pch=19,cex.lab="1.5",xlab="year",ylab="change in body size")
abline(h=0,lty=2,col="darkgray")
dev.off()
pred<-(cvGscanBody1[[1]][,3]+cvGscanBody2[[1]][,3])/2
## three values were off by orders of magnitude, here replace mean with more reasonable value
pred[which.min(cvGscanBody1[[1]][,3])]<-cvGscanBody2[[1]][which.min(cvGscanBody1[[1]][,3]),3]
pred[which.max(cvGscanBody1[[1]][,3])]<-cvGscanBody2[[1]][which.max(cvGscanBody1[[1]][,3]),3]
pred[which.max(cvGscanBody2[[1]][,3])]<-cvGscanBody1[[1]][which.min(cvGscanBody2[[1]][,3]),3]
o<-lm(dat$y[-c(1:2)]~pred[-c(1:2)])
pdf("dp-GscanBodyObsVpredict.pdf",width=5,height=5)
par(mar=c(5,5,0.5,0.5))
plot(pred,dp,pch=19,xlab="predicted change",ylab="observed change",cex.lab=1.4)
abline(o$coefficients,col="darkgray",lwd=1.5)
dev.off()
# Estimate Std. Error t value Pr(>|t|)
#(Intercept) -0.005377 0.033181 -0.162 0.872
#pred[-c(1:2)] -0.159569 0.111927 -1.426 0.163
#Residual standard error: 0.199 on 35 degrees of freedom
#Multiple R-squared: 0.05488, Adjusted R-squared: 0.02788
#F-statistic: 2.032 on 1 and 35 DF, p-value: 0.1628
forecastGscanBodyComb<-(forecastGscanBody1[[1]]+forecastGscanBody2[[1]])/2
apply(forecastGscanBodyComb,2,median)
# 0.3951535 0.6263869 -0.2028047 0.9260241
# median r^2, r, lb r, ub r
pdf("dp-GscanBodyNyrPredict.pdf",width=5,height=5)
par(mar=c(5,5,.5,.5))
plot(drop,forecastGscanBodyComb[,1],type='b',xlab="years predicted",ylab=expression(paste("predictive ",r^2)),cex.lab=1.4,ylim=c(0,1),pch=19)
dev.off()
pdf("dp-GscanBodyNyrPredictCorCI.pdf",width=5,height=5)
par(mar=c(5,5,.5,.5))
plot(drop,forecastGscanBodyComb[,2],type='n',xlab="years predicted",ylab="predictive correlation",cex.lab=1.4,ylim=c(-1,1),pch=19)
polygon(c(drop[-1],rev(drop[-1])),c(forecastGscanBodyComb[-1,3],rev(forecastGscanBodyComb[-1,4])),col="gray",density=50)
points(drop,forecastGscanBodyComb[,2],type='b',pch=19)
abline(h=0,lty=2,col="darkblue")
dev.off()
save(list=ls(),file="predictGscanBodySize.rdat")
######################### G. scandens, beak pc1 #########################
dat<-read.csv("scandens.csv")
y1<-dat$PC1.body[1:40] ## pc1 body size
y2<-dat$PC1.beak[1:40] ## pc1 beak size
y3<-dat$PC2.beak[1:40] ## pc2 beak size
p<-y2
## dp
dp<-p[2:40]-p[1:39]
dat<-list(y=dp,N=39)
## determine best arma model based on DIC
dicGscanBeakA1<-fitArmaModels(dat=dat,niter=100000,thin=50,nchains=10)
dicGscanBeakA2<-fitArmaModels(dat=dat,niter=100000,thin=50,nchains=10)
## best model was 6, ARMA(1,2)
bestMod<-"dp-arma12model"
theta<-1
phi<-2
cvGscanBeakA1<-crossValidation(dat=dat,niter=100000,thin=50,burn=50000,nchains=10,mymod=bestMod,init.theta=theta,init.phi=phi)
cvGscanBeakA2<-crossValidation(dat=dat,niter=100000,thin=50,burn=50000,nchains=10,mymod=bestMod,init.theta=theta,init.phi=phi)
## forecasting
drop<-3:10
forecastGscanBeakA1<-forecastChange(dat=dat,niter=100000,thin=50,burn=50000,nchains=10,mymod=bestMod,init.theta=theta,init.phi=phi,drop=drop)
forecastGscanBeakA2<-forecastChange(dat=dat,niter=100000,thin=50,burn=50000,nchains=10,mymod=bestMod,init.theta=theta,init.phi=phi,drop=drop)
## make plots
yr<-1973:2012
pdf("p-GscanBeak1TimeSeriesRaw.pdf",width=6,height=5)
par(mar=c(5,5,0.5,0.5))
plot(yr,p,type='b',pch=19,cex.lab="1.5",xlab="year",ylab="beak PC1")
dev.off()
yr<-1973:2011
pdf("dp-GscanBeak1TimeSeriesRaw.pdf",width=6,height=5)
par(mar=c(5,5,0.5,0.5))
plot(yr,dp,type='b',pch=19,cex.lab="1.5",xlab="year",ylab="change in beak PC1")
abline(h=0,lty=2,col="darkgray")
dev.off()
pred<-(cvGscanBeakA1[[1]][,3]+cvGscanBeakA2[[1]][,3])/2
## one value was off by orders of magnitude, here replace mean with more reasonable value
pred[which.min(cvGscanBeakA1[[1]][,3])]<-cvGscanBeakA2[[1]][which.min(cvGscanBeakA1[[1]][,3]),3]
o<-lm(dat$y[-c(1:2)]~pred[-c(1:2)])
pdf("dp-GscanBeak1ObsVpredict.pdf",width=5,height=5)
par(mar=c(5,5,0.5,0.5))
plot(pred,dp,pch=19,xlab="predicted change",ylab="observed change",cex.lab=1.4)
abline(o$coefficients,col="darkgray",lwd=1.5)
dev.off()
# Estimate Std. Error t value Pr(>|t|)
#(Intercept) 0.02175 0.02391 0.909 0.369
#pred[-c(1:2)] -0.22193 0.16008 -1.386 0.174
#Residual standard error: 0.1425 on 35 degrees of freedom
#Multiple R-squared: 0.05206, Adjusted R-squared: 0.02497
#F-statistic: 1.922 on 1 and 35 DF, p-value: 0.1744
forecastGscanBeakAComb<-(forecastGscanBeakA1[[1]]+forecastGscanBeakA2[[1]])/2
apply(forecastGscanBeakAComb,2,median)
# 0.1395220 0.2741161 -0.7107950 0.8694696
# median r^2, r, lb r, ub r
pdf("dp-GscanBeak1NyrPredict.pdf",width=5,height=5)
par(mar=c(5,5,.5,.5))
plot(drop,forecastGscanBeakAComb[,1],type='b',xlab="years predicted",ylab=expression(paste("predictive ",r^2)),cex.lab=1.4,ylim=c(0,1),pch=19)
dev.off()
pdf("dp-GscanBeak1NyrPredictCorCI.pdf",width=5,height=5)
par(mar=c(5,5,.5,.5))
plot(drop,forecastGscanBeakAComb[,2],type='n',xlab="years predicted",ylab="predictive correlation",cex.lab=1.4,ylim=c(-1,1),pch=19)
polygon(c(drop[-1],rev(drop[-1])),c(forecastGscanBeakAComb[-1,3],rev(forecastGscanBeakAComb[-1,4])),col="gray",density=50)
points(drop,forecastGscanBeakAComb[,2],type='b',pch=19)
abline(h=0,lty=2,col="darkblue")
dev.off()
save(list=ls(),file="predictGscanBeakSizeA.rdat")
######################### G. scandens, beak pc2 #########################
dat<-read.csv("scandens.csv")
y1<-dat$PC1.body[1:40] ## pc1 body size
y2<-dat$PC1.beak[1:40] ## pc1 beak size
y3<-dat$PC2.beak[1:40] ## pc2 beak size
p<-y3
## dp
dp<-p[2:40]-p[1:39]
dat<-list(y=dp,N=39)
## determine best arma model based on DIC
dicGscanBeakB1<-fitArmaModels(dat=dat,niter=100000,thin=50,nchains=10)
dicGscanBeakB2<-fitArmaModels(dat=dat,niter=100000,thin=50,nchains=10)
## best model was 6, ARMA(1,2)
bestMod<-"dp-arma12model"
theta<-1
phi<-2
cvGscanBeakB1<-crossValidation(dat=dat,niter=100000,thin=50,burn=50000,nchains=10,mymod=bestMod,init.theta=theta,init.phi=phi)
cvGscanBeakB2<-crossValidation(dat=dat,niter=100000,thin=50,burn=50000,nchains=10,mymod=bestMod,init.theta=theta,init.phi=phi)
## forecasting
drop<-3:10
forecastGscanBeakB1<-forecastChange(dat=dat,niter=100000,thin=50,burn=50000,nchains=10,mymod=bestMod,init.theta=theta,init.phi=phi,drop=drop)
forecastGscanBeakB2<-forecastChange(dat=dat,niter=100000,thin=50,burn=50000,nchains=10,mymod=bestMod,init.theta=theta,init.phi=phi,drop=drop)
## make plots
yr<-1973:2012
pdf("p-GscanBeak2TimeSeriesRaw.pdf",width=6,height=5)
par(mar=c(5,5,0.5,0.5))
plot(yr,p,type='b',pch=19,cex.lab="1.5",xlab="year",ylab="beak PC2")
dev.off()
yr<-1973:2011
pdf("dp-GscanBeak2TimeSeriesRaw.pdf",width=6,height=5)
par(mar=c(5,5,0.5,0.5))
plot(yr,dp,type='b',pch=19,cex.lab="1.5",xlab="year",ylab="change in beak PC2")
abline(h=0,lty=2,col="darkgray")
dev.off()
pred<-(cvGscanBeakB1[[1]][,3]+cvGscanBeakB2[[1]][,3])/2
## as with others, replacing absurd values with the estimate from the analysis that gave a reasonable result
x<-which(abs(pred) > 10)
## 18 25 29 ## note 18 is absurdly high in both cases, so conservatively has been replaced with NA
pred[x[1]]<-NA
pred[x[2]]<-cvGscanBeakB2[[1]][x[2],3]
pred[x[3]]<-cvGscanBeakB1[[1]][x[3],3]
o<-lm(dat$y[-c(1:2)]~pred[-c(1:2)])
pdf("dp-GscanBeak2ObsVpredict.pdf",width=5,height=5)
par(mar=c(5,5,0.5,0.5))
plot(pred,dp,pch=19,xlab="predicted change",ylab="observed change",cex.lab=1.4)
abline(o$coefficients,col="darkgray",lwd=1.5)
dev.off()
# Estimate Std. Error t value Pr(>|t|)
#(Intercept) -0.002938 0.012453 -0.236 0.815
#pred[-c(1:2)] 0.546308 0.371567 1.470 0.151
#Residual standard error: 0.06135 on 34 degrees of freedom
# (1 observation deleted due to missingness)
#Multiple R-squared: 0.05978, Adjusted R-squared: 0.03213
#F-statistic: 2.162 on 1 and 34 DF, p-value: 0.1507
forecastGscanBeakBComb<-(forecastGscanBeakB1[[1]]+forecastGscanBeakB2[[1]])/2
apply(forecastGscanBeakBComb,2,median)
# 0.18602100 0.08118622 -0.76748883 0.63518650
# median r^2, r, lb r, ub r
pdf("dp-GscanBeak2NyrPredict.pdf",width=5,height=5)
par(mar=c(5,5,.5,.5))
plot(drop,forecastGscanBeakBComb[,1],type='b',xlab="years predicted",ylab=expression(paste("predictive ",r^2)),cex.lab=1.4,ylim=c(0,1),pch=19)
dev.off()
pdf("dp-GscanBeak2NyrPredictCorCI.pdf",width=5,height=5)
par(mar=c(5,5,.5,.5))
plot(drop,forecastGscanBeakBComb[,2],type='n',xlab="years predicted",ylab="predictive correlation",cex.lab=1.4,ylim=c(-1,1),pch=19)
polygon(c(drop[-1],rev(drop[-1])),c(forecastGscanBeakBComb[-1,3],rev(forecastGscanBeakBComb[-1,4])),col="gray",density=50)
points(drop,forecastGscanBeakBComb[,2],type='b',pch=19)
abline(h=0,lty=2,col="darkblue")
dev.off()
save(list=ls(),file="predictGscanBeakSizeB.rdat")
######################### Panaxia dominula; fisher/wright medionigra #########################
dat<-read.csv("data_medionigra.csv")
p<-dat$medionigraAlleleCnt/dat$sampleSize2N
p<-p[1:39] ## conseq. years 1940-1978
## dp
dp<-p[2:39]-p[1:38]
dat<-list(y=dp,N=38)
## determine best arma model based on DIC
dicMedio1<-fitArmaModels(dat=dat,niter=100000,thin=50,nchains=10)
dicModio2<-fitArmaModels(dat=dat,niter=100000,thin=50,nchains=10)
## best model was 4, ARMA(1,1)
bestMod<-"dp-arma11model"
theta<-1
phi<-1
cvMedio1<-crossValidation(dat=dat,niter=100000,thin=50,burn=50000,nchains=10,mymod=bestMod,init.theta=theta,init.phi=phi)
cvMedio2<-crossValidation(dat=dat,niter=100000,thin=50,burn=50000,nchains=10,mymod=bestMod,init.theta=theta,init.phi=phi)
## forecasting
drop<-3:10
forecastMedio1<-forecastChange(dat=dat,niter=100000,thin=50,burn=50000,nchains=10,mymod=bestMod,init.theta=theta,init.phi=phi,drop=drop)
forecastMedio2<-forecastChange(dat=dat,niter=100000,thin=50,burn=50000,nchains=10,mymod=bestMod,init.theta=theta,init.phi=phi,drop=drop)
## make plots
yr<-1940:1978
pdf("p-MedioTimeSeriesRaw.pdf",width=6,height=5)
par(mar=c(5,5,0.5,0.5))
plot(yr,p,type='b',pch=19,cex.lab="1.5",xlab="year",ylab="morph frequency")
dev.off()
yr<-1940:1977
pdf("dp-MedioTimeSeriesRaw.pdf",width=6,height=5)
par(mar=c(5,5,0.5,0.5))
plot(yr,dp,type='b',pch=19,cex.lab="1.5",xlab="year",ylab="change in morph frequency")
abline(h=0,lty=2,col="darkgray")
dev.off()
pred<-(cvMedio1[[1]][,3]+cvMedio2[[1]][,3])/2
o<-lm(dat$y[-c(1:2)]~pred[-c(1:2)])
pdf("dp-MedioObsVpredict.pdf",width=5,height=5)
par(mar=c(5,5,0.5,0.5))
plot(pred,dp,pch=19,xlab="predicted change",ylab="observed change",cex.lab=1.4)
abline(o$coefficients,col="darkgray",lwd=1.5)
dev.off()
# Estimate Std. Error t value Pr(>|t|)
#(Intercept) -0.0005844 0.0021393 -0.273 0.786
#pred[-c(1:2)] 0.4419179 0.5766178 0.766 0.449
#Residual standard error: 0.01207 on 34 degrees of freedom
#Multiple R-squared: 0.01698, Adjusted R-squared: -0.01193
#F-statistic: 0.5874 on 1 and 34 DF, p-value: 0.4487
forecastMedioComb<-(forecastMedio1[[1]]+forecastMedio2[[1]])/2
apply(forecastMedioComb,2,median)
# 0.180800958 0.000594893 -0.804753467 0.709522671
# median r^2, r, lb r, ub r
pdf("dp-MedioNyrPredict.pdf",width=5,height=5)
par(mar=c(5,5,.5,.5))
plot(drop,forecastMedioComb[,1],type='b',xlab="years predicted",ylab=expression(paste("predictive ",r^2)),cex.lab=1.4,ylim=c(0,1),pch=19)
dev.off()
pdf("dp-MedioNyrPredictCorCI.pdf",width=5,height=5)
par(mar=c(5,5,.5,.5))
plot(drop,forecastMedioComb[,2],type='n',xlab="years predicted",ylab="predictive correlation",cex.lab=1.4,ylim=c(-1,1),pch=19)
polygon(c(drop[-1],rev(drop[-1])),c(forecastMedioComb[-1,3],rev(forecastMedioComb[-1,4])),col="gray",density=50)
points(drop,forecastMedioComb[,2],type='b',pch=19)
abline(h=0,lty=2,col="darkblue")
dev.off()
save(list=ls(),file="predictMedionigra.rdat")
######################### Biston betularia; peppered moth #########################
dat<-read.csv("mothdat.csv")
p<-dat$bb_mel/(dat$bb_mel+dat$bb_nonmel)
## dp
dp<-p[2:34]-p[1:33]
n<-(dat$bb_mel+dat$bb_nonmel)
nbar<-(n[2:34]+n[1:33])/2
## 29, 30 and 31 all have nbar < 10; ding 29:33 (the end)
dp<-dp[-c(29:33)]
p<-p[-c(30:34)]
dat<-list(y=dp,N=28)
## determine best arma model based on DIC
dicPepperMoth1<-fitArmaModels(dat=dat,niter=100000,thin=50,nchains=10)
dicPepperMoth2<-fitArmaModels(dat=dat,niter=100000,thin=50,nchains=10)
## best model was 1, then 2, ARMA(1,0) ... 3 and 4 similar
bestMod<-"dp-arma10model"
theta<-1
phi<-0
cvPepperMoth1<-crossValidation(dat=dat,niter=100000,thin=50,burn=50000,nchains=10,mymod=bestMod,init.theta=theta,init.phi=phi)
cvPepperMoth2<-crossValidation(dat=dat,niter=100000,thin=50,burn=50000,nchains=10,mymod=bestMod,init.theta=theta,init.phi=phi)
## forecasting
drop<-3:10
forecastPepperMoth1<-forecastChange(dat=dat,niter=100000,thin=50,burn=50000,nchains=10,mymod=bestMod,init.theta=theta,init.phi=phi,drop=drop)
forecastPepperMoth2<-forecastChange(dat=dat,niter=100000,thin=50,burn=50000,nchains=10,mymod=bestMod,init.theta=theta,init.phi=phi,drop=drop)
## make plots
yr<-1967:1995
pdf("p-PepperTimeSeriesRaw.pdf",width=6,height=5)
par(mar=c(5,5,0.5,0.5))
plot(yr,p,type='b',pch=19,cex.lab="1.5",xlab="year",ylab="morph frequency")
dev.off()
yr<-1967:1994
pdf("dp-PepperTimeSeriesRaw.pdf",width=6,height=5)
par(mar=c(5,5,0.5,0.5))
plot(yr,dp,type='b',pch=19,cex.lab="1.5",xlab="year",ylab="change in morph frequency")
abline(h=0,lty=2,col="darkgray")
dev.off()
pred<-(cvPepperMoth1[[1]][,3]+cvPepperMoth2[[1]][,3])/2
o<-lm(dat$y[-c(1:2)]~pred[-c(1:2)])
pdf("dp-PepperObsVpredict.pdf",width=5,height=5)
par(mar=c(5,5,0.5,0.5))
plot(pred,dp,pch=19,xlab="predicted change",ylab="observed change",cex.lab=1.4)
abline(o$coefficients,col="darkgray",lwd=1.5)
dev.off()
# Estimate Std. Error t value Pr(>|t|)
#(Intercept) -0.05174 0.01186 -4.362 0.000335 ***
#pred[-c(1:2)] -1.80925 0.29897 -6.052 8.05e-06 ***
#Residual standard error: 0.05128 on 19 degrees of freedom
# (5 observations deleted due to missingness)
#Multiple R-squared: 0.6584, Adjusted R-squared: 0.6404
#F-statistic: 36.62 on 1 and 19 DF, p-value: 8.049e-06
forecastPepperComb<-(forecastPepperMoth1[[1]]+forecastPepperMoth2[[1]])/2
apply(forecastPepperComb,2,median)
# 0.02692756 0.15945025 -0.84112676 0.91261486
# median r^2, r, lb r, ub r
pdf("dp-PepperNyrPredict.pdf",width=5,height=5)
par(mar=c(5,5,.5,.5))
plot(drop,forecastPepperComb[,1],type='b',xlab="years predicted",ylab=expression(paste("predictive ",r^2)),cex.lab=1.4,ylim=c(0,1),pch=19)
dev.off()
pdf("dp-PepperNyrPredictCorCI.pdf",width=5,height=5)
par(mar=c(5,5,.5,.5))
plot(drop,forecastPepperComb[,2],type='n',xlab="years predicted",ylab="predictive correlation",cex.lab=1.4,ylim=c(-1,1),pch=19)
polygon(c(drop[-1],rev(drop[-1])),c(forecastPepperComb[-1,3],rev(forecastPepperComb[-1,4])),col="gray",density=50)
points(drop,forecastPepperComb[,2],type='b',pch=19)
abline(h=0,lty=2,col="darkblue")
dev.off()
save(list=ls(),file="predictPepperMoth.rdat")