Uncertainty estimation and plotting in dsm

hi_len.R

# for LJT

# this wasn't as simple as I thought, it was the summary() method that does the smart thing
# so there is a bit of split() nonsense to deal with...

library(Distance)
library(dsm)
library(ggplot2)


## fit model, do normal stuff...

# load the Gulf of Mexico dolphin data (see ?mexdolphins)
data(mexdolphins)

# fit a detection function and look at the summary
hr.model <- ds(distdata, max(distdata$distance),
               key = "hn", adjustment = NULL)

# fit a simple smooth of x and y
mod1 <- dsm(count~ s(depth), hr.model, segdata, obsdata)

# setup data
preddata$width <- preddata$height <- sqrt(preddata$area)
pred.new <- split(preddata, 1:nrow(preddata))


## var gam
# calculate
mod1.var <- dsm.var.gam(mod1, pred.new, off.set=preddata$area)

preddata$Nhat <- unlist(mod1.var$pred)
preddata$CV <- sqrt(mod1.var$pred.var)/unlist(mod1.var$pred)

# plot abundance (ggplot2)
p <- ggplot(preddata) +
  geom_tile(aes(x=x, y=y, fill=Nhat), width=preddata$width, height=preddata$height) +
  labs(fill="Density") +
  theme_minimal() +
  coord_equal()
print(p)
# plot uncertainty (using ggplot2)
p <- ggplot(preddata) +
  geom_tile(aes(x=x, y=y, fill=CV), width=preddata$width, height=preddata$height) +
  labs(fill="CV") +
  theme_minimal() +
  coord_equal()
print(p)


## var gam
# calculate
mod1.var <- dsm.var.prop(mod1, pred.new, off.set=preddata$area)

preddata$Nhat <- unlist(mod1.var$pred)
preddata$CV <- sqrt(mod1.var$pred.var)/unlist(mod1.var$pred)

# plot abundance (ggplot2)
p <- ggplot(preddata) +
  geom_tile(aes(x=x, y=y, fill=Nhat), width=preddata$width, height=preddata$height) +
  labs(fill="Density") +
  theme_minimal() +
  coord_equal()
print(p)
# plot uncertainty (using ggplot2)
p <- ggplot(preddata) +
  geom_tile(aes(x=x, y=y, fill=CV), width=preddata$width, height=preddata$height) +
  labs(fill="CV") +
  theme_minimal() +
  coord_equal()
print(p)


## moving block

mod1.var <- dsm.var.movblk(mod1, preddata, off.set=preddata$area, n.boot=100, block.size=3)
n <- mod1.var$n.boot
cell.se <- sqrt((mod1.var$short.var$sumx.sq/(n - 1)) - ((mod1.var$short.var$sumx/n)^2 *(n/(n - 1))))
preddata$Nhat <- mod1.var$short.var$sumx/n
preddata$CV <- cell.se/preddata$Nhat

# plot abundance (ggplot2)
p <- ggplot(preddata) +
  geom_tile(aes(x=x, y=y, fill=Nhat), width=preddata$width, height=preddata$height) +
  labs(fill="Density") +
  theme_minimal() +
  coord_equal()
print(p)
# plot uncertainty (using ggplot2)
p <- ggplot(preddata) +
  geom_tile(aes(x=x, y=y, fill=CV), width=preddata$width, height=preddata$height) +
  labs(fill="CV") +
  theme_minimal() +
  coord_equal()
print(p)