Make GAM term plots that extend outside the range of the data to show how bad extrapolation is for Laura M

extra_plot.R

load("best_model.Rdata")
library(dsm)
library(raster)

# lazily get the plot data for the rug plot
plotdat <- plot(M)

# load the raster and mudge it into the format I want
dists <- stack("NA_Shore_Dist_10km_mean_10km.img")
dists <- as.data.frame(dists)
dists <- dists$NA_Shore_Dist_10km_mean_10km
dists <- dists[!is.na(dists)]

# make a data frame that has the max to min value of dist to shore
dat <- data.frame(DistToShore = seq(min(dists), max(dists), len=1000),
                  area_km2    = 10000*10000)

## make our cool plot

# predict extrapolated values
extr_term <- predict(M ,dat, type="terms", se=TRUE)

# make predictions for the data "inside"
inside_dat <- dat
inside_dat <- inside_dat[inside_dat$DistToShore > 0 &
                         inside_dat$DistToShore < 6e5,]
inside_term <- predict(M, inside_dat, type="terms", se=TRUE)

# plot the "bad" extrapolation in red
plot(dat$DistToShore, extr_term$fit, type="l", col="red",
     xlab="DistToShore", ylab="s(DistToShore)")
# overlay in green the "good" bit where the data is
lines(inside_dat$DistToShore, inside_term$fit, col="green")

# what about the standard errors?
lines(dat$DistToShore, extr_term$fit+2*extr_term$se.fit, lty=2)
lines(dat$DistToShore, extr_term$fit-2*extr_term$se.fit, lty=2)
rug(plotdat[[1]]$raw)