#StackBounty: #r #bayesian #mixed-model #circular-statistics #stan Interpretation of coefficients in mixed-effects model with circular …

Bounty: 50

I have a dataset from an experiment where wild ants were surveyed continuously for 24 hours under a number of temperature treatments (chambers). Whenever an ant was observed, the species of the ant and the time, rounded to the nearest hour, was recorded. This is circular data because the observations cover the entire 24-hour period (at least some ants are active at any time of day or night). I calculated the circular median time within each species and chamber. The null hypothesis is that an individual species does not change its median time with change in temperature.

I fit a mixed-effects model with the R package brms (a wrapper for Stan software) using a von Mises distribution (with default link functions) for the response, with temperature as a fixed effect and species as a random effect (each species has both random slope and random intercept). I had to transform the hour values to radians such that 0:00 maps to $-pi$, 12:00 maps to 0, and 24:00 maps to $pi$.

I am confused about how to interpret the species-level coefficients. I see the highest coefficient on a species that basically shows no change in response to temperature treatment but where the median time crosses midnight. I am concerned that I set up the model wrong or that I am interpreting the coefficients wrong.

data

library(circular)
library(brms)

dat <- structure(list(sp = c("apla", "apla", "apla", "apla", "apla", 
"apla", "apla", "apru", "apru", "apru", "apru", "apru", "apru", 
"apru", "apru", "apru", "apru", "apru", "apru", "caca", "caca", 
"caca", "caca", "caca", "caca", "caca", "caca", "caca", "caca", 
"caca", "cape", "cape", "cape", "cape", "cape", "cape", "cape", 
"cape", "cape", "cape", "cape", "cape", "crli", "crli", "crli", 
"crli", "crli", "crli", "crli", "crli", "crli", "crli", "crli", 
"crli", "fosu", "fosu", "fosu", "fosu", "fosu", "fosu", "fosu", 
"fosu", "fosu", "fosu", "fosu", "prim", "prim", "prim", "prim", 
"prim", "prim", "prim", "prim", "prim", "prim", "prim", "prim"
), chamber = c(1, 2, 4, 5, 6, 7, 8, 1, 2, 3, 4, 5, 6, 7, 8, 9, 
10, 11, 12, 1, 2, 3, 4, 5, 6, 7, 8, 9, 11, 12, 1, 2, 3, 4, 5, 
6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 
1, 2, 3, 4, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6, 7, 8, 9, 
10, 11, 12), temperature = c(3.5, 0, 2, 0, 1.5, 3, 5, 3.5, 0, 
4.5, 2, 0, 1.5, 3, 5, 5.5, 2.5, 0, 4, 3.5, 0, 4.5, 2, 0, 1.5, 
3, 5, 5.5, 0, 4, 3.5, 0, 4.5, 2, 0, 1.5, 3, 5, 5.5, 2.5, 0, 4, 
3.5, 0, 4.5, 2, 0, 1.5, 3, 5, 5.5, 2.5, 0, 4, 3.5, 0, 4.5, 2, 
1.5, 3, 5, 5.5, 2.5, 0, 4, 3.5, 0, 4.5, 2, 0, 1.5, 3, 5, 5.5, 
2.5, 0, 4), median_time = structure(c(11, 8, 14, 17.5, 16, 9, 
8, 20, 9, 13, 11, 9, 7, 9, 14, 6, 22, 7, 19, 23, 1, 23, 23, 2, 
0, 1, 23, 2, 1, 2, 15, 19.508716014162, 21, 20, 3, 12, 22, 21, 
1, 23, 0.999999999999999, 12, 23, 0.999999999999999, 0.999999999999999, 
17, 2, 3, 17, 0.999999999999999, 0.999999999999999, 16, 14, 0, 
12.3324823150422, 14, 13, 12, 10, 12, 18, 15, 9.65973937593219, 
15, 13, 0.999999999999999, 23, 0.999999999999999, 6, 21, 17, 
4, 0.999999999999999, 4, 4, 2, 3), medians = 11, circularp = list(
    type = "angles", units = "hours", template = "none", modulo = "2pi", 
    zero = 0, rotation = "counter"), class = c("circular", "numeric"
)), median_time_radians = c(-0.26179938779915, -1.0471975511966, 
0.523598775598299, 1.43989663289532, 1.0471975511966, -0.785398163397448, 
-1.0471975511966, 2.0943951023932, -0.785398163397448, 0.261799387799149, 
-0.26179938779915, -0.785398163397448, -1.30899693899575, -0.785398163397448, 
0.523598775598299, -1.5707963267949, 2.61799387799149, -1.30899693899575, 
1.83259571459405, 2.87979326579064, -2.87979326579064, 2.87979326579064, 
2.87979326579064, -2.61799387799149, -3.14159265358979, -2.87979326579064, 
2.87979326579064, -2.61799387799149, -2.87979326579064, -2.61799387799149, 
0.785398163397447, 1.96577725566528, 2.35619449019234, 2.0943951023932, 
-2.35619449019234, 0, 2.61799387799149, 2.35619449019234, -2.87979326579064, 
2.87979326579064, -2.87979326579064, 0, 2.87979326579064, -2.87979326579064, 
-2.87979326579064, 1.30899693899575, -2.61799387799149, -2.35619449019234, 
1.30899693899575, -2.87979326579064, -2.87979326579064, 1.0471975511966, 
0.523598775598298, -3.14159265358979, 0.0870436665320824, 0.523598775598299, 
0.261799387799149, 0, -0.523598775598299, -4.44089209850063e-16, 
1.5707963267949, 0.785398163397448, -0.612678798671407, 0.785398163397447, 
0.261799387799149, -2.87979326579064, 2.87979326579064, -2.87979326579064, 
-1.5707963267949, 2.35619449019234, 1.30899693899575, -2.09439510239319, 
-2.87979326579064, -2.0943951023932, -2.0943951023932, -2.61799387799149, 
-2.35619449019234)), class = "data.frame", row.names = c(NA, 
-77L))

model

priors <- prior_string("student_t(3, 0, 5)", class = "sd")

fit <- brm(median_time_radians ~ temperature + (temperature | sp), 
                            family = von_mises(), 
                            prior = priors,
                            data = median_times,
                            control = list(adapt_delta = 0.9),
                            chains = 2, iter = 7500, warmup = 5000, seed = 12345)

species level coefficients

coef(fit)$sp[,,'temperature']

      Estimate  Est.Error          Q2.5      Q97.5
apla -0.3153341 0.23798523  -0.892426917  0.0289234
apru  0.2865710 0.27866258   0.002069992  0.8184251
caca -6.5935606 3.15748526 -14.064381326 -2.5290273
cape  3.0701637 2.21674069  -0.253182098  7.5921491
crli  3.2702919 1.82584857   1.068027298  7.8987657
fosu  0.0571131 0.08858313  -0.101666321  0.2462271
prim -3.3404271 1.61870242  -7.440654851 -1.3915963

I’m confused why the species caca has the highest absolute value of its coefficient even though its median time barely changes — all its median values are between 23:00 and 2:00, so its trend crosses midnight but the times do not change much. I would appreciate any help interpreting these coefficients, or coefficients from a mixed-effects model with circular response more generally.


Get this bounty!!!

Leave a Reply

This site uses Akismet to reduce spam. Learn how your comment data is processed.