use GAMM models for full circular wind direction data

using a GAM modeling approach enables to include circular data as fixed effects. I wish to analyze the effect of wind on the flight behaviour of birds. In order to do so I want to add wind speed and wind direction as two distinct factors to my model. While wind speed is linear, I would use a smooth term for the wind direction taking its circularity into account. However, as far as I understand it, this transforms the data to a linear scale, but the model does not take into account that 360° is similar to 0°, is that right? Interestingly, I am not able to find evidence for this approach being actually flawed or wrong, still I also don't see anyone using it in my field of research. I am aware that one possibility to overcome the problem is combining wind speed and direction and calculate a GLM model including a northward and eastward wind component. But I am interested in whether the GAM approach would also give accurate results or is problematic in any way and if so, why.

• If the aim is to model the effect on behaviour of birds using features of the wind, then you don't have a circular response variable. Can you clarify what you are trying to model with what? A cyclic smooth s(wind_direction, bs = "cc") with knots = list(wind_direction = c(0,360)) also passed to gam() will handle the fact that the predictor wind_direction is a circular variable with equivalence at 0 and 360 degrees. Oct 14 '21 at 9:48
• the model looks like: does_the_bird_fly ~ wind_speed+s(wind_direction, bs= "cc"), family = binomial, knots= list(wind_direction = c(0,360))
– Some
Oct 15 '21 at 7:21
• Then you don't have "circular data as response variables." The response is binary. You could edit your question to fix this as different approaches would be required to handle a circular response. Circular responses require us to specify a conditional distribution for the response that meets the requirement of treating values close to 0 and 360 as being, well, close. An example of such a distribution would be the von Mises en.wikipedia.org/wiki/Von_Mises_distribution . Oct 15 '21 at 8:53
• In your case the constraints on the response are that observations are individual binary events, suggesting a Bernoulli distribution, or binary events summed over some number of observations, suggesting a binomial distribution, the former being a special case of the latter. Oct 15 '21 at 8:53

There is nothing per se wrong or flawed about representing wind direction, a circular variable, via a cyclic cubic regression spline. That this approach is not as common (or used at all) in your field could simply be due to people in your field not being as familiar with modern GAMs, splines, or cyclic splines in particular. They may be familiar with the decomposition using cosine and sine functions and hence just following the crowd or familiar.

More generally, GAMs are just GLMs. GAMs are certainly fancy GLMs but they are still GLMs. If I were to generate the model matrix that the GAM is creating and plug it into glm() (via say glm.fit()), I'd be fitting a GAM. I wouldn't be doing any smoothness selection if I use glm.fit(), but I would be fitting a GAM as a GLM. What sets modern GAMs apart, is the way we choose the wiggliness of the fitted functions; in modern GAMs we set some upper limit on the expected wiggliness (via the k argument in {mgcv}) — there are other ways if we're fitting modern GAMs in a fully Bayesian setting — and then we apply a wiggliness penalty to the fit (the log likelihood), which penalises overly complex ≡ wiggly estimated functions. This is a form of penalization or regularization, similar to the ridge or lasso penalties that you might have encountered in the GLM world.

If we take a step back and think about what both approaches (cosine-sine decomposition and cyclic splines) are doing, they are both representing the information in the data in some way that respects the circular nature of the original data. These representations don't suffer from the circularity feature of the original data and hence are more amenable to be used as covariates in a model.

Both the cosine and sine decomposition and the cyclic CRS approach can be viewed as basis expansions of the original data, in the same way that adding polynomials of a covariate $$x$$ ($$x^0 = 1$$, $$x^1 = x$$, $$x^2$$, $$x^3$$, etc) to a model can represent non-linear relationships. The Cyclic CRS does the same thing as the cosine, sine decompsition qualitatively; it provides continuous variables that encode in some way features of the original data, via a basis expansion.

The Cyclic CRS is a richer decomposition but it isn't as immediately intuitive as the northlyness and easterlyness interpretation say of the cosine sine decomposition. We can get around that by plotting the estimated smooth function (perhaps using polar coordinates so you really see the circular nature of the effect of the covariate on the response).

The main issues I have encountered with using cyclic CRS are

• that you really have to have data that cover the full range of (in this case) wind direction, and
• that you must specify the end points where values of the covariate map onto one another (0 and 360 in the case of wind direction).

If you don't have data that cover the full range and don't specify the join point, the spline is going to try to make the data join at incorrect places. If you have data covering the interval 20-270 degrees, the spline will be constrained to make 20 == 270, thus introducing a bias.

