Spatially explicit models for distance sampling data
Practical “solution” 4: Model checking
David L Miller
Aims
By the end of this practical, you should feel comfortable:
- Looking at how changing
kchanges smooths - Plotting residuals vs. covariates
- Fitting models to the residuals
- Using
gam.checkandrqgam.check - Using
obs_exp
Load data and packages
library(Distance)## Loading required package: mrds## This is mrds 2.2.3
## Built: R 4.0.2; ; 2020-08-01 10:33:56 UTC; unix##
## Attaching package: 'Distance'## The following object is masked from 'package:mrds':
##
## create.binslibrary(dsm)## Loading required package: mgcv## Loading required package: nlme## This is mgcv 1.8-33. For overview type 'help("mgcv-package")'.## Loading required package: numDeriv## This is dsm 2.3.0
## Built: R 4.0.2; ; 2020-07-16 23:56:50 UTC; unixlibrary(ggplot2)
library(knitr)Load the data and the fitted dsm objects from the previous exercises:
load("spermwhale.RData")
load("dsms.RData")
load("df-models.RData")We’re going to check the DSMs that we fitted in the previous practical, then save those that we think are good!
Pre-model fitting
As we did in the previous exercise, we must remove the observations from the spatial data that we excluded when we fitted the detection function, i.e. those observations at distances greater than the truncation.
obs <- obs[obs$distance <= df_hn$ddf$meta.data$width, ]Looking at how changing k changes smooths
First checking that k is big enough, we should really do this during model fitting, but we’ve separated this up for the practical exercises.
First look at the text output of gam.check, are the values of k' for your models close to the edf in the outputted table. Here’s a silly example where I’ve deliberately set k too small:
dsm_k_check_eg <- dsm(count ~ s(Depth, k=4),
df_hn, segs, obs,
family=tw())
gam.check(dsm_k_check_eg)##
## Method: REML Optimizer: outer newton
## full convergence after 8 iterations.
## Gradient range [-1.66079e-07,1.171305e-07]
## (score 392.2823 & scale 5.216272).
## Hessian positive definite, eigenvalue range [0.8702242,302.9225].
## Model rank = 4 / 4
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(Depth) 3.0 2.9 0.78 0.17Generally if the EDF is close to the value of k you supplied, it is worth doubling k and refitting to see what happens. You can always switch back to the smaller k if there is little difference. The ?choose.k manual page can offer some guidance.
Continuing with that example, if we double k:
dsm_k_check_eg <- dsm(count ~ s(Depth, k=8),
df_hn, segs, obs,
family=tw())
gam.check(dsm_k_check_eg)##
## Method: REML Optimizer: outer newton
## full convergence after 8 iterations.
## Gradient range [-5.675949e-08,4.752562e-08]
## (score 392.2128 & scale 5.216315).
## Hessian positive definite, eigenvalue range [1.313321,300.6316].
## Model rank = 8 / 8
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(Depth) 7.00 4.47 0.78 0.26We get something much more reasonable. Doubling again
dsm_k_check_eg <- dsm(count ~ s(Depth, k=16),
df_hn, segs, obs,
family=tw())
gam.check(dsm_k_check_eg)##
## Method: REML Optimizer: outer newton
## full convergence after 7 iterations.
## Gradient range [-4.245788e-07,3.205948e-07]
## (score 392.2614 & scale 5.213115).
## Hessian positive definite, eigenvalue range [1.440887,300.5744].
## Model rank = 16 / 16
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(Depth) 15.00 4.83 0.78 0.17We see almost no improvement, so might refit the second model here with k=8.
Check to see if any of your models need this treatment.
Checking k size
I saved the following models last time: dsm_nb_xy_ms_p1_4, dsm_tw_xy_ms, dsm_nb_noxy_ms_p1_3, dsm_nb_xy, dsm_tw_xy, dsm_nb_xy_ms_ae and dsm_tw_xy_ms_ae.
I’ll run gam.check on each. Note that I’ve used the fig.show='hide' “chunk” option in RMarkdown to suppress the plots here and save space
gam.check(dsm_nb_xy_ms_p1_4)##
## Method: REML Optimizer: outer newton
## full convergence after 6 iterations.
## Gradient range [-1.246014e-07,3.24631e-05]
## (score 382.759 & scale 1).
## Hessian positive definite, eigenvalue range [0.3417744,30.6163].
## Model rank = 48 / 48
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(x,y) 29.000 1.864 0.55 <2e-16 ***
## s(Depth) 9.000 3.418 0.69 0.045 *
## s(EKE) 9.000 0.856 0.74 0.680
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1gam.check(dsm_tw_xy_ms)##
## Method: REML Optimizer: outer newton
## full convergence after 11 iterations.
## Gradient range [-7.566865e-07,4.88331e-07]
## (score 385.0946 & scale 4.573293).
## Hessian positive definite, eigenvalue range [0.3357644,301.6658].
## Model rank = 48 / 48
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(x,y) 29.000 1.897 0.62 <2e-16 ***
## s(Depth) 9.000 3.695 0.81 0.48
## s(EKE) 9.000 0.811 0.82 0.52
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1gam.check(dsm_nb_noxy_ms_p1_3)##
## Method: REML Optimizer: outer newton
## full convergence after 6 iterations.
## Gradient range [-2.657814e-08,5.557278e-06]
## (score 383.0449 & scale 1).
## Hessian positive definite, eigenvalue range [0.3565081,30.73282].
## Model rank = 28 / 28
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(Depth) 9.000 3.415 0.69 0.025 *
## s(EKE) 9.000 0.866 0.74 0.700
## s(NPP) 9.000 2.814 0.68 0.025 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1gam.check(dsm_nb_xy)##
## Method: REML Optimizer: outer newton
## full convergence after 4 iterations.
## Gradient range [2.05065e-09,4.38541e-08]
## (score 394.963 & scale 1).
## Hessian positive definite, eigenvalue range [2.166966,29.85537].
## Model rank = 25 / 25
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(x,y) 24.0 13.9 0.51 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1gam.check(dsm_tw_xy)##
## Method: REML Optimizer: outer newton
## full convergence after 7 iterations.
## Gradient range [-3.976573e-07,8.201047e-08]
## (score 396.9528 & scale 4.966045).
## Hessian positive definite, eigenvalue range [2.604446,301.2817].
## Model rank = 25 / 25
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(x,y) 24.0 14.5 0.6 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1gam.check(dsm_nb_xy_ms_ae)##
## Method: REML Optimizer: outer newton
## full convergence after 7 iterations.
## Gradient range [1.666612e-10,2.633464e-09]
## (score 505.0342 & scale 1).
## Hessian positive definite, eigenvalue range [0.2941817,44.37103].
## Model rank = 48 / 48
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(x,y) 29.000 1.335 0.46 <2e-16 ***
## s(Depth) 9.000 3.503 0.55 0.025 *
## s(EKE) 9.000 0.772 0.60 0.795
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1gam.check(dsm_tw_xy_ms_ae)##
## Method: REML Optimizer: outer newton
## full convergence after 11 iterations.
## Gradient range [-2.828872e-06,1.348677e-06]
## (score 455.5013 & scale 8.662789).
## Hessian positive definite, eigenvalue range [0.3430708,293.2559].
## Model rank = 48 / 48
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(x,y) 29.000 1.900 0.66 <2e-16 ***
## s(Depth) 9.000 3.772 0.84 0.59
## s(EKE) 9.000 0.826 0.84 0.60
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1(That might have taken a while to run!)
Sifting through these results, most of the cases were the \(p\)-value is significant (at any level beyond 0.1) it’s the case that the EDF is far from \(k^\prime\). The only models I am worried about are dsm_nb_xy and dsm_tw_xy where the \(p\)-value is singificant and \(k^\prime\) is “near” the EDF (it’s not that near but is is over half, so let’s check it out).
For the Tweedie model we set k=25 and got an EDF of 14.5, so we’ll double k to 50:
dsm_tw_xy <- dsm(count~s(x,y, bs="ts", k=50),
df_hn, segs, obs,
family=tw())
summary(dsm_tw_xy)##
## Family: Tweedie(p=1.283)
## Link function: log
##
## Formula:
## count ~ s(x, y, bs = "ts", k = 50) + offset(off.set)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -20.6281 0.2518 -81.92 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(x,y) 20.34 49 2.068 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.125 Deviance explained = 37.5%
## -REML = 406.92 Scale est. = 4.6374 n = 949gam.check(dsm_tw_xy)##
## Method: REML Optimizer: outer newton
## full convergence after 6 iterations.
## Gradient range [-0.0001153251,2.507349e-05]
## (score 406.9169 & scale 4.637444).
## Hessian positive definite, eigenvalue range [1.816522,297.6612].
## Model rank = 50 / 50
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(x,y) 49.0 20.3 0.62 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Okay that’s a moderate jump there! We can see now the EDF is less than half of \(k^\prime\) but as we can see the \(p\)-value is still significant. It’s not likely that increasing \(k\) will do anything here. As I said in the lecture, we need to look at \(k^\prime\), the EDF and the \(p\)-value here to see what’s going on.
We can do the same thing with the negative binomial model, there the EDF was 13.9 when k=25.
dsm_nb_xy <- dsm(count~s(x,y, bs="ts", k=50),
df_hn, segs, obs,
family=nb())
summary(dsm_nb_xy)##
## Family: Negative Binomial(0.108)
## Link function: log
##
## Formula:
## count ~ s(x, y, bs = "ts", k = 50) + offset(off.set)
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -20.6787 0.2423 -85.35 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(x,y) 19.75 49 81.49 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.102 Deviance explained = 41.6%
## -REML = 405.47 Scale est. = 1 n = 949gam.check(dsm_nb_xy)##
## Method: REML Optimizer: outer newton
## full convergence after 4 iterations.
## Gradient range [-0.0001709692,0.0001083025]
## (score 405.4683 & scale 1).
## Hessian positive definite, eigenvalue range [1.705091,27.92868].
## Model rank = 50 / 50
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(x,y) 49.0 19.8 0.54 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Again we see a wee jump in the EDF but the \(p\)-value is still significant. For the same reasons as above, we’re not going to worry too much about that.
Plotting residuals vs. covariates
Now following from the slides, let’s try plotting residuals vs. covariates. We can use these diagnostics (along with those in the following section) to assess how well the model has accounted for the structure in the data. This involves a little data manipulation and a little thought about how to chop up the covariate values.
I’ll show an example here using the not very good model we fitted above, so we can see what happens when things go wrong:
dsm_k_low <- dsm(count ~ s(Depth, k=4),
df_hn, segs, obs,
family=tw())
# get the data used to fit the model
resid_dat <- dsm_k_low$data
# get the residuals
resid_dat$resids <- residuals(dsm_k_low)
# now let's chop x into 20 chunks
resid_dat$x_chop <- cut(resid_dat$x,
seq(min(resid_dat$x), max(resid_dat$x), len=20))
# and plot the boxplot
boxplot(resids~x_chop, data=resid_dat, ylab="Deviance residuals", xlab="x")
This seems to show that something is going on in the far east of the data – there’s a wide range of residuals there. Note that for our data the most easterly points are also North (since our data is somewhat diagonal in shape, with respect to the compass). So let’s see if this is an issue with x alone (East-West) or y (North-Southness) is an issue too.
# now let's chop y into 20 chunks
resid_dat$y_chop <- cut(resid_dat$y,
seq(min(resid_dat$y), max(resid_dat$y), len=20))
# and plot the boxplot
boxplot(resids~y_chop, data=resid_dat, ylab="Deviance residuals", xlab="y")
We can also do this chopping with Depth (the only covariate in this model):
# now let's chop Depth into 20 chunks
resid_dat$Depth_chop <- cut(resid_dat$Depth,
seq(min(resid_dat$Depth), max(resid_dat$Depth), len=20))
boxplot(resids~Depth_chop, data=resid_dat, ylab="Deviance residuals", xlab="Depth")
Now we see something more pathological, there is pattern in these residuals.
Let’s make the same plot for a model with an spatial smooth in it:
dsm_better <- dsm(count ~ s(x,y) + s(Depth),
df_hn, segs, obs,
family=tw())
# get the data from the larger k model
resid_dat <- dsm_better$data
# get the residuals
resid_dat$resids <- residuals(dsm_better)
# now let's chop Depth into 20 chunks
resid_dat$Depth_chop <- cut(resid_dat$Depth,
seq(min(resid_dat$Depth), max(resid_dat$Depth), len=20))
boxplot(resids~Depth_chop, data=resid_dat, ylab="Deviance residuals", xlab="Depth")
This looks better (note the change in vertical axis scale).
Try this out with your models, using the covariates in the model and perhaps x and y too. You may want to change the number of chunks (via the len argument to cut()).
Looking for patterns in the residuals
Trying this out with the models without smooths of space in them (the habitat models):
# get the data from the larger k model
resid_dat <- dsm_nb_noxy_ms_p1_3$data
# get the residuals
resid_dat$resids <- residuals(dsm_nb_noxy_ms_p1_3)
# now let's chop x into 20 chunks
resid_dat$x_chop <- cut(resid_dat$x,
seq(min(resid_dat$x), max(resid_dat$x), len=20))
# now let's chop y into 20 chunks
resid_dat$y_chop <- cut(resid_dat$y,
seq(min(resid_dat$y), max(resid_dat$y), len=20))
par(mfrow=c(1,2))
boxplot(resids~x_chop, data=resid_dat,
ylab="Deviance residuals", xlab="x", las=2)
boxplot(resids~y_chop, data=resid_dat,
ylab="Deviance residuals", xlab="y", las=2)
Comparing this to the model with space in it:
# get the data from the larger k model
resid_dat <- dsm_nb_xy_ms_p1_4$data
# get the residuals
resid_dat$resids <- residuals(dsm_nb_xy_ms_p1_4)
# now let's chop x into 20 chunks
resid_dat$x_chop <- cut(resid_dat$x,
seq(min(resid_dat$x), max(resid_dat$x), len=20))
# now let's chop y into 20 chunks
resid_dat$y_chop <- cut(resid_dat$y,
seq(min(resid_dat$y), max(resid_dat$y), len=20))
par(mfrow=c(1,2))
boxplot(resids~x_chop, data=resid_dat,
ylab="Deviance residuals", xlab="x", las=2)
boxplot(resids~y_chop, data=resid_dat,
ylab="Deviance residuals", xlab="y", las=2)
Neither are perfect, there is some pattern in them but we don’t see a really tight pattern as we do in the pathological versions above.
Fitting models to the residuals
We can also fit a model to the residuals to see if there is residual structure left in the model. The idea here is that we know that the model won’t fit perfectly everywhere, but is the difference between model and data (the residuals) systematic. Going back to our “bad” model above, we can see what happens when things go wrong:
# get the data used to fit the model
resid_dat <- dsm_k_low$data
# get the residuals
resid_dat$resids <- residuals(dsm_k_low)
# fit the model with the same predictors, but larger k
resid_model_depth <- gam(resids ~ s(Depth, k=40),
data=resid_dat,
family=gaussian(), method="REML")Note that we used the gam function from mgcv (which is what happens inside dsm anyway), to fit our residual model. The important things to note are:
gam()in place ofdsm().- The response is
resids, the column we added above. - Theoretically we expect the deviance residuals to be distributed normally, so we set
family=gaussian(). - Need to set
method="REML"which happens automatically indsm(), and makes sure the best method is used for fitting the smooths.
We want to look at the EDFs associated with the terms in the model:
summary(resid_model_depth)##
## Family: gaussian
## Link function: identity
##
## Formula:
## resids ~ s(Depth, k = 40)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.54782 0.03492 -15.69 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(Depth) 3.304 4.14 8.315 1.34e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0344 Deviance explained = 3.77%
## -REML = 1422 Scale est. = 1.1571 n = 949We can see that there’s still a bit of pattern there, but not much (you can also try using plot() to investigate this).
Now trying that with the better model:
# get the data used to fit the model
resid_dat_better <- dsm_better$data
# get the residuals
resid_dat_better$resids <- residuals(dsm_better)
# fit the model with the same predictors, but larger k
resid_model_better <- gam(resids ~ s(Depth, bs="ts", k=40) + s(x,y, bs="ts", k=100),
data=resid_dat_better,
family=gaussian(), method="REML")Now again looking at the EDFs:
summary(resid_model_better)##
## Family: gaussian
## Link function: identity
##
## Formula:
## resids ~ s(Depth, bs = "ts", k = 40) + s(x, y, bs = "ts", k = 100)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.45018 0.03071 -14.66 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(Depth) 2.435 39 0.550 6.79e-06 ***
## s(x,y) 1.028 99 0.045 0.0255 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0296 Deviance explained = 3.32%
## -REML = 1302.3 Scale est. = 0.8952 n = 949We see that there’s not much going on in this model. Again you can use plot() here, remembering to compare the vertical axis limits and look at the width of the confidence bands.
Try this out with your model(s). What do you think is missing?
Fitting to the residuals
Trying this out again with dsm_nb_noxy_ms_p1_3 (habitat model) and dsm_nb_xy_ms_p1_4 (model with space included).
I’ll start with the habitat model, where we didn’t include s(x,y) from the start. Let’s fit a model to the residuals that includes all of the possible covariates, to see if we’re missing anything…
# make a template data frame from the fitted data
resid_dat_hab <- dsm_nb_noxy_ms_p1_3$data
# get the residuals
resid_dat_hab$resids <- residuals(dsm_nb_noxy_ms_p1_3)
# fit the model with all predictors, but larger k
resid_model_hab <- gam(resids ~ s(x,y, bs="ts", k=60) +
s(Depth, bs="ts", k=20) +
s(DistToCAS, bs="ts", k=20) +
s(SST, bs="ts", k=20) +
s(EKE, bs="ts", k=20) +
s(NPP, bs="ts", k=20),
data=resid_dat_hab,
family=gaussian(), method="REML")
summary(resid_model_hab)##
## Family: gaussian
## Link function: identity
##
## Formula:
## resids ~ s(x, y, bs = "ts", k = 60) + s(Depth, bs = "ts", k = 20) +
## s(DistToCAS, bs = "ts", k = 20) + s(SST, bs = "ts", k = 20) +
## s(EKE, bs = "ts", k = 20) + s(NPP, bs = "ts", k = 20)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.24747 0.01272 -19.46 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(x,y) 9.335e-01 59 0.089 0.0114 *
## s(Depth) 2.988e+00 19 2.735 <2e-16 ***
## s(DistToCAS) 2.041e-05 19 0.000 0.8395
## s(SST) 2.313e-04 19 0.000 0.6214
## s(EKE) 7.794e-01 19 0.183 0.0332 *
## s(NPP) 4.388e-05 19 0.000 0.7063
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0696 Deviance explained = 7.42%
## -REML = 468.19 Scale est. = 0.15353 n = 949Most of these EDFs are very very small, that’s good! It means we’ve explained most of what’s going on in the data using this model. We can see this a little better by plotting the model terms:
plot(resid_model_hab, scheme=2, pages=1)
We can see that the EDF of around 1 for s(x,y) here is giving a slight gradient in the spatial smooth of the residuals, but no complex pattern (it’s a plane through the space). Depth might cause use a bit more concern.
Let’s try refitting the model, increasing the k for s(Depth) and s(EKE):
dsm_nb_noxy_ms_refit <- dsm(count~s(Depth, bs="ts", k=20) +
#s(DistToCAS, bs="ts") + # 1
#s(SST, bs="ts") + # 2
s(EKE, bs="ts", k=20) +
s(NPP, bs="ts"),
df_hn, segs, obs,
family=nb())
summary(dsm_nb_noxy_ms_refit)##
## Family: Negative Binomial(0.115)
## Link function: log
##
## Formula:
## count ~ s(Depth, bs = "ts", k = 20) + s(EKE, bs = "ts", k = 20) +
## s(NPP, bs = "ts") + offset(off.set)
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -20.7750 0.2281 -91.09 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(Depth) 3.4649 19 37.088 < 2e-16 ***
## s(EKE) 0.8657 19 5.532 0.0106 *
## s(NPP) 2.8152 9 26.063 1.73e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0966 Deviance explained = 40.7%
## -REML = 383.13 Scale est. = 1 n = 949We don’t see any change here, so it’s not that k is the issue. More likely there’s just not enough information in the data to account for the variation that’s there. For example NPP is a rough guide to net primary productivity (it’s a measure of epipelagic micronekton, much smaller than the squid we usually see sperm whales eat), so this may not be “close enough” to the prey that sperm whales are interested in.
Trying the same thing for the model that has a smooth of space in it…
# make a template data frame from the fitted data
resid_dat_xy <- dsm_nb_xy_ms_p1_4$data
# get the residuals
resid_dat_xy$resids <- residuals(dsm_nb_xy_ms_p1_4)
# fit the model with the same predictors, but larger k
resid_model_xy <- gam(resids ~ s(x,y, bs="ts", k=60) +
s(Depth, bs="ts", k=20) +
s(DistToCAS, bs="ts", k=20) +
s(SST, bs="ts", k=20) +
s(EKE, bs="ts", k=20) +
s(NPP, bs="ts", k=20),
data=resid_dat_xy,
family=gaussian(), method="REML")
summary(resid_model_xy)##
## Family: gaussian
## Link function: identity
##
## Formula:
## resids ~ s(x, y, bs = "ts", k = 60) + s(Depth, bs = "ts", k = 20) +
## s(DistToCAS, bs = "ts", k = 20) + s(SST, bs = "ts", k = 20) +
## s(EKE, bs = "ts", k = 20) + s(NPP, bs = "ts", k = 20)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.24892 0.01261 -19.74 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(x,y) 1.266e+01 59 0.569 9.83e-05 ***
## s(Depth) 2.510e+00 19 2.194 < 2e-16 ***
## s(DistToCAS) 5.518e-01 19 0.064 0.0478 *
## s(SST) 1.349e-04 19 0.000 0.9364
## s(EKE) 6.698e-01 19 0.106 0.0652 .
## s(NPP) 3.077e-04 19 0.000 0.9658
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.103 Deviance explained = 11.9%
## -REML = 480.85 Scale est. = 0.15087 n = 949Again we see that space and EKE have small effects (close to EDF of 1) and Depth has a slightly larger one. Looking at a plot of these effects:
plot(resid_model_xy, scheme=2, pages=1)
We see very very similar results to the above. Let’s try refitting that model, but increasing k for s(x,y), s(Depth) and s(EKE):
dsm_nb_xy_ms_refit <- dsm(count~ #s(x,y, bs="ts", k=60) + # 1
s(Depth, bs="ts", k=20) +
#s(DistToCAS, bs="ts") + # 2
#s(SST, bs="ts") + # 3
s(EKE, bs="ts", k=20) +
s(NPP, bs="ts"),
df_hn, segs, obs,
family=nb())
summary(dsm_nb_xy_ms_refit)##
## Family: Negative Binomial(0.115)
## Link function: log
##
## Formula:
## count ~ s(Depth, bs = "ts", k = 20) + s(EKE, bs = "ts", k = 20) +
## s(NPP, bs = "ts") + offset(off.set)
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -20.7750 0.2281 -91.09 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(Depth) 3.4649 19 37.088 < 2e-16 ***
## s(EKE) 0.8657 19 5.532 0.0106 *
## s(NPP) 2.8152 9 26.063 1.73e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0966 Deviance explained = 40.7%
## -REML = 383.13 Scale est. = 1 n = 949So! We’ve now ended-up with the same model as dsm_nb_noxy_ms_refit!
This again shows that we need to look at multiple outputs, adjust and see what’s going on.
Let’s do the residual plot again for the terms in the model:
# get the data from the larger k model
resid_dat <- dsm_nb_xy_ms_refit$data
# get the residuals
resid_dat$resids <- residuals(dsm_nb_xy_ms_refit)
# now let's chop Depth into 20 chunks
resid_dat$Depth_chop <- cut(resid_dat$Depth,
seq(min(resid_dat$Depth), max(resid_dat$Depth), len=20))
# now let's chop EKE into 20 chunks
resid_dat$EKE_chop <- cut(resid_dat$EKE,
seq(min(resid_dat$EKE), max(resid_dat$EKE), len=20))
# now let's chop NPP into 20 chunks
resid_dat$NPP_chop <- cut(resid_dat$NPP,
seq(min(resid_dat$NPP), max(resid_dat$NPP), len=20))
par(mfrow=c(1,3))
boxplot(resids~Depth_chop, data=resid_dat,
ylab="Deviance residuals", xlab="Depth", las=2)
boxplot(resids~EKE_chop, data=resid_dat,
ylab="Deviance residuals", xlab="EKE", las=2)
boxplot(resids~NPP_chop, data=resid_dat,
ylab="Deviance residuals", xlab="NPP", las=2)
We can see there’s still something going on with Depth, but as I said above, this may be down to still not having a good covariate for prey.
Finally we’ll do the same thing with the Tweedie model with space in it dsm_tw_xy_ms. First a reminder of what’s in that model:
summary(dsm_tw_xy_ms)##
## Family: Tweedie(p=1.279)
## Link function: log
##
## Formula:
## count ~ s(x, y, bs = "ts") + s(Depth, bs = "ts") + s(EKE, bs = "ts") +
## offset(off.set)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -20.673 0.234 -88.33 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(x,y) 1.8969 29 0.707 1.84e-05 ***
## s(Depth) 3.6949 9 5.024 < 2e-16 ***
## s(EKE) 0.8106 9 0.470 0.0216 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.105 Deviance explained = 34.9%
## -REML = 385.09 Scale est. = 4.5733 n = 949Now doing the same thing again and refitting to the residuals:
# make a template data frame from the fitted data
resid_dat_xy <- dsm_tw_xy_ms$data
# get the residuals
resid_dat_xy$resids <- residuals(dsm_tw_xy_ms)
# fit the model with all predictors, but larger k
resid_model_xy <- gam(resids ~ s(x,y, bs="ts", k=60) +
s(Depth, bs="ts", k=20) +
s(DistToCAS, bs="ts", k=20) +
s(SST, bs="ts", k=20) +
s(EKE, bs="ts", k=20) +
s(NPP, bs="ts", k=20),
data=resid_dat_xy,
family=gaussian(), method="REML")
summary(resid_model_xy)##
## Family: gaussian
## Link function: identity
##
## Formula:
## resids ~ s(x, y, bs = "ts", k = 60) + s(Depth, bs = "ts", k = 20) +
## s(DistToCAS, bs = "ts", k = 20) + s(SST, bs = "ts", k = 20) +
## s(EKE, bs = "ts", k = 20) + s(NPP, bs = "ts", k = 20)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.49222 0.03223 -15.27 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(x,y) 1.076e+01 59 0.440 0.000874 ***
## s(Depth) 1.789e+00 19 0.904 5.42e-06 ***
## s(DistToCAS) 6.035e-01 19 0.078 0.042294 *
## s(SST) 2.494e-04 19 0.000 0.762203
## s(EKE) 9.816e-03 19 0.000 0.304668
## s(NPP) 2.082e-04 19 0.000 0.979887
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0494 Deviance explained = 6.26%
## -REML = 1366 Scale est. = 0.98562 n = 949Again, most of these EDFs are very very small. Let’s again plot the model terms:
plot(resid_model_xy, scheme=2, pages=1)
Looking at the summary results and the plot, I think we need to increase k for s(x,y), s(Depth) and s(DistToCAS) (which wasn’t in our model!):
dsm_tw_xy_ms_refit <- dsm(count~s(x,y, bs="ts", k=60)+
s(Depth, bs="ts", k=20) +
#s(DistToCAS, bs="ts", k=20) + # 1
#s(SST, bs="ts") + # 2
s(EKE, bs="ts"),# +
#s(NPP, bs="ts"),
df_hn, segs, obs,
family=tw())
summary(dsm_tw_xy_ms_refit)##
## Family: Tweedie(p=1.279)
## Link function: log
##
## Formula:
## count ~ s(x, y, bs = "ts", k = 60) + s(Depth, bs = "ts", k = 20) +
## s(EKE, bs = "ts") + offset(off.set)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -20.673 0.234 -88.35 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(x,y) 1.8967 59 0.348 1.82e-05 ***
## s(Depth) 3.7784 19 2.385 < 2e-16 ***
## s(EKE) 0.8105 9 0.470 0.0216 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.106 Deviance explained = 34.9%
## -REML = 385.18 Scale est. = 4.572 n = 949So we’re back to where we started! In this case increasing k for DistToCAS didn’t have any impact on our model selection in the end.
Using rqgam.check
We can use rqgam.check() to look at the residual check plots for this model. Here we’re looking for pattern in the plots that would indicate that there is unmodelled structure. Ideally the plot should look like a messy blob of points.
rqgam.check(dsm_better)## Loading required namespace: tweedie
Try this out with your models and see if there is any pattern in the randomised quantile residuals.
Checking randomized quantile residuals
Let’s compare between negative binomial and Tweedie models.
First negative binomial, as refitted above:
rqgam.check(dsm_nb_noxy_ms_refit)
Then Tweedie:
rqgam.check(dsm_tw_xy_ms)
Both of these seem fine, not pathological like those in the lecture slides.
Using obs_exp
To compare the observed vs. expected counts from the model, we need to aggregate the data at some level. We can use obs_exp() to do this. For continuous covariates in the detection function or spatial moel, we just need to specify the cutpoints. For example:
obs_exp(dsm_better, "Depth", c(0, 1000, 2000, 3000, 4000, 6000))## (0,1e+03] (1e+03,2e+03] (2e+03,3e+03] (3e+03,4e+03] (4e+03,6e+03]
## Observed 4.000000 52.53333 139.1667 35.00000 8.000000
## Expected 5.163824 49.92310 126.8073 38.65673 8.134937Try this out for your models.
Looking at observed vs. expected counts
Again trying this for our negative binomial habitat model using the Depth covariate:
obs_exp(dsm_nb_noxy_ms_refit, "Depth", c(0, 1000, 2000, 3000, 4000, 6000))## (0,1e+03] (1e+03,2e+03] (2e+03,3e+03] (3e+03,4e+03] (4e+03,6e+03]
## Observed 4.000000 52.53333 139.1667 35.00000 8.000000
## Expected 4.973135 39.26691 122.2847 51.45393 7.425874This looks pretty good! We can try with EKE and NPP too:
# for EKE the majority of the values are between 0 and
# 0.05, so we'll put more bins there and make bins larger for larger values -- see hist(segs$EKE)
obs_exp(dsm_nb_noxy_ms_refit, "EKE", c(0, 0.01, 0.025, 0.05, 0.1, 0.4, 0.7))## (0,0.01] (0.01,0.025] (0.025,0.05] (0.05,0.1] (0.1,0.4] (0.4,0.7]
## Observed 161.9500 37.7500 17.00000 15.00000 7.000000 0.00000000
## Expected 134.0855 42.0134 22.50145 18.88298 7.920128 0.00107712# we can get the default breaks we have from hist()
npp_breaks <- hist(segs$NPP, plot=FALSE)$breaks
obs_exp(dsm_nb_noxy_ms_refit, "NPP", npp_breaks)## (0,200] (200,400] (400,600] (600,800] (800,1e+03] (1e+03,1.2e+03]
## Observed 11.450000 7.50000 98.38333 81.36667 28.00000 11.000000
## Expected 9.531847 13.22469 82.90801 88.24576 24.38344 6.143246
## (1.2e+03,1.4e+03] (1.4e+03,1.6e+03] (1.6e+03,1.8e+03] (1.8e+03,2e+03]
## Observed 1.0000000 0.0000000 0.000000000 0.000000000
## Expected 0.6519371 0.3115906 0.002891405 0.001170438Again, these seem pretty reasonable!
Doing the same thing for the Tweedie model dsm_tw_xy_ms, again using the default histogram breaks for x and y:
x_breaks <- hist(segs$x, plot=FALSE)$breaks
obs_exp(dsm_tw_xy_ms, "x", x_breaks)## (-8e+05,-7e+05] (-7e+05,-6e+05] (-6e+05,-5e+05] (-5e+05,-4e+05]
## Observed 0.00000000 0.0000000 1.000000 0.0000000
## Expected 0.00316796 0.3395217 0.561756 0.8881759
## (-4e+05,-3e+05] (-3e+05,-2e+05] (-2e+05,-1e+05] (-1e+05,0] (0,1e+05]
## Observed 0.00000 7.450000 43.33333 50.91667 54.00000
## Expected 3.27297 8.245829 25.72470 43.54315 49.57802
## (1e+05,2e+05] (2e+05,3e+05] (3e+05,4e+05] (4e+05,5e+05] (5e+05,6e+05]
## Observed 34.00000 21.00000 13.0000 1.000000 4.000000
## Expected 48.27922 22.73238 11.3278 9.773835 7.411754
## (6e+05,7e+05]
## Observed 9.000000
## Expected 1.645921y_breaks <- hist(segs$y, plot=FALSE)$breaks
obs_exp(dsm_tw_xy_ms, "y", y_breaks)## (-7e+05,-6e+05] (-6e+05,-5e+05] (-5e+05,-4e+05] (-4e+05,-3e+05]
## Observed 0.00000000 0.00000000 1.0000000 0.0000
## Expected 0.06294614 0.05824977 0.2168808 2.8786
## (-3e+05,-2e+05] (-2e+05,-1e+05] (-1e+05,0] (0,1e+05] (1e+05,2e+05]
## Observed 7.450000 0.000000 0.000000 4.000000 2.00000
## Expected 2.376874 3.428002 3.866229 8.603532 11.44869
## (2e+05,3e+05] (3e+05,4e+05] (4e+05,5e+05] (5e+05,6e+05] (6e+05,7e+05]
## Observed 41.33333 52.91667 29.00000 57.00000 37.00000
## Expected 20.73613 19.88760 41.55495 58.52634 51.50198
## (7e+05,8e+05]
## Observed 7.000000
## Expected 8.181202obs_exp(dsm_tw_xy_ms, "Depth", c(0, 1000, 2000, 3000, 4000, 6000))## (0,1e+03] (1e+03,2e+03] (2e+03,3e+03] (3e+03,4e+03] (4e+03,6e+03]
## Observed 4.00000 52.53333 139.1667 35.00000 8.000000
## Expected 5.80883 47.76153 131.9859 39.23381 8.538091# again need to be careful with EKE
obs_exp(dsm_tw_xy_ms, "EKE", c(0, 0.01, 0.025, 0.05, 0.1, 0.4, 0.7))## (0,0.01] (0.01,0.025] (0.025,0.05] (0.05,0.1] (0.1,0.4] (0.4,0.7]
## Observed 161.9500 37.75000 17.00000 15.00000 7.000000 0.000000000
## Expected 151.2765 34.55754 19.68545 18.86281 8.944266 0.001596079The EKE and Depth results look good, but the spatial terms look less good. That’s just because we’re taking marginals of a 2 dimensional distribution. We might do better looking at “areas of interest” or other covariates here. For example we can look at DistToCAS even though it’s not in the model:
# again need to inspect hist(segs$DistToCAS) first
obs_exp(dsm_tw_xy_ms, "DistToCAS", c(-1e-10, 1000, 2000, 5000, 10000, 350000))## (-1e-10,1e+03] (1e+03,2e+03] (2e+03,5e+03] (5e+03,1e+04]
## Observed 44.05000 43.83333 28.50000 24.2000
## Expected 38.94865 31.87421 29.49987 31.5549
## (1e+04,3.5e+05]
## Observed 98.11667
## Expected 101.45057This looks pretty good!
Saving models
Now save the models that you’d like to use to predict with later (that have good check results!): I recommend saving a few models so you can compare the results later.
# add your models here
save(dsm_nb_noxy_ms_refit, dsm_tw_xy_ms,
file="dsms-checked.RData")