Angular Regression Workflow

This vignette gives a compact package workflow using only datasets distributed with CircularRegression. It focuses on the package interface rather than on a full ecological analysis.

The response y.dir and angular covariates are in radians. The modifier z.meadow is non-negative and enters through the weighted term x.meadow:z.meadow.

Select a Reference Direction

The homogeneous angular model fixes one plain angular term as the reference. The package can select the reference direction using the empirical mean-cosine criterion.

ref <- select_reference_angle(
  y.dir ~ y.prec + x.meadow:z.meadow,
  data = d
)
ref
#> Reference angle selection (Rivest et al., 2016 criterion)
#> Selected reference: y.prec 
#> 
#>   angle mean_cos
#>  y.prec   0.9818

Fit Homogeneous and Consensus Models

The homogeneous model can be fitted directly once the reference is known.

fit_hom <- angular(
  y.dir ~ y.prec + x.meadow:z.meadow,
  data = d,
  reference = c("name", ref$ref_name)
)

summary(fit_hom)
#> Call:
#> angular(formula = y.dir ~ y.prec + x.meadow:z.meadow, data = d, 
#>     reference = c("name", ref$ref_name))
#> 
#> Reference angle: y.prec 
#> Maximum cosine: 0.9832 
#> Concentration parameter: 30.03 
#> 
#> Residuals:
#>       Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
#> -8.645e-01 -6.029e-02  3.576e-05  6.479e-03  7.892e-02  7.836e-01 
#> 
#> Parameters:
#>                     estimate     stderr    z value P(|z|>.)
#> x.meadow:z.meadow 0.00099798 0.00031944 3.12410811   0.0018
#> 
#> Number of observations: 164

The consensus model estimates one concentration-like coefficient for each directional term.

fit_cons <- consensus(
  y.dir ~ y.prec + x.meadow:z.meadow,
  data = d
)

summary(fit_cons)
#> Call:
#> consensus(formula = y.dir ~ y.prec + x.meadow:z.meadow, data = d)
#> 
#> Maximum log-likelihood: -68.139 
#> AIC: 140.28 
#> BIC: 146.48 
#> 
#> Residuals:
#>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
#> -0.55486 -0.09113  0.03415  0.05666  0.17710  0.72394 
#> 
#> Kappa coefficients:
#>                    estimate Std. Error   z value P(|z|>.)
#> y.prec            1.4160378  0.3360335 4.2139777        0
#> x.meadow:z.meadow 0.0179803  0.0032064 5.6076674        0
#> 
#> Number of observations: 164

The default high-level workflow fits the consensus model first, then fits the homogeneous model with the selected reference direction.

fit_two <- angular_two_step(
  y.dir ~ y.prec + x.meadow:z.meadow,
  data = d
)

coef(fit_two)
#> [1] 0.0009996682

The coefficient for x.meadow:z.meadow is interpreted relative to the selected reference direction. In this subset, its sign and magnitude describe how the weighted meadow direction contributes to the fitted mean resultant vector after the reference direction has coefficient one.

Predictions and Residuals

Predictions are angles in radians. Residuals are wrapped circular differences.

new_d <- d[seq_len(5), c("y.prec", "x.meadow", "z.meadow")]

predict(fit_hom, newdata = new_d)
#> [1] -2.251274 -2.251129 -2.251053 -2.250919 -2.248901
head(fitted(fit_hom))
#> [1] -2.251274 -2.251129 -2.251053 -2.250919 -2.248901 -2.247789
head(residuals(fit_hom))
#> [1] -0.01745609 -0.01765817 -0.01766949 -0.01478973 -0.01523846 -0.01635804

A lightweight base-R residual plot is often sufficient for a first diagnostic check.

plot(
  fitted(fit_hom),
  residuals(fit_hom),
  xlab = "Fitted direction",
  ylab = "Circular residual",
  pch = 19,
  col = "gray30"
)
abline(h = 0, lty = 2)

Special-Case Wrappers

The special-case wrappers keep the two-step fitting machinery but return natural-parameter summaries for model forms used in the circular-regression literature.

mean_fit <- meanDirectionModel(
  data = d["y.dir"],
  response = "y.dir"
)

dec_fit <- decentredPredictorModel(
  data = d[c("y.dir", "y.prec")],
  response = "y.dir",
  w = "y.prec"
)

pres_fit <- presnellModel(
  data = na.omit(bison[500:579, c("y.dir", "z.gap")]),
  response = "y.dir",
  w = "z.gap"
)

jam_fit <- jammalamadakaModel(
  data = d[c("y.dir", "y.prec")],
  response = "y.dir",
  w = "y.prec"
)

mean_fit$natural_parameters
#>      estimate se_model se_robust
#> mu -0.1110785       NA 0.0585229
dec_fit$natural_parameters
#>             estimate   se_model se_robust
#> r_hat     0.02301754 0.03437148 0.0269269
#> alpha_hat 1.99883664 0.92772363 0.7881173
#> beta_hat  1.99839768 0.91811141 0.7871155
pres_fit$natural_parameters
#>        estimate  se_model se_robust
#> beta2 -1.678551  93827.63  77463.39
#> beta3  5.217216 724208.19 598340.81
#> beta4 -5.213951 723756.55 597981.78
jam_fit$natural_parameters
#>          estimate se_model    se_robust
#> gamma2    8840299       NA 4.440917e+13
#> gamma3 -176048747       NA 9.132929e+14
#> gamma4    4669233       NA 2.356235e+13
#> gamma5 -163036174       NA 8.471331e+14
#> gamma6    2929010       NA 1.464703e+13

Random-Intercept Fit

The random-intercept model is separate from circular_regression() in this version. It can be used when clustered circular observations are available. The following small example uses the Sandhopper data and converts degree variables to radians before modeling.

data(Sandhopper)

sh <- Sandhopper[seq_len(24), ]
sh$y <- sh$LN1 * pi / 180
sh$ref <- sh$Azimuth * pi / 180
sh$wind <- sh$DirW * pi / 180
sh$wind_speed <- sh$SpeedW / max(sh$SpeedW, na.rm = TRUE)

fit_re <- angular_re(
  y ~ ref + wind:wind_speed,
  data = sh,
  cluster = sh$Anim,
  control = list(maxit = 50, reltol = 1e-8)
)

coef(fit_re)
#> wind:wind_speed         kappa_e         kappa_a 
#>    2.368792e-02    1.529995e+04    1.204056e+00
head(ranef.angular_re(fit_re))
#>          1          2          3          4          5          6 
#>  1.5589175  0.8046885  0.1060834 -1.2383562 -0.2396195  0.4169630
head(predict(fit_re))
#> [1] 2.193462 2.214674 2.250103 2.285628 2.316580 2.358079

This example is intentionally small so that it remains suitable for package checks and pkgdown builds.