2
$\begingroup$

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.

$\endgroup$
2
+50
$\begingroup$

Your estimates of coefficient values in circular coordinates seem to be suffering from a common problem with periodic data, aliasing. You model the data with a von Mises distribution:

$$ f(x\mid\mu,\kappa)=\frac{e^{\kappa\cos(x-\mu)}}{2\pi I_0(\kappa)}$$

where $\mu$ is the measure of location, $\kappa$ is the measure of concentration, and $I_0(\kappa)$ is the modified Bessel function of order 0. You are modeling $\mu$, in radians, as a further function of species and temperature (covering a range from 0 to 5.5 in these data). You are particularly interested in how the slope of the relationship between $\mu$ and temperature differs among species.

Your estimates of these slopes are expressed in radians per unit temperature. Note that these 7 coefficient estimates are all near multiples of $\pi$: $-2\pi, -\pi, 0,\pi$. Furthermore, the average of the 7 slope coefficients for the species, the average change of $\mu$ per unit temperature, is close to a change of $-\pi$ radians per unit change in temperature! That doesn't make sense for these biological data.

Here's what I suspect is the problem. The cosine function used to model $\mu$ (from the observed values $x$ and the associated species and temperature values) returns to 0 at intervals of $\pi$. So it's possible that a model with very rapid changes of $\mu$ with temperature could fit the data as well as the much more modest changes that you know, based on the subject matter, are much more likely to be the case.

This is similar to aliasing while sampling periodic signals, in which you can't distinguish a low frequency from higher multiples of that frequency. One solution in digital signal processing is to impose a low-pass filter. I see two ways to impose such a filter here.

For one, although I don't do much Bayesian modeling, I suspect that your prior on the temperature coefficient was much too wide. Your prior knowledge suggests a very narrow range of possible values for these slopes, I would guess on the order of $\pm 0.5$ radians per unit temperature or less in these data. Set your prior accordingly.

As an alternative, pre-center the data for each species at 0 radians, keeping track of the associated offset for each species. Then just do a standard linear mixed model, so there is no periodicity beyond what's in your initial coding of the data. That will model intercepts as differences around the individual species offsets, so to get actual intercept estimates you will have to add back those offsets. But the slopes should be handled just fine, they won't suffer from aliasing, and they should be similar to what a narrow-prior Bayesian model would provide, given the similarity between the von Mises and normal distributions.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ I'm confused -- would it not be possible for a species to have a strongly negative relationship between median time and temperature? Or is there no meaningful distinction between positive and negative relationships in the circular model? $\endgroup$ – qdread Jul 16 at 15:17
  • 1
    $\begingroup$ @qdread in my answer yesterday I confused coef() for ranef(). I should have known better. The coef() values "are the sum of population-level effects and corresponding group- level effects." My answer erroneously took the values to represent the random effects instead. I'll need to look into this some more. Your result might have to do with the break at 0000/2400. I will keep a slightly edited version of this answer up for now, but should have something better in a couple of days when I have time to look at this more closely. $\endgroup$ – EdM Jul 16 at 15:56
  • $\begingroup$ I appreciate the help. I will go ahead and award the bounty since I have already spent almost all my reputation on it :-O and it must be awarded within the next 24 hours $\endgroup$ – qdread Jul 16 at 16:54
  • 1
    $\begingroup$ @qdread think I found the real problem, finally. $\endgroup$ – EdM Jul 17 at 13:42
  • $\begingroup$ This is a very thoughtful answer that makes a lot of sense. The issues you highlight would also explain why I had a lot of trouble with convergence. $\endgroup$ – qdread Jul 18 at 0:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.