```
library(GLMMcosinor)
library(dplyr)
#>
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#>
#> filter, lag
#> The following objects are masked from 'package:base':
#>
#> intersect, setdiff, setequal, union
library(ggplot2)
```

GLMMcosinor allows specification of mixed models accounting for fixed and/or random effects. Mixed model specification follows the lme4 format. See their vignette, Fitting Linear Mixed-Effects Models Using lme4, for details about how to specify mixed models.

### Data with subject-level differences

To illustrate an example of using a model with random effects on the
cosinor components, we will first simulate some data with
`id`

-level differences in amplitude and acrophase.

```
f_sample_id <- function(id_num,
n = 30,
mesor,
amp,
acro,
family = "gaussian",
sd = 0.2,
period,
n_components,
beta.group = TRUE) {
data <- simulate_cosinor(
n = n,
mesor = mesor,
amp = amp,
acro = acro,
family = family,
sd = sd,
period = period,
n_components = n_components
)
data$subject <- id_num
data
}
dat_mixed <- do.call(
"rbind",
lapply(1:30, function(x) {
f_sample_id(
id_num = x,
mesor = rnorm(1, mean = 0, sd = 1),
amp = rnorm(1, mean = 3, sd = 0.5),
acro = rnorm(1, mean = 1.5, sd = 0.2),
period = 24,
n_components = 1
)
})
) |>
mutate(subject = as.factor(subject))
```

A quick graph shows how there are individual differences in terms of MESOR, amplitude and phase.

```
dat_mixed |>
ggplot(aes(times, Y, col = subject)) +
geom_point() +
geom_line() +
theme_bw()
```

### A single component model with random effects

For the model, we should include a random effect for the MESOR, amplitude and acrophase as these are clustered within individuals.

In the model formula, we can use the special `amp_acro[n]`

which represents the n^{th} cosinor component. In this case, we
only have one component so we use `amp_acro1`

. Following the
lme4-style mixed model formula, we add our random effect
for this component and the intercept term (MESOR) clustered within
subjects by using `(1 + amp_acro1 | subject)`

. The code below
fits this model

```
mixed_mod <- cglmm(
Y ~ amp_acro(times,
n_components = 1,
period = 24
) + (1 + amp_acro1 | subject),
data = dat_mixed
)
```

This works by replacing the amp_acro1 with the relevant cosinor
components when the data is rearranged and the formula created. The
formula created can be accessed using `.$formula`

, and shows
the `amp_acro1`

is replaced by the `main_rrr1`

and
`main_sss1`

(the cosine and sine components of time that also
appear in the fixed effects).

```
mixed_mod$formula
#> Y ~ main_rrr1 + main_sss1 + (1 + main_rrr1 + main_sss1 | subject)
#> <environment: 0x559311317d10>
```

The mixed model can also be plotted using `autoplot`

, but
some of the plotting features that are available for fixed-effects
models may not be available for mixed-effect models.

`autoplot(mixed_mod, superimpose.data = TRUE)`

The summary of the model shows that the input means for MESOR, amplitude and acrophase are similar to what we specified in the simulation (0, 3, and 1.5, respectively).

```
summary(mixed_mod)
#>
#> Conditional Model
#> Raw model coefficients:
#> estimate standard.error lower.CI upper.CI p.value
#> (Intercept) -0.03231528 0.18976318 -0.40424427 0.33961 0.86477998
#> main_rrr1 0.27614529 0.07553744 0.12809463 0.42420 0.00025644 ***
#> main_sss1 2.91288462 0.07833708 2.75934677 3.06642 < 2.22e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Transformed coefficients:
#> estimate standard.error lower.CI upper.CI p.value
#> (Intercept) -0.03231528 0.18976318 -0.40424427 0.33961 0.86478
#> amp1 2.92594481 0.07824558 2.77258630 3.07930 < 2e-16 ***
#> acr1 1.47627749 0.02580137 1.42570774 1.52685 < 2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

We can see that the predicted values from the model closely resemble the patterns we see in the input data.

```
data.frame(pred = predict(mixed_mod)) |>
cbind(dat_mixed) |>
ggplot() +
geom_point(aes(x = times, y = Y, col = subject)) +
geom_line(aes(x = times, y = pred, col = subject))
```

This looks like a good model fit for these data. We can highlight the
importance of using a mixed model in this situation rather than a fixed
effects only model by creating that (bad) model and comparing the two by
using the Akaike information criterion using `AIC()`

.

```
fixed_effects_mod <- cglmm(
Y ~ amp_acro(times, n_components = 1, period = 24),
data = dat_mixed
)
AIC(fixed_effects_mod$fit)
#> [1] 2765.136
AIC(mixed_mod$fit)
#> [1] 82.13652
```

Aside from not being able to be useful to see the differences between subjects from the model, we end up with much worse model fit and likely biased and/or imprecise estimates of our fixed effects that we are interested in!