gavinsimpson / gratia

ggplot-based graphics and useful functions for GAMs fitted using the mgcv package
https://gavinsimpson.github.io/gratia/
Other
202 stars 28 forks source link

Error with AR1 model : 'x' must be a factor or numeric vector. Supplied <logical> #124

Closed gabdevocean closed 2 years ago

gabdevocean commented 2 years ago

Hi Gavin!

I've been loving gratia, but came across an issue when I tried to use draw() on an AR1 model.

(valRho <- acf(resid(df), plot=FALSE)$acf[2])

mod = bam(sensorvalue ~ 
                             s(azimuth, bs="cc") +
                             s(julianday) +
                             s(f_id, bs="re", k=17) + 
                            asfactor_tide_phase + 
                            s(distance),
                          data = df2,
                          AR.start=df2$start.event, rho=valRho, 
                          method="REML") 

draw(depth_mod_rho)

Gives the following error: Error: 'x' must be a factor or numeric vector. Supplied <logical>

draw() functions perfectly well when I take out the AR.start and rho, or when I remove asfactor_tide_phase (factor, no smooth). I'm a little new to gams, so this may very well be an issue with something else, but have you had previous issues with AR1 models that contain factors?

Thank you!

Best, Gabby

gavinsimpson commented 2 years ago

Thanks @gabdevocean --- I'll take a look. I should probably be ignoring variables like this that aren't related to the smooths, but hadn't considered AR models in bam() when writing the code that pulls out the variables needed for each smooth or the model as a whole.

This should be a reasonably easy fix. Will look this evening.

gabdevocean commented 2 years ago

Hi Gavin! draw() is now working with bam AR model, however, I can no longer plot my parametric effects (in any model), even with parametric=TRUE. Draw() output is only showing smooths.

gavinsimpson commented 2 years ago

Yep; this is an oversight on my part; I haven't implemented a new alternative to evaluate_parametric_terms() that works with the new smooth_estimates() functionality. Thanks to a potential covid exposure my weekend has suddenly become less busy so I'll see about implementing this as a priority

gabdevocean commented 2 years ago

No worries -- just wanted to let you know in case (slash to see if my code was the problem)! And oh no... Wishing you the best! 🤞

gavinsimpson commented 2 years ago

Reproducible example for the AR(1) with bam() - purely to develop a test:

#-- An AR(1) example using bam() with factor by --------------------------------
# from ?magic
## simulate truth
set.seed(1)
n <- 400
sig <- 2
x <- 0:(n-1)/(n-1)
## produce scaled covariance matrix for AR1 errors...
rho <- 0.6
V <- corMatrix(Initialize(corAR1(rho), data.frame(x = x)))
Cv <- chol(V)  # t(Cv) %*% Cv=V
## Simulate AR1 errors ...
e1 <- t(Cv) %*% rnorm(n, 0, sig) # so cov(e) = V * sig^2
e2 <- t(Cv) %*% rnorm(n, 0, sig) # so cov(e) = V * sig^2
## Observe truth + AR1 errors
f1 <- 0.2*x^11*(10*(1-x))^6+10*(10*x)^3*(1-x)^10
f2 <- (1280 * x^4) * (1- x)^4
df <- data.frame(x = rep(x, 2), f = c(f1, f2), y = c(f1 + e1, f2 + e2),
                 series = as.factor(rep(c("A", "B"), each = n)))
#rm(x, f1, f2, e1, e2, V, Cv)
AR.start <- rep(FALSE, n*2)
AR.start[c(1, n+1)] <- TRUE
## fit GAM using `bam()` with known correlation
## first just to a single series
m_ar1 <- bam(y ~ s(x, k = 20), data = df[seq_len(n), ], rho = rho,
             AR.start = NULL)
## now as a factor by smooth to model both series
m_ar1_by <- bam(y ~ series + s(x, k = 20, by = series), data = df, rho = rho,
                AR.start = AR.start)

This is working: I think the issue reported was unrelated to the AR aspects but rather related to variables being identified for the smooths, and OP was using a version prior to that being fixed.

The above example does show the parametric terms either so that needs to be fixed, but I'll do that in a separate issue.