Introduction to Circular Regression

CircularRegression models a circular response, such as an angle or movement direction, using circular explanatory variables. Angles are assumed to be in radians. The package follows the angular regression framework of Rivest, Duchesne, Nicosia and Fortin (2016), where the conditional mean direction is obtained from the argument of a resultant vector.

For a homogeneous angular model with reference direction \(x_1\), a simple two-term mean direction can be written as

\[ \mu_i(\beta) = \operatorname{atan2}\{ \sin(x_{1i}) + \beta z_i \sin(x_{2i}), \cos(x_{1i}) + \beta z_i \cos(x_{2i}) \}. \]

The formula syntax mirrors this construction. A plain term such as x1 represents an angular covariate. A term such as x2:z2 represents an angular covariate x2 weighted by a finite non-negative modifier z2.

Simulated example

library(CircularRegression)

wrap_angle <- function(x) atan2(sin(x), cos(x))

set.seed(123)
n <- 150
x1 <- runif(n, -pi, pi)
x2 <- runif(n, -pi, pi)
z2 <- runif(n, 0.2, 1.8)
beta <- 0.35

mu <- atan2(
  sin(x1) + beta * z2 * sin(x2),
  cos(x1) + beta * z2 * cos(x2)
)
y <- wrap_angle(mu + rnorm(n, sd = 0.12))

dat <- data.frame(y = y, x1 = x1, x2 = x2, z2 = z2)

The recommended high-level interface is circular_regression(). Its default method is "two_step": it first fits the consensus model, then fits the homogeneous angular model using the selected reference direction.

fit <- circular_regression(y ~ x1 + x2:z2, data = dat)
fit
#> Call:
#> circular_regression(formula = y ~ x1 + x2:z2, data = dat)
#> 
#> Method: two_step 
#> Number of observations: 150 
#> 
#> Call:
#> angular(formula = formula, data = data, reference = c("name", 
#>     ref_name), initbeta = initbeta, control = control_angular, 
#>     na.action = na.action)
#> 
#> Reference angle: x1 
#> Maximum cosine: 0.9928 
#> Concentration parameter: 69.26 
#> 
#> Parameters (non-reference terms):
#>       Estimate Std. Error z value  P(|z|>.)    
#> x2:z2 0.362624   0.011194  32.394 < 2.2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The usual S3 generics are available.

coef(fit)
#> [1] 0.3626242
head(fitted(fit))
#> [1] -1.6985278  1.7346547 -0.6174325  2.9131913  3.0755212 -2.5073521
head(residuals(fit))
#> [1]  0.02935248 -0.08183534 -0.08501247  0.08440631 -0.13025813  0.22079331
head(predict(fit))
#> [1] -1.6985278  1.7346547 -0.6174325  2.9131913  3.0755212 -2.5073521

Predictions for new angular covariates are returned as angles in radians.

new_dat <- dat[1:5, c("x1", "x2", "z2")]
predict(fit, newdata = new_dat)
#> [1] -1.6985278  1.7346547 -0.6174325  2.9131913  3.0755212

Model-specific functions

The original model-specific interfaces remain available. Use angular() when you want direct control of the homogeneous model and the reference direction.

fit_hom <- angular(
  y ~ x1 + x2:z2,
  data = dat,
  reference = c("name", "x1")
)

summary(fit_hom)
#> Call:
#> angular(formula = y ~ x1 + x2:z2, data = dat, reference = c("name", 
#>     "x1"))
#> 
#> Reference angle: x1 
#> Maximum cosine: 0.9928 
#> Concentration parameter: 69.26 
#> 
#> Residuals:
#>      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
#> -0.308781 -0.063432  0.011776  0.009991  0.089255  0.311168 
#> 
#> Parameters:
#>        estimate    stderr   z value P(|z|>.)
#> x2:z2  0.362624  0.011194 32.394221        0
#> 
#> Number of observations: 150

Use consensus() to fit the consensus formulation directly.

fit_cons <- consensus(y ~ x1 + x2:z2, data = dat)
summary(fit_cons)
#> Call:
#> consensus(formula = y ~ x1 + x2:z2, data = dat)
#> 
#> Maximum log-likelihood: 61.994 
#> AIC: -119.99 
#> BIC: -113.97 
#> 
#> Residuals:
#>      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
#> -0.330792 -0.075215 -0.005747  0.006829  0.103917  0.333623 
#> 
#> Kappa coefficients:
#>       estimate Std. Error  z value P(|z|>.)
#> x1    21.52345    2.45392  8.77105        0
#> x2:z2  6.10845    0.84379  7.23931        0
#> 
#> Number of observations: 150

Practical notes

All angles should be supplied in radians. Missing values are handled through the na.action argument passed to model.frame(). Modifiers in x:z terms must be finite and non-negative, because they act as vector weights in the mean-direction construction.

Reference:

Rivest, L.-P., Duchesne, T., Nicosia, A. and Fortin, D. (2016). A general angular regression model for the analysis of data on animal movement in ecology. Journal of the Royal Statistical Society: Series C (Applied Statistics), 65(3), 445-463.