## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 6,
  fig.height = 4
)

## -----------------------------------------------------------------------------
library(CircularRegression)

data(bison)

d <- bison[seq_len(180), c("y.dir", "y.prec", "x.meadow", "z.meadow")]
d <- na.omit(d)

## -----------------------------------------------------------------------------
ref <- select_reference_angle(
  y.dir ~ y.prec + x.meadow:z.meadow,
  data = d
)
ref

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

summary(fit_hom)

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

summary(fit_cons)

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

coef(fit_two)

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

predict(fit_hom, newdata = new_d)
head(fitted(fit_hom))
head(residuals(fit_hom))

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

## -----------------------------------------------------------------------------
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
dec_fit$natural_parameters
pres_fit$natural_parameters
jam_fit$natural_parameters

## -----------------------------------------------------------------------------
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)
head(ranef.angular_re(fit_re))
head(predict(fit_re))