In statguy/STREM: Winter track counts
library(knitr)
single_width <- 2.75591 * 2
double_width <- single_width*2 + 0.5 * 2
single_height <- single_width
double_height <- double_width
#opts_knit$set(base.dir="~/.pandoc")
opts_chunk$set(echo=FALSE, warning=FALSE, message=FALSE, tidy=FALSE, fig.width=single_width, fig.height=single_height, cache=TRUE, cache.path="~/tmp/cache/")
base_size <- 16
library(parallel)
library(doMC)
registerDoMC(cores=detectCores())
library(STREM)
source("~/git/STREM/setup/WTC-Boot.R")
library(sp)
library(ggplot2)
library(rasterVis)
library(reshape2)
library(plyr)
context <- Context$new(resultDataDirectory=wd.data.results, processedDataDirectory=wd.data.processed, rawDataDirectory=wd.data.raw, scratchDirectory=wd.scratch, figuresDirectory=wd.figures)
study <- FinlandWTCStudy$new(context=context, response="canis.lupus")
study$studyArea$loadBoundary(thin=F, tolerance=0.001)
boundaryDF <- ggplot2::fortify(study$studyArea$boundary)
responses <- c("canis.lupus", "lynx.lynx", "rangifer.tarandus.fennicus")
Figure 1: Survey routes
surveyRoutes <- study$loadSurveyRoutes(findLengths=FALSE)
surveyRoutesSP <- SpatialLinesDataFrame(surveyRoutes$surveyRoutes, data=data.frame(x=1:length(surveyRoutes$surveyRoutes)), match.ID=FALSE)
surveyRoutesDF <- ggplot2::fortify(surveyRoutesSP)
p <- ggplot(boundaryDF, aes(long, lat, group=group)) + geom_polygon(colour="black", fill=NA) +
geom_polygon(data=surveyRoutesDF, aes(long, lat, group=group), colour="blue", fill=NA) +
coord_equal() + theme_raster()
print(p)
#saveFigure(p, filename="SurveyRoutes.svg", bg="transparent")
Number of counted survey routes
intersections <- study$loadIntersections(predictDistances=FALSE)
x <- ddply(as.data.frame(intersections$intersections), .(year), nrow)
s <- sd(x$V1)
m <- mean(x$V1)
cat("mean = ", m, "\n")
cat("sd = ", s, "\n")
cat("range = ", m+c(-s,s), "\n")
% of intersections D > 1
intersections <- FinlandWTCIntersections$new(study=study, maxDuration=1)
intersections$saveIntersections()
Figure 2: Crossing rates (D > 1 removed)
getIntersectionsRate <- function(responses, context) {
library(plyr)
x <- ldply(responses, function(response, context) {
study <- FinlandWTCStudy$new(context=context, response=response)
intersections <- study$loadIntersections(predictDistances=FALSE)
x <- ddply(intersections$intersections@data, .(year), function(x, response) {
data.frame(intersections=sum(x$intersections / (x$length/1000) / x$duration) / nrow(x), response=study$getPrettyResponse(response))
}, response=response)
return(x)
}, context=context)
return(x)
}
intersectionsRate <- getIntersectionsRate(responses, context)
p <- ggplot(intersectionsRate, aes(year, intersections)) + geom_line(size=1) + facet_grid(~response, scales="free_y") +
xlab("Year") + ylab("Crossings rate") + theme_presentation(base_size)
print(p)
#saveFigure(p, filename="IntersectionsRate.svg", bg="transparent")
GPS figures
getTracks <- function(responses, context) {
library(plyr)
library(ggplot2)
library(adehabitatLT)
l_ply(responses, function(response, context) {
study <- FinlandWTCStudy$new(context=context, response=response)
tracks <- FinlandWTCTracks$new(study=study)
tracks$loadTracks()
message("Response = ", response)
message("Number of individuals in the cleaned data = ", length(unique(adehabitatLT::id(tracks$tracks))))
d <- as.POSIXlt(ld(tracks$tracks)$date)
message("Years in the cleaned data = ", paste(sort(unique(d$year+1900)), collapse=","), " (", length(unique(d$year)), ")")
message("Number of samples in the cleaned data = ", length(d))
message("Number of days in the cleaned data = ", length(unique(paste(d$year, d$yday))))
}, context=context)
}
getTracks(responses=responses, context=context)
Figure 3: GPS movements, juvenile is no more than 1 years old individual
library(plyr)
response <- "canis.lupus"
study <- FinlandWTCStudy$new(context=context, response=response)
tracks <- FinlandWTCTracks$new(study=study)
tracks$loadTracks()
tracks$loadMetadata()
xtabs(sex~born, tracks$metadata)
library(reshape2)
dcast(tracks$metadata, .~sex)
dcast(tracks$metadata, .~born)
dcast(tracks$metadata, .~dispersing)$"TRUE"
dcast(tracks$metadata, .~homerange)$"TRUE"
dcast(tracks$metadata, .~alpha)
x <- ld(tracks$tracks)$dt/3600
sort(x[x<=24])
sort(x[x<=24])*3600
tracksSP <- tracks$getSpatialLines(variables=.(id,year,burst))
x <- tracksSP@data
x$juvenile <- if ("born" %in% colnames(x)) factor(ifelse(x$year - x$born <= 1, "juvenile", "adult")) else factor("unknown")
x$juvenile <- factor(x$juvenile, levels=c(levels(x$juvenile), "unknown"))
x$juvenile[is.na(x$juvenile)] <- "unknown"
dcast(x, .~juvenile)
getTracks <- function(responses, context) {
library(plyr)
library(ggplot2)
x <- ldply(responses, function(response, context) {
study <- FinlandWTCStudy$new(context=context, response=response)
tracks <- FinlandWTCTracks$new(study=study)
tracks$loadTracks()
tracksSP <- tracks$getSpatialLines(variables=.(id,year,burst))
x <- tracksSP@data
if (!"sex" %in% colnames(x)) x$sex <- "unknown"
x$juvenile <- if ("born" %in% colnames(x)) factor(ifelse(x$year - x$born <= 1, "juvenile", "adult")) else factor("unknown")
x$juvenile <- factor(x$juvenile, levels=c(levels(x$juvenile), "unknown"))
x$juvenile[is.na(x$juvenile)] <- "unknown"
x$id <- sapply(tracksSP@lines, function(x) x@ID)
x$sexjuv <- factor(paste(x$sex, x$juvenile, sep=", "))
tracksDF <- ggplot2::fortify(tracksSP)
tracksDF <- plyr::join(tracksDF, x, by="id")
tracksDF$response <- study$getPrettyResponse(response)
return(tracksDF)
}, context=context)
return(x)
}
tracks <- getTracks(responses=responses, context=context)
#levels(tracks$sexjuv)
p <- ggplot(tracks, aes(long, lat, group=group)) +
geom_polygon(data=boundaryDF, aes(long, lat, group=group), colour="darkgray", fill=NA) +
geom_path(aes(colour=sexjuv, linetype=sexjuv), alpha=.7) +
scale_color_manual("Sex, age", values=c("darkred","red","red","darkblue","blue","blue","black")) + scale_linetype_manual("Sex, age", values=c("solid","dotted","solid","solid","dotted","solid","solid")) +
facet_grid(~response) + theme_raster(base_size, aspect.ratio=1, legend.position=c(0.5,-0.1), legend.direction="horizontal", legend.position="bottom") +
guides(colour=guide_legend(nrow=2)) + coord_cartesian(ylim=range(tracks$lat)+c(-1,1)*50000)
print(p)
#saveFigure(p, filename="GPS.svg", bg="transparent")
Figure 4: Straight-line-distance-error
ss <- SimulationStudy$new(response="Intensive")$setup(context=context)
tracks <- ss$loadTracks(iteration=as.integer(1))
tracks$tracks <- subset(tracks$tracks, id == 1 & year == 2001 & month == 1 & day %in% c(2,3))
intervals <- MovementSampleIntervals$new(study=ss)
thinnedTracks <- intervals$getThinnedTracksSampleIntervals(tracks=tracks)
plotTracks <- data.frame()
for (i in seq(1, 11, by=5)) {
x <- thinnedTracks$tracksList[[i]]$tracks
x$thin <- i
x <- subset(x, date >= as.POSIXct("2001-01-03 00:00:00") & date <= as.POSIXct("2001-01-03 03:40:00"))
plotTracks <- rbind(plotTracks, x)
}
plotTracks$thin <- factor(plotTracks$thin)
plotTracks$dt <- factor(plotTracks$dt / 60)
segment <- aes(x=x, y=y, xend=c(tail(x, n=-1), NA), yend=c(tail(y, n=-1), NA), group=dt, colour=dt)
ending <- arrow(length=unit(0.7, "cm"))
p <- ggplot(plotTracks, mapping=segment) +
theme_raster(base_size, legend.position=c(0.5,0.85), legend.background=element_rect(color="grey")) +
#guides(colour=guide_legend(title=expression(paste(Delta, "t (min)")))) +
guides(colour=guide_legend(title="Sampling interval (min)")) +
scale_colour_manual(values=c("steelblue1","steelblue3","steelblue4")) +
#ggtitle("Straight-line distance error") +
geom_path(size=1, arrow=ending) + geom_point(size=2)
#geom_segment(size=2, arrow=ending) # doesn't work with grouping
print(p)
#saveFigure(p, filename="StraightLineDistanceError.svg", bg="transparent")
Raw GPS sample intervals
intervals <- llply(responses, function(response, context) {
study <- FinlandWTCStudy$new(context=context, response=response)
tracks <- FinlandWTCTracks$new(study=study)
tracks$loadTracks()
intervals <- MovementSampleIntervals$new(study=study)
intervals$loadSampleIntervals()
return(table(intervals$intervals$intervalH))
}, context=context)
intervals
Figure 5: Power law fits for distance - sampling interval
distances <- ldply(responses, function(response, context) {
study <- FinlandWTCStudy$new(context=context, response=response)
intervals <- study$loadSampleIntervals()
x <- intervals$intervals
x$response <- study$getPrettyResponse(response=response)
return(x)
}, context=context)
intervals <- ldply(responses, function(response, context) {
study <- FinlandWTCStudy$new(context=context, response=response)
intervals <- study$loadSampleIntervals()
distances <- intervals$getDistanceCurve()
distances$response <- study$getPrettyResponse(response=response)
return(distances)
}, context=context)
p <- ggplot(distances, aes(intervalH, distanceKm)) + geom_point(alpha=.3) +
ylab("Distance / day (km)") + xlab("Sampling interval (h)") +
geom_line(data=intervals, aes(intervalH, distanceKm), size=1) +
geom_smooth(data=intervals, aes(ymin=distanceKm25, ymax=distanceKm975), stat="identity") +
ylim(c(0, max(distances$distanceKm))) +
facet_grid(~response) + theme_presentation(base_size)
print(p)
#saveFigure(p, filename=paste("DistanceCorrection-", response, ".svg", sep=""), bg="transparent")
Table 1: Observed distances - estimated distances
observedDistances <- ldply(responses, function(response, context) {
study <- FinlandWTCStudy$new(context=context, response=response)
tracks <- study$loadTracks()
intervals <- tracks$getSampleIntervals()
return(data.frame(distance=intervals$intervals$dist, logDistance=log(intervals$intervals$dist), response=study$getPrettyResponse(response)))
}, context=context)
predictedDistances <- ldply(responses, function(response, context) {
study <- FinlandWTCStudy$new(context=context, response=response, distanceCovariatesModel=~populationDensity+rrday+snow+tday-1, trackSampleInterval=2)
intersections <- study$loadIntersections()
intersections$predictDistances()
distances <- intersections$intersections$distance / 1000
return(data.frame(distance=distances, logDistance=log(distances), response=study$getPrettyResponse(response)))
}, context=context)
cat("Observed distances (km/day):\n")
ddply(observedDistances, .(response), summarize, mean=mean(distance), sd=sd(distance), min=min(distance), max=max(distance))
cat("Predicted distances (km/day):\n")
ddply(predictedDistances, .(response), summarize, mean=mean(distance), sd=sd(distance), min=min(distance), max=max(distance))
cat("Prediction interval = 2h")
Distances from literature
- Wolf: 7-25 km/day
- Lynx:
- FFR:
Table 2: Power law fit estimates for distance - sampling interval
- fixed_effect = human population density
- fixed_effect1 = precipitation
- fixed_effect2 = snow cover
- fixed_effect3 = temperature
# Fixed effects: populationDensity + rrday + snow + tday
intervals <- l_ply(responses, function(response, context) {
study <- FinlandWTCStudy$new(context=context, response=response)
intervals <- study$loadSampleIntervals()
cat("response = ", study$getPrettyResponse(response=response), "\n")
intervals$estimatedValuesSummary()
cat("\n\n")
}, context=context)
Table 3: Crude distance estimates
responses <- c("canis.lupus","lynx.lynx","rangifer.tarandus.fennicus")
for (response in responses) {
context <- Context$new(resultDataDirectory=wd.data.results, processedDataDirectory=wd.data.processed, rawDataDirectory=wd.data.raw, scratchDirectory=wd.scratch, figuresDirectory=wd.figures)
study <- FinlandWTCStudy$new(context=context, response=response, distanceCovariatesModel=~populationDensity+rrday+snow+tday-1, trackSampleInterval=2)
populationSize <- FinlandPopulationSize$new(study=study)
populationSize$loadValidationData()
intersections <- study$loadIntersections()
A <- 338424 # km^2
M <- ddply(intersections$getData(), .(year), function(x) sum(x$length))$V1 / 1000
N <- populationSize$sizeData$Validation
x <- ddply(intersections$getData(), .(year), function(x) sum(x$intersections))$V1
C <- 2/pi * A / M
l <- C * x / N
cat(response, ": crude distance mean = ", mean(l, na.rm=T), " median = ", median(l, na.rm=T), " sd = ", sd(l, na.rm=T), " km/d\n")
h <- study$loadHabitatWeightsRaster()
h.masked <- mask(h, study$studyArea$boundary)
n.cell <- length(h.masked) - freq(h.masked, value=NA)
A.discrete <- n.cell * prod(res(h.masked)) / 1000^2
h.cell <- cellStats(h.masked, sum)
A.h <- h.cell * prod(res(h.masked)) / 1000^2 * A / A.discrete # A.discrete is a little bit smaller than A due to the discretization, so apply a correction
C <- 2/pi * A.h / M
l <- C * x / N
cat(response, ": crude weighted distance mean = ", mean(l, na.rm=T), " median = ", median(l, na.rm=T), " sd = ", sd(l, na.rm=T), " km/d\n")
cat(response, ": observed number of intersections:\n")
print(1/C * l)
cat("Area correction = ", A/A.discrete, "\n")
}
Figure 6: Habitat weights
getHabitatWeights <- function(responses, context) {
library(plyr)
x <- ldply(responses, function(response, context) {
study <- FinlandWTCStudy$new(context=context, response=response)
weights <- study$loadHabitatWeights()$getHabitatSelectionWeights()
n <- names(weights)
print(n)
weights <- data.frame(habitat=factor(n, levels=n), weights)
weights$response <- study$getPrettyResponse(response)
return(weights)
}, context=context)
return(x)
}
weights <- getHabitatWeights(responses=responses, context=context)
p <- ggplot(weights, aes(habitat, weights, fill=habitat)) +
geom_bar(stat="identity") + facet_grid(~response) + scale_fill_manual(values=c("#beaed4","#ffff99","#7fc97f","#fdc086","#386cb0")) +
xlab("") + ylab("Weight") + theme_presentation(base_size, axis.text.x=element_text(angle=90, hjust=1))
print(p)
#saveFigure(p, filename="HabitatWeights.svg", bg="transparent")
Population density estimates
plotObservations <- function(i, params) {
subdata <- params$data[params$data$year == min(params$data$year) + i - 1,]
subdata$Observed <- subdata$intersections>0
geom_point(mapping=aes(coords.x1, coords.x2, colour=Observed), data=subdata)
}
getPopulationDensity <- function(responses, spatialModels, timeModels, context, withHabitatWeights=FALSE, saveDensityPlots=TRUE) {
populationDensity <- list()
for (i in 1:length(responses)) {
response <- responses[i]
timeModel <- timeModels[i]
spatialModel <- spatialModels[i]
cat("Response = ", response, ", time model = ", timeModel, "\n")
study <- FinlandWTCStudy$new(context=context, response=response, distanceCovariatesModel=~populationDensity+rrday+snow+tday-1, trackSampleInterval=2)
estimates <- if (spatialModel) FinlandSmoothModelSpatioTemporal(study=study)$setModelName("nbinomial", paste("matern", timeModel, sep="-"))$loadEstimates()
else FinlandSmoothModelTemporal(study=study)$setModelName("nbinomial", timeModel)$loadEstimates()
estimates$offsetScale <- 1000^2 # TODO: quickfix
estimates$collectEstimates()
habitatWeightsRaster <- study$loadHabitatWeightsRaster()
data <- cbind(estimates$getUnscaledObservationCoordinates(), estimates$data)
populationDensity[[i]] <- estimates$getPopulationDensity(templateRaster=habitatWeightsRaster, getSD=FALSE)
populationDensity[[i]]$mean$animate(name="PopulationDensity-mean", legend=expression(bold(Individuals / km^2)),
ggfun=plotObservations, ggfunParams=list(data=data))
populationDensity[[i]]$mean$weight(habitatWeightsRaster)
populationDensity[[i]]$mean$animate(name="WeightedPopulationDensity-mean", legend=expression(bold(Individuals / km^2)),
ggfun=plotObservations, ggfunParams=list(data=data))
#populationDensity[[i]] <- study$getPopulationDensity(model=estimates, withHabitatWeights=withHabitatWeights, saveDensityPlots=saveDensityPlots)
#populationDensity[[i]] <- study$getPopulationDensity2(withHabitatWeights=withHabitatWeights)
}
return(populationDensity)
}
timeModels <- c("ar1", "ar1", "rw2")
spatialModels <- c(T,T,F)
weightedPopulationDensity <- getPopulationDensity(responses=responses, spatialModels=spatialModels, timeModels=timeModels, context=context, withHabitatWeights=TRUE, saveDensityPlots=TRUE)
Table 4: Adjusted total population size factors
populationSize <- data.frame()
for (i in 1:length(responses)) {
response <- responses[i]
study <- FinlandWTCStudy$new(context=context, response=response)
x <- weightedPopulationDensity[[i]]$mean$integrate(volume=FinlandPopulationSize$new(study=study))
x$loadValidationData()
y <- x$sizeData
y$response <- study$getPrettyResponse(response)
cat("Response = ", response, ", coef to match the validation data:\n")
print(coef(x$match()))
if (response == "canis.lupus") y$"Adjusted estimated" <- y$Estimated * .75
else if (response == "lynx.lynx") {
y$"Adjusted estimated" <- y$Estimated * .3
y$Estimated <- y$Estimated / 2
}
else if (response == "rangifer.tarandus.fennicus") {
y$"Adjusted estimated" <- y$Estimated / 2
#y$Validation <- with(y, approx(Year, Validation, n=nrow(y)))$y
}
populationSize <- rbind(populationSize, y)
}
Figure 7: Total population size estimates
l <- levels(populationSize$Year)
s <- seq(1, length(l), by=2)
breaks <- l[s]
labels <- l[s]
x <- melt(populationSize, id.vars=c("Year","response"), variable.name="Variable")
p <- ggplot(x, aes(Year, value, group=Variable, colour=Variable, fill=Variable)) + facet_wrap(~response, scales="free_y") +
geom_bar(subset=.(Variable=="Validation"), stat="identity") + geom_line(subset=.(Variable!="Validation"), size=2) +
scale_colour_manual(values=c("steelblue","violetred1","palegreen3")) +
scale_fill_manual(values=c(NA,NA,"palegreen3")) +
scale_x_discrete(breaks=breaks, labels=labels) +
xlab("Year") + ylab("Population size") +
theme_presentation(base_size, axis.text.x=element_text(angle=90, hjust=1)) + theme(legend.position="bottom", strip.text.x=element_blank())
#x <- melt(populationSize, id.vars=c("Year","response"), variable.name="Variable")
#p <- ggplot(x, aes(Year, value, group=Variable, colour=Variable)) + geom_line(size=1) + facet_wrap(~response, scales="free_y") +
# scale_colour_manual(values=c("steelblue","violetred1","steelblue1")) +
# xlab("Year") + ylab("Population size") +
# theme_presentation(base_size, axis.text.x=element_text(angle=90, hjust=1)) + theme(legend.position="bottom")
print(p)
#saveFigure(p, filename="PopulationSize.svg", bg="transparent")
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(knitr) single_width <- 2.75591 * 2 double_width <- single_width*2 + 0.5 * 2 single_height <- single_width double_height <- double_width #opts_knit$set(base.dir="~/.pandoc") opts_chunk$set(echo=FALSE, warning=FALSE, message=FALSE, tidy=FALSE, fig.width=single_width, fig.height=single_height, cache=TRUE, cache.path="~/tmp/cache/") base_size <- 16
library(parallel) library(doMC) registerDoMC(cores=detectCores()) library(STREM) source("~/git/STREM/setup/WTC-Boot.R") library(sp) library(ggplot2) library(rasterVis) library(reshape2) library(plyr) context <- Context$new(resultDataDirectory=wd.data.results, processedDataDirectory=wd.data.processed, rawDataDirectory=wd.data.raw, scratchDirectory=wd.scratch, figuresDirectory=wd.figures) study <- FinlandWTCStudy$new(context=context, response="canis.lupus") study$studyArea$loadBoundary(thin=F, tolerance=0.001) boundaryDF <- ggplot2::fortify(study$studyArea$boundary) responses <- c("canis.lupus", "lynx.lynx", "rangifer.tarandus.fennicus")
Figure 1: Survey routes
surveyRoutes <- study$loadSurveyRoutes(findLengths=FALSE) surveyRoutesSP <- SpatialLinesDataFrame(surveyRoutes$surveyRoutes, data=data.frame(x=1:length(surveyRoutes$surveyRoutes)), match.ID=FALSE) surveyRoutesDF <- ggplot2::fortify(surveyRoutesSP) p <- ggplot(boundaryDF, aes(long, lat, group=group)) + geom_polygon(colour="black", fill=NA) + geom_polygon(data=surveyRoutesDF, aes(long, lat, group=group), colour="blue", fill=NA) + coord_equal() + theme_raster() print(p) #saveFigure(p, filename="SurveyRoutes.svg", bg="transparent")
Number of counted survey routes
intersections <- study$loadIntersections(predictDistances=FALSE) x <- ddply(as.data.frame(intersections$intersections), .(year), nrow) s <- sd(x$V1) m <- mean(x$V1) cat("mean = ", m, "\n") cat("sd = ", s, "\n") cat("range = ", m+c(-s,s), "\n")
% of intersections D > 1
intersections <- FinlandWTCIntersections$new(study=study, maxDuration=1) intersections$saveIntersections()
Figure 2: Crossing rates (D > 1 removed)
getIntersectionsRate <- function(responses, context) { library(plyr) x <- ldply(responses, function(response, context) { study <- FinlandWTCStudy$new(context=context, response=response) intersections <- study$loadIntersections(predictDistances=FALSE) x <- ddply(intersections$intersections@data, .(year), function(x, response) { data.frame(intersections=sum(x$intersections / (x$length/1000) / x$duration) / nrow(x), response=study$getPrettyResponse(response)) }, response=response) return(x) }, context=context) return(x) } intersectionsRate <- getIntersectionsRate(responses, context) p <- ggplot(intersectionsRate, aes(year, intersections)) + geom_line(size=1) + facet_grid(~response, scales="free_y") + xlab("Year") + ylab("Crossings rate") + theme_presentation(base_size) print(p) #saveFigure(p, filename="IntersectionsRate.svg", bg="transparent")
GPS figures
getTracks <- function(responses, context) { library(plyr) library(ggplot2) library(adehabitatLT) l_ply(responses, function(response, context) { study <- FinlandWTCStudy$new(context=context, response=response) tracks <- FinlandWTCTracks$new(study=study) tracks$loadTracks() message("Response = ", response) message("Number of individuals in the cleaned data = ", length(unique(adehabitatLT::id(tracks$tracks)))) d <- as.POSIXlt(ld(tracks$tracks)$date) message("Years in the cleaned data = ", paste(sort(unique(d$year+1900)), collapse=","), " (", length(unique(d$year)), ")") message("Number of samples in the cleaned data = ", length(d)) message("Number of days in the cleaned data = ", length(unique(paste(d$year, d$yday)))) }, context=context) } getTracks(responses=responses, context=context)
Figure 3: GPS movements, juvenile is no more than 1 years old individual
library(plyr) response <- "canis.lupus" study <- FinlandWTCStudy$new(context=context, response=response) tracks <- FinlandWTCTracks$new(study=study) tracks$loadTracks() tracks$loadMetadata() xtabs(sex~born, tracks$metadata) library(reshape2) dcast(tracks$metadata, .~sex) dcast(tracks$metadata, .~born) dcast(tracks$metadata, .~dispersing)$"TRUE" dcast(tracks$metadata, .~homerange)$"TRUE" dcast(tracks$metadata, .~alpha) x <- ld(tracks$tracks)$dt/3600 sort(x[x<=24]) sort(x[x<=24])*3600 tracksSP <- tracks$getSpatialLines(variables=.(id,year,burst)) x <- tracksSP@data x$juvenile <- if ("born" %in% colnames(x)) factor(ifelse(x$year - x$born <= 1, "juvenile", "adult")) else factor("unknown") x$juvenile <- factor(x$juvenile, levels=c(levels(x$juvenile), "unknown")) x$juvenile[is.na(x$juvenile)] <- "unknown" dcast(x, .~juvenile)
getTracks <- function(responses, context) { library(plyr) library(ggplot2) x <- ldply(responses, function(response, context) { study <- FinlandWTCStudy$new(context=context, response=response) tracks <- FinlandWTCTracks$new(study=study) tracks$loadTracks() tracksSP <- tracks$getSpatialLines(variables=.(id,year,burst)) x <- tracksSP@data if (!"sex" %in% colnames(x)) x$sex <- "unknown" x$juvenile <- if ("born" %in% colnames(x)) factor(ifelse(x$year - x$born <= 1, "juvenile", "adult")) else factor("unknown") x$juvenile <- factor(x$juvenile, levels=c(levels(x$juvenile), "unknown")) x$juvenile[is.na(x$juvenile)] <- "unknown" x$id <- sapply(tracksSP@lines, function(x) x@ID) x$sexjuv <- factor(paste(x$sex, x$juvenile, sep=", ")) tracksDF <- ggplot2::fortify(tracksSP) tracksDF <- plyr::join(tracksDF, x, by="id") tracksDF$response <- study$getPrettyResponse(response) return(tracksDF) }, context=context) return(x) } tracks <- getTracks(responses=responses, context=context) #levels(tracks$sexjuv) p <- ggplot(tracks, aes(long, lat, group=group)) + geom_polygon(data=boundaryDF, aes(long, lat, group=group), colour="darkgray", fill=NA) + geom_path(aes(colour=sexjuv, linetype=sexjuv), alpha=.7) + scale_color_manual("Sex, age", values=c("darkred","red","red","darkblue","blue","blue","black")) + scale_linetype_manual("Sex, age", values=c("solid","dotted","solid","solid","dotted","solid","solid")) + facet_grid(~response) + theme_raster(base_size, aspect.ratio=1, legend.position=c(0.5,-0.1), legend.direction="horizontal", legend.position="bottom") + guides(colour=guide_legend(nrow=2)) + coord_cartesian(ylim=range(tracks$lat)+c(-1,1)*50000) print(p) #saveFigure(p, filename="GPS.svg", bg="transparent")
Figure 4: Straight-line-distance-error
ss <- SimulationStudy$new(response="Intensive")$setup(context=context) tracks <- ss$loadTracks(iteration=as.integer(1)) tracks$tracks <- subset(tracks$tracks, id == 1 & year == 2001 & month == 1 & day %in% c(2,3)) intervals <- MovementSampleIntervals$new(study=ss) thinnedTracks <- intervals$getThinnedTracksSampleIntervals(tracks=tracks) plotTracks <- data.frame() for (i in seq(1, 11, by=5)) { x <- thinnedTracks$tracksList[[i]]$tracks x$thin <- i x <- subset(x, date >= as.POSIXct("2001-01-03 00:00:00") & date <= as.POSIXct("2001-01-03 03:40:00")) plotTracks <- rbind(plotTracks, x) } plotTracks$thin <- factor(plotTracks$thin) plotTracks$dt <- factor(plotTracks$dt / 60)
segment <- aes(x=x, y=y, xend=c(tail(x, n=-1), NA), yend=c(tail(y, n=-1), NA), group=dt, colour=dt) ending <- arrow(length=unit(0.7, "cm")) p <- ggplot(plotTracks, mapping=segment) + theme_raster(base_size, legend.position=c(0.5,0.85), legend.background=element_rect(color="grey")) + #guides(colour=guide_legend(title=expression(paste(Delta, "t (min)")))) + guides(colour=guide_legend(title="Sampling interval (min)")) + scale_colour_manual(values=c("steelblue1","steelblue3","steelblue4")) + #ggtitle("Straight-line distance error") + geom_path(size=1, arrow=ending) + geom_point(size=2) #geom_segment(size=2, arrow=ending) # doesn't work with grouping print(p) #saveFigure(p, filename="StraightLineDistanceError.svg", bg="transparent")
Raw GPS sample intervals
intervals <- llply(responses, function(response, context) { study <- FinlandWTCStudy$new(context=context, response=response) tracks <- FinlandWTCTracks$new(study=study) tracks$loadTracks() intervals <- MovementSampleIntervals$new(study=study) intervals$loadSampleIntervals() return(table(intervals$intervals$intervalH)) }, context=context) intervals
Figure 5: Power law fits for distance - sampling interval
distances <- ldply(responses, function(response, context) { study <- FinlandWTCStudy$new(context=context, response=response) intervals <- study$loadSampleIntervals() x <- intervals$intervals x$response <- study$getPrettyResponse(response=response) return(x) }, context=context) intervals <- ldply(responses, function(response, context) { study <- FinlandWTCStudy$new(context=context, response=response) intervals <- study$loadSampleIntervals() distances <- intervals$getDistanceCurve() distances$response <- study$getPrettyResponse(response=response) return(distances) }, context=context) p <- ggplot(distances, aes(intervalH, distanceKm)) + geom_point(alpha=.3) + ylab("Distance / day (km)") + xlab("Sampling interval (h)") + geom_line(data=intervals, aes(intervalH, distanceKm), size=1) + geom_smooth(data=intervals, aes(ymin=distanceKm25, ymax=distanceKm975), stat="identity") + ylim(c(0, max(distances$distanceKm))) + facet_grid(~response) + theme_presentation(base_size) print(p) #saveFigure(p, filename=paste("DistanceCorrection-", response, ".svg", sep=""), bg="transparent")
Table 1: Observed distances - estimated distances
observedDistances <- ldply(responses, function(response, context) { study <- FinlandWTCStudy$new(context=context, response=response) tracks <- study$loadTracks() intervals <- tracks$getSampleIntervals() return(data.frame(distance=intervals$intervals$dist, logDistance=log(intervals$intervals$dist), response=study$getPrettyResponse(response))) }, context=context) predictedDistances <- ldply(responses, function(response, context) { study <- FinlandWTCStudy$new(context=context, response=response, distanceCovariatesModel=~populationDensity+rrday+snow+tday-1, trackSampleInterval=2) intersections <- study$loadIntersections() intersections$predictDistances() distances <- intersections$intersections$distance / 1000 return(data.frame(distance=distances, logDistance=log(distances), response=study$getPrettyResponse(response))) }, context=context)
cat("Observed distances (km/day):\n") ddply(observedDistances, .(response), summarize, mean=mean(distance), sd=sd(distance), min=min(distance), max=max(distance)) cat("Predicted distances (km/day):\n") ddply(predictedDistances, .(response), summarize, mean=mean(distance), sd=sd(distance), min=min(distance), max=max(distance)) cat("Prediction interval = 2h")
Distances from literature
- Wolf: 7-25 km/day
- Lynx:
- FFR:
Table 2: Power law fit estimates for distance - sampling interval
- fixed_effect = human population density
- fixed_effect1 = precipitation
- fixed_effect2 = snow cover
- fixed_effect3 = temperature
# Fixed effects: populationDensity + rrday + snow + tday intervals <- l_ply(responses, function(response, context) { study <- FinlandWTCStudy$new(context=context, response=response) intervals <- study$loadSampleIntervals() cat("response = ", study$getPrettyResponse(response=response), "\n") intervals$estimatedValuesSummary() cat("\n\n") }, context=context)
Table 3: Crude distance estimates
responses <- c("canis.lupus","lynx.lynx","rangifer.tarandus.fennicus") for (response in responses) { context <- Context$new(resultDataDirectory=wd.data.results, processedDataDirectory=wd.data.processed, rawDataDirectory=wd.data.raw, scratchDirectory=wd.scratch, figuresDirectory=wd.figures) study <- FinlandWTCStudy$new(context=context, response=response, distanceCovariatesModel=~populationDensity+rrday+snow+tday-1, trackSampleInterval=2) populationSize <- FinlandPopulationSize$new(study=study) populationSize$loadValidationData() intersections <- study$loadIntersections() A <- 338424 # km^2 M <- ddply(intersections$getData(), .(year), function(x) sum(x$length))$V1 / 1000 N <- populationSize$sizeData$Validation x <- ddply(intersections$getData(), .(year), function(x) sum(x$intersections))$V1 C <- 2/pi * A / M l <- C * x / N cat(response, ": crude distance mean = ", mean(l, na.rm=T), " median = ", median(l, na.rm=T), " sd = ", sd(l, na.rm=T), " km/d\n") h <- study$loadHabitatWeightsRaster() h.masked <- mask(h, study$studyArea$boundary) n.cell <- length(h.masked) - freq(h.masked, value=NA) A.discrete <- n.cell * prod(res(h.masked)) / 1000^2 h.cell <- cellStats(h.masked, sum) A.h <- h.cell * prod(res(h.masked)) / 1000^2 * A / A.discrete # A.discrete is a little bit smaller than A due to the discretization, so apply a correction C <- 2/pi * A.h / M l <- C * x / N cat(response, ": crude weighted distance mean = ", mean(l, na.rm=T), " median = ", median(l, na.rm=T), " sd = ", sd(l, na.rm=T), " km/d\n") cat(response, ": observed number of intersections:\n") print(1/C * l) cat("Area correction = ", A/A.discrete, "\n") }
Figure 6: Habitat weights
getHabitatWeights <- function(responses, context) { library(plyr) x <- ldply(responses, function(response, context) { study <- FinlandWTCStudy$new(context=context, response=response) weights <- study$loadHabitatWeights()$getHabitatSelectionWeights() n <- names(weights) print(n) weights <- data.frame(habitat=factor(n, levels=n), weights) weights$response <- study$getPrettyResponse(response) return(weights) }, context=context) return(x) } weights <- getHabitatWeights(responses=responses, context=context) p <- ggplot(weights, aes(habitat, weights, fill=habitat)) + geom_bar(stat="identity") + facet_grid(~response) + scale_fill_manual(values=c("#beaed4","#ffff99","#7fc97f","#fdc086","#386cb0")) + xlab("") + ylab("Weight") + theme_presentation(base_size, axis.text.x=element_text(angle=90, hjust=1)) print(p) #saveFigure(p, filename="HabitatWeights.svg", bg="transparent")
Population density estimates
plotObservations <- function(i, params) { subdata <- params$data[params$data$year == min(params$data$year) + i - 1,] subdata$Observed <- subdata$intersections>0 geom_point(mapping=aes(coords.x1, coords.x2, colour=Observed), data=subdata) } getPopulationDensity <- function(responses, spatialModels, timeModels, context, withHabitatWeights=FALSE, saveDensityPlots=TRUE) { populationDensity <- list() for (i in 1:length(responses)) { response <- responses[i] timeModel <- timeModels[i] spatialModel <- spatialModels[i] cat("Response = ", response, ", time model = ", timeModel, "\n") study <- FinlandWTCStudy$new(context=context, response=response, distanceCovariatesModel=~populationDensity+rrday+snow+tday-1, trackSampleInterval=2) estimates <- if (spatialModel) FinlandSmoothModelSpatioTemporal(study=study)$setModelName("nbinomial", paste("matern", timeModel, sep="-"))$loadEstimates() else FinlandSmoothModelTemporal(study=study)$setModelName("nbinomial", timeModel)$loadEstimates() estimates$offsetScale <- 1000^2 # TODO: quickfix estimates$collectEstimates() habitatWeightsRaster <- study$loadHabitatWeightsRaster() data <- cbind(estimates$getUnscaledObservationCoordinates(), estimates$data) populationDensity[[i]] <- estimates$getPopulationDensity(templateRaster=habitatWeightsRaster, getSD=FALSE) populationDensity[[i]]$mean$animate(name="PopulationDensity-mean", legend=expression(bold(Individuals / km^2)), ggfun=plotObservations, ggfunParams=list(data=data)) populationDensity[[i]]$mean$weight(habitatWeightsRaster) populationDensity[[i]]$mean$animate(name="WeightedPopulationDensity-mean", legend=expression(bold(Individuals / km^2)), ggfun=plotObservations, ggfunParams=list(data=data)) #populationDensity[[i]] <- study$getPopulationDensity(model=estimates, withHabitatWeights=withHabitatWeights, saveDensityPlots=saveDensityPlots) #populationDensity[[i]] <- study$getPopulationDensity2(withHabitatWeights=withHabitatWeights) } return(populationDensity) } timeModels <- c("ar1", "ar1", "rw2") spatialModels <- c(T,T,F) weightedPopulationDensity <- getPopulationDensity(responses=responses, spatialModels=spatialModels, timeModels=timeModels, context=context, withHabitatWeights=TRUE, saveDensityPlots=TRUE)
Table 4: Adjusted total population size factors
populationSize <- data.frame() for (i in 1:length(responses)) { response <- responses[i] study <- FinlandWTCStudy$new(context=context, response=response) x <- weightedPopulationDensity[[i]]$mean$integrate(volume=FinlandPopulationSize$new(study=study)) x$loadValidationData() y <- x$sizeData y$response <- study$getPrettyResponse(response) cat("Response = ", response, ", coef to match the validation data:\n") print(coef(x$match())) if (response == "canis.lupus") y$"Adjusted estimated" <- y$Estimated * .75 else if (response == "lynx.lynx") { y$"Adjusted estimated" <- y$Estimated * .3 y$Estimated <- y$Estimated / 2 } else if (response == "rangifer.tarandus.fennicus") { y$"Adjusted estimated" <- y$Estimated / 2 #y$Validation <- with(y, approx(Year, Validation, n=nrow(y)))$y } populationSize <- rbind(populationSize, y) }
Figure 7: Total population size estimates
l <- levels(populationSize$Year) s <- seq(1, length(l), by=2) breaks <- l[s] labels <- l[s] x <- melt(populationSize, id.vars=c("Year","response"), variable.name="Variable") p <- ggplot(x, aes(Year, value, group=Variable, colour=Variable, fill=Variable)) + facet_wrap(~response, scales="free_y") + geom_bar(subset=.(Variable=="Validation"), stat="identity") + geom_line(subset=.(Variable!="Validation"), size=2) + scale_colour_manual(values=c("steelblue","violetred1","palegreen3")) + scale_fill_manual(values=c(NA,NA,"palegreen3")) + scale_x_discrete(breaks=breaks, labels=labels) + xlab("Year") + ylab("Population size") + theme_presentation(base_size, axis.text.x=element_text(angle=90, hjust=1)) + theme(legend.position="bottom", strip.text.x=element_blank()) #x <- melt(populationSize, id.vars=c("Year","response"), variable.name="Variable") #p <- ggplot(x, aes(Year, value, group=Variable, colour=Variable)) + geom_line(size=1) + facet_wrap(~response, scales="free_y") + # scale_colour_manual(values=c("steelblue","violetred1","steelblue1")) + # xlab("Year") + ylab("Population size") + # theme_presentation(base_size, axis.text.x=element_text(angle=90, hjust=1)) + theme(legend.position="bottom") print(p) #saveFigure(p, filename="PopulationSize.svg", bg="transparent")
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.