Skip to contents

mixed-models

Oliver Jayasinghe and Rex Parsons

Source: vignettes/mixed-models.Rmd
mixed-models.Rmd
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 nth 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!