#read in the data - first column is years, second column months:: sessions = read.table("sessions.txt") #Random walk innovation errors: err.gamma.mu = rep(0.4, 12) err.gamma.sigma = rep(0.075, 12) err.epsilon.mu = rep(0.4, 28) err.epsilon.sigma = rep(0.075, 28) err.beta.mu = 0.3 err.beta.sigma = 0.08 #Now run the random effects model for site/seasonal effects (see "fitRC.R"):: R = RE(n = 50000, hyper = c(0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0, 0, 100, 100), sessions = sessions) #seasonal.aprobs2 = matrix(ncol=2, nrow=12) #for(i in 1:12){ # seasonal.aprobs2[i,] = c(mean(R$APROB.GAMMA.MU[,i], na.rm=TRUE), mean(R$APROB.GAMMA.SIGMA[,i], na.rm=TRUE))} #site.aprobs2 = matrix(ncol=2, nrow=28) #for(i in 1:28){ # site.aprobs2[i,] = c(mean(R$APROB.EPSILON.MU[,i], na.rm=TRUE), mean(R$APROB.EPSILON.SIGMA[,i], #na.rm=TRUE))} #trend.aprobs2 = vector("numeric", length = 2) #trend.aprobs2 = c(mean(R$APROB.BETA.MU, na.rm=TRUE), mean(R$APROB.BETA.SIGMA, na.rm=TRUE)) #par(mfrow=c(3,2)) #plot(seasonal.aprobs1[,1], seasonal.aprobs2[,1]) #abline(0,1) #plot(seasonal.aprobs1[,2], seasonal.aprobs2[,2]) #abline(0,1) #plot(site.aprobs1[,1], site.aprobs2[,1]) #abline(0,1) #plot(site.aprobs1[,2], site.aprobs2[,2]) #abline(0,1) #Extract and plot seasonal and site effects. This assumes the MCMC results are stored in "R" n_m = dim(R$GAMMA.MU)[2] #No. of seasons n_s = dim(R$EPSILON.MU)[2] #No. of sites B = 1000 #Burn-in period #First, seasonal effects: seasonal.summaries = matrix(ncol = 6, nrow = n_m) for(i in 1:n_m){ seasonal.summaries[i, 2] = mean(R$GAMMA.MU[(B+1):dim(R$GAMMA.MU)[1],i]) seasonal.summaries[i, 1] = quantile(R$GAMMA.MU[(B+1):dim(R$GAMMA.MU)[1],i], 0.2) seasonal.summaries[i, 3] = quantile(R$GAMMA.MU[(B+1):dim(R$GAMMA.MU)[1],i], 0.8) seasonal.summaries[i, 5] = mean(R$GAMMA.SIGMA[(B+1):dim(R$GAMMA.SIGMA)[1],i]) seasonal.summaries[i, 4] = quantile(R$GAMMA.SIGMA[(B+1):dim(R$GAMMA.SIGMA)[1],i], 0.2) seasonal.summaries[i, 6] = quantile(R$GAMMA.SIGMA[(B+1):dim(R$GAMMA.SIGMA)[1],i], 0.8)} pdf("re2.pdf") par(mfrow=c(4,2)) plot(seasonal.summaries[,2], pch=16, lwd=3, ylim=c(min(seasonal.summaries[,1:3]), max(seasonal.summaries[,1:3])), xlab="season", ylab="gamma.mu",main="Seasonal effects for mu") for(i in 1:n_m){ lines(c(i,i), c(seasonal.summaries[i,1], seasonal.summaries[i,3]))} plot(seasonal.summaries[,5], pch=16, lwd=3, xlab="season", ylab="gamma.sigma", ylim=c(min(seasonal.summaries[,4:6]), max(seasonal.summaries[,4:6])),main="Seasonal effects for sigma") for(i in 1:n_m){ lines(c(i,i), c(seasonal.summaries[i,4], seasonal.summaries[i,6]))} #Now, zone effects: lat = c(55.21991, 55.81432, 53.26850, 51.27338, 51.27036, 54.72890, 55.78273, 54.31434, 51.58770, 51.26730, 52.20911, 52.21002, 52.66499, 51.74698, 50.95870, 50.95654, 53.13957, 53.13955, 52.87349, 52.06541, 51.15721, 51.15672, 50.97787, 50.97629, 53.13721, 53.13883, 52.30624, 52.30775) long = c(-3.41047, -4.06255, -2.80313, 0.03689, 0.03996, -1.52275, -4.03449, -2.63503, -2.83749, -0.61361, -0.94353, -0.94355, -1.96741, -4.06537, -1.44842, -1.44853, -2.33755, -2.33567, -2.16146, -2.15724, -1.34629, -1.34409, -3.14783, -3.14681, -1.33318, -1.32842, -1.12281, -1.12087) site.summaries = matrix(ncol = 6, nrow = n_s) for(i in 1:n_s){ site.summaries[i, 2] = mean(R$EPSILON.MU[(B+1):dim(R$EPSILON.MU)[1],i]) site.summaries[i, 1] = quantile(R$EPSILON.MU[(B+1):dim(R$EPSILON.MU)[1],i], 0.05) site.summaries[i, 3] = quantile(R$EPSILON.MU[(B+1):dim(R$EPSILON.MU)[1],i], 0.95) site.summaries[i, 5] = mean(R$EPSILON.SIGMA[(B+1):dim(R$EPSILON.SIGMA)[1],i]) site.summaries[i, 4] = quantile(R$EPSILON.SIGMA[(B+1):dim(R$EPSILON.SIGMA)[1],i], 0.05) site.summaries[i, 6] = quantile(R$EPSILON.SIGMA[(B+1):dim(R$EPSILON.SIGMA)[1],i], 0.95)} plot(long,site.summaries[,2], pch=16, lwd=3, ylim=c(min(site.summaries[,1:3]), max(site.summaries[,1:3])), xlab="longitude", ylab="epsilon.mu", main="Site effects for mu") for(i in 1:n_s){ lines(c(long[i],long[i]), c(site.summaries[i,1], site.summaries[i,3]))} fit.loess = loess(site.summaries[,2]~long) sm = cbind(long, fit.loess$fitted) sm2 = sm[order(sm[,1]),] lines(sm2[,1], sm2[,2], col="red") plot(long, site.summaries[,5], pch=16, lwd=3, ylim=c(min(site.summaries[,4:6]), max(site.summaries[,4:6])), xlab="longitude", ylab="epsilon.sigma",main="Site effects for sigma") for(i in 1:n_s){ lines(c(long[i],long[i]), c(site.summaries[i,4], site.summaries[i,6]))} fit.loess = loess(site.summaries[,5]~long) sm = cbind(long, fit.loess$fitted) sm2 = sm[order(sm[,1]),] lines(sm2[,1], sm2[,2], col="red") plot(lat,site.summaries[,2], pch=16, lwd=3, ylim=c(min(site.summaries[,1:3]), max(site.summaries[,1:3])), xlab="latitude", ylab="epsilon.mu", main="Site effects for mu") for(i in 1:n_s){ lines(c(lat[i],lat[i]), c(site.summaries[i,1], site.summaries[i,3]))} fit.loess = loess(site.summaries[,2]~lat) sm = cbind(lat, fit.loess$fitted) sm2 = sm[order(sm[,1]),] lines(sm2[,1], sm2[,2], col="red") plot(lat, site.summaries[,5], pch=16, lwd=3, ylim=c(min(site.summaries[,4:6]), max(site.summaries[,4:6])), xlab="latitude", ylab="epsilon.sigma", main="Site effects for sigma") for(i in 1:n_s){ lines(c(lat[i],lat[i]), c(site.summaries[i,4], site.summaries[i,6]))} fit.loess = loess(site.summaries[,5]~lat) sm = cbind(lat, fit.loess$fitted) sm2 = sm[order(sm[,1]),] lines(sm2[,1], sm2[,2], col="red") plot(density(R$BETA.MU[(B+1):length(R$BETA.MU)], bw=0.06), main="Trend effect for mu", xlab="beta_mu", ylab="Posterior density") plot(density(R$BETA.SIGMA[(B+1):length(R$BETA.MU)],bw=0.1), main="Trend effect for sigma", xlab="beta_sigma", ylab="Posterior density") #Now plot map:: library(maps) library(mapdata) map('worldHires', 'UK', xlim=c(-6,1.8), ylim=c(50, 59), xlab="longitude") points(long, lat, col="red", pch=20) #ll = seq(-87, -79) #for(i in 1:length(ll)){ # abline(v = ll[i], lty=2, col="gray")} #box() #axis(1) #plot(coord[,1], site.summaries[,2], ylim=c(min(site.summaries), max(site.summaries)), pch=20, xlab="longitude", ylab="epsilon") #for(i in 1:n_s){ # lines(c(coord[i,1], coord[i,1]), c(site.summaries[i,1], site.summaries[i,3]))} #fit.loess = loess(site.summaries[,2]~coord[,1]) #sm = cbind(coord[,1], fit.loess$fitted) #sm2 = sm[order(sm[,1]),] #lines(sm2[,1], sm2[,2], col="red")