####################################### # OUTPUT SUMMARIES FOR GEOSPATIAL MODEL ####################################### # Produces summaries of model parameters, does some posterior predictive # model checks and produces some figures # Assumes variables yrs, nsites.active, effort.positiveEffort, xyEffort, lmnc # and geomodel.res are defined # Final plot assumes file RefugeCoords.csv is in current directory outName <- paste("GeospatialModel", format(Sys.Date(), "%Y%m%d"), " - ", sep="") # Save entire bugs() output as .RData file (can read-in later with 'load' command) save(geomodel.res, file = paste(outName, "bugsOut.RData", sep = "")) # Also save a data matrix, for potential later re-use d2 <- data.frame(Site = as.numeric(rownames(effort.positiveEffort))) # Add x and y location of monitoring locations d2 <- cbind(d2, xyEffort) # Add mean click count tmp.names <- names(d2); d2 <- cbind(d2, exp(lmnc) - 1); names(d2) <- c(tmp.names, paste0("mnc1", 1:nyrs)) # Add log mean click count tmp.names <- names(d2); d2 <- cbind(d2, lmnc); names(d2) <- c(tmp.names, paste0("lmnc1", 1:nyrs)) # Add sampling effort tmp.names <- names(d2); d2 <- cbind(d2, effort.positiveEffort); names(d2) <- c(tmp.names, paste0("n1", 1:nyrs)) save(d2, file = "d2.RData") #Convert geomodel.res to an mcmc object, to facilitate checks and summaries geomodel.res.mcmc <- as.mcmc(geomodel.res) # Geometric mean population change, as an example 2*(1 - pnorm(abs(geweke.diag(geomodel.res.mcmc$sims.matrix[, 'geochng'])$z))) heidel.diag(geomodel.res.mcmc$sims.matrix[,'geochng']) bm(geomodel.res.mcmc$sims.matrix[, 'geochng']) summary(geomodel.res.mcmc$sims.matrix[, 'geochng']) tiff('GeospatialModel-plot-posteriorGeoChngRate.tiff', width = 1200, height = 800) par(mar = c(6, 6, 0.2, 0.2)) plot(density(geomodel.res.mcmc$sims.matrix[, 'geochng'], bw = 0.04), ylab = 'Posterior Density', main = '', xlab = 'Annual Proportional Change', cex.axis = 2, cex.lab = 3, lwd = 5) dev.off() quantile(geomodel.res.mcmc$sims.matrix[, 'geochng'], p = c(0.025, 0.975)) write.csv(geomodel.res.mcmc$sims.matrix[, 'geochng'], 'geochng.csv', row.names = FALSE) # Examine lambda values from each year; also write to csv for subsequent use # in model averaging mean(geomodel.res.mcmc$sims.matrix[, 'chng[1]']) quantile(geomodel.res.mcmc$sims.matrix[, 'chng[1]'], p = c(0.025, .5, 0.975)) write.csv(geomodel.res.mcmc$sims.matrix[, 'chng[1]'], 'geochng11to12.csv', row.names = FALSE) mean(geomodel.res.mcmc$sims.matrix[, 'chng[2]']) quantile(geomodel.res.mcmc$sims.matrix[, 'chng[2]'], p = c(0.025, .5, 0.975)) write.csv(geomodel.res.mcmc$sims.matrix[, 'chng[2]'], 'geochng12to13.csv', row.names = FALSE) mean(geomodel.res.mcmc$sims.matrix[, 'chng[3]']) quantile(geomodel.res.mcmc$sims.matrix[, 'chng[3]'], p = c(0.025, .5, 0.975)) write.csv(geomodel.res.mcmc$sims.matrix[, 'chng[3]'], 'geochng13to14.csv', row.names = FALSE) mean(geomodel.res.mcmc$sims.matrix[, 'chng[4]']) quantile(geomodel.res.mcmc$sims.matrix[, 'chng[4]'], p = c(0.025, .5, 0.975)) write.csv(geomodel.res.mcmc$sims.matrix[, 'chng[4]'], 'geochng14to15.csv', row.names = FALSE) mean(geomodel.res.mcmc$sims.matrix[, 'chng[5]']) quantile(geomodel.res.mcmc$sims.matrix[, 'chng[5]'], p = c(0.025, .5, 0.975)) write.csv(geomodel.res.mcmc$sims.matrix[, 'chng[5]'], 'geochng15to16.csv', row.names = FALSE) mean(geomodel.res.mcmc$sims.matrix[, 'chng[6]']) quantile(geomodel.res.mcmc$sims.matrix[, 'chng[6]'], p = c(0.025, .5, 0.975)) write.csv(geomodel.res.mcmc$sims.matrix[, 'chng[6]'], 'geochng16to17.csv', row.names = FALSE) mean(geomodel.res.mcmc$sims.matrix[, 'chng[7]']) quantile(geomodel.res.mcmc$sims.matrix[, 'chng[7]'], p = c(0.025, .5, 0.975)) write.csv(geomodel.res.mcmc$sims.matrix[, 'chng[7]'], 'geochng17to18.csv', row.names = FALSE) # Create a data frame with smoothed predicted counts, based on d2 data frame d3 <- d2 for(i in 1:nyrs){ explam.names <- paste0(paste0("explam[", i, ","), 1:nsites.active, "]") explam <- geomodel.res.mcmc$sims.matrix[, explam.names] spSmooth <- apply(explam, 2, median) tmp.names <- names(d3) d3 <- cbind(d3, spSmooth) names(d3) <- c(tmp.names, paste0("smo1", i)) } # Perform posterior predictive model checks for each monitoring station and year # One set of checks are based on the chi-sq discrepancies of Gelman, Meng and # Stern (1996, equation 8), with stations passing the check if their discrepancy # is less than the 99th percentile critical value from a ChiSq, so we expect one # or two failures # The other set are posterior simulations, from which p-values are calculated as # proportion of simulations greater than or equal to the observed value - referred to by Gelman # et al (Bayesian Data Analysis 2014) as marginal predictive checks. q.val <- qchisq(0.99, df = 1) marginal.pred.p <- matrix(NA, nyrs, nsites.active) for(i in 1:nyrs){ lmnc.yr <- d2[, paste0("lmnc1", i)] ind <- !is.na(lmnc.yr) which.ind <- which(ind) lmnc.yr <- lmnc.yr[ind] explam.names <- paste0(paste0("explam[", i, ","), 1:nsites.active, "]") d2MCMC <- log(geomodel.res.mcmc$sims.matrix[, explam.names[ind]] + 1) postPredCheck <- rep(NA, times = sum(ind)) n1 <- d3[ind , paste0("n1", i)] for(j in 1:sum(ind)) { #Chi-sq check Rit <- (d2MCMC[, j] - lmnc.yr[j]) ^ 2 / (geomodel.res.mcmc$sims.matrix[, "sigmaEps"] ^2 / n1[j]) postPredCheck[j] <- mean(Rit) < q.val #Marginal predictive check post.pred <- rnorm(nrow(d2MCMC), d2MCMC[, j], geomodel.res.mcmc$sims.matrix[, "sigmaEps"]/sqrt(n1[j])) marginal.pred.p[i, which.ind[j]] <- sum(post.pred >= lmnc.yr[j]) / length(post.pred) } print(i) print(postPredCheck) } #Plot spatial predictions by year d3SPDF <- d3 coordinates(d3SPDF) <- c("x","y") max.clicks <- 1200 #Create refuge coordinates polygon refuge.coords <- read.csv(file = "RefugeCoords.csv") Sr1 <- Polygon(refuge.coords) Sp1 <- Polygons(list(Sr1), "Refuge") spRefuge <- SpatialPolygons(list(Sp1)) #Define utility function for colors getcolor<-function(clicks, max.clicks, c.ramp, log.scale){ crampfunc<-colorRamp(c.ramp) if(log.scale) { logval <- log10(clicks + 1) logval[logval < 0] <- 0 scaledlogval <- logval / log10(max.clicks + 1) return(rgb(crampfunc(scaledlogval), maxColorValue = 255)) } else { return(rgb(crampfunc(clicks/max.clicks), maxColorValue = 255)) } } #Plot function add.plot<-function(year, refuge.outline, coords, effort, is.no.effort, clicks, max.clicks, c.ramp = c("white", "darkblue"), log.scale = TRUE, highlight.outliers = FALSE){ #Makes a single plot cexstart <- 0.5 cexmult <- 1.4 plot(refuge.outline) text(4 - 2.577163, 50 - 7.61422, year, cex = 4, pos = 4) n4cex = as.integer(cut(effort, breaks = c(0,1,15,30,45,63))) n4cex[n4cex == 0] = 1 points(coords, col = getcolor(clicks, max.clicks, c.ramp, log.scale = log.scale), pch = 19, cex = cexstart + cexmult * n4cex) points(coords, pch = 1, cex = cexstart + cexmult * n4cex) if(highlight.outliers){ #draw a red boundary around those that are within 5% of zero or max.clicks outliers <- (clicks <= max.clicks*0.05) | (clicks >= max.clicks*0.95) points(coords[outliers, 1], coords[outliers, 2], pch = 1, cex = cexstart + cexmult * n4cex[outliers], lwd = 7, col = "red") } if(sum(is.no.effort) > 0){ points(coords[is.no.effort, 1], coords[is.no.effort, 2], pch = 1, cex = 6) points(coords[is.no.effort, 1], coords[is.no.effort, 2], pch = 4, cex = 4) } } #Note - this makes a very large (publication quality) plot file! tiff(file = 'GeospatialModel-plot-SpSmooth.tiff', width=7012, height=11250, res=720) par(mar = c(0, 0, 0, 0)) layout(matrix(c(1, 2, 3, 4, 5, 6, 7, 8, 9), nrow = 3, byrow = TRUE)) add.plot("2011", spRefuge, coordinates(d3SPDF), d3SPDF@data$n11, is.na(d3SPDF@data$mnc11), d3SPDF@data$smo11,max.clicks) add.plot("2012", spRefuge, coordinates(d3SPDF), d3SPDF@data$n12, is.na(d3SPDF@data$mnc12), d3SPDF@data$smo12,max.clicks) add.plot("2013", spRefuge, coordinates(d3SPDF), d3SPDF@data$n13, is.na(d3SPDF@data$mnc13), d3SPDF@data$smo13,max.clicks) add.plot("2014", spRefuge, coordinates(d3SPDF), d3SPDF@data$n14, is.na(d3SPDF@data$mnc14), d3SPDF@data$smo14,max.clicks) add.plot("2015", spRefuge, coordinates(d3SPDF), d3SPDF@data$n15, is.na(d3SPDF@data$mnc15), d3SPDF@data$smo15,max.clicks) add.plot("2016", spRefuge, coordinates(d3SPDF), d3SPDF@data$n16, is.na(d3SPDF@data$mnc16), d3SPDF@data$smo16,max.clicks) add.plot("2017", spRefuge, coordinates(d3SPDF), d3SPDF@data$n17, is.na(d3SPDF@data$mnc17), d3SPDF@data$smo17,max.clicks) add.plot("2018", spRefuge, coordinates(d3SPDF), d3SPDF@data$n18, is.na(d3SPDF@data$mnc18), d3SPDF@data$smo18,max.clicks) #Add legends plot(c(0, 1), c(0, 1), type = "n", xlab = "", ylab = "", bty = "n", xaxt = "n", yaxt = "n") #Clicks/day legend n.click.levels <- 1000 click.levels <- 10 ^ (seq(log10(max.clicks + 1), 0, length = n.click.levels)) - 1 legend.colors<-getcolor(click.levels, max.clicks, log.scale = TRUE, c.ramp = c("white", "darkblue")) rasterImage(as.raster(legend.colors), 0.05, 0.05, 0.2, 0.85) rect(0.05, 0.05, 0.2, 0.85) #Work out where the labels go legend.labels <- c(0, 1, 10, 100, max.clicks) posn<-numeric(length(legend.labels)) for(i in 1:length(legend.labels)) { posn[i] <- which.min(click.levels > legend.labels[i]) } posn <- (n.click.levels - posn) / n.click.levels for(i in 1:length(legend.labels)){ text(0.2, (0.85 - 0.05) * posn[i] + 0.05, legend.labels[i], cex = 2.5, pos = 4) } text(0.18,0.95, labels = 'clicks/day', cex = 2.5) #Sample days legend cexstart <- 0.5 cexmult <- 1.4 points(c(0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7), c(0.1667, 0.3, 0.4333, 0.5667, 0.7, 0.7, 0.7), pch = c(1, 1, 1, 1, 1, 1, 4), cex = c(cexstart + cexmult * 5, cexstart + cexmult * 4, cexstart + cexmult * 3, cexstart + cexmult * 2, cexstart + cexmult, 6, 4), lwd = 2) text(0.8, 0.83, labels = 'days of\nsampling', cex = 2.5) text(c(0.7, 0.7, 0.7, 0.7, 0.7), c(0.1667, 0.3, 0.4333, 0.5667, 0.7), pos = 4, labels = c('46-63', '31-45', '16-30', '1-15', '0'), offset = 1.7, cex = 2.5) dev.off() #Marginal predictive check p-values c.ramp <- c("yellow", "blue") tiff(file = 'GeospatialModel-plot-post-pred-p.tiff', width=7012, height=11250, res=720) par(mar = c(0, 0, 0, 0)) layout(matrix(c(1, 2, 3, 4, 5, 6, 7, 8, 9), nrow = 3, byrow = TRUE)) tmp.marginal.p <- marginal.pred.p[1, ] tmp.marginal.p [is.na(tmp.marginal.p)] <- 0.5 add.plot("2011", spRefuge, coordinates(d3SPDF), d3SPDF@data$n11, is.na(d3SPDF@data$mnc11), tmp.marginal.p, 1, c.ramp = c.ramp, log.scale = FALSE, highlight.outliers = TRUE) tmp.marginal.p <- marginal.pred.p[2, ] tmp.marginal.p [is.na(tmp.marginal.p)] <- 0.5 add.plot("2012", spRefuge, coordinates(d3SPDF), d3SPDF@data$n12, is.na(d3SPDF@data$mnc12), tmp.marginal.p, 1, c.ramp = c.ramp, log.scale = FALSE, highlight.outliers = TRUE) tmp.marginal.p <- marginal.pred.p[3, ] tmp.marginal.p [is.na(tmp.marginal.p)] <- 0.5 add.plot("2013", spRefuge, coordinates(d3SPDF), d3SPDF@data$n13, is.na(d3SPDF@data$mnc13), tmp.marginal.p, 1, c.ramp = c.ramp, log.scale = FALSE, highlight.outliers = TRUE) tmp.marginal.p <- marginal.pred.p[4, ] tmp.marginal.p [is.na(tmp.marginal.p)] <- 0.5 add.plot("2014", spRefuge, coordinates(d3SPDF), d3SPDF@data$n14, is.na(d3SPDF@data$mnc14), tmp.marginal.p, 1, c.ramp = c.ramp, log.scale = FALSE, highlight.outliers = TRUE) tmp.marginal.p <- marginal.pred.p[5, ] tmp.marginal.p [is.na(tmp.marginal.p)] <- 0.5 add.plot("2015", spRefuge, coordinates(d3SPDF), d3SPDF@data$n15, is.na(d3SPDF@data$mnc15), tmp.marginal.p, 1, c.ramp = c.ramp, log.scale = FALSE, highlight.outliers = TRUE) tmp.marginal.p <- marginal.pred.p[6, ] tmp.marginal.p [is.na(tmp.marginal.p)] <- 0.5 add.plot("2016", spRefuge, coordinates(d3SPDF), d3SPDF@data$n16, is.na(d3SPDF@data$mnc16), tmp.marginal.p, 1, c.ramp = c.ramp, log.scale = FALSE, highlight.outliers = TRUE) tmp.marginal.p <- marginal.pred.p[7, ] tmp.marginal.p [is.na(tmp.marginal.p)] <- 0.5 add.plot("2017", spRefuge, coordinates(d3SPDF), d3SPDF@data$n17, is.na(d3SPDF@data$mnc17), tmp.marginal.p, 1, c.ramp = c.ramp, log.scale = FALSE, highlight.outliers = TRUE) tmp.marginal.p <- marginal.pred.p[8, ] tmp.marginal.p [is.na(tmp.marginal.p)] <- 0.5 add.plot("2018", spRefuge, coordinates(d3SPDF), d3SPDF@data$n18, is.na(d3SPDF@data$mnc18), tmp.marginal.p, 1, c.ramp = c.ramp, log.scale = FALSE, highlight.outliers = TRUE) #Add legend plot(c(0,1), c(0,1), type = "n", xlab = "", ylab = "", bty = "n", xaxt = "n", yaxt = "n") n.p.levels <- 2 p.levels<-seq(0, 1, length = n.p.levels) legend.colors<-getcolor(p.levels, 1, c.ramp = rev(c.ramp), log.scale = FALSE) rasterImage(as.raster(legend.colors), 0.05, 0.05, 0.2, 0.85) rect(0.05, 0.05, 0.2, 0.85) #Work out where the labels go legend.labels<-seq(0, 1, length = 5) for(i in 1:length(legend.labels)){ text(0.2, (0.85 - 0.05) * legend.labels[i] + 0.05, legend.labels[i], cex = 2.5, pos = 4) } text(0.18, 0.95, labels = 'p-value', cex = 2.5) dev.off() #Add a plot of the mixture model predictive check p-value results #Marginal predictive check p-values c.ramp <- c("yellow", "blue") tiff(file = 'MixtureModel-plot-post-pred-p.tiff', width=7012, height=11250, res=720) par(mar = c(0, 0, 0, 0)) layout(matrix(c(1, 2, 3, 4, 5, 6, 7, 8, 9), nrow = 3, byrow = TRUE)) tmp.marginal.p <- marginal.pred.p.byyear.bysite[1, ] tmp.marginal.p [is.na(tmp.marginal.p)] <- 0.5 add.plot("2011",spRefuge,coordinates(d3SPDF), d3SPDF@data$n11, is.na(d3SPDF@data$mnc11), tmp.marginal.p, 1, c.ramp = c.ramp, log.scale = FALSE, highlight.outliers = TRUE) tmp.marginal.p <- marginal.pred.p.byyear.bysite[2, ] tmp.marginal.p [is.na(tmp.marginal.p)] <- 0.5 add.plot("2012", spRefuge,coordinates(d3SPDF), d3SPDF@data$n12, is.na(d3SPDF@data$mnc12), tmp.marginal.p, 1, c.ramp = c.ramp, log.scale = FALSE, highlight.outliers = TRUE) tmp.marginal.p <- marginal.pred.p.byyear.bysite[3, ] tmp.marginal.p [is.na(tmp.marginal.p)] <- 0.5 add.plot("2013", spRefuge,coordinates(d3SPDF), d3SPDF@data$n13, is.na(d3SPDF@data$mnc13), tmp.marginal.p, 1, c.ramp = c.ramp, log.scale = FALSE, highlight.outliers = TRUE) tmp.marginal.p <- marginal.pred.p.byyear.bysite[4, ] tmp.marginal.p [is.na(tmp.marginal.p)] <- 0.5 add.plot("2014", spRefuge,coordinates(d3SPDF), d3SPDF@data$n14, is.na(d3SPDF@data$mnc14), tmp.marginal.p, 1, c.ramp = c.ramp, log.scale = FALSE, highlight.outliers = TRUE) tmp.marginal.p <- marginal.pred.p.byyear.bysite[5, ] tmp.marginal.p [is.na(tmp.marginal.p)] <- 0.5 add.plot("2015", spRefuge,coordinates(d3SPDF), d3SPDF@data$n15, is.na(d3SPDF@data$mnc15), tmp.marginal.p, 1, c.ramp = c.ramp, log.scale = FALSE, highlight.outliers = TRUE) tmp.marginal.p <- marginal.pred.p.byyear.bysite[6, ] tmp.marginal.p [is.na(tmp.marginal.p)] <- 0.5 add.plot("2016", spRefuge,coordinates(d3SPDF), d3SPDF@data$n16, is.na(d3SPDF@data$mnc16), tmp.marginal.p, 1, c.ramp = c.ramp, log.scale = FALSE, highlight.outliers = TRUE) tmp.marginal.p <- marginal.pred.p.byyear.bysite[7, ] tmp.marginal.p [is.na(tmp.marginal.p)] <- 0.5 add.plot("2017", spRefuge,coordinates(d3SPDF), d3SPDF@data$n17, is.na(d3SPDF@data$mnc17), tmp.marginal.p, 1, c.ramp = c.ramp, log.scale = FALSE, highlight.outliers = TRUE) tmp.marginal.p <- marginal.pred.p.byyear.bysite[8, ] tmp.marginal.p [is.na(tmp.marginal.p)] <- 0.5 add.plot("2018", spRefuge,coordinates(d3SPDF), d3SPDF@data$n18, is.na(d3SPDF@data$mnc18), tmp.marginal.p, 1, c.ramp = c.ramp, log.scale = FALSE, highlight.outliers = TRUE) #Add legend plot(c(0, 1), c(0, 1), type = "n", xlab = "", ylab = "", bty = "n", xaxt = "n", yaxt = "n") n.p.levels <- 2 p.levels<-seq(0, 1, length = n.p.levels) legend.colors<-getcolor(p.levels,1, c.ramp = rev(c.ramp), log.scale = FALSE) rasterImage(as.raster(legend.colors), 0.05, 0.05, 0.2, 0.85) rect(0.05, 0.05, 0.2, 0.85) #Work out where the labels go legend.labels<-seq(0, 1, length = 5) for(i in 1:length(legend.labels)){ text(0.2, (0.85 - 0.05) * legend.labels[i] + 0.05, legend.labels[i], cex = 2.5, pos = 4) } text(0.18, 0.95, labels = 'p-value', cex = 2.5) dev.off()