In statguy/STREM: Winter track counts
#library(devtools); install_github("statguy/STREM")
library(parallel)
library(doMC)
registerDoMC(cores=detectCores())
library(STREM)
source("~/git/STREM/setup/WTC-Boot.R")
library(ggplot2)
library(scales)
library(grid)
library(plyr)
library(gridExtra)
context <- Context$new(resultDataDirectory=wd.data.results, processedDataDirectory=wd.data.processed, rawDataDirectory=wd.data.raw, scratchDirectory=wd.scratch, figuresDirectory=wd.figures, figuresHost=wd.figuresHost)
scenarios <- c("A","B","C","D","E","F")
allScenarios <- c(scenarios, paste0(scenarios, "combined"), paste0(scenarios, "10days"))
modelNames <- c("SmoothModel-nbinomial-matern-ar1", "SmoothModel-nbinomial-ar1", "SmoothModelMean-nbinomial-ar1", "FMPModel", "SmoothModelMean-nbinomial-ar1-priors1")
iterations <- as.integer(c(1,1,1,1,3,1))
#populationSizeOverEstimate <- 1e999
prettyScenario <- function(scenario) substr(scenario, 1, 1)
prettyModelName <- function(x) {
lookup <- data.frame(x=modelNames, y=c("ST", "VT", "T1", "FMP", "T2"), stringsAsFactors=F)
lookup$y[match(x, lookup$x)]
}
#studyArea <- FinlandStudyArea$new(context=context)$setup(tolerance=0.001)
boundaryDF <- ggplot2::fortify(FinlandWTCStudy(context=context, response="canis.lupus")$studyArea$boundary)
cutoff <- c(1,1000)
study <- getMSS(scenario=scenarios[1])$study
study$studyArea$boundary@polygons[[1]]@area / 1000^2 # Study area area
200 / (study$studyArea$boundary@polygons[[1]]@area / 1000^2) # Initial population density in Finland
intersections <- study$loadIntersections(iteration=as.integer(4))
ddply(intersections$intersections@data, .(response, year), function(x) sum(x$length) / 1000) # Total sampling effort per year (km)
#ddply(intersections$intersections@data, .(response, year), function(x) sum(x$intersections))
# Mean encounter rates
#aggrates <- data.frame()
aggrates <- ldply(allScenarios, function(scenario) {
study <- getMSS(scenario=scenario)$study
rates <- data.frame()
for (iteration in 1:50) {
intersections <- try(study$loadIntersections(iteration=as.integer(iteration)))
if (inherits(intersections, "try-error")) next
x <- ddply(intersections$intersections@data, .(response, year), function(x) data.frame(rate=sum(x$intersections) / (sum(x$length)/1000) / sum(x$duration)))
rates <- rbind(rates, x)
}
y <- ddply(rates, .(year), function(x, scenario) data.frame(scenario=scenario, rate.mean=mean(x$rate), rate.sd=sd(x$rate)), scenario=scenario)
#aggrates <- rbind(aggrates, y)
return(y)
}, .parallel=T)
subset(aggrates, year==2020)
###
library(raster)
mss <- getMSS("E")
study <- mss$study
surveyRoutes <- study$surveyRoutes
#habitatTypes <- raster::extract(habitat, surveyRoutes$surveyRoutes[1:3])
study$studyArea$readRasterIntoMemory()
habitat <- study$studyArea$habitat
habitatTypes <- raster::extract(habitat, surveyRoutes$surveyRoutes)
habitatWeights <- CORINEHabitatWeights$new()
forestProportions <- c()
for (i in seq_along(habitatTypes)) {
classifiedHabitatTypes <- habitatWeights$classify(habitatTypes[[i]]) == 3 # Forest
forestProportions <- c(forestProportions, sum(classifiedHabitatTypes) / length(classifiedHabitatTypes))
}
#mean(forestProportion)
sum(forestProportion * surveyRoutes$lengths) / sum(surveyRoutes$lengths) # = 0.9404947
###
study$studyArea$boundary@polygons[[1]]@area
habitatCut <- raster::mask(habitat, study$studyArea$boundary)
habitatWeights <- CORINEHabitatWeights$new(study=study)
habitatFreq <- freq(habitatCut, progress="text")
x <- habitatFreq
p <- c()
for (i in 1:50) {
habitatPreferences <- HabitatSelection$new(study=study, iteration=as.integer(i))
fileName <- habitatPreferences$getHabitatSelectionFileName()
habitatSelection <- try(habitatPreferences$loadHabitatSelection())
if (inherits(habitatSelection, "try-error")) next
habitatWeights$setHabitatSelectionWeights(habitatSelection)
y <- subset(habitatWeights$weights, habitat %in% x[,1])
effectiveArea <- sum(y$weight * x[,2]) * prod(res(habitat))
fullArea <- sum(x[-nrow(x),]) * prod(res(habitat))
p <- c(p, effectiveArea / fullArea)
}
mean(1/p) # = 1.286095
p <- c()
for (i in 1:50) {
estimates <- study$getModel(modelName=modelNames[4], iteration=as.integer(i))
popsize <- try(study$getPopulationSize(estimates=estimates, save=FALSE))
if (inherits(popsize, "try-error")) next
p <- c(p, popsize$sizeData$Estimated)
}
i <- 1
trackData <- data.frame()
surveyRoutesData <- data.frame()
for (scenario in scenarios) {
mss <- getMSS(scenario, nSurveyRoutes=500)
study <- mss$study
tracks <- study$loadTracks(iteration=iterations[i])
intersections <- tracks$countIntersections(surveyRoutes=study$surveyRoutes, save=F)
tracks$tracks <- subset(tracks$tracks, year==2001 & yday<59)
x <- tracks$toGGDF()
x$scenario <- scenario
x$observation <- F
intersections$tracks$tracks <- subset(intersections$tracks$tracks, year==2001)
y <- intersections$tracks$toGGDF()
y$scenario <- scenario
y$observation <- T
levels(y$group) <- paste(levels(y$group), "2")
trackData <- rbind(trackData, x, y)
x <- intersections$surveyRoutes$toGGDF()
x$scenario <- scenario
observedAt <- intersections$intersections$intersections[intersections$intersections$intersections$intersections > 0 & intersections$intersections$intersections$year == 2001,]$surveyRoute
y <- data.frame(id=unique(x$id), observation=FALSE)
if (length(observedAt) > 0) y[y$id %in% observedAt,]$observation <- TRUE
x <- plyr::join(x, y, by="id")
surveyRoutesData <- rbind(surveyRoutesData, x)
i <- i + 1
}
trackData$scenario <- factor(trackData$scenario)
#window <- c(3.375, 3.4, 6.9, 6.925)*1e6
window <- (c(3.425, 3.45, 6.95, 6.975)+c(-.01,.01,-.01,.01)*9)*1e6
habitatRaster <- crop(study$studyArea$habitat, extent(window))
z <- sampleRegular(habitatRaster, 10000, asRaster=TRUE)
habitatData <- as.data.frame(rasterToPoints(z))
colnames(habitatData) <- c("x","y","z")
habitatData$scenario <- factor("E", levels=scenarios)
habitatData2 <- habitatData
habitatData2$scenario <- factor("F", levels=scenarios)
habitatData <- rbind(habitatData, habitatData2)
habitatData$z <- factor(habitatData$z, levels=0:255)
col <- study$studyArea$habitat@legend@colortable
names(col) <- 0:255
habitatData$x <- habitatData$x / 1000
habitatData$y <- habitatData$y / 1000
trackData$long <- trackData$long / 1000
trackData$lat <- trackData$lat / 1000
surveyRoutesData$long <- surveyRoutesData$long / 1000
surveyRoutesData$lat <- surveyRoutesData$lat / 1000
window <- window / 1000
p <- ggplot(habitatData) + geom_raster(aes(x, y, fill=z), alpha=.4) + scale_fill_manual(values=col) +
geom_path(data=trackData, aes(long, lat, group=group, colour=observation), alpha=1, size=.3) +#, alpha=observation)) +
scale_colour_manual(values=c("gray55","black")) + # scale_alpha_manual(values=c(.5,1)) + # gray40
geom_path(data=surveyRoutesData, aes(long, lat, group=group), size=.3, colour="blue", alpha=.5) +
coord_fixed(ratio=1, xlim=window[1:2], ylim=window[3:4]) +
facet_wrap(~scenario) + #coord_equal()
theme_raster() + theme(panel.border=element_rect(fill="transparent")) +
xlab("X-coordinate (km)") + ylab("Y-coordinate (km)") +
theme(legend.position="none", panel.margin=unit(0, "lines"), panel.grid.minor=element_blank(), panel.grid.major=element_blank(), plot.margin=unit(c(0,0.1,0,0), "lines")) # top, right, bottom, and left
saveFigure(p, "manuscript-figure-1.png", width=17, height=15, dpi=200) # draft quality
saveFigure(p, "manuscript-figure-1.pdf", width=17, height=15, dpi=300)
# scale on the map in km
c((window[2]-window[1]), (window[4]-window[3]))/1000
getSurveyRouteLocations = function(scenarios, modelName, iterations, years) {
ddply(data.frame(scenario=scenarios, iteration=iterations, stringsAsFactors=FALSE), .(scenario, iteration), function(x, modelName, years) {
message("Scenario = ", x$scenario, ", iteration = ", x$iteration)
study <- getStudy(scenario=x$scenario)
estimates <- study$getModel(modelName=modelName, iteration=x$iteration)
estimates$modelName <- modelName
estimates <- study$loadEstimates(estimates)
z <- cbind(estimates$getUnscaledObservationCoordinates(), estimates$data)
surveyRoutes <- subset(z, year %in% years, select=c(x,y,intersections))
surveyRoutes$intersections <- surveyRoutes$intersections > 0
surveyRoutes$Scenario <- x$scenario
surveyRoutes$Iteration <- x$iteration
return(surveyRoutes)
}, modelName=modelName, years=years)
}
#x<-data.frame(scenario="C", iteration=as.integer(3), stringsAsFactors=F)
getAgentLocations = function(scenarios, modelName, iterations, years) {
ddply(data.frame(scenario=scenarios, iteration=iterations, stringsAsFactors=FALSE), .(scenario, iteration), function(x, modelName, years) {
study <- getStudy(scenario=x$scenario)
tracks <- study$loadTracks(iteration=x$iteration)
tracks$tracks <- subset(tracks$tracks, year %in% years & yday<59)
intersections <- study$loadIntersections(iteration=x$iteration)
# The observation day has not been recorded, ouch... just randomize a day from the 60 days track
observedBursts <- aaply(intersections$intersectionsMatrix, 2, function(x) if (any(x>0)) TRUE else FALSE)
#observationTracks <- tracks$randomizeObservationDayTracks(days=1)
observationTracks <- tracks
# Return mean location of the unobserved agents
a <- ddply(observationTracks$tracks, .(id, year, burst), function(x, observedBursts) {
name <- with(x[1,], paste(id, year, burst))
if (observedBursts[name] == FALSE) return(data.frame(x=mean(x$x), y=mean(x$y), year=x$year[1]))
#if (observedBursts[name] == FALSE) return(data.frame(x=x$x, y=x$y, year=x$year))
NULL
}, observedBursts=observedBursts)[,c("x","y","year")]
a$Observed <- F
# But for the observed agents we know the location, which is the transect
index <- intersections$intersections$intersections > 0 & intersections$intersections$year %in% years
b <- if (all(index == FALSE)) data.frame() else {
x <- intersections$intersections[index,]
y <- data.frame(coordinates(x), year=x$year)
y$Observed <- T
y
}
agentLocations <- rbind(a, b)
agentLocations$Scenario <- x$scenario
agentLocations$Model <- x$modelName
return(agentLocations)
}, modelName=modelName, years=years)
}
getEstimatedDensity = function(scenarios, modelName, iterations, year, sd=FALSE, weightHabitats=F) {
ddply(data.frame(scenario=scenarios, year=rep(year, length(scenarios)), iteration=iterations, stringsAsFactors=FALSE), .(scenario, iteration), function(x, modelName, sd) {
message("Scenario = ", x$scenario, ", iteration = ", x$iteration)
withHabitatWeights <- if (weightHabitats && substr(x$scenario, 1, 1) %in% c("E","F")) T else F
study <- getStudy(scenario=x$scenario, withHabitatWeights=withHabitatWeights)
estimates <- study$getModel(modelName=modelName, iteration=x$iteration)
estimates$modelName <- modelName
estimates <- study$loadEstimates(estimates)
estimates$collectEstimates()
estimates$data <- subset(estimates$data, year == x$year)
density <- estimates$getPopulationDensity(.parallel=T)
density <- if (withHabitatWeights) {
habitatWeights <- study$getHabitatWeights(iteration=x$iteration)
habitatWeightsRaster <- habitatWeights$getWeightsRaster(save=F)
density$weight(habitatWeightsRaster)
} else wdensity <- density$copy()
if (is.null(density) || maxValue(density$rasterStack[[1]]) > 1) stop("Estimation failed.")
density <- as.data.frame(density$rasterStack[[1]], xy=T)
density$Scenario <- x$scenario
density$Iteration <- x$iteration
density$Model <- modelName
density$Year <- x$year
return(density)
}, modelName=modelName, sd=sd)
}
study <- getStudy(scenario="A")
year <- 2010
yearLayer <- paste0("X", year)
#iterations <- as.integer(c(1,4,1,4,5,38))
#iterations <- as.integer(c(1,4,1,4,5,16))
iterations <- as.integer(c(1,4,1,4,5,21))
densities <- getEstimatedDensity(scenarios=scenarios, modelName=modelNames[1], iterations=iterations, year=year)
wdensities <- getEstimatedDensity(scenarios=scenarios, modelName=modelNames[1], iterations=iterations, year=year, weightHabitats=T)
surveyRoutes <- getSurveyRouteLocations(scenarios=scenarios, modelName=modelNames[1], iterations=iterations, years=year)
agents <- getAgentLocations(scenarios=scenarios, modelName=modelNames[1], iterations=iterations, years=year)
limits <- extent(study$studyArea$boundary)
# Scale density to 1/km^2
scaleDensity <- function(x, yearLayer, weight) {
x[,yearLayer] <- x[,yearLayer] * weight
return(x)
}
wscale <- 1 / (prod(res(study$getTemplateRaster())) / 1000^2)
densitiesScaled <- ddply(densities, .(scenario), scaleDensity, yearLayer=yearLayer, weight=wscale)
wdensitiesScaled <- ddply(wdensities, .(scenario), scaleDensity, yearLayer=yearLayer, weight=wscale)
gridTemplate <- raster(extent(study$studyArea$boundary), nrows=12, ncols=6, crs=study$studyArea$proj4string)
gridTemplate@extent@xmax <- gridTemplate@extent@xmin + dim(gridTemplate)[2] * 100 * 1e3
gridTemplate@extent@ymax <- gridTemplate@extent@ymin + dim(gridTemplate)[1] * 100 * 1e3
boundaryDF <- ggplot2::fortify(FinlandWTCStudy(context=context, response="canis.lupus")$studyArea$boundary)
wdensitiesScaled$x <- wdensitiesScaled$x / 1000 - 3000
wdensitiesScaled$y <- wdensitiesScaled$y / 1000 - 6000
densitiesScaled$x <- densitiesScaled$x / 1000 - 3000
densitiesScaled$y <- densitiesScaled$y / 1000 - 6000
boundaryDF$lat <- boundaryDF$lat / 1000 - 6000
boundaryDF$long <- boundaryDF$long / 1000 - 3000
agents$x <- agents$x / 1000 - 3000
agents$y <- agents$y / 1000 - 6000
limits@xmin <- limits@xmin / 1000 - 3000
limits@ymin <- limits@ymin / 1000 - 6000
limits@xmax <- limits@xmax / 1000 - 3000
limits@ymax <- limits@ymax / 1000 - 6000
surveyRoutes$x <- surveyRoutes$x / 1000 - 3000
surveyRoutes$y <- surveyRoutes$y / 1000 - 6000
gridTemplate@extent@xmin <- round(gridTemplate@extent@xmin / 1000) - 3000
gridTemplate@extent@ymin <- round(gridTemplate@extent@ymin / 1000) - 6000
gridTemplate@extent@xmax <- round(gridTemplate@extent@xmax / 1000) - 3000
gridTemplate@extent@ymax <- round(gridTemplate@extent@ymax / 1000) - 6000
p1 <- ggplot(wdensitiesScaled) +
geom_raster(aes_string(x="x", y="y", fill=yearLayer)) +
stat_contour(data=densitiesScaled, aes_string(x="x", y="y", z=yearLayer, colour=yearLayer), size=0.1) +
scale_fill_gradientn(colours=c(rep("#00A600FF", 60), terrain.colors(40)), na.value="transparent",
guide=guide_legend(title=expression(paste("Population density ", (individuals/km^2)))), trans="log10",
breaks=trans_breaks("log10", function(x) 10^x),
labels=trans_format("log10", math_format(10^.x))) +
geom_path(data=boundaryDF, aes(long, lat, group=group), colour="black", size=0.5) +
geom_point(data=agents, aes(x, y, colour=Observed), size=1.3, colour="white") +
geom_point(data=agents, aes(x, y, colour=Observed), size=1) + scale_color_manual(values=c("darkgray", "black"), guide=FALSE) +
scale_x_continuous(limits=c(limits@xmin, limits@xmax), expand=c(0.01,0.01), breaks=pretty_breaks(n=3)) +
scale_y_continuous(limits=c(limits@ymin, limits@ymax), expand=c(0.01,0.01)) +
facet_grid(~scenario) + coord_equal() + theme_raster() +
theme(legend.position="bottom", legend.margin=unit(-1, units="line"), legend.title=element_text(size=10), legend.text=element_text(size=10)) + #top, right, bottom, and left margins
theme(strip.background=element_blank(), panel.background=element_rect(fill="transparent", colour=NA), panel.margin=unit(0, "lines"), panel.grid.minor=element_blank(), panel.grid.major=element_blank()) +
xlab("X-coordinate (km)") + ylab("Y-coordinate (km)")
g1 <- arrangeGrob(p1, main=textGrob("a)", x=0.01, y=0.01, just="left", gp=gpar(fontsize=10)))
saveFigure(g1, "manuscript-figure-2a.png", width=17, height=8.07, dpi=200) # draft quality
saveFigure(g1, "manuscript-figure-2a.pdf", width=17, height=8.07, dpi=300)
xlabels <- seq(gridTemplate@extent@xmin, gridTemplate@extent@xmax, 1e5/1000)
xlabels[c(2,4,6,7)] <- ""
p2 <- ggplot(boundaryDF) + geom_path(aes(long, lat, group=group), colour="grey", size=0.5) +
geom_point(data=surveyRoutes, aes(x, y, colour=intersections), size=1) + scale_color_manual(values=c("darkgrey", "black")) +
scale_x_continuous(breaks=seq(gridTemplate@extent@xmin, gridTemplate@extent@xmax, 1e5/1000), expand=c(0.01,0.01), labels=xlabels) +
scale_y_continuous(breaks=seq(gridTemplate@extent@ymin, gridTemplate@extent@ymax, 1e5/1000), expand=c(0.01,0.01)) +
facet_grid(~scenario) + coord_equal() + theme_raster() +
theme(panel.grid.major=element_line(colour="grey")) + #top, right, bottom, and left margins
theme(strip.background=element_blank(), panel.background=element_rect(fill="transparent", colour=NA), panel.margin=unit(0, "lines"), panel.grid.minor=element_blank(), legend.position="none", axis.ticks=element_line(colour="grey", size=0.2)) +
xlab("X-coordinate (km)") + ylab("Y-coordinate (km)")
g2 <- arrangeGrob(p2, main=textGrob("b)", x=0.01, y=0.01, just="left", gp=gpar(fontsize=10)))
saveFigure(g2, "manuscript-figure-2b.png", width=17, height=7.2, dpi=200) # draft quality
saveFigure(g2, "manuscript-figure-2b.pdf", width=17, height=7.2, dpi=300)
p <- arrangeGrob(g1, g2, heights=c(8.07, 7.2)/(8.07+7.2))
saveFigure(p, "manuscript-figure-2.png", width=17, height=8.07+7.2, dpi=200) # draft quality
saveFigure(p, "manuscript-figure-2.pdf", width=17, height=8.07+7.2, dpi=300)
cutoffPopulationSizes <- function(populationSizes, cutoff=c(0,Inf), models) {
x <- subset(populationSizes, modelName %in% models & cutoff[1] <= Estimated & Estimated <= cutoff[2])
y <- subset(populationSizes, !(modelName %in% models))
rbind(x, y)
}
getIntersections <- function(scenarios) {
study <- getStudy(scenario="A")
validation <- Validation(study=study)
intersections <- validation$getValidationTemporalIntersections(scenarios=scenarios)
intersections$Year <- intersections$year
intersections
}
getPopulationSizes <- function(scenarios, modelNames, populationSizeCutoff=Inf) {
study <- getStudy(scenario="A") # Get study object, scenario doesn't matter
validation <- Validation(study=study, populationSizeCutoff=populationSizeCutoff)
populationSizes <- validation$getValidationTemporalPopulationSizes(scenarios=scenarios, modelNames=modelNames)
}
getMeanPopulationSizes <- function(populationSizes) {
library(plyr)
meanPopulationSizes <- ddply(populationSizes, .(scenario, modelName), function(x) {
x <- data.frame(True=mean(x$Observed), Estimated=mean(x$Estimated), Scenario=x$scenario[1])
if (is.infinite(x$True)) x$True <- NA
if (is.infinite(x$Estimated)) x$Estimated <- NA
return(x)
})
return(meanPopulationSizes)
}
getLogPopulationSizes <- function(populationSizes) {
logPopulationSizes <- populationSizes
logPopulationSizes$Observed <- log10(logPopulationSizes$Observed + 1)
logPopulationSizes$Estimated <- log10(logPopulationSizes$Estimated + 1)
return(logPopulationSizes)
}
populationSizes <- getPopulationSizes(scenarios=scenarios, modelNames=modelNames[c(1,3,4,5)])
meanPopulationSizes <- getMeanPopulationSizes(cutoffPopulationSizes(populationSizes, cutoff, models=modelNames[c(1,3,5)]))
logPopulationSizes <- getLogPopulationSizes(cutoffPopulationSizes(populationSizes, cutoff, models=modelNames[c(1,3,5)]))
logMeanPopulationSizes <- getMeanPopulationSizes(logPopulationSizes)
getPopulationSizesTable <- function(populationSizes, meanPopulationSizes, overEstimate, models, file, ...) {
library(plyr)
x <- ddply(populationSizes, .(Scenario, modelName), function(x, overEstimate, models) {
xrows <- nrow(x)
old <- x
x <- cutoffPopulationSizes(x, overEstimate, models=models)
if (nrow(x) == 0) {
print(old)
stop("ERROR")
}
if (nrow(x) == 0) return(NULL)
failed <- (1 - nrow(x) / xrows) * 100
y <- lm(Estimated~Observed, x)
#prop <- x$Estimated/x$Observed
prop <- (x$Estimated+1)/(x$Observed+1)
prop <- prop[!is.infinite(prop) & prop > 0]
prop.mean <- mean(log(prop))
prop.sd <- sd(log(prop))
prop.min <- quantile(log(prop), .025)
prop.max <- quantile(log(prop), .975)
slope.sd <- coef(summary(y))[2,2]
#print(prop)
#message("prop = ", mean(prop), " log mean = ", prop.mean)
df <- data.frame(Intercept=coef(y)[1], Slope=coef(y)[2], SlopeSD=slope.sd, R2=summary(y)$r.squared, EstTrue=prop.mean, EstTrueSD=prop.sd, EstTrueMin=prop.min, EstTrueMax=prop.max, Failed=failed, N=nrow(x))
#print(df)
return(df)
}, overEstimate=overEstimate, models=models, .inform=T)
y <- unique(meanPopulationSizes)
y$Scenario <- y$scenario
z <- plyr::join(x, y)
z$Model <- prettyModelName(z$modelName)
z$Scenario <- prettyScenario(z$Scenario)
return(z)
}
table1 <- getPopulationSizesTable(populationSizes, meanPopulationSizes, overEstimate=cutoff, models=modelNames[c(1,3,5)])
logTable1 <- getPopulationSizesTable(logPopulationSizes, logMeanPopulationSizes, overEstimate=log10(cutoff+1), models=modelNames[c(1,3,5)])
logTable1$Failed <- table1$Failed
#subset(populationSizes, Scenario=="C" & modelName=="FMPModel")
#subset(populationSizes, Scenario=="B" & modelName=="FMPModel")
writeTable <- function(df, filename, context) {
library(ReporteRs)
doc <- docx()
flextable <- FlexTable(df)
addFlexTable(doc, flextable)
dir <- normalizePath(context$figuresDirectory)
writeDoc(doc, file=file.path(dir, filename))
}
facet_wrap_labeller <- function(gg.plot,labels=NULL) {
#works with R 3.0.1 and ggplot2 0.9.3.1
require(gridExtra)
g <- ggplotGrob(gg.plot)
gg <- g$grobs
strips <- grep("strip_t", names(gg))
for(ii in seq_along(labels)) {
modgrob <- getGrob(gg[[strips[ii]]], "strip.text",
grep=TRUE, global=TRUE)
gg[[strips[ii]]]$children[[modgrob$name]] <- editGrob(modgrob,label=labels[ii])
}
g$grobs <- gg
class(g) = c("arrange", "ggplot",class(g))
g
}
plotPopulationSizesTable <- function(logTable) {
library(reshape2)
library(dplyr)
library(gtable)
library(gridExtra)
x <- melt(logTable, id.vars=c("Scenario", "Model"), measure.vars=c("EstTrue", "Slope", "R2", "Failed"))
bestFit <- data.frame(variable=unique(x$variable), y=c(0, 1, 1, 0))
errorBars <- x
errorBars$value <- NA
errorBars$value[errorBars$variable=="EstTrue"] <- logTable$EstTrue
errorBars$qmin[errorBars$variable=="EstTrue"] <- logTable$EstTrueMin
errorBars$qmax[errorBars$variable=="EstTrue"] <- logTable$EstTrueMax
margin <- unit(c(0.1, 0.1, -.75, 0), units="line")
p1 <- x %>% filter(variable=="EstTrue") %>%
ggplot(aes(Scenario, value, group=Model)) +
geom_hline(data=data.frame(y=0), aes(yintercept=y), color="darkgray") +
geom_errorbar(data=errorBars, aes(Scenario, group=Model, ymax=qmax, ymin=qmin), position=position_dodge(width=.7), colour="darkgray") +
geom_point(aes(shape=Model), position=position_dodge(width=.7)) +
xlab("") + ylab(expression(alpha)) +
theme_ms() + theme(panel.border=element_rect(fill="transparent"),
legend.title=element_text(size=10), legend.position="none",
strip.text.x=element_text(size=10), strip.text.y=element_text(size=10),
plot.margin=margin) #top, right, bottom, and left margins
g1 <- arrangeGrob(p1, main=textGrob("a)", x=0.01, y=0.01, just="left", gp=gpar(fontsize=10)))
p2 <- x %>% filter(variable=="Slope") %>%
ggplot(aes(Scenario, value, group=Model)) +
geom_hline(data=data.frame(y=1), aes(yintercept=y), color="darkgray") +
geom_point(aes(shape=Model), position=position_dodge(width=.7)) +
xlab("") + ylab(expression(beta)) +
theme_ms() + theme(panel.border=element_rect(fill="transparent"),
legend.title=element_text(size=10), legend.position="none",
strip.text.x=element_text(size=10), strip.text.y=element_text(size=10),
plot.margin=margin)
g2 <- arrangeGrob(p2, main=textGrob("b)", x=0.01, y=0.01, just="left", gp=gpar(fontsize=10)))
p3 <- x %>% filter(variable=="R2") %>%
ggplot(aes(Scenario, value, group=Model)) +
geom_hline(data=data.frame(y=1), aes(yintercept=y), color="darkgray") +
geom_point(aes(shape=Model), position=position_dodge(width=.7)) +
xlab("") + ylab(expression(R^2)) +
theme_ms() + theme(panel.border=element_rect(fill="transparent"),
legend.title=element_text(size=10), legend.position="none",
strip.text.x=element_text(size=10), strip.text.y=element_text(size=10),
plot.margin=margin)
g3 <- arrangeGrob(p3, main=textGrob("c)", x=0.01, y=0.01, just="left", gp=gpar(fontsize=10)))
p4 <- x %>% filter(variable=="Failed") %>%
ggplot(aes(Scenario, value, group=Model)) +
geom_hline(data=data.frame(y=0), aes(yintercept=y), color="darkgray") +
geom_point(aes(shape=Model), position=position_dodge(width=.7)) +
xlab("") + ylab("Failed (%)") + ylim(0,100) +
theme_ms() + theme(panel.border=element_rect(fill="transparent"),
legend.title=element_text(size=10), legend.position="none",
strip.text.x=element_text(size=10), strip.text.y=element_text(size=10),
plot.margin=margin)
g4 <- arrangeGrob(p4, main=textGrob("d)", x=0.01, y=0.01, just="left", gp=gpar(fontsize=10)))
legend <- gtable_filter(ggplot_gtable(ggplot_build(p1 + theme(legend.position="bottom", legend.margin=unit(-1, units="line")))), "guide-box")
lheight <- sum(legend$height)
theight <- unit(10, "points")
#p <- arrangeGrob(g1, g2, g3, g4, ncol=2)
p1 <- arrangeGrob(g1, g3, nrow=2)
p2 <- arrangeGrob(g2, g4, nrow=2)
p1 <- arrangeGrob(p1, textGrob("Scenario", x=0.58, gp=gpar(fontsize=10)), heights=unit.c(unit(1, "npc") - theight, theight))
p2 <- arrangeGrob(p2, textGrob("Scenario", x=0.57, gp=gpar(fontsize=10)), heights=unit.c(unit(1, "npc") - theight, theight))
p <- arrangeGrob(p1, p2, ncol=2)
#p <- arrangeGrob(p, textGrob("Scenario", gp=gpar(fontsize=10)), heights=unit.c(unit(1, "npc") - theight, theight))
p <- arrangeGrob(p, legend, heights=unit.c(unit(1, "npc") - lheight, lheight), ncol=1)
return(p)
}
logTable <- logTable1
p <- plotPopulationSizesTable(logTable)
saveFigure(p, "manuscript-figure-3.png", width=17, dpi=200) # draft quality
saveFigure(p, "manuscript-figure-3.pdf", width=17, dpi=300)
logTable <- logTable1
logTablePretty <- logTable
logTablePretty$Slope <- round(logTablePretty$Slope, 2)
logTablePretty$R2 <- round(logTablePretty$R2, 2)
logTablePretty$EstTrue <- round(logTablePretty$EstTrue, 2)
logTablePretty$EstTrueSD <- round(logTablePretty$EstTrueSD, 2)
logTablePretty$Failed <- round(logTablePretty$Failed, 0)
logTablePretty <- logTablePretty[,c("Scenario","Model","Slope","R2","EstTrue","EstTrueSD","Failed")]
writeTable(subset(logTablePretty, Model != "VT"), "table-1.docx", context)
logPopulationSizesSample <- subset(logPopulationSizes, scenario == "A")
logPopulationSizesSample <- cutoffPopulationSizes(logPopulationSizesSample, cutoff, models=modelNames[c(1,3,5)])
logPopulationSizesSample$modelName <- prettyModelName(logPopulationSizesSample$modelName)
logPopulationSizesSample <- rbind(
subset(logPopulationSizesSample, modelName=="ST"),
subset(logPopulationSizesSample, modelName=="T1"),
subset(logPopulationSizesSample, modelName=="T2"),
subset(logPopulationSizesSample, modelName=="FMP"))
lm(Estimated~Observed, subset(logPopulationSizes, modelName=="FMPModel" & Scenario=="A"))
lm(Estimated~Observed, subset(logPopulationSizesSample, modelName=="FMP" & Scenario=="A"))
plotLogPopulationSizeSample <- function(size, modelName, case, xlim, ylim) {
library(ggplot2)
library(scales)
size$Scenario <- prettyScenario(size$Scenario)
meanSize <- getMeanPopulationSizes(size)
meanSize$Scenario <- prettyScenario(meanSize$Scenario)
p <- ggplot(size, aes(Observed, Estimated, colour=modelName)) +
geom_point(size=1, colour="gray40", alpha=.5) +
geom_abline(colour="black", linetype="longdash") +
geom_point(data=meanSize, aes(True, Estimated), size=2, shape=15, colour="black")
#if (!any(abs(size$Estimated) == Inf)) p <- p + geom_smooth(method="lm", se=F, fullrange=T, colour="black")
ls <- ddply(size, .(modelName), function(x) if (!any(abs(x$Estimated) == Inf)) coef(lm(Estimated~Observed, x)) else c(NA,NA))
print(ls)
p <- p + geom_abline(data=ls, aes(intercept=`(Intercept)`, slope=Observed))
p <- p + facet_wrap(~modelName) +
scale_x_continuous(limits=xlim, labels=trans_format("identity", math_format(10^.x))) +
scale_y_continuous(limits=ylim, labels=trans_format("identity", math_format(10^.x))) +
xlab("True population size") + ylab("Estimated population size") +
theme_ms2() + theme(panel.border=element_rect(fill="transparent"),
plot.title=element_text(size=10),
legend.margin=unit(0, units="line"))
return(p)
}
p <- plotLogPopulationSizeSample(subset(logPopulationSizesSample, modelName %in% c("ST","T1","T2","FMP")), modelName=c("ST","T1","T2","FMP"), "a", xlim=log10(c(1, 10^3)), ylim=log10(c(1,10^3)))
saveFigure(p, "manuscript-figure-3sample.png", height=7.87, width=17) # draft quality
saveFigure(p, "manuscript-figure-3sample.pdf", height=7.87, width=17, dpi=300)
plotLogPopulationSize <- function(size, modelName, case, xlim, ylim) {
library(ggplot2)
library(scales)
library(grid)
size$Scenario <- prettyScenario(size$Scenario)
meanSize <- getMeanPopulationSizes(size)
meanSize$Scenario <- prettyScenario(meanSize$Scenario)
p <- ggplot(size, aes(Observed, Estimated, colour=modelName)) +
geom_point(size=1, colour="gray40", alpha=.5) +
geom_abline(colour="black", linetype="longdash") +
geom_point(data=meanSize, aes(True, Estimated), size=2, shape=15, colour="black")
ls <- ddply(size, .(Scenario), function(x) if (!any(abs(x$Estimated) == Inf)) coef(lm(Estimated~Observed, x)) else c(NA,NA))
p <- p + geom_abline(data=ls, aes(intercept=`(Intercept)`, slope=Observed))
p <- p + facet_wrap(~Scenario) +
scale_x_continuous(limits=xlim, labels=trans_format("identity", math_format(.x))) +
scale_y_continuous(limits=ylim, labels=trans_format("identity", math_format(.x))) +
xlab(expression(paste(Log[10], " true population size"))) + ylab(expression(paste(Log[10], " estimated population size"))) + ggtitle(paste("Model:", prettyModelName(modelName[1]))) +
theme_ms2() + theme(panel.border=element_rect(fill="transparent"),
legend.margin=unit(0, units="line"),
panel.margin=unit(0, units="line"))
return(p)
}
#modelNames <- modelNames[c(2:4,1)]
#xlim <- log10(c(1,1200)); ylim <- xlim
xlim <- log10(c(1,2000)); ylim <- xlim
p <- plotLogPopulationSize(subset(logPopulationSizes, modelName==modelNames[1]), modelName=modelNames[1], "1", xlim=xlim, ylim=ylim)
saveFigure(p, "manuscript-figure-3a1.png", width=17, dpi=200)
saveFigure(p, "manuscript-figure-3a1.pdf", width=17, dpi=300)
p <- plotLogPopulationSize(subset(logPopulationSizes, modelName==modelNames[3]), modelName=modelNames[3], "1", xlim=xlim, ylim=ylim)
saveFigure(p, "manuscript-figure-3a2.png", width=17, dpi=200)
saveFigure(p, "manuscript-figure-3a2.pdf", width=17, dpi=300)
p <- plotLogPopulationSize(subset(logPopulationSizes, modelName==modelNames[4]), modelName=modelNames[4], "1", xlim=xlim, ylim=ylim)
saveFigure(p, "manuscript-figure-3a3.png", width=17, dpi=200)
saveFigure(p, "manuscript-figure-3a3.pdf", width=17, dpi=300)
p <- plotLogPopulationSize(subset(logPopulationSizes, modelName==modelNames[5]), modelName=modelNames[5], "1", xlim=xlim, ylim=ylim)
saveFigure(p, "manuscript-figure-3a4.png", width=17, dpi=200)
saveFigure(p, "manuscript-figure-3a4.pdf", width=17, dpi=300)
ddply(logTable, .(modelName, Case, Scenario), function(x) {
data.frame(p=mean(exp(x$EstTrue)))
})
ddply(logTable, .(modelName, Case), function(x) {
data.frame(p=mean(exp(x$EstTrue)))
})
ddply(logTable, .(modelName, Scenario), function(x) {
data.frame(p=mean(exp(x$EstTrue)))
})
p <- ddply(subset(logTable, Scenario %in% c("A","E")), .(modelName, Case), function(x) {
p <- 1 - exp(subset(x, Scenario %in% "A")$EstTrue) / exp(subset(x, Scenario %in% "E")$EstTrue)
data.frame(p=p)
})
mean(p$p)
p <- ddply(subset(logTable, Scenario %in% c("A","F")), .(modelName, Case), function(x) {
p <- 1 - exp(subset(x, Scenario %in% "A")$EstTrue) / exp(subset(x, Scenario %in% "F")$EstTrue)
data.frame(p=p)
})
plotIntersectionPopulationSizeTable <- function(ratioTable) {
library(reshape2)
library(grid)
x <- melt(ratioTable, id.vars=c("Scenario", "Case"), measure.vars=c("SizeCount"))
errorBars <- x
errorBars$value <- NA
errorBars$value[errorBars$variable=="SizeCount"] <- ratioTable$SizeCount
errorBars$qmin[errorBars$variable=="SizeCount"] <- ratioTable$SizeCountMin
errorBars$qmax[errorBars$variable=="SizeCount"] <- ratioTable$SizeCountMax
p <- ggplot(x, aes(Scenario, value, group=Case)) +
geom_errorbar(data=errorBars, aes(Scenario, group=Case, ymax=qmax, ymin=qmin), position=position_dodge(width=.7), colour="darkgray") +
geom_point(aes(shape=Case), position=position_dodge(width=0.7)) +
ylab("Ratio") + #guides(shape=guide_legend()) +
theme_minimal(base_size=10) + theme(panel.border=element_rect(fill="transparent"), legend.title=element_text(size=10),
strip.text.y=element_text(size=10), legend.position="right",
#legend.background=element_rect(fill="transparent"),
plot.margin=unit(c(0.2, 0, 0, 0), units="line"),
legend.margin=unit(-1, units="line"),
panel.margin=unit(0, units="line"))# , legend.key=element_rect(fill="transparent"))
return(p)
}
intersections <- getIntersections(scenarios)
intersections2000 <- getIntersections(scenarios=paste0(scenarios, "combined"))
intersections10days <- getIntersections(scenarios=paste0(scenarios, "10days"))
intersections10days$intersections <- intersections10days$intersections / 10
intersectionsx <- plyr::join(populationSizes, intersections)
intersections2000x <- plyr::join(populationSizes2000, intersections2000)
intersections10daysx <- plyr::join(populationSizes10Days, intersections10days)
rer.stat <- function(x) data.frame(SizeCount=mean(x$intersections/x$Observed),
SizeCountMin=quantile(x$intersections/x$Observed, 0.025),
SizeCountMax=quantile(x$intersections/x$Observed, 0.975))
i <- ddply(intersectionsx, .(Scenario), rer.stat); i$Scenario=prettyScenario(i$Scenario); i$Case <- "1"
i2000 <- ddply(intersections2000x, .(Scenario), rer.stat); i2000$Scenario=prettyScenario(i2000$Scenario); i2000$Case <- "2"
i10days <- ddply(intersections10daysx, .(Scenario), rer.stat); i10days$Scenario=prettyScenario(i10days$Scenario); i10days$Case <- "3"
iall <- rbind(i, i2000, i10days)
iall
p <- plotIntersectionPopulationSizeTable(iall)
saveFigure(p, "manuscript-figure-4x.png", height=7.87/2) # draft quality
saveFigure(p, "manuscript-figure-4x.pdf", height=7.87/2, dpi=300)
#writeTable(i, "table-2a.docx", context)
#writeTable(i2000, "table-2b.docx", context)
#writeTable(i10days, "table-2c.docx", context)
getPopulationCorrelations <- function(scenarios, modelNames, cutoff=c(-Inf,Inf)) {
study <- getStudy(scenario="A")
validation <- Validation(study=study, populationSizeCutoff=cutoff)
populationCorrelations <- validation$getValidationSpatialPopulationSizes(scenarios=scenarios, modelNames=modelNames)
return(populationCorrelations)
}
plotPopulationCorrelations <- function(populationCorrelations, case) {
populationCorrelations$Scenario <- prettyScenario(populationCorrelations$Scenario)
p <- ggplot(populationCorrelations, aes(True, Correlation)) +
geom_hline(data=data.frame(y=0), aes(yintercept=y), color="darkgray") +
geom_point(size=1, colour="black") +
facet_wrap(~Scenario, ncol=2) + xlab("True population size") +
scale_y_continuous(limits=c(-1,1)) +
#ggtitle(paste("Case:", case)) +
theme_ms() + theme(panel.border=element_rect(fill="transparent"),
panel.margin=unit(0, units="line"))#, axis.text.x=element_text(angle=45, hjust=.5, vjust=1))
return(p)
}
pc1a <- getPopulationCorrelations(scenarios="A", modelNames=modelNames[1], cutoff=cutoff)
pc1b <- getPopulationCorrelations(scenarios="B", modelNames=modelNames[1], cutoff=cutoff)
pc1c <- getPopulationCorrelations(scenarios="C", modelNames=modelNames[1], cutoff=cutoff)
pc1d <- getPopulationCorrelations(scenarios="D", modelNames=modelNames[1], cutoff=cutoff)
pc1e <- getPopulationCorrelations(scenarios="E", modelNames=modelNames[1], cutoff=cutoff)
pc1f <- getPopulationCorrelations(scenarios="F", modelNames=modelNames[1], cutoff=cutoff)
pc1 <- rbind(pc1a,pc1b,pc1c,pc1d,pc1e,pc1f)
p1 <- plotPopulationCorrelations(pc1, "Case 1")
saveFigure(p1, "manuscript-figure-4a.png")
saveFigure(p1, "manuscript-figure-4a.pdf", dpi=300)
getCIProportions <- function(scenarios, modelNames, panel, scenarioNames, cutoff=Inf) {
study <- getStudy(scenario="A")
validation <- Validation(study=study, populationSizeCutoff=cutoff)
credibilityIntervals <- list()
if ("FMPModel" %in% modelNames) {
for (scenario in scenarios) {
study <- getStudy(scenario=scenario)
validation <- Validation(study=study, populationSizeCutoff=cutoff)
iterations <- validation$getCredibilityIntervalsValidationIterations("FMPModel")
p95 <- c()
p50 <- c()
for (i in iterations) {
load(file=validation$getCredibilityIntervalsValidationFileName(modelName="FMPModel", iteration=i))
p95 <- c(p95, bootCI$p95)
p50 <- c(p50, bootCI$p50)
}
a <- data.frame(scenario=scenario, Proportion=mean(p95, na.rm=T), n=NA, True=NA, Estimated=NA, Estimated.q1=NA, Estimated.q2=NA, modelName="FMPModel", Interval="95%")
b <- data.frame(scenario=scenario, Proportion=mean(p50, na.rm=T), n=NA, True=NA, Estimated=NA, Estimated.q1=NA, Estimated.q2=NA, modelName="FMPModel", Interval="50%")
credibilityIntervals <- rbind(credibilityIntervals, rbind(a, b))
}
}
a <- validation$getValidatedCredibilityIntervalsProportions(scenarios=scenarios, modelNames=modelNames[!modelNames %in% "FMPModel"])
b <- validation$getValidatedCredibilityIntervalsProportions(scenarios=scenarios, modelNames=modelNames[!modelNames %in% "FMPModel"], probs=c(.25, .75), probsName="50%")
credibilityIntervals <- rbind(credibilityIntervals, rbind(a, b))
credibilityIntervals$panel <- panel
return(credibilityIntervals)
}
cip1 <- getCIProportions(scenarios=scenarios, modelNames=modelNames[c(1,3,4,5)], panel="Case 1", scenarioNames=scenarios, cutoff=cutoff)
plotCIProportions <- function(credibilityIntervals) {
credibilityIntervals$Model <- prettyModelName(credibilityIntervals$modelName)
credibilityIntervals$CrI <- credibilityIntervals$Interval
p <- ggplot(credibilityIntervals, aes(scenario, Proportion, shape=Model, colour=CrI)) +
geom_hline(aes(yintercept=.5), color="darkgrey") +
geom_hline(aes(yintercept=.95), color="darkgrey") +
geom_point(aes(shape=Model, group=Model), position=position_dodge(width=.7)) +
scale_colour_manual("CrI / CI", values=c("black", "grey"), guide=guide_legend(override.aes=list(size=3.5, shape=15))) +
theme_ms() +
theme(panel.border=element_rect(fill="transparent")) +
xlab("Scenario") + ylab("Proportion") + ylim(0, 1) +
theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(),
panel.margin=unit(0, units="line"),
legend.direction="vertical",
legend.box="horizontal",
legend.key.height=unit(0.7, "line"))
#guides(scale_colour_manual=guide_colourbar())
library(gridExtra)
g_legend <- function(a.gplot) {
tmp <- ggplot_gtable(ggplot_build(a.gplot))
leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
legend <- tmp$grobs[[leg]]
legend
}
legend <- g_legend(p)
p <- arrangeGrob(p + theme(legend.position="none"), legend, heights=c(6/10,4/10))
return(p)
}
cip <- cip1
p <- plotCIProportions(cip)
saveFigure(p, "manuscript-figure-5.png", height=7.87*2/3)
saveFigure(p, "manuscript-figure-5.pdf", height=7.87*2/3, dpi=300)
statguy/STREM documentation built on May 30, 2019, 9:43 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
#library(devtools); install_github("statguy/STREM") library(parallel) library(doMC) registerDoMC(cores=detectCores()) library(STREM) source("~/git/STREM/setup/WTC-Boot.R") library(ggplot2) library(scales) library(grid) library(plyr) library(gridExtra) context <- Context$new(resultDataDirectory=wd.data.results, processedDataDirectory=wd.data.processed, rawDataDirectory=wd.data.raw, scratchDirectory=wd.scratch, figuresDirectory=wd.figures, figuresHost=wd.figuresHost) scenarios <- c("A","B","C","D","E","F") allScenarios <- c(scenarios, paste0(scenarios, "combined"), paste0(scenarios, "10days")) modelNames <- c("SmoothModel-nbinomial-matern-ar1", "SmoothModel-nbinomial-ar1", "SmoothModelMean-nbinomial-ar1", "FMPModel", "SmoothModelMean-nbinomial-ar1-priors1") iterations <- as.integer(c(1,1,1,1,3,1)) #populationSizeOverEstimate <- 1e999 prettyScenario <- function(scenario) substr(scenario, 1, 1) prettyModelName <- function(x) { lookup <- data.frame(x=modelNames, y=c("ST", "VT", "T1", "FMP", "T2"), stringsAsFactors=F) lookup$y[match(x, lookup$x)] } #studyArea <- FinlandStudyArea$new(context=context)$setup(tolerance=0.001) boundaryDF <- ggplot2::fortify(FinlandWTCStudy(context=context, response="canis.lupus")$studyArea$boundary) cutoff <- c(1,1000)
study <- getMSS(scenario=scenarios[1])$study study$studyArea$boundary@polygons[[1]]@area / 1000^2 # Study area area 200 / (study$studyArea$boundary@polygons[[1]]@area / 1000^2) # Initial population density in Finland intersections <- study$loadIntersections(iteration=as.integer(4)) ddply(intersections$intersections@data, .(response, year), function(x) sum(x$length) / 1000) # Total sampling effort per year (km) #ddply(intersections$intersections@data, .(response, year), function(x) sum(x$intersections)) # Mean encounter rates #aggrates <- data.frame() aggrates <- ldply(allScenarios, function(scenario) { study <- getMSS(scenario=scenario)$study rates <- data.frame() for (iteration in 1:50) { intersections <- try(study$loadIntersections(iteration=as.integer(iteration))) if (inherits(intersections, "try-error")) next x <- ddply(intersections$intersections@data, .(response, year), function(x) data.frame(rate=sum(x$intersections) / (sum(x$length)/1000) / sum(x$duration))) rates <- rbind(rates, x) } y <- ddply(rates, .(year), function(x, scenario) data.frame(scenario=scenario, rate.mean=mean(x$rate), rate.sd=sd(x$rate)), scenario=scenario) #aggrates <- rbind(aggrates, y) return(y) }, .parallel=T) subset(aggrates, year==2020) ### library(raster) mss <- getMSS("E") study <- mss$study surveyRoutes <- study$surveyRoutes #habitatTypes <- raster::extract(habitat, surveyRoutes$surveyRoutes[1:3]) study$studyArea$readRasterIntoMemory() habitat <- study$studyArea$habitat habitatTypes <- raster::extract(habitat, surveyRoutes$surveyRoutes) habitatWeights <- CORINEHabitatWeights$new() forestProportions <- c() for (i in seq_along(habitatTypes)) { classifiedHabitatTypes <- habitatWeights$classify(habitatTypes[[i]]) == 3 # Forest forestProportions <- c(forestProportions, sum(classifiedHabitatTypes) / length(classifiedHabitatTypes)) } #mean(forestProportion) sum(forestProportion * surveyRoutes$lengths) / sum(surveyRoutes$lengths) # = 0.9404947 ### study$studyArea$boundary@polygons[[1]]@area habitatCut <- raster::mask(habitat, study$studyArea$boundary) habitatWeights <- CORINEHabitatWeights$new(study=study) habitatFreq <- freq(habitatCut, progress="text") x <- habitatFreq p <- c() for (i in 1:50) { habitatPreferences <- HabitatSelection$new(study=study, iteration=as.integer(i)) fileName <- habitatPreferences$getHabitatSelectionFileName() habitatSelection <- try(habitatPreferences$loadHabitatSelection()) if (inherits(habitatSelection, "try-error")) next habitatWeights$setHabitatSelectionWeights(habitatSelection) y <- subset(habitatWeights$weights, habitat %in% x[,1]) effectiveArea <- sum(y$weight * x[,2]) * prod(res(habitat)) fullArea <- sum(x[-nrow(x),]) * prod(res(habitat)) p <- c(p, effectiveArea / fullArea) } mean(1/p) # = 1.286095 p <- c() for (i in 1:50) { estimates <- study$getModel(modelName=modelNames[4], iteration=as.integer(i)) popsize <- try(study$getPopulationSize(estimates=estimates, save=FALSE)) if (inherits(popsize, "try-error")) next p <- c(p, popsize$sizeData$Estimated) }
i <- 1 trackData <- data.frame() surveyRoutesData <- data.frame() for (scenario in scenarios) { mss <- getMSS(scenario, nSurveyRoutes=500) study <- mss$study tracks <- study$loadTracks(iteration=iterations[i]) intersections <- tracks$countIntersections(surveyRoutes=study$surveyRoutes, save=F) tracks$tracks <- subset(tracks$tracks, year==2001 & yday<59) x <- tracks$toGGDF() x$scenario <- scenario x$observation <- F intersections$tracks$tracks <- subset(intersections$tracks$tracks, year==2001) y <- intersections$tracks$toGGDF() y$scenario <- scenario y$observation <- T levels(y$group) <- paste(levels(y$group), "2") trackData <- rbind(trackData, x, y) x <- intersections$surveyRoutes$toGGDF() x$scenario <- scenario observedAt <- intersections$intersections$intersections[intersections$intersections$intersections$intersections > 0 & intersections$intersections$intersections$year == 2001,]$surveyRoute y <- data.frame(id=unique(x$id), observation=FALSE) if (length(observedAt) > 0) y[y$id %in% observedAt,]$observation <- TRUE x <- plyr::join(x, y, by="id") surveyRoutesData <- rbind(surveyRoutesData, x) i <- i + 1 } trackData$scenario <- factor(trackData$scenario) #window <- c(3.375, 3.4, 6.9, 6.925)*1e6 window <- (c(3.425, 3.45, 6.95, 6.975)+c(-.01,.01,-.01,.01)*9)*1e6 habitatRaster <- crop(study$studyArea$habitat, extent(window)) z <- sampleRegular(habitatRaster, 10000, asRaster=TRUE) habitatData <- as.data.frame(rasterToPoints(z)) colnames(habitatData) <- c("x","y","z") habitatData$scenario <- factor("E", levels=scenarios) habitatData2 <- habitatData habitatData2$scenario <- factor("F", levels=scenarios) habitatData <- rbind(habitatData, habitatData2) habitatData$z <- factor(habitatData$z, levels=0:255) col <- study$studyArea$habitat@legend@colortable names(col) <- 0:255 habitatData$x <- habitatData$x / 1000 habitatData$y <- habitatData$y / 1000 trackData$long <- trackData$long / 1000 trackData$lat <- trackData$lat / 1000 surveyRoutesData$long <- surveyRoutesData$long / 1000 surveyRoutesData$lat <- surveyRoutesData$lat / 1000 window <- window / 1000 p <- ggplot(habitatData) + geom_raster(aes(x, y, fill=z), alpha=.4) + scale_fill_manual(values=col) + geom_path(data=trackData, aes(long, lat, group=group, colour=observation), alpha=1, size=.3) +#, alpha=observation)) + scale_colour_manual(values=c("gray55","black")) + # scale_alpha_manual(values=c(.5,1)) + # gray40 geom_path(data=surveyRoutesData, aes(long, lat, group=group), size=.3, colour="blue", alpha=.5) + coord_fixed(ratio=1, xlim=window[1:2], ylim=window[3:4]) + facet_wrap(~scenario) + #coord_equal() theme_raster() + theme(panel.border=element_rect(fill="transparent")) + xlab("X-coordinate (km)") + ylab("Y-coordinate (km)") + theme(legend.position="none", panel.margin=unit(0, "lines"), panel.grid.minor=element_blank(), panel.grid.major=element_blank(), plot.margin=unit(c(0,0.1,0,0), "lines")) # top, right, bottom, and left saveFigure(p, "manuscript-figure-1.png", width=17, height=15, dpi=200) # draft quality saveFigure(p, "manuscript-figure-1.pdf", width=17, height=15, dpi=300) # scale on the map in km c((window[2]-window[1]), (window[4]-window[3]))/1000
getSurveyRouteLocations = function(scenarios, modelName, iterations, years) { ddply(data.frame(scenario=scenarios, iteration=iterations, stringsAsFactors=FALSE), .(scenario, iteration), function(x, modelName, years) { message("Scenario = ", x$scenario, ", iteration = ", x$iteration) study <- getStudy(scenario=x$scenario) estimates <- study$getModel(modelName=modelName, iteration=x$iteration) estimates$modelName <- modelName estimates <- study$loadEstimates(estimates) z <- cbind(estimates$getUnscaledObservationCoordinates(), estimates$data) surveyRoutes <- subset(z, year %in% years, select=c(x,y,intersections)) surveyRoutes$intersections <- surveyRoutes$intersections > 0 surveyRoutes$Scenario <- x$scenario surveyRoutes$Iteration <- x$iteration return(surveyRoutes) }, modelName=modelName, years=years) } #x<-data.frame(scenario="C", iteration=as.integer(3), stringsAsFactors=F) getAgentLocations = function(scenarios, modelName, iterations, years) { ddply(data.frame(scenario=scenarios, iteration=iterations, stringsAsFactors=FALSE), .(scenario, iteration), function(x, modelName, years) { study <- getStudy(scenario=x$scenario) tracks <- study$loadTracks(iteration=x$iteration) tracks$tracks <- subset(tracks$tracks, year %in% years & yday<59) intersections <- study$loadIntersections(iteration=x$iteration) # The observation day has not been recorded, ouch... just randomize a day from the 60 days track observedBursts <- aaply(intersections$intersectionsMatrix, 2, function(x) if (any(x>0)) TRUE else FALSE) #observationTracks <- tracks$randomizeObservationDayTracks(days=1) observationTracks <- tracks # Return mean location of the unobserved agents a <- ddply(observationTracks$tracks, .(id, year, burst), function(x, observedBursts) { name <- with(x[1,], paste(id, year, burst)) if (observedBursts[name] == FALSE) return(data.frame(x=mean(x$x), y=mean(x$y), year=x$year[1])) #if (observedBursts[name] == FALSE) return(data.frame(x=x$x, y=x$y, year=x$year)) NULL }, observedBursts=observedBursts)[,c("x","y","year")] a$Observed <- F # But for the observed agents we know the location, which is the transect index <- intersections$intersections$intersections > 0 & intersections$intersections$year %in% years b <- if (all(index == FALSE)) data.frame() else { x <- intersections$intersections[index,] y <- data.frame(coordinates(x), year=x$year) y$Observed <- T y } agentLocations <- rbind(a, b) agentLocations$Scenario <- x$scenario agentLocations$Model <- x$modelName return(agentLocations) }, modelName=modelName, years=years) } getEstimatedDensity = function(scenarios, modelName, iterations, year, sd=FALSE, weightHabitats=F) { ddply(data.frame(scenario=scenarios, year=rep(year, length(scenarios)), iteration=iterations, stringsAsFactors=FALSE), .(scenario, iteration), function(x, modelName, sd) { message("Scenario = ", x$scenario, ", iteration = ", x$iteration) withHabitatWeights <- if (weightHabitats && substr(x$scenario, 1, 1) %in% c("E","F")) T else F study <- getStudy(scenario=x$scenario, withHabitatWeights=withHabitatWeights) estimates <- study$getModel(modelName=modelName, iteration=x$iteration) estimates$modelName <- modelName estimates <- study$loadEstimates(estimates) estimates$collectEstimates() estimates$data <- subset(estimates$data, year == x$year) density <- estimates$getPopulationDensity(.parallel=T) density <- if (withHabitatWeights) { habitatWeights <- study$getHabitatWeights(iteration=x$iteration) habitatWeightsRaster <- habitatWeights$getWeightsRaster(save=F) density$weight(habitatWeightsRaster) } else wdensity <- density$copy() if (is.null(density) || maxValue(density$rasterStack[[1]]) > 1) stop("Estimation failed.") density <- as.data.frame(density$rasterStack[[1]], xy=T) density$Scenario <- x$scenario density$Iteration <- x$iteration density$Model <- modelName density$Year <- x$year return(density) }, modelName=modelName, sd=sd) } study <- getStudy(scenario="A") year <- 2010 yearLayer <- paste0("X", year) #iterations <- as.integer(c(1,4,1,4,5,38)) #iterations <- as.integer(c(1,4,1,4,5,16)) iterations <- as.integer(c(1,4,1,4,5,21)) densities <- getEstimatedDensity(scenarios=scenarios, modelName=modelNames[1], iterations=iterations, year=year) wdensities <- getEstimatedDensity(scenarios=scenarios, modelName=modelNames[1], iterations=iterations, year=year, weightHabitats=T) surveyRoutes <- getSurveyRouteLocations(scenarios=scenarios, modelName=modelNames[1], iterations=iterations, years=year) agents <- getAgentLocations(scenarios=scenarios, modelName=modelNames[1], iterations=iterations, years=year) limits <- extent(study$studyArea$boundary) # Scale density to 1/km^2 scaleDensity <- function(x, yearLayer, weight) { x[,yearLayer] <- x[,yearLayer] * weight return(x) } wscale <- 1 / (prod(res(study$getTemplateRaster())) / 1000^2) densitiesScaled <- ddply(densities, .(scenario), scaleDensity, yearLayer=yearLayer, weight=wscale) wdensitiesScaled <- ddply(wdensities, .(scenario), scaleDensity, yearLayer=yearLayer, weight=wscale) gridTemplate <- raster(extent(study$studyArea$boundary), nrows=12, ncols=6, crs=study$studyArea$proj4string) gridTemplate@extent@xmax <- gridTemplate@extent@xmin + dim(gridTemplate)[2] * 100 * 1e3 gridTemplate@extent@ymax <- gridTemplate@extent@ymin + dim(gridTemplate)[1] * 100 * 1e3 boundaryDF <- ggplot2::fortify(FinlandWTCStudy(context=context, response="canis.lupus")$studyArea$boundary) wdensitiesScaled$x <- wdensitiesScaled$x / 1000 - 3000 wdensitiesScaled$y <- wdensitiesScaled$y / 1000 - 6000 densitiesScaled$x <- densitiesScaled$x / 1000 - 3000 densitiesScaled$y <- densitiesScaled$y / 1000 - 6000 boundaryDF$lat <- boundaryDF$lat / 1000 - 6000 boundaryDF$long <- boundaryDF$long / 1000 - 3000 agents$x <- agents$x / 1000 - 3000 agents$y <- agents$y / 1000 - 6000 limits@xmin <- limits@xmin / 1000 - 3000 limits@ymin <- limits@ymin / 1000 - 6000 limits@xmax <- limits@xmax / 1000 - 3000 limits@ymax <- limits@ymax / 1000 - 6000 surveyRoutes$x <- surveyRoutes$x / 1000 - 3000 surveyRoutes$y <- surveyRoutes$y / 1000 - 6000 gridTemplate@extent@xmin <- round(gridTemplate@extent@xmin / 1000) - 3000 gridTemplate@extent@ymin <- round(gridTemplate@extent@ymin / 1000) - 6000 gridTemplate@extent@xmax <- round(gridTemplate@extent@xmax / 1000) - 3000 gridTemplate@extent@ymax <- round(gridTemplate@extent@ymax / 1000) - 6000 p1 <- ggplot(wdensitiesScaled) + geom_raster(aes_string(x="x", y="y", fill=yearLayer)) + stat_contour(data=densitiesScaled, aes_string(x="x", y="y", z=yearLayer, colour=yearLayer), size=0.1) + scale_fill_gradientn(colours=c(rep("#00A600FF", 60), terrain.colors(40)), na.value="transparent", guide=guide_legend(title=expression(paste("Population density ", (individuals/km^2)))), trans="log10", breaks=trans_breaks("log10", function(x) 10^x), labels=trans_format("log10", math_format(10^.x))) + geom_path(data=boundaryDF, aes(long, lat, group=group), colour="black", size=0.5) + geom_point(data=agents, aes(x, y, colour=Observed), size=1.3, colour="white") + geom_point(data=agents, aes(x, y, colour=Observed), size=1) + scale_color_manual(values=c("darkgray", "black"), guide=FALSE) + scale_x_continuous(limits=c(limits@xmin, limits@xmax), expand=c(0.01,0.01), breaks=pretty_breaks(n=3)) + scale_y_continuous(limits=c(limits@ymin, limits@ymax), expand=c(0.01,0.01)) + facet_grid(~scenario) + coord_equal() + theme_raster() + theme(legend.position="bottom", legend.margin=unit(-1, units="line"), legend.title=element_text(size=10), legend.text=element_text(size=10)) + #top, right, bottom, and left margins theme(strip.background=element_blank(), panel.background=element_rect(fill="transparent", colour=NA), panel.margin=unit(0, "lines"), panel.grid.minor=element_blank(), panel.grid.major=element_blank()) + xlab("X-coordinate (km)") + ylab("Y-coordinate (km)") g1 <- arrangeGrob(p1, main=textGrob("a)", x=0.01, y=0.01, just="left", gp=gpar(fontsize=10))) saveFigure(g1, "manuscript-figure-2a.png", width=17, height=8.07, dpi=200) # draft quality saveFigure(g1, "manuscript-figure-2a.pdf", width=17, height=8.07, dpi=300) xlabels <- seq(gridTemplate@extent@xmin, gridTemplate@extent@xmax, 1e5/1000) xlabels[c(2,4,6,7)] <- "" p2 <- ggplot(boundaryDF) + geom_path(aes(long, lat, group=group), colour="grey", size=0.5) + geom_point(data=surveyRoutes, aes(x, y, colour=intersections), size=1) + scale_color_manual(values=c("darkgrey", "black")) + scale_x_continuous(breaks=seq(gridTemplate@extent@xmin, gridTemplate@extent@xmax, 1e5/1000), expand=c(0.01,0.01), labels=xlabels) + scale_y_continuous(breaks=seq(gridTemplate@extent@ymin, gridTemplate@extent@ymax, 1e5/1000), expand=c(0.01,0.01)) + facet_grid(~scenario) + coord_equal() + theme_raster() + theme(panel.grid.major=element_line(colour="grey")) + #top, right, bottom, and left margins theme(strip.background=element_blank(), panel.background=element_rect(fill="transparent", colour=NA), panel.margin=unit(0, "lines"), panel.grid.minor=element_blank(), legend.position="none", axis.ticks=element_line(colour="grey", size=0.2)) + xlab("X-coordinate (km)") + ylab("Y-coordinate (km)") g2 <- arrangeGrob(p2, main=textGrob("b)", x=0.01, y=0.01, just="left", gp=gpar(fontsize=10))) saveFigure(g2, "manuscript-figure-2b.png", width=17, height=7.2, dpi=200) # draft quality saveFigure(g2, "manuscript-figure-2b.pdf", width=17, height=7.2, dpi=300) p <- arrangeGrob(g1, g2, heights=c(8.07, 7.2)/(8.07+7.2)) saveFigure(p, "manuscript-figure-2.png", width=17, height=8.07+7.2, dpi=200) # draft quality saveFigure(p, "manuscript-figure-2.pdf", width=17, height=8.07+7.2, dpi=300)
cutoffPopulationSizes <- function(populationSizes, cutoff=c(0,Inf), models) { x <- subset(populationSizes, modelName %in% models & cutoff[1] <= Estimated & Estimated <= cutoff[2]) y <- subset(populationSizes, !(modelName %in% models)) rbind(x, y) } getIntersections <- function(scenarios) { study <- getStudy(scenario="A") validation <- Validation(study=study) intersections <- validation$getValidationTemporalIntersections(scenarios=scenarios) intersections$Year <- intersections$year intersections } getPopulationSizes <- function(scenarios, modelNames, populationSizeCutoff=Inf) { study <- getStudy(scenario="A") # Get study object, scenario doesn't matter validation <- Validation(study=study, populationSizeCutoff=populationSizeCutoff) populationSizes <- validation$getValidationTemporalPopulationSizes(scenarios=scenarios, modelNames=modelNames) } getMeanPopulationSizes <- function(populationSizes) { library(plyr) meanPopulationSizes <- ddply(populationSizes, .(scenario, modelName), function(x) { x <- data.frame(True=mean(x$Observed), Estimated=mean(x$Estimated), Scenario=x$scenario[1]) if (is.infinite(x$True)) x$True <- NA if (is.infinite(x$Estimated)) x$Estimated <- NA return(x) }) return(meanPopulationSizes) } getLogPopulationSizes <- function(populationSizes) { logPopulationSizes <- populationSizes logPopulationSizes$Observed <- log10(logPopulationSizes$Observed + 1) logPopulationSizes$Estimated <- log10(logPopulationSizes$Estimated + 1) return(logPopulationSizes) } populationSizes <- getPopulationSizes(scenarios=scenarios, modelNames=modelNames[c(1,3,4,5)]) meanPopulationSizes <- getMeanPopulationSizes(cutoffPopulationSizes(populationSizes, cutoff, models=modelNames[c(1,3,5)])) logPopulationSizes <- getLogPopulationSizes(cutoffPopulationSizes(populationSizes, cutoff, models=modelNames[c(1,3,5)])) logMeanPopulationSizes <- getMeanPopulationSizes(logPopulationSizes) getPopulationSizesTable <- function(populationSizes, meanPopulationSizes, overEstimate, models, file, ...) { library(plyr) x <- ddply(populationSizes, .(Scenario, modelName), function(x, overEstimate, models) { xrows <- nrow(x) old <- x x <- cutoffPopulationSizes(x, overEstimate, models=models) if (nrow(x) == 0) { print(old) stop("ERROR") } if (nrow(x) == 0) return(NULL) failed <- (1 - nrow(x) / xrows) * 100 y <- lm(Estimated~Observed, x) #prop <- x$Estimated/x$Observed prop <- (x$Estimated+1)/(x$Observed+1) prop <- prop[!is.infinite(prop) & prop > 0] prop.mean <- mean(log(prop)) prop.sd <- sd(log(prop)) prop.min <- quantile(log(prop), .025) prop.max <- quantile(log(prop), .975) slope.sd <- coef(summary(y))[2,2] #print(prop) #message("prop = ", mean(prop), " log mean = ", prop.mean) df <- data.frame(Intercept=coef(y)[1], Slope=coef(y)[2], SlopeSD=slope.sd, R2=summary(y)$r.squared, EstTrue=prop.mean, EstTrueSD=prop.sd, EstTrueMin=prop.min, EstTrueMax=prop.max, Failed=failed, N=nrow(x)) #print(df) return(df) }, overEstimate=overEstimate, models=models, .inform=T) y <- unique(meanPopulationSizes) y$Scenario <- y$scenario z <- plyr::join(x, y) z$Model <- prettyModelName(z$modelName) z$Scenario <- prettyScenario(z$Scenario) return(z) } table1 <- getPopulationSizesTable(populationSizes, meanPopulationSizes, overEstimate=cutoff, models=modelNames[c(1,3,5)]) logTable1 <- getPopulationSizesTable(logPopulationSizes, logMeanPopulationSizes, overEstimate=log10(cutoff+1), models=modelNames[c(1,3,5)]) logTable1$Failed <- table1$Failed #subset(populationSizes, Scenario=="C" & modelName=="FMPModel") #subset(populationSizes, Scenario=="B" & modelName=="FMPModel") writeTable <- function(df, filename, context) { library(ReporteRs) doc <- docx() flextable <- FlexTable(df) addFlexTable(doc, flextable) dir <- normalizePath(context$figuresDirectory) writeDoc(doc, file=file.path(dir, filename)) } facet_wrap_labeller <- function(gg.plot,labels=NULL) { #works with R 3.0.1 and ggplot2 0.9.3.1 require(gridExtra) g <- ggplotGrob(gg.plot) gg <- g$grobs strips <- grep("strip_t", names(gg)) for(ii in seq_along(labels)) { modgrob <- getGrob(gg[[strips[ii]]], "strip.text", grep=TRUE, global=TRUE) gg[[strips[ii]]]$children[[modgrob$name]] <- editGrob(modgrob,label=labels[ii]) } g$grobs <- gg class(g) = c("arrange", "ggplot",class(g)) g } plotPopulationSizesTable <- function(logTable) { library(reshape2) library(dplyr) library(gtable) library(gridExtra) x <- melt(logTable, id.vars=c("Scenario", "Model"), measure.vars=c("EstTrue", "Slope", "R2", "Failed")) bestFit <- data.frame(variable=unique(x$variable), y=c(0, 1, 1, 0)) errorBars <- x errorBars$value <- NA errorBars$value[errorBars$variable=="EstTrue"] <- logTable$EstTrue errorBars$qmin[errorBars$variable=="EstTrue"] <- logTable$EstTrueMin errorBars$qmax[errorBars$variable=="EstTrue"] <- logTable$EstTrueMax margin <- unit(c(0.1, 0.1, -.75, 0), units="line") p1 <- x %>% filter(variable=="EstTrue") %>% ggplot(aes(Scenario, value, group=Model)) + geom_hline(data=data.frame(y=0), aes(yintercept=y), color="darkgray") + geom_errorbar(data=errorBars, aes(Scenario, group=Model, ymax=qmax, ymin=qmin), position=position_dodge(width=.7), colour="darkgray") + geom_point(aes(shape=Model), position=position_dodge(width=.7)) + xlab("") + ylab(expression(alpha)) + theme_ms() + theme(panel.border=element_rect(fill="transparent"), legend.title=element_text(size=10), legend.position="none", strip.text.x=element_text(size=10), strip.text.y=element_text(size=10), plot.margin=margin) #top, right, bottom, and left margins g1 <- arrangeGrob(p1, main=textGrob("a)", x=0.01, y=0.01, just="left", gp=gpar(fontsize=10))) p2 <- x %>% filter(variable=="Slope") %>% ggplot(aes(Scenario, value, group=Model)) + geom_hline(data=data.frame(y=1), aes(yintercept=y), color="darkgray") + geom_point(aes(shape=Model), position=position_dodge(width=.7)) + xlab("") + ylab(expression(beta)) + theme_ms() + theme(panel.border=element_rect(fill="transparent"), legend.title=element_text(size=10), legend.position="none", strip.text.x=element_text(size=10), strip.text.y=element_text(size=10), plot.margin=margin) g2 <- arrangeGrob(p2, main=textGrob("b)", x=0.01, y=0.01, just="left", gp=gpar(fontsize=10))) p3 <- x %>% filter(variable=="R2") %>% ggplot(aes(Scenario, value, group=Model)) + geom_hline(data=data.frame(y=1), aes(yintercept=y), color="darkgray") + geom_point(aes(shape=Model), position=position_dodge(width=.7)) + xlab("") + ylab(expression(R^2)) + theme_ms() + theme(panel.border=element_rect(fill="transparent"), legend.title=element_text(size=10), legend.position="none", strip.text.x=element_text(size=10), strip.text.y=element_text(size=10), plot.margin=margin) g3 <- arrangeGrob(p3, main=textGrob("c)", x=0.01, y=0.01, just="left", gp=gpar(fontsize=10))) p4 <- x %>% filter(variable=="Failed") %>% ggplot(aes(Scenario, value, group=Model)) + geom_hline(data=data.frame(y=0), aes(yintercept=y), color="darkgray") + geom_point(aes(shape=Model), position=position_dodge(width=.7)) + xlab("") + ylab("Failed (%)") + ylim(0,100) + theme_ms() + theme(panel.border=element_rect(fill="transparent"), legend.title=element_text(size=10), legend.position="none", strip.text.x=element_text(size=10), strip.text.y=element_text(size=10), plot.margin=margin) g4 <- arrangeGrob(p4, main=textGrob("d)", x=0.01, y=0.01, just="left", gp=gpar(fontsize=10))) legend <- gtable_filter(ggplot_gtable(ggplot_build(p1 + theme(legend.position="bottom", legend.margin=unit(-1, units="line")))), "guide-box") lheight <- sum(legend$height) theight <- unit(10, "points") #p <- arrangeGrob(g1, g2, g3, g4, ncol=2) p1 <- arrangeGrob(g1, g3, nrow=2) p2 <- arrangeGrob(g2, g4, nrow=2) p1 <- arrangeGrob(p1, textGrob("Scenario", x=0.58, gp=gpar(fontsize=10)), heights=unit.c(unit(1, "npc") - theight, theight)) p2 <- arrangeGrob(p2, textGrob("Scenario", x=0.57, gp=gpar(fontsize=10)), heights=unit.c(unit(1, "npc") - theight, theight)) p <- arrangeGrob(p1, p2, ncol=2) #p <- arrangeGrob(p, textGrob("Scenario", gp=gpar(fontsize=10)), heights=unit.c(unit(1, "npc") - theight, theight)) p <- arrangeGrob(p, legend, heights=unit.c(unit(1, "npc") - lheight, lheight), ncol=1) return(p) } logTable <- logTable1 p <- plotPopulationSizesTable(logTable) saveFigure(p, "manuscript-figure-3.png", width=17, dpi=200) # draft quality saveFigure(p, "manuscript-figure-3.pdf", width=17, dpi=300) logTable <- logTable1 logTablePretty <- logTable logTablePretty$Slope <- round(logTablePretty$Slope, 2) logTablePretty$R2 <- round(logTablePretty$R2, 2) logTablePretty$EstTrue <- round(logTablePretty$EstTrue, 2) logTablePretty$EstTrueSD <- round(logTablePretty$EstTrueSD, 2) logTablePretty$Failed <- round(logTablePretty$Failed, 0) logTablePretty <- logTablePretty[,c("Scenario","Model","Slope","R2","EstTrue","EstTrueSD","Failed")] writeTable(subset(logTablePretty, Model != "VT"), "table-1.docx", context)
logPopulationSizesSample <- subset(logPopulationSizes, scenario == "A") logPopulationSizesSample <- cutoffPopulationSizes(logPopulationSizesSample, cutoff, models=modelNames[c(1,3,5)]) logPopulationSizesSample$modelName <- prettyModelName(logPopulationSizesSample$modelName) logPopulationSizesSample <- rbind( subset(logPopulationSizesSample, modelName=="ST"), subset(logPopulationSizesSample, modelName=="T1"), subset(logPopulationSizesSample, modelName=="T2"), subset(logPopulationSizesSample, modelName=="FMP")) lm(Estimated~Observed, subset(logPopulationSizes, modelName=="FMPModel" & Scenario=="A")) lm(Estimated~Observed, subset(logPopulationSizesSample, modelName=="FMP" & Scenario=="A")) plotLogPopulationSizeSample <- function(size, modelName, case, xlim, ylim) { library(ggplot2) library(scales) size$Scenario <- prettyScenario(size$Scenario) meanSize <- getMeanPopulationSizes(size) meanSize$Scenario <- prettyScenario(meanSize$Scenario) p <- ggplot(size, aes(Observed, Estimated, colour=modelName)) + geom_point(size=1, colour="gray40", alpha=.5) + geom_abline(colour="black", linetype="longdash") + geom_point(data=meanSize, aes(True, Estimated), size=2, shape=15, colour="black") #if (!any(abs(size$Estimated) == Inf)) p <- p + geom_smooth(method="lm", se=F, fullrange=T, colour="black") ls <- ddply(size, .(modelName), function(x) if (!any(abs(x$Estimated) == Inf)) coef(lm(Estimated~Observed, x)) else c(NA,NA)) print(ls) p <- p + geom_abline(data=ls, aes(intercept=`(Intercept)`, slope=Observed)) p <- p + facet_wrap(~modelName) + scale_x_continuous(limits=xlim, labels=trans_format("identity", math_format(10^.x))) + scale_y_continuous(limits=ylim, labels=trans_format("identity", math_format(10^.x))) + xlab("True population size") + ylab("Estimated population size") + theme_ms2() + theme(panel.border=element_rect(fill="transparent"), plot.title=element_text(size=10), legend.margin=unit(0, units="line")) return(p) } p <- plotLogPopulationSizeSample(subset(logPopulationSizesSample, modelName %in% c("ST","T1","T2","FMP")), modelName=c("ST","T1","T2","FMP"), "a", xlim=log10(c(1, 10^3)), ylim=log10(c(1,10^3))) saveFigure(p, "manuscript-figure-3sample.png", height=7.87, width=17) # draft quality saveFigure(p, "manuscript-figure-3sample.pdf", height=7.87, width=17, dpi=300)
plotLogPopulationSize <- function(size, modelName, case, xlim, ylim) { library(ggplot2) library(scales) library(grid) size$Scenario <- prettyScenario(size$Scenario) meanSize <- getMeanPopulationSizes(size) meanSize$Scenario <- prettyScenario(meanSize$Scenario) p <- ggplot(size, aes(Observed, Estimated, colour=modelName)) + geom_point(size=1, colour="gray40", alpha=.5) + geom_abline(colour="black", linetype="longdash") + geom_point(data=meanSize, aes(True, Estimated), size=2, shape=15, colour="black") ls <- ddply(size, .(Scenario), function(x) if (!any(abs(x$Estimated) == Inf)) coef(lm(Estimated~Observed, x)) else c(NA,NA)) p <- p + geom_abline(data=ls, aes(intercept=`(Intercept)`, slope=Observed)) p <- p + facet_wrap(~Scenario) + scale_x_continuous(limits=xlim, labels=trans_format("identity", math_format(.x))) + scale_y_continuous(limits=ylim, labels=trans_format("identity", math_format(.x))) + xlab(expression(paste(Log[10], " true population size"))) + ylab(expression(paste(Log[10], " estimated population size"))) + ggtitle(paste("Model:", prettyModelName(modelName[1]))) + theme_ms2() + theme(panel.border=element_rect(fill="transparent"), legend.margin=unit(0, units="line"), panel.margin=unit(0, units="line")) return(p) } #modelNames <- modelNames[c(2:4,1)] #xlim <- log10(c(1,1200)); ylim <- xlim xlim <- log10(c(1,2000)); ylim <- xlim p <- plotLogPopulationSize(subset(logPopulationSizes, modelName==modelNames[1]), modelName=modelNames[1], "1", xlim=xlim, ylim=ylim) saveFigure(p, "manuscript-figure-3a1.png", width=17, dpi=200) saveFigure(p, "manuscript-figure-3a1.pdf", width=17, dpi=300) p <- plotLogPopulationSize(subset(logPopulationSizes, modelName==modelNames[3]), modelName=modelNames[3], "1", xlim=xlim, ylim=ylim) saveFigure(p, "manuscript-figure-3a2.png", width=17, dpi=200) saveFigure(p, "manuscript-figure-3a2.pdf", width=17, dpi=300) p <- plotLogPopulationSize(subset(logPopulationSizes, modelName==modelNames[4]), modelName=modelNames[4], "1", xlim=xlim, ylim=ylim) saveFigure(p, "manuscript-figure-3a3.png", width=17, dpi=200) saveFigure(p, "manuscript-figure-3a3.pdf", width=17, dpi=300) p <- plotLogPopulationSize(subset(logPopulationSizes, modelName==modelNames[5]), modelName=modelNames[5], "1", xlim=xlim, ylim=ylim) saveFigure(p, "manuscript-figure-3a4.png", width=17, dpi=200) saveFigure(p, "manuscript-figure-3a4.pdf", width=17, dpi=300)
ddply(logTable, .(modelName, Case, Scenario), function(x) { data.frame(p=mean(exp(x$EstTrue))) }) ddply(logTable, .(modelName, Case), function(x) { data.frame(p=mean(exp(x$EstTrue))) }) ddply(logTable, .(modelName, Scenario), function(x) { data.frame(p=mean(exp(x$EstTrue))) }) p <- ddply(subset(logTable, Scenario %in% c("A","E")), .(modelName, Case), function(x) { p <- 1 - exp(subset(x, Scenario %in% "A")$EstTrue) / exp(subset(x, Scenario %in% "E")$EstTrue) data.frame(p=p) }) mean(p$p) p <- ddply(subset(logTable, Scenario %in% c("A","F")), .(modelName, Case), function(x) { p <- 1 - exp(subset(x, Scenario %in% "A")$EstTrue) / exp(subset(x, Scenario %in% "F")$EstTrue) data.frame(p=p) })
plotIntersectionPopulationSizeTable <- function(ratioTable) { library(reshape2) library(grid) x <- melt(ratioTable, id.vars=c("Scenario", "Case"), measure.vars=c("SizeCount")) errorBars <- x errorBars$value <- NA errorBars$value[errorBars$variable=="SizeCount"] <- ratioTable$SizeCount errorBars$qmin[errorBars$variable=="SizeCount"] <- ratioTable$SizeCountMin errorBars$qmax[errorBars$variable=="SizeCount"] <- ratioTable$SizeCountMax p <- ggplot(x, aes(Scenario, value, group=Case)) + geom_errorbar(data=errorBars, aes(Scenario, group=Case, ymax=qmax, ymin=qmin), position=position_dodge(width=.7), colour="darkgray") + geom_point(aes(shape=Case), position=position_dodge(width=0.7)) + ylab("Ratio") + #guides(shape=guide_legend()) + theme_minimal(base_size=10) + theme(panel.border=element_rect(fill="transparent"), legend.title=element_text(size=10), strip.text.y=element_text(size=10), legend.position="right", #legend.background=element_rect(fill="transparent"), plot.margin=unit(c(0.2, 0, 0, 0), units="line"), legend.margin=unit(-1, units="line"), panel.margin=unit(0, units="line"))# , legend.key=element_rect(fill="transparent")) return(p) } intersections <- getIntersections(scenarios) intersections2000 <- getIntersections(scenarios=paste0(scenarios, "combined")) intersections10days <- getIntersections(scenarios=paste0(scenarios, "10days")) intersections10days$intersections <- intersections10days$intersections / 10 intersectionsx <- plyr::join(populationSizes, intersections) intersections2000x <- plyr::join(populationSizes2000, intersections2000) intersections10daysx <- plyr::join(populationSizes10Days, intersections10days) rer.stat <- function(x) data.frame(SizeCount=mean(x$intersections/x$Observed), SizeCountMin=quantile(x$intersections/x$Observed, 0.025), SizeCountMax=quantile(x$intersections/x$Observed, 0.975)) i <- ddply(intersectionsx, .(Scenario), rer.stat); i$Scenario=prettyScenario(i$Scenario); i$Case <- "1" i2000 <- ddply(intersections2000x, .(Scenario), rer.stat); i2000$Scenario=prettyScenario(i2000$Scenario); i2000$Case <- "2" i10days <- ddply(intersections10daysx, .(Scenario), rer.stat); i10days$Scenario=prettyScenario(i10days$Scenario); i10days$Case <- "3" iall <- rbind(i, i2000, i10days) iall p <- plotIntersectionPopulationSizeTable(iall) saveFigure(p, "manuscript-figure-4x.png", height=7.87/2) # draft quality saveFigure(p, "manuscript-figure-4x.pdf", height=7.87/2, dpi=300) #writeTable(i, "table-2a.docx", context) #writeTable(i2000, "table-2b.docx", context) #writeTable(i10days, "table-2c.docx", context)
getPopulationCorrelations <- function(scenarios, modelNames, cutoff=c(-Inf,Inf)) { study <- getStudy(scenario="A") validation <- Validation(study=study, populationSizeCutoff=cutoff) populationCorrelations <- validation$getValidationSpatialPopulationSizes(scenarios=scenarios, modelNames=modelNames) return(populationCorrelations) } plotPopulationCorrelations <- function(populationCorrelations, case) { populationCorrelations$Scenario <- prettyScenario(populationCorrelations$Scenario) p <- ggplot(populationCorrelations, aes(True, Correlation)) + geom_hline(data=data.frame(y=0), aes(yintercept=y), color="darkgray") + geom_point(size=1, colour="black") + facet_wrap(~Scenario, ncol=2) + xlab("True population size") + scale_y_continuous(limits=c(-1,1)) + #ggtitle(paste("Case:", case)) + theme_ms() + theme(panel.border=element_rect(fill="transparent"), panel.margin=unit(0, units="line"))#, axis.text.x=element_text(angle=45, hjust=.5, vjust=1)) return(p) } pc1a <- getPopulationCorrelations(scenarios="A", modelNames=modelNames[1], cutoff=cutoff) pc1b <- getPopulationCorrelations(scenarios="B", modelNames=modelNames[1], cutoff=cutoff) pc1c <- getPopulationCorrelations(scenarios="C", modelNames=modelNames[1], cutoff=cutoff) pc1d <- getPopulationCorrelations(scenarios="D", modelNames=modelNames[1], cutoff=cutoff) pc1e <- getPopulationCorrelations(scenarios="E", modelNames=modelNames[1], cutoff=cutoff) pc1f <- getPopulationCorrelations(scenarios="F", modelNames=modelNames[1], cutoff=cutoff) pc1 <- rbind(pc1a,pc1b,pc1c,pc1d,pc1e,pc1f) p1 <- plotPopulationCorrelations(pc1, "Case 1") saveFigure(p1, "manuscript-figure-4a.png") saveFigure(p1, "manuscript-figure-4a.pdf", dpi=300)
getCIProportions <- function(scenarios, modelNames, panel, scenarioNames, cutoff=Inf) { study <- getStudy(scenario="A") validation <- Validation(study=study, populationSizeCutoff=cutoff) credibilityIntervals <- list() if ("FMPModel" %in% modelNames) { for (scenario in scenarios) { study <- getStudy(scenario=scenario) validation <- Validation(study=study, populationSizeCutoff=cutoff) iterations <- validation$getCredibilityIntervalsValidationIterations("FMPModel") p95 <- c() p50 <- c() for (i in iterations) { load(file=validation$getCredibilityIntervalsValidationFileName(modelName="FMPModel", iteration=i)) p95 <- c(p95, bootCI$p95) p50 <- c(p50, bootCI$p50) } a <- data.frame(scenario=scenario, Proportion=mean(p95, na.rm=T), n=NA, True=NA, Estimated=NA, Estimated.q1=NA, Estimated.q2=NA, modelName="FMPModel", Interval="95%") b <- data.frame(scenario=scenario, Proportion=mean(p50, na.rm=T), n=NA, True=NA, Estimated=NA, Estimated.q1=NA, Estimated.q2=NA, modelName="FMPModel", Interval="50%") credibilityIntervals <- rbind(credibilityIntervals, rbind(a, b)) } } a <- validation$getValidatedCredibilityIntervalsProportions(scenarios=scenarios, modelNames=modelNames[!modelNames %in% "FMPModel"]) b <- validation$getValidatedCredibilityIntervalsProportions(scenarios=scenarios, modelNames=modelNames[!modelNames %in% "FMPModel"], probs=c(.25, .75), probsName="50%") credibilityIntervals <- rbind(credibilityIntervals, rbind(a, b)) credibilityIntervals$panel <- panel return(credibilityIntervals) } cip1 <- getCIProportions(scenarios=scenarios, modelNames=modelNames[c(1,3,4,5)], panel="Case 1", scenarioNames=scenarios, cutoff=cutoff) plotCIProportions <- function(credibilityIntervals) { credibilityIntervals$Model <- prettyModelName(credibilityIntervals$modelName) credibilityIntervals$CrI <- credibilityIntervals$Interval p <- ggplot(credibilityIntervals, aes(scenario, Proportion, shape=Model, colour=CrI)) + geom_hline(aes(yintercept=.5), color="darkgrey") + geom_hline(aes(yintercept=.95), color="darkgrey") + geom_point(aes(shape=Model, group=Model), position=position_dodge(width=.7)) + scale_colour_manual("CrI / CI", values=c("black", "grey"), guide=guide_legend(override.aes=list(size=3.5, shape=15))) + theme_ms() + theme(panel.border=element_rect(fill="transparent")) + xlab("Scenario") + ylab("Proportion") + ylim(0, 1) + theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(), panel.margin=unit(0, units="line"), legend.direction="vertical", legend.box="horizontal", legend.key.height=unit(0.7, "line")) #guides(scale_colour_manual=guide_colourbar()) library(gridExtra) g_legend <- function(a.gplot) { tmp <- ggplot_gtable(ggplot_build(a.gplot)) leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box") legend <- tmp$grobs[[leg]] legend } legend <- g_legend(p) p <- arrangeGrob(p + theme(legend.position="none"), legend, heights=c(6/10,4/10)) return(p) } cip <- cip1 p <- plotCIProportions(cip) saveFigure(p, "manuscript-figure-5.png", height=7.87*2/3) saveFigure(p, "manuscript-figure-5.pdf", height=7.87*2/3, dpi=300)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.