# How to construct a GAM for 1/0 response of whale behaviours with environmental variables?

Firstly, I am novice at GAM construction, as I am sure will become clear very quickly! I have 2 questions related to what I have come across. I want to use a GAM in MGCV to figure which temporal, spatial and environmental variables have a relationship, and what that relationship looks like, on the behaviour of humpback whales.

The data were collected from a fixed location, where all variables were recorded and a 'random' pod of humpbacks were watched for 2 minutes. Each behaviour we were interested in was marked with a '1' or '0' if it was shown during those 2 minutes. These data span from DoY 152 - 314 (~June - Nov seasonally) each year from 2011-2020.

    str(hws)
$$DATE: POSIXct, format: "2011-07-14" "2011-07-18" ...$$ DOY : num  195 199 202 202 202 202 206 206 206 206 ...
$$YR : Factor w/ 10 levels "2011","2012",..: 1 1 1 1 1 1 1 1 1 1 ...$$ W (WEATHER)   : Factor w/ 5 levels "C","O","OC","PC",..: 4 1 1 1 1 1 1 1 1 1 ...
$$SS (SEASTATE) : Factor w/ 5 levels "1","2","3","4",..: 4 2 2 2 2 2 2 2 2 2 ...$$ TIME: num  0.688 0.689 0.662 0.667 0.688 ...
$$ED (EST_DIST) : num 5 4 9 6 8 7 7 7 7 7 ...$$ TD (TRAV_DIR)  : Factor w/ 4 levels "E","N","S","W": NA 3 2 2 2 1 2 2 2 2 ...
$$POD (POD_SIZE) : num 1 1 3 3 1 1 3 3 2 1 ...$$ BL  : Factor w/ 2 levels "0","1": 1 2 1 1 1 1 2 2 2 1 ...
$$BA : Factor w/ 2 levels "0","1": 1 2 1 1 2 1 1 1 1 1 ...$$ BR  : Factor w/ 2 levels "0","1": 2 1 2 2 1 2 2 2 2 2 ...
$$FS : Factor w/ 2 levels "0","1": 1 1 2 1 1 1 2 2 2 2 ...$$ TS  : Factor w/ 2 levels "0","1": 1 1 2 1 1 1 1 1 1 1 ...
$$PS : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...$$ FUD : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
\$ CALF: num  0 0 0 0 0 0 0 0 0 0 ...


I can't find a way to return TIME (time of day) as hh:mm:ss. But anyway, I started by using BR (breach) as the response and, from my limited mgcv argument knowledge, try to fit a simple model to kick it off...

br <- gam(BR ~ s(DOY, k = 100) + YR + s(W, bs = "re") + SS + s(ED, k = 7) +
s(POD, k = 3) + TD, data = hws, method = "REML",
family = binomial, select = TRUE)
summary(br)

Family: binomial

Formula:
BR ~ s(DOY, k = 100) + YR + s(W, bs = "re") + SS + s(ED,
k = 7) + s(POD, k = 3) + TD

Parametric coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.26349    0.23082  -5.474 4.40e-08 ***
YR2012      -0.86836    0.15212  -5.708 1.14e-08 ***
YR2013      -0.85179    0.16504  -5.161 2.45e-07 ***
YR2014      -0.65376    0.14903  -4.387 1.15e-05 ***
YR2015      -0.99120    0.14619  -6.780 1.20e-11 ***
YR2016      -1.04035    0.24481  -4.250 2.14e-05 ***
YR2017      -0.64429    0.16962  -3.798 0.000146 ***
YR2018      -0.52502    0.14828  -3.541 0.000399 ***
YR2019      -1.84776    0.18141 -10.186  < 2e-16 ***
YR2020      -0.25963    0.30593  -0.849 0.396064
SS2          0.35263    0.07311   4.823 1.41e-06 ***
SS3          0.81142    0.09892   8.202 2.36e-16 ***
SS4         -0.31774    0.29623  -1.073 0.283437
SS5          1.66526    0.37249   4.471 7.80e-06 ***
TDN          0.25232    0.18230   1.384 0.166349
TDS          0.19895    0.18415   1.080 0.279973
TDW         -0.41442    0.42258  -0.981 0.326747
---
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(DOY) 35.8742     99 248.997  < 2e-16 ***
s(W)    1.6572      4   4.653   0.0515 .
s(ED)   0.9869      6  71.903  < 2e-16 ***
s(POD)  0.9833      2  58.170 1.56e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.0742   Deviance explained = 8.09%
-REML = 3739.7  Scale est. = 1         n = 8229

gam.check(br)

Method: REML   Optimizer: outer newton
full convergence after 14 iterations.
(score 3739.729 & scale 1).
Hessian positive definite, eigenvalue range [2.726934e-07,2.117211].
Model rank =  129 / 129

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(DOY) 99.000 35.874    0.91  <2e-16 ***
s(W)    5.000  1.657      NA      NA
s(ED)   6.000  0.987    0.95   0.335
s(POD)  2.000  0.983    0.94   0.085 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1



So from what I can interpret, I ought to probably drop travel direction (TD) and increase k in s(DOY, k=100), but sea state 4 (SS4) and year 2020 also don't have significant P-values.

Q1. What can I do about that, seeing as they're levels within factors? Must I just live with them, and interpret these levels as not being able to be explained? Or should I be telling the model that these levels are ordered?

Q2. The Q-Q and residual plots are mental, is this because I am using 1/0 data? Or I have used totally the wrong family? Please help!!