---
title: "Multicomponent cosinor modeling"
author: "Oliver Jayasinghe and Rex Parsons"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{multiple-components}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r, echo=F}
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:

```{r, eval=F}
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
)
```


```{r, message=F, warning=F}
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.10205`represents the estimate for amplitude of `group 0` for the first component

-   `[group=1]:amp1 = 1.99964`represents the estimate for amplitude of `group 1` for the first component

-   `[group=0]:amp2 = 0.40175`represents the estimate for amplitude of `group 0` for the second component

-   `[group=1]:amp2 = 1.00057`represents the estimate for amplitude of `group 1` for the second component

-   *Similarly for acrophase estimates*

```{r}
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.

```{r, eval=F}
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
)
```

```{r, message=F, warning=F}
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.39795`and 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.09273`radians

-   The estimate for the acrophase of the first component for `group 1` is `[group=1]:acr1 = 0.99984`radians

-   The estimate for the acrophase of the second component is `acr2 = 1.50512`radians and is the same for both `group 0` and `group 1`

```{r}
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:

```{r}
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:

```{r, message=F, warning=F}
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:

```{r, eval=F}
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
)
```


```{r, message=F, warning=F}
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.