## ----eval=FALSE---------------------------------------------------------- ## install.packages(c("dsm", "Distance", "knitr", "captioner", "ggplot2", "rgdal", ## "maptools", "plyr", "tweedie")) ## ----echo=FALSE, warning=FALSE------------------------------------------- library(knitr) library(captioner) # make a numbering for figures fig_nums <- captioner(prefix="Figure") fig_c <- function(p) fig_nums(p, display="cite") ## ----loadlibraries, comment=NA------------------------------------------- library(dsm) library(ggplot2) # plotting options gg.opts <- theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(), panel.background=element_blank()) # make the results reproducible set.seed(11123) ## ----loaddata------------------------------------------------------------ data(mexdolphins) attach(mexdolphins) ## ----head-segdata-------------------------------------------------------- kable(segdata[1:3, ], digits=2) ## ----head-distdata------------------------------------------------------- kable(distdata[1:3, ], digits=2) ## ----head-obsdata-------------------------------------------------------- kable(obsdata[1:3, ], digits=2) ## ----head-preddata------------------------------------------------------- kable(preddata[1:3, ], digits=2) ## ----projectsurvey, results='hide', warning=FALSE------------------------ library(rgdal, quietly = TRUE) library(maptools) # tell R that the survey.area object is currently in lat/long proj4string(survey.area) <- CRS("+proj=longlat +datum=WGS84") # proj 4 string # using http://spatialreference.org/ref/esri/north-america-lambert-conformal-conic/ lcc_proj4 <- CRS("+proj=lcc +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs ") # project using LCC survey.area <- spTransform(survey.area, CRSobj=lcc_proj4) # simplify the object survey.area <- data.frame(survey.area@polygons[[1]]@Polygons[[1]]@coords) names(survey.area) <- c("x", "y") ## ----echo=FALSE---------------------------------------------------------- fig_nums(name="areawithtransects", caption="The survey area with transect lines.", display=FALSE) ## ----areawithtransects, fig.cap="", fig.height=4------------------------- p <- qplot(data=survey.area, x=x, y=y, geom="polygon",fill=I("lightblue"), ylab="y", xlab="x", alpha=I(0.7)) p <- p + coord_equal() p <- p + geom_line(aes(x,y,group=Transect.Label),data=segdata) p <- p + gg.opts print(p) ## ----echo=FALSE---------------------------------------------------------- fig_nums(name="projection-compare", caption="Comparison between unprojected latitude and longitude (left) and the prediction grid projected using the North American Lambert Conformal Conic projection (right)", display=FALSE) ## ----projection-compare, fig.cap="", fig.height=3------------------------ par(mfrow=c(1,2)) # put pred.polys into lat/long pred_latlong <- spTransform(pred.polys,CRSobj=CRS("+proj=longlat +datum=WGS84")) # plot latlong plot(pred_latlong, xlab="Longitude", ylab="Latitude") axis(1); axis(2); box() # plot as projected plot(pred.polys, xlab="Northing", ylab="Easting") axis(1); axis(2); box() ## ----ggpoly-------------------------------------------------------------- # given the argument fill (the covariate vector to use as the fill) and a name, # return a geom_polygon object # fill must be in the same order as the polygon data library(plyr) grid_plot_obj <- function(fill, name, sp){ # what was the data supplied? names(fill) <- NULL row.names(fill) <- NULL data <- data.frame(fill) names(data) <- name spdf <- SpatialPolygonsDataFrame(sp, data) spdf@data$id <- rownames(spdf@data) spdf.points <- fortify(spdf, region="id") spdf.df <- join(spdf.points, spdf@data, by="id") # seems to store the x/y even when projected as labelled as "long" and "lat" spdf.df$x <- spdf.df$long spdf.df$y <- spdf.df$lat geom_polygon(aes_string(x="x",y="y",fill=name, group="group"), data=spdf.df) } ## ----echo=FALSE---------------------------------------------------------- fig_nums(name="EDA-plots", caption="Exploratory plot of the distance sampling data. Top row, left to right: histograms of distance and cluster size; bottom row: plot of distance against cluster size and plot of distances against Beaufort sea state.", display=FALSE) ## ----EDA-plots, fig.height=7, fig.width=7, fig.cap="", results="hide"---- par(mfrow=c(2,2)) # histograms hist(distdata$distance,main="",xlab="Distance (m)") hist(distdata$size,main="",xlab="Cluster size") # plots of distance vs. cluster size plot(distdata$distance, distdata$size, main="", xlab="Distance (m)", ylab="Group size", pch=19, cex=0.5, col=gray(0.7)) # lm fit l.dat <- data.frame(distance=seq(0,8000,len=1000)) lo <- lm(size~distance, data=distdata) lines(l.dat$distance, as.vector(predict(lo,l.dat))) plot(distdata$distance,distdata$beaufort, main="", xlab="Distance (m)", ylab="Beaufort sea state", pch=19, cex=0.5, col=gray(0.7)) ## ----echo=FALSE---------------------------------------------------------- fig_nums(name="spatialEDA", caption="Plot of depth values over the survey area with transects and observations overlaid. Point size is proportional to the group size for each observation.", display=FALSE) ## ----spatialEDA, fig.cap=""---------------------------------------------- p <- ggplot() + grid_plot_obj(preddata$depth, "Depth", pred.polys) + coord_equal() p <- p + labs(fill="Depth",x="x",y="y",size="Group size") p <- p + geom_line(aes(x, y, group=Transect.Label), data=segdata) p <- p + geom_point(aes(x, y, size=size), data=distdata, colour="red",alpha=I(0.7)) p <- p + gg.opts print(p) ## ----loadDistance-------------------------------------------------------- library(Distance) ## ----hrmodel, comment=NA------------------------------------------------- detfc.hr.null<-ds(distdata, max(distdata$distance), key="hr", adjustment=NULL) ## ----hrmodelsummary, comment=NA------------------------------------------ summary(detfc.hr.null) ## ----echo=FALSE---------------------------------------------------------- fig_nums(name="hr-detfct", caption="Plot of the fitted detection function (left) and goodness of fit plot (right) for the hazard-rate model.", display=FALSE) ## ----hr-detfct, fig.cap="", fig.width=9, fig.height=6, results="hide"---- par(mfrow=c(1,2)) plot(detfc.hr.null, showpoints=FALSE, pl.den=0, lwd=2) ddf.gof(detfc.hr.null$ddf) ## ----hrcovdf, message=FALSE, cache=TRUE, warning=FALSE------------------- detfc.hr.beau<-ds(distdata, max(distdata$distance), formula=~as.factor(beaufort), key="hr", adjustment=NULL) ## ----hrcovdfsummary, comment=NA------------------------------------------ summary(detfc.hr.beau) ## ------------------------------------------------------------------------ dsm.xy <- dsm(count~s(x,y), detfc.hr.null, segdata, obsdata, method="REML") ## ---- comment=NA--------------------------------------------------------- summary(dsm.xy) ## ----echo=FALSE---------------------------------------------------------- fig_nums(name="visgam1", caption="Plot of the spatial smooth in `dsm.xy`, values are relative abundances. White/yellow indicates high values, red low indicates low values.", display=FALSE) ## ----visgam1, fig.cap=""------------------------------------------------- vis.gam(dsm.xy, plot.type="contour", view=c("x","y"), asp=1, type="response", contour.col="black", n.grid=100) ## ----depthmodel, comment=NA---------------------------------------------- dsm.xy.depth <- dsm(count~s(x,y,k=10) + s(depth,k=20), detfc.hr.null, segdata, obsdata, method="REML") summary(dsm.xy.depth) ## ----echo=FALSE---------------------------------------------------------- fig_nums(name="dsm-xy-depth-depth", caption="Plot of the smooth of depth in `dsm.xy.depth`. This \"hockey stick\" shaped smooth indicates that there is little predictive value in depth beyond 500m.", display=FALSE) ## ----dsm-xy-depth-depth, fig.cap=""-------------------------------------- plot(dsm.xy.depth, select=2) ## ----modelcomp, cache=FALSE---------------------------------------------- # make a data.frame to print out mod_results <- data.frame("Model name" = c("`dsm.xy`", "`dsm.xy.depth`"), "Description" = c("Bivariate smooth of location, quasipoisson", "Bivariate smooth of location, smooth of depth, quasipoisson"), "Deviance explained" = c(unlist(lapply(list(dsm.xy, dsm.xy.depth), function(x){paste0(round(summary(x)$dev.expl*100,2),"%")})))) ## ----results-table, results='asis'--------------------------------------- kable(mod_results, col.names=c("Model name", "Description", "Deviance explained")) ## ------------------------------------------------------------------------ dsm.xy.pred <- predict(dsm.xy, preddata, preddata$area) ## ----echo=FALSE---------------------------------------------------------- fig_nums(name="dsm-xy-preds", caption="Predicted density surface for `dsm.xy`.", display=FALSE) ## ----dsm-xy-preds, fig.cap=""-------------------------------------------- p <- ggplot() + grid_plot_obj(dsm.xy.pred, "Abundance", pred.polys) + coord_equal() +gg.opts p <- p + geom_path(aes(x=x, y=y),data=survey.area) p <- p + labs(fill="Abundance") print(p) ## ----dsm.xy-abund, comment=NA-------------------------------------------- sum(dsm.xy.pred) ## ----echo=FALSE---------------------------------------------------------- fig_nums(name="dsm-xy-depth-preds", caption="Predicted density surface for `dsm.xy.depth`.", display=FALSE) ## ----dsm-xy-depth-preds, fig.cap=""-------------------------------------- dsm.xy.depth.pred <- predict(dsm.xy.depth, preddata, preddata$area) p <- ggplot() + grid_plot_obj(dsm.xy.depth.pred, "Abundance", pred.polys) + coord_equal() +gg.opts p <- p + geom_path(aes(x=x, y=y), data=survey.area) p <- p + labs(fill="Abundance") print(p) ## ----dsm.xy.depth-abund, comment=NA-------------------------------------- sum(dsm.xy.depth.pred) ## ----dsm.xy-varprop, cache=FALSE----------------------------------------- preddata.varprop <- split(preddata, 1:nrow(preddata)) dsm.xy.varprop <- dsm.var.prop(dsm.xy, pred.data=preddata.varprop, off.set=preddata$area) ## ---- comment=NA--------------------------------------------------------- summary(dsm.xy.varprop) ## ----echo=FALSE---------------------------------------------------------- fig_nums(name="dsm-xyvarplot", caption="Plot of the coefficient of variation for the study area with transect lines and observations overlaid. Note the increase in CV away from the transect lines.", display=FALSE) ## ----dsm-xyvarplot, fig.cap=""------------------------------------------- p <- ggplot() + grid_plot_obj(sqrt(dsm.xy.varprop$pred.var)/unlist(dsm.xy.varprop$pred), "CV", pred.polys) + coord_equal() +gg.opts p <- p + geom_path(aes(x=x, y=y),data=survey.area) p <- p + geom_line(aes(x, y, group=Transect.Label), data=segdata) print(p) ## ----dsm.xy.depth-varprop, cache=FALSE, comment=NA----------------------- dsm.xy.depth.varprop <- dsm.var.prop(dsm.xy.depth, pred.data=preddata.varprop, off.set=preddata$area) summary(dsm.xy.depth.varprop)