I'm new to GAMs and am trying to figure out how to parameterize and interpret models where the response variable is percent of salmon migrating through a given route and I think there might be autocorrelation. Apologies in advance for the long post but you can skip to my question at the end if it's helpful.
The data contains percent (my response variable), Julian date (JulianDate), Year, and stock (GroupName). There is also some biological reason to believe that at least some stocks differ in their percents depending upon which "cycle" the year of migration is in (CycleLine, e.g. percents may be lower for some stocks for Cycle3 years).
BACKGROUND MODEL PARAMETERIZATION
I think the percent should be explained by JulianDate, GroupName, and maybe CycleLine. I started with:
gam1 <- bam(Percent ~ GroupName +
s(JulianDate) + s(JulianDate, by = GroupName),
data = data,
method = "fREML",
family = betar(),
discrete = T,
select = T,
# select = T will place extra penalties to help with variable selection
nthreads = 2)
based off of this question where @GavinSimpson suggested estimating a global or average smooth effect of x on y (the s(JulianDate) term) plus a smooth difference term (the second s(JulianDate, by = GroupName) term).
summary(gam1)
Produces:
Family: Beta regression(3.059)
Link function: logit
Formula:
Percent ~ GroupName + s(JulianDate) + s(JulianDate, by = GroupName)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.16096 0.04339 -3.709 0.000210 ***
GroupNameStock2 0.18920 0.04975 3.803 0.000145 ***
GroupNameStock4 -0.51227 0.05981 -8.565 < 2e-16 ***
GroupNameStock3 0.49874 0.05098 9.783 < 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(JulianDate) 6.8299873 9 59.983 < 2e-16 ***
s(JulianDate):GroupNameStock1 0.0001348 9 0.000 0.394002
s(JulianDate):GroupNameStock2 0.2757218 9 0.036 0.278520
s(JulianDate):GroupNameStock4 5.7796322 9 2.487 0.000215 ***
s(JulianDate):GroupNameStock3 3.5601339 9 5.724 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.397 Deviance explained = 44.1%
fREML = 2947.5 Scale est. = 1 n = 4385
gratia::draw(gam1)
produces:
But the R-sq is only 0.397 and when I check forecast::auto.arima(residuals(gam1)) I get
Series: residuals(gam1)
ARIMA(2,1,2)
Coefficients:
ar1 ar2 ma1 ma2
1.4314 -0.5199 -1.482 0.4889
s.e. 0.0723 0.0618 0.076 0.0754
sigma^2 = 0.1825: log likelihood = -2490.74
AIC=4991.48 AICc=4991.5 BIC=5023.41
As I mentioned above, I have reason to believe that CycleLine could improve the fit and maybe deal with the autocorrelation:
gam2 <- bam(Percent ~ GroupName + CycleLine +
s(JulianDate) +
s(JulianDate, by = GroupName) +
s(JulianDate, by = CycleLine),
data = data,
method = "fREML",
family = betar(),
discrete = T,
select = T,
nthreads = 2)
The summary(gam2)$r.sq is still only 0.426 though and forecast::auto.arima(residuals(gam2)) produces:
Series: residuals(gam2)
ARIMA(2,1,2)
Coefficients:
ar1 ar2 ma1 ma2
1.3935 -0.4912 -1.4397 0.4464
s.e. 0.0709 0.0605 0.0744 0.0739
sigma^2 = 0.1908: log likelihood = -2587.87
AIC=5185.73 AICc=5185.75 BIC=5217.66
Even when I scrap CycleLine and just look at the interaction between JulianDate and Year the model actually has a slightly higher r-squared of 0.523 but still has un-accounted for autocorrelation:
gam3 <- bam(Percent ~ GroupName +
te(JulianDate, Year, by = GroupName),
data = data,
method = "fREML",
family = betar(),
discrete = T,
select = T,
nthreads = 2)
auto.arima(residuals(gam3)) produces:
Series: residuals(gam3)
ARIMA(2,0,5) with zero mean
Overall, all the auto.arima() checks and acf() checks show there is unaccounted for autocorrelation:
par(mfrow = c(1, 4), cex = 1.1)
acf(residuals(gam1))
acf(residuals(gam2))
acf(residuals(gam3))
AUTOCORRELATION QUESTION
Given that background, I'm struggling with how to correctly account for the autocorrelation. It seems like adding explanatory variables doesn't capture it so I was advised here to look into setting the rho and AR.start variables or increasing my k value but I'm not sure that I'm parameterizing those correctly. I'm not sure what to set my k value to and when I try to add rho and AR.start I get strange results. For example, using a variation of gam2:
# First mark the start of each time series as TRUE,
# and all other data points as FALSE
data %<>%
arrange(GroupName, Year, CycleLine, JulianDate) %>%
as.data.frame()
simdat <- itsadug::start_event(data, column = "JulianDate",
event = c("GroupName", "CycleLine"),
label.event = "Event")
r1 <- itsadug::start_value_rho(gam2, plot = TRUE) # value = 0.879586
gam4 <- bam(Percent ~ GroupName + CycleLine +
s(JulianDate) +
s(JulianDate, by = GroupName) +
s(JulianDate, by = CycleLine),
data = simdat,
rho = r1,
AR.start = simdat$start.event,
method = "fREML",
family = betar(),
discrete = T,
select = T,
nthreads = 2)
```
but I think something must be wrong since the `summary(gam4)$r.sq` is -0.099 and when I look at the uncorrected vs. 'corrected' residuals I'm still seeing autocorrelation:
```r
# Uncorrected versus corrected residuals:
par(mfrow = c(1, 2), cex = 1.1)
itsadug::acf_resid(gam2)
itsadug::acf_resid(gam4)
How can I better account for the autocorrelation? My understanding that rho is just for an AR1 but the results of auto.arima suggest that the process is more complicated? And the r1 value doesn't match the output of forecast::auto.arima(residuals(gam2))$coef (see values above).
So, to summarize, any feedback would be greatly appreciated on how I can better account for autocorrelation. Even when I have JulianDate and CycleLine/Year I see autocorrelation in the residuals AND when I set rho and AR.start I still see autocorrelation. Am I miss-specifying rho and AR.start or do I need a more complex way of accounting for the autocorrelation (and what might that be)?
