Multicomponent cosinor modeling"
In GLMMcosinor: Fit a Cosinor Model Using a Generalized Mixed Modeling Framework

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(GLMMcosinor)
withr::with_seed(
  50,
  {
    testdata_two_components <- simulate_cosinor(1000,
      n_period = 10,
      mesor = 7,
      amp = c(0.1, 0.4),
      acro = c(1, 1.5),
      beta.mesor = 4.4,
      beta.amp = c(2, 1),
      beta.acro = c(1, -1.5),
      family = "poisson",
      period = c(12, 6),
      n_components = 2,
      beta.group = TRUE
    )

    testdata_two_components_grouped <- simulate_cosinor(1000,
      n_period = 5,
      mesor = 3.7,
      amp = c(0.1, 0.4),
      acro = c(1, 1.5),
      beta.mesor = 4,
      beta.amp = c(2, 0.4),
      beta.acro = c(1, 1.5),
      family = "poisson",
      period = c(12, 6),
      n_components = 2,
      beta.group = TRUE
    )

    testdata_three_components <- simulate_cosinor(1000,
      n_period = 2,
      mesor = 7,
      amp = c(0.1, 0.4, 0.5),
      acro = c(1, 1.5, 0.1),
      beta.mesor = 4.4,
      beta.amp = c(2, 1, 0.4),
      beta.acro = c(1, -1.5, -1),
      family = "poisson",
      period = c(12, 6, 12),
      n_components = 3,
      beta.group = TRUE
    )
  }
)

{GLMMcosinor} allows specification of multi-component cosinor models. This is useful if there are multiple explanatory variables with known periods affecting the response variable.

Generating a two-component model

To generate a multi-component model, set n_components in the amp_acro() part of the formula to the desired number of components. Then, optionally assign groups to each component in the group argument. If only one group entry is supplied but n_components is greater than 1, then the single group entry will be matched to each component.

The period argument must also match the length of n_components, where the order of the periods corresponds to their assigned component. For example, if n_components = 2, and period = c(12,6), then the first component has a period of 12 and the second a period of 6. Similarly to the group argument, if only one period is supplied despite n_components being greater than 1, then this period will be matched to each component.

For example:

library(GLMMcosinor)
testdata_two_components <- simulate_cosinor(
  1000,
  n_period = 10,
  mesor = 7,
  amp = c(0.1, 0.4),
  acro = c(1, 1.5),
  beta.mesor = 4.4,
  beta.amp = c(2, 1),
  beta.acro = c(1, -1.5),
  family = "poisson",
  period = c(12, 6),
  n_components = 2,
  beta.group = TRUE
)

object <- cglmm(
  Y ~ group + amp_acro(
    time_col = times,
    n_components = 2,
    period = c(12, 6),
    group = c("group", "group")
  ),
  data = testdata_two_components,
  family = poisson()
)
object

In the output, the suffix on the estimates for amplitude and acrophase represents its component:

[group=0]:amp1 = 0.10205represents the estimate for amplitude of group 0 for the first component
[group=1]:amp1 = 1.99964represents the estimate for amplitude of group 1 for the first component
[group=0]:amp2 = 0.40175represents the estimate for amplitude of group 0 for the second component
[group=1]:amp2 = 1.00057represents the estimate for amplitude of group 1 for the second component
Similarly for acrophase estimates

autoplot(object)

If a multicomponent model has one component that is grouped with other components that aren't, the vector input for group must still be the same length as n_components but have the non-grouped components represented as group = NA.

For example, if wanted only the first component to have a grouped component, we would specify the group argument as group = c("group", NA)) . Here, the first component is grouped by group, and the second component is not grouped. The data was simulated such that the second component was the same for both groups.

testdata_two_components_grouped <- simulate_cosinor(
  1000,
  n_period = 5,
  mesor = 3.7,
  amp = c(0.1, 0.4),
  acro = c(1, 1.5),
  beta.mesor = 4,
  beta.amp = c(2, 0.4),
  beta.acro = c(1, 1.5),
  family = "poisson",
  period = c(12, 6),
  n_components = 2,
  beta.group = TRUE
)

object <- cglmm(
  Y ~ group + amp_acro(
    time_col = times,
    n_components = 2,
    period = c(12, 6),
    group = c("group", NA)
  ),
  data = testdata_two_components_grouped,
  family = poisson()
)
object

We would interpret the output the transformed coefficients as follows:

MESOR for group 0 is 3.69558.
MESOR difference to group 0 for group 1 is [group=1] = 0.31184
The estimate for the amplitude of the first component for group 0 is[group=0]:amp1 = 0.10752
The estimate for the amplitude of the first component for group 1 is [group=1]:amp1 = 1.99248
The estimate for the amplitude of the second component is amp2 = 0.39795and the same for both group 0 and group 1
The estimate for the acrophase of the first component for group 0 is [group=0]:acr1 = 1.09273radians
The estimate for the acrophase of the first component for group 1 is [group=1]:acr1 = 0.99984radians
The estimate for the acrophase of the second component is acr2 = 1.50512radians and is the same for both group 0 and group 1

autoplot(object, superimpose.data = TRUE)

In this example, it is not strictly necessary to specify group = c("group", NA)) since specifying group = c("group","group")still yields accurate estimates:

object <- cglmm(
  Y ~ group + amp_acro(
    time_col = times,
    n_components = 2,
    period = c(12, 6),
    group = c("group", "group")
  ),
  data = testdata_two_components_grouped,
  family = poisson()
)
object

If a multicomponent model is specified (n_components > 1) but the length of group or period is 1, then it will be assumed that the one group and/or period values specified apply to [all]{.underline} components. For example, if n_components = 2 , but group = "group", then the one element in this group vector will be replicated to produce group = c("group","group")which now has a length that matches n_components. The same applies for period.

For instance, the following two cglmm() calls fit the same models:

cglmm(
  Y ~ group + amp_acro(times,
    n_components = 2,
    period = 12,
    group = "group"
  ),
  data = testdata_two_components,
  family = poisson()
)


cglmm(
  Y ~ group + amp_acro(times,
    n_components = 2,
    period = c(12, 12),
    group = c("group", "group")
  ),
  data = testdata_two_components,
  family = poisson()
)

Generating a three-component model

The plot below shows a 3-component model with the simulated data overlayed:

testdata_three_components <- simulate_cosinor(
  1000,
  n_period = 2,
  mesor = 7,
  amp = c(0.1, 0.4, 0.5),
  acro = c(1, 1.5, 0.1),
  beta.mesor = 4.4,
  beta.amp = c(2, 1, 0.4),
  beta.acro = c(1, -1.5, -1),
  family = "poisson",
  period = c(12, 6, 12),
  n_components = 3,
  beta.group = TRUE
)

object <- cglmm(
  Y ~ group + amp_acro(times,
    n_components = 3,
    period = c(12, 6, 12),
    group = "group"
  ),
  data = testdata_three_components,
  family = poisson()
)
autoplot(object,
  superimpose.data = TRUE,
  x_str = "group",
  predict.ribbon = FALSE,
  data_opacity = 0.08
)

Generating models with n-components

Generating a model with n components simply involves setting n_components to be the number of desired components and ensuring that the period argument is a vector where each element corresponds the period of its respective component.