If you do specify the end points where values of the covariate wrap around to map onto one another again, but you don't have data that cover the full range of the variable, then you may introduce uncertainty and noise into the model estimates. In the case of these cyclic splines, if you observed data values in the range 90-270 and represented this directional variable as cyclical, you are going to have the effect on the response for a lot of the range of this variable totally unsupported by data; this would involve extrapolation beyond the observed data and this could affect the fit to the data where you have observations. It may well be better to ignore the circularity in wind direction if you only observe values over a restricted part of the range of values that could be observed.

You also have to be careful about choosing the end points as it can make a difference to fitted and predicted values. For wind direction we have a clear point at which the data values wrap around: 0 ≡ 360! However, for a variable related to seasonality, this equivalence point is not always as immediately obvious.

If you had monthly data for example with months measured as 1, 2, ..., 12, you shouldn't set the end points to be c(1, 12) because December isn't exactly the same as January, unlike 0 ≡ 360. So what do you do? In some limited testing, IIRC, setting the knots to be c(0.5, 12.5) seemed to produce fits with lowest error against simulated data (compared with say c(1, 13)), perhaps because you are only pushing the end points only a little bit beyond each end of the observed data rather than a bit more at one end only. I don't know if this holds more broadly, and we stopped looking into it as the differences in RMSEP were small, so caveat emptor.

• (+1) I would not mess myself with cyclic splines unless and until I was clear that cosines and sines didn't work well enough. Indeed, I want to emphasise that you are not restricted to a single pair of sine and cosine terms. With environmental (including ecological) data it is often a very good first approximation that one direction promotes something most and the opposite direction promotes it least, so a single pair of sine and cosine terms often work well. Your set-up may be more complicated, e.g. winds from different directions depending on time of day or time of year. Oct 15 '21 at 10:40
• It's a bit odd that using trigonometric predictors in regression-like models is not always discussed at length even in texts on directional statistics. I see far more problems in which directional variables are predictors than problems in which they are responses, but that may be personal. I tried my hand at a tutorial in stata-journal.com/article.html?article=st0116 The Stata content should not rule out this being of some use or interest to people using other software. Oct 15 '21 at 10:43
• Coming back to this late. I think this is personal preference; you're using one particular basis expansion of the covariate, I'm using another basis expansion, & there are likely other approaches we could take. I've used Fourier basis expansions before & in a number of cases, using a single pair of sin/cos functions didn't capture the asymmetry of data. So I'd add more functions but now I'm in a situation where I need to do model selection etc. Typically these terms are in models with many other terms including penalised splines so a penalised cyclic spline was a better/expedient solution. Nov 3 '21 at 8:45
• Splines are Marmite: I can sense which way you jump. Wind direction can be simple or a mix depending on season and location, e.g up and down valley winds, onshore and offshore winds. Nov 3 '21 at 12:31

If the effect of the cyclical component on the response can possibly change smoothly over time, an alternative is to use a smooth modulation model (like the one proposed by Eilers et al. (2008) and Marx et al. (2010)).

The only difference with the standard cosine-sine decomposition is that we can allow the coefficients associated to the periodic terms to vary in a smooth way over, for example, time. The model itself is a special case of GAM (also known as varying-coefficient model).

To make this suggestion a bit more concrete, I add below a small simulated example (I hope it is clear enough).

# Packages
library(splines)
library(Matrix)
library(mgcv)
set.seed(123)

# Simulate some data
N  = 360
x  = seq(1, N, by = 1)
sx = sin(x/7) * (0.15 * x)
cx = cos(x/7) * (-0.5 * x)
tr = 2 + 1.5 * x
ys = tr + cx + sx
y  = rnorm(N, ys, 20)

# Create sin-cos independent variables for gam model
sv = sin(x/7)
cv = cos(x/7)

# Estimate gam
mod = gam(y ~ s(x, bs = 'ps') + s(x, by = sv, bs = 'ps') + s(x, by = cv, bs = 'ps'))
fit = predict(mod, type = 'response')
trm = predict(mod, type = 'terms')

# Plot estimates
par(mfrow =c(3, 1), mar= c(2,2,2,2))
plot(x, y, main = 'Data and fit')
lines(x, fit, col = 2)

plot(x, tr, main = 'Linear trend and its estimate')
lines(x,  trm[, 1]  + attr(trm, 'constant'), col = 2)

plot(x, sx + cx, main = 'Modulated sine-cosine term and fit')
lines(x,  trm[, 2] + trm[, 3] , col = 2) 