The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.
Circular Statistics and Distribution
Overlays
library(radiatR)
library(ggplot2)
Overview
The main radiatR vignette covers the path from raw tracking
files to a plotted set of headings. This vignette picks up from a
heading data frame and shows the analysis layer:
- Dispersion summaries —
circ_dispersion(), sector_summary()
- Parametric fits —
vonmises_fit(),
wrappedcauchy_fit()
- Hypothesis tests —
test_uniformity(),
test_mean_directions(), test_concentration(),
all with multiple-comparison correction
- Correlation —
circ_cor()
(circular-linear and circular-circular)
- Distribution overlays —
add_angle_rose(), add_vonmises_density(),
add_wrappedcauchy_density(),
add_circular_kde()
- Significance geometry —
add_critical_r() (Rayleigh / V-test circle) and
add_critical_v_line() (V-test boundary)
Every statistics function takes a data frame with a heading column in
radians (default name "heading") and returns a tidy data
frame, so the results drop straight into dplyr,
knitr::kable(), or further plotting.
A heading data frame
We derive one heading per trial from the bundled
cpunctatus dataset with the ring-crossing rule, then attach
each trial’s target half-width (arc) for grouping. This is
the same construction used in the main vignette.
data(cpunctatus)
hd <- derive_headings(cpunctatus, rule = "crossing",
circ0 = 0.2, circ1 = 0.4,
coords = "relative")
#> Warning: derive_headings(rule = 'crossing'): 25 of 251 headings (10.0%), from
#> 25 of 235 trials, are missing and are excluded from circular statistics.
#> Rule-based failures are often non-random and can bias results; inspect attr(x,
#> "missing_ids").
names(hd)[names(hd) == "id"] <- "trial_id"
# attach target half-width (arc) from the dataset, by trial
arc_map <- unique(cpunctatus@data[, c("trial_id", "arc")])
hd <- merge(hd, arc_map, by = "trial_id")
hd$arc <- factor(hd$arc)
# keep trials with a defined crossing heading
hd <- hd[is.finite(hd$heading), , drop = FALSE]
head(hd[, c("trial_id", "arc", "heading")])
#> trial_id arc heading
#> 1 10_1_1 10 0.2777794
#> 2 10_10_1 10 3.0376774
#> 3 10_11_1 10 6.0710737
#> 4 10_12_1 10 2.2074837
#> 5 10_13_1 10 0.9322361
#> 6 10_14_1 10 4.7955587
The heading column is in radians, reference-relative (0
= toward the target). Everything below operates on that column.
If you only need a summary and not the heading frame itself,
circ_summary_headings() derives headings from the tracks
with a chosen rule and returns the circular summary in one call.
group_by = NULL pools all trials;
group_by = "id" (the default) gives one row per trial:
circ_summary_headings(cpunctatus, rule = "distal", group_by = NULL)
#> mean_dir resultant_R kappa n group
#> 1 1.814938 0.04215183 0.08437866 235 all
Dispersion summaries
circ_dispersion() returns the mean direction, resultant
length R, and circular standard deviation. Grouping by
arc gives one row per condition.
circ_dispersion(hd, group_col = "arc")
#> arc mean_dir resultant_R circ_sd n
#> 1 10 1.980062802 0.25662879 1.6493178 35
#> 2 15 0.002875648 0.18589914 1.8344214 27
#> 3 20 2.960834000 0.07603959 2.2700225 24
#> 4 30 5.830185048 0.06016175 2.3709570 34
#> 5 40 0.075648794 0.64058845 0.9437882 19
#> 6 5 0.394217868 0.17579270 1.8646446 30
#> 7 50 6.109431857 0.44404605 1.2742268 24
#> 8 0 2.276557489 0.08000698 2.2475059 33
R runs from 0 (uniformly scattered) to 1 (all headings
identical); the circular SD moves the opposite way. For dense per-frame
heading series — gaze direction from a tethered subject, say —
sector_summary() bins the angles and reports dwell
proportions per sector:
sector_summary(hd, sectors = 8L)
#> sector mid_angle count proportion
#> 1 -158degrees -2.7488936 27 0.11946903
#> 2 -112degrees -1.9634954 18 0.07964602
#> 3 -68degrees -1.1780972 20 0.08849558
#> 4 -22degrees -0.3926991 42 0.18584071
#> 5 22degrees 0.3926991 45 0.19911504
#> 6 68degrees 1.1780972 13 0.05752212
#> 7 112degrees 1.9634954 35 0.15486726
#> 8 158degrees 2.7488936 26 0.11504425
Parametric fits
vonmises_fit() estimates the mean direction \(\mu\) and concentration \(\kappa\) by maximum likelihood, with
asymptotic standard errors and a confidence interval on \(\mu\):
vonmises_fit(hd, group_col = "arc")[, c("arc", "mu_deg", "kappa", "n")]
#> arc mu_deg kappa n
#> 1 10 113.4492417 0.5310863 35
#> 2 15 0.1647625 0.3784077 27
#> 3 20 169.6432921 0.1525210 24
#> 4 30 -25.9550030 0.1205419 34
#> 5 40 4.3343566 1.6868180 19
#> 6 5 22.5870200 0.3571578 30
#> 7 50 -9.9553394 0.9900343 24
#> 8 0 130.4371359 0.1605288 33
wrappedcauchy_fit() is the heavier-tailed alternative —
more robust when the data have outliers or weak directionality. Its
concentration \(\rho\) is bounded to
\([0, 1)\):
wrappedcauchy_fit(hd, group_col = "arc")[, c("arc", "mu_deg", "rho", "n")]
#> arc mu_deg rho n
#> 1 10 118.338958 0.25480084 35
#> 2 15 4.745037 0.27375349 27
#> 3 20 159.559101 0.09056095 24
#> 4 30 339.103222 0.08571126 34
#> 5 40 6.248322 0.75905693 19
#> 6 5 25.796263 0.15937039 30
#> 7 50 351.847489 0.68966189 24
#> 8 0 133.599471 0.08467292 33
Both return a row per group, so a quick merge() puts the
two concentration estimates side by side for comparison.
Hypothesis tests
Equal mean directions
test_mean_directions() is the Watson-Williams test — the
circular analogue of a one-way ANOVA on the mean angle. The omnibus form
asks whether any group differs:
test_mean_directions(hd, group_col = "arc")
#> n_groups statistic df1 df2 p_value test
#> 1 8 8.230315 7 218 6.70879e-09 Watson-Williams
Set pairwise = TRUE for all pairwise comparisons;
p_adjust is strongly recommended here because the number of
comparisons grows quickly:
pw <- test_mean_directions(hd, group_col = "arc",
pairwise = TRUE, p_adjust = "holm")
head(pw[order(pw$p_value_adj), ])
#> group1 group2 statistic df1 df2 p_value test p_value_adj
#> 14 20 30 68.71016 1 56 2.590906e-11 Watson-Williams 7.254536e-10
#> 22 30 0 45.35938 1 65 5.086470e-09 Watson-Williams 1.373347e-07
#> 6 10 50 31.92401 1 57 5.341214e-07 Watson-Williams 1.388716e-05
#> 4 10 40 24.34873 1 52 8.658990e-06 Watson-Williams 2.164747e-04
#> 1 10 15 16.68434 1 60 1.328071e-04 Watson-Williams 3.187369e-03
#> 13 15 0 16.32869 1 58 1.586887e-04 Watson-Williams 3.649841e-03
Equal concentrations
test_concentration() checks whether the groups are
equally concentrated (the circular analogue of a test for equal
variances), a key assumption behind the Watson-Williams test above:
test_concentration(hd, group_col = "arc")
#> statistic df p_value test
#> 1 14.23259 7 0.04719587 equal.kappa
Model selection
circ_model_select() ranks five candidate distributions
by AICc – uniform, unimodal,
axial, unimodal_uniform (a directed mode over
a uniform background), and bimodal (an asymmetric two-mode
mixture) – and reports Akaike weights. It answers whether a heading
sample is best described as undirected, singly directed, axially
oriented, directed-with-scatter, or genuinely bimodal.
hd_ms <- derive_headings(cpunctatus, rule = "distal")
arc_map <- unique(cpunctatus@data[, c("trial_id", "arc")])
names(arc_map)[1] <- "id"
hd_ms <- merge(as.data.frame(hd_ms), arc_map, by = "id")
circ_model_select(hd_ms, group_col = "arc")
#> arc model n k logLik AIC AICc BIC dAICc
#> 1 10 uniform 29 0 -53.29843 106.59687 106.59687 106.59687 0.0000000
#> 3 10 axial 29 2 -51.50436 107.00871 107.47025 109.74330 0.8733814
#> 4 10 unimodal_uniform 29 3 -51.04300 108.08600 109.04600 112.18789 2.4491293
#> 2 10 unimodal 29 2 -53.15387 110.30775 110.76929 113.04234 4.1724154
#> 5 10 bimodal 29 5 -49.34376 108.68752 111.29622 115.52400 4.6993497
#> 11 15 uniform 30 0 -55.13631 110.27262 110.27262 110.27262 0.0000000
#> 41 15 unimodal_uniform 30 3 -51.82459 109.64918 110.57226 113.85277 0.2996343
#> 21 15 unimodal 30 2 -53.31094 110.62187 111.06632 113.42427 0.7936929
#> 31 15 axial 30 2 -53.79386 111.58773 112.03217 114.39012 1.7595485
#> 51 15 bimodal 30 5 -50.88927 111.77855 114.27855 118.78454 4.0059258
#> 12 20 uniform 27 0 -49.62268 99.24536 99.24536 99.24536 0.0000000
#> 32 20 axial 27 2 -47.74062 99.48124 99.98124 102.07291 0.7358736
#> 52 20 bimodal 27 5 -43.99541 97.99081 100.84796 104.47000 1.6025959
#> 22 20 unimodal 27 2 -49.47757 102.95513 103.45513 105.54681 4.2097721
#> 42 20 unimodal_uniform 27 3 -49.47757 104.95513 105.99861 108.84264 6.7532506
#> 13 30 uniform 37 0 -68.00145 136.00290 136.00290 136.00290 0.0000000
#> 33 30 axial 37 2 -66.94814 137.89628 138.24923 141.11812 2.2463222
#> 23 30 unimodal 37 2 -67.37960 138.75921 139.11215 141.98105 3.1092481
#> 53 30 bimodal 37 5 -63.72366 137.44733 139.38281 145.50192 3.3799106
#> 43 30 unimodal_uniform 37 3 -67.37961 140.75921 141.48648 145.59196 5.4835806
#> 44 40 unimodal_uniform 20 3 -27.62695 61.25391 62.75391 64.24111 0.0000000
#> 54 40 bimodal 20 5 -27.19870 64.39739 68.68311 69.37605 5.9291975
#> 14 40 uniform 20 0 -36.75754 73.51508 73.51508 73.51508 10.7611738
#> 34 40 axial 20 2 -36.25323 76.50646 77.21235 78.49793 14.4584376
#> 24 40 unimodal 20 2 -36.67200 77.34400 78.04988 79.33547 15.2959744
#> 15 5 uniform 31 0 -56.97419 113.94838 113.94838 113.94838 0.0000000
#> 55 5 bimodal 31 5 -51.14723 112.29447 114.69447 119.46440 0.7460906
#> 25 5 unimodal 31 2 -56.81408 117.62817 118.05674 120.49614 4.1083591
#> 35 5 axial 31 2 -56.97248 117.94497 118.37354 120.81294 4.4251616
#> 45 5 unimodal_uniform 31 3 -56.81349 119.62698 120.51587 123.92894 6.5674914
#> 16 50 uniform 26 0 -47.78480 95.56961 95.56961 95.56961 0.0000000
#> 46 50 unimodal_uniform 26 3 -44.93562 95.87123 96.96214 99.64552 1.3925323
#> 26 50 unimodal 26 2 -47.13245 98.26490 98.78664 100.78109 3.2170284
#> 36 50 axial 26 2 -47.50728 99.01456 99.53630 101.53075 3.9666892
#> 56 50 bimodal 26 5 -44.86639 99.73279 102.73279 106.02327 7.1631803
#> 37 0 axial 35 2 -61.26559 126.53118 126.90618 129.64188 0.0000000
#> 47 0 unimodal_uniform 35 3 -60.53792 127.07584 127.85004 131.74189 0.9438571
#> 17 0 uniform 35 0 -64.32570 128.65139 128.65139 128.65139 1.7452142
#> 27 0 unimodal 35 2 -62.35996 128.71991 129.09491 131.83061 2.1887313
#> 57 0 bimodal 35 5 -59.49074 128.98148 131.05045 136.75822 4.1442667
#> weight
#> 1 0.4630456847
#> 3 0.2992068080
#> 4 0.1360824640
#> 2 0.0574904028
#> 5 0.0441746405
#> 11 0.3243473868
#> 41 0.2792194402
#> 21 0.2181032681
#> 31 0.1345641252
#> 51 0.0437657797
#> 12 0.4353638826
#> 32 0.3013418647
#> 52 0.1953678628
#> 22 0.0530532503
#> 42 0.0148731396
#> 13 0.5600667577
#> 33 0.1821617229
#> 23 0.1183246276
#> 53 0.1033478731
#> 43 0.0360990187
#> 44 0.9457263081
#> 54 0.0487816619
#> 14 0.0043551815
#> 34 0.0006857339
#> 24 0.0004511146
#> 15 0.5092330373
#> 55 0.3506752001
#> 25 0.0652826414
#> 35 0.0557192092
#> 45 0.0190899120
#> 16 0.5364612848
#> 46 0.2673953487
#> 26 0.1073914095
#> 36 0.0738214844
#> 56 0.0149304725
#> 37 0.3996279357
#> 47 0.2492871407
#> 17 0.1669888184
#> 27 0.1337762051
#> 57 0.0503199001
The five models mirror the Schnute-Groot family implemented by the
CircMLE package (Fitak & Johnsen 2017); radiatR fits them natively
(the two mixtures by numerical maximum likelihood) and uses
human-readable labels rather than the M1-M5
codes.
Circular correlation
circ_cor() measures the association between headings and
a covariate. With x_type = "linear" (the default) it
computes the circular-linear correlation — here, whether heading
direction is associated with the numeric target half-width:
hd$arc_num <- as.numeric(as.character(hd$arc))
circ_cor(hd, x_col = "arc_num", angle_col = "heading", x_type = "linear")
#> r n type statistic df p_value
#> 1 0.2175965 226 circular-linear 10.7007 2 0.00474648
The returned r is unsigned (association strength, 0–1);
the test statistic \(n r^2\) is
approximately \(\chi^2_2\). For two
angular variables, pass x_type = "circular" to get Fisher’s
\(\rho \in [-1, 1]\).
Circular regression
Where circ_cor() measures the strength of an
association, circ_regression() models it: it fits
the Fisher–Lee circular–linear regression of a heading on one or more
linear covariates, via a formula interface over
circular::lm.circular(). To show that the fit recovers a
known effect, we simulate trajectories whose mean heading shifts with a
predictor — the per-condition mean_slope controls that
shift — and then fit the model back:
cond <- data.frame(condition = "demo", n_trials = 120, ref_mean = 0,
concentration_base = 12, mean_slope = 0.6,
predictor_mean = 0, predictor_sd = 1)
s <- simulate_tracks(conditions = cond, n_points = 8, seed = 1)
reg <- s[!duplicated(s$trial_id), c("predictor", "final_heading")]
names(reg)[2] <- "heading"
fit <- circ_regression(reg, heading ~ predictor)
summary(fit)
#> term estimate se statistic p_value conf.low conf.high
#> 1 predictor 0.3281699 0.02007503 16.34717 0 0.2888236 0.3675163
The coefficient table recovers a positive,
significant slope for predictor, confirming the
simulated effect. Note the magnitude: the Fisher–Lee model links the
mean heading to the linear predictor through a \(2\arctan(\cdot)\) function, so the fitted
slope is on that link scale and is attenuated relative
to the plain mean_slope used to simulate. Interpret the
sign and significance of the coefficients directly, and
use predict() to read the effect back on the angular
(heading) scale:
predict(fit, data.frame(predictor = c(-1, 0, 1)))
#> [1] 5.6202847 6.2544774 0.6054847
The predicted mean heading tracks the predictor — rotating with its
value — exactly as the simulation set up.
The fitted relationship can be drawn straight onto the circular
panel: each covariate value becomes a mean-direction arrow,
colour-graded by the predictor, so you see the mean heading sweep as the
covariate increases.
radiate(headings_frame(hd, heading, units = "radians")) +
add_circ_mean(fitted_directions(fit, at = seq(-2, 2, length.out = 7)),
colour_col = "predictor")
Distribution overlays
The overlay layers draw an angular distribution in the same Cartesian
unit-disc space as radiate(), so they compose onto a bare
circular canvas with +. We build that canvas once:
canvas <- ggplot() + coord_fixed() +
add_circ(radius = 1) +
theme_void()
A rose diagram bins the headings into wedges; the
parametric and non-parametric density curves overlay on top. Giving
every layer the same scale aligns their radii, so the
shapes can be compared directly:
vm <- vonmises_fit(hd) # pooled fit (group_col = NULL)
wc <- wrappedcauchy_fit(hd)
canvas +
add_angle_rose(hd, bins = 12, scale = 0.8, fill = "grey80") +
add_circular_kde(hd, scale = 0.8, colour = "tomato") +
add_vonmises_density(vm, scale = 0.8, colour = "steelblue") +
add_wrappedcauchy_density(wc, scale = 0.8, colour = "darkorange")
The grey wedges are the empirical rose; the
steelblue curve is the von Mises fit,
darkorange the wrapped Cauchy, and
tomato a circular kernel density estimate. Where the
parametric curves track the rose closely the fit is good; a systematic
gap (especially in the tails) is the cue to prefer the heavier-tailed
wrapped Cauchy.
To show the individual headings rather than a
smoothed distribution, add_stacked_headings() fans
coincident headings into radial columns of dots, reducing the
overplotting of points stacked at the same angle. It needs a
headings_frame (the object derive_headings()
returns), so we use a fresh one here — hd above was demoted
to a plain data frame by the merge()/subset:
hd_frame <- derive_headings(cpunctatus, rule = "distal")
radiate(hd_frame, show_markers = FALSE) + add_stacked_headings(hd_frame)
step sets the gap between dots and
start_sep the inward offset of the first dot from the rim;
pass axial = TRUE to draw each datum at both poles.
Per-group overlays
Heading-summary overlays subdivide by a column to draw one element
per group, coloured by it – handy for comparing conditions on one plot.
The mean-direction arrow
(add_heading_arrow()) takes a colour_col; its
bootstrap confidence interval
(add_heading_interval()) takes a group_col.
Here both are split by target half-width (arc):
set.seed(1) # the interval is bootstrapped
canvas +
add_heading_arrow(hd, colour_col = "arc") +
add_heading_interval(hd, group_col = "arc", stat = "bootstrap_ci")
Each arc group gets its own mean vector and CI arc in a
matching colour. The grouping is independent of any faceting, so you can
(for example) colour by cohort while faceting radiate() by
treatment. radiate()’s own built-in arrow follows the same
idea via arrow_colour_col, drawing one arrow per group:
radiate(cpunctatus, group_col = "trial_id",
show_arrow = TRUE, arrow_colour_col = "arc")
Significance geometry
The remaining two helpers draw the decision boundary of a
significance test directly in resultant-length space (radius 0–1), so it
can be read against the observed mean vector.
add_critical_r() draws the Rayleigh critical
circle: the mean resultant length needed for significance at
level alpha, namely \(r_\text{crit} = \sqrt{-\log(\alpha)/n}\).
If the mean vector reaches past the circle, the headings are
significantly non-uniform.
disp <- circ_dispersion(hd)
mean_vec <- data.frame(x = disp$resultant_R * cos(disp$mean_dir),
y = disp$resultant_R * sin(disp$mean_dir))
canvas +
add_critical_r(hd, alpha = 0.05, test = "rayleigh") +
geom_segment(data = mean_vec,
aes(x = 0, y = 0, xend = x, yend = y),
arrow = arrow(length = unit(0.15, "inches")),
linewidth = 1)
test_uniformity(hd, test = "rayleigh")
#> statistic p_value n test
#> 1 0.1298181 0.02217652 226 rayleigh
The arrow tip lies well outside the dashed critical circle, matching
the tiny Rayleigh p-value above.
add_critical_v_line() is the corresponding boundary for
the V-test, which tests uniformity against a
specified direction mu0 (here 0, i.e. toward the
target). Unlike the Rayleigh circle, the V-test boundary is a straight
line perpendicular to mu0: significance requires the mean
vector’s projection onto mu0 to exceed \(c = z_\alpha / \sqrt{2n}\). Set
show_region = TRUE to shade the rejection side.
canvas +
add_critical_v_line(hd, mu0 = 0, alpha = 0.05, show_region = TRUE) +
geom_segment(data = mean_vec,
aes(x = 0, y = 0, xend = x, yend = y),
arrow = arrow(length = unit(0.15, "inches")),
linewidth = 1)
When the experiment has an a priori expected direction (the
target), the V-test is more powerful than the omnibus Rayleigh test, and
the line makes the one-sided nature of the decision visible.
Periodic data: time of day or year
Headings are not the only angles. Any periodic
variable — the time of day an event occurs (chronobiology), the day of
the year (phenology), or a generic cycle — maps onto the circle, and
then every tool above applies. as_angle() does the mapping
in the same clock/calendar convention as the circumference labels, so a
06:00 event lands on the "6" o’clock
label.
set.seed(1)
# 200 events: a strong morning cluster (~09:00) and a weaker evening one (~21:00)
times <- as.POSIXct("2020-06-01", tz = "UTC") +
c(rnorm(150, 9 * 3600, 1.2 * 3600), rnorm(50, 21 * 3600, 3600))
events <- headings_frame(data.frame(heading = as_angle(times, "day")),
col = heading, units = "radians")
# is the timing non-uniform over the day?
test_uniformity(events, test = "rayleigh")
#> statistic p_value n test
#> 1 0.4789747 1.183343e-20 200 rayleigh
radiate(events, angle_labels = "hours")
The dots cluster at their clock positions, and the Rayleigh test
confirms the day is not uniform. Use period = "year" for
Date/POSIXct day-of-year data (pair it with
angle_labels = "months"), or a positive number for a custom
cycle length.
Where next
- The main vignette covers building these heading
data frames and the trajectory/heading-overlay plotting layers.
- The Loaders vignette covers reading tracking-tool
exports into the
Tracks objects these analyses start
from.
- Every function above is documented individually under
Reference, grouped by role (summaries, parametric fitting,
correlation, hypothesis tests, and distribution overlays). ```
These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.