R/ovary/ExploreMixedModels.R

In premodelling/mixed: Modelling Multilevel & Longitudinal Problems

# Title     : ExploreMixedModels.R
# Objective : Explore a variety of generalised mixed models in relation to the ovary data
# Created by: greyhypotheses
# Created on: 03/04/2022


#'
#' @param data: The data frame of ovary data
#'
ExploreMixedModels <- function (data) {


  # GLM
  model <- lme(fixed = follicles ~ cos(2*pi*Time) + sin(2*pi*Time), data = data,  method = 'ML')
  models <- list('1' = model)


  # GLMM
  model <- lme(fixed = follicles ~ cos(2*pi*Time) + sin(2*pi*Time), random = ~1 | Mare,
                 data = data,  method = 'ML')
  models <- append(models, list('2' = model))

  model <- lme(fixed = follicles ~ cos(2*pi*Time) + sin(2*pi*Time),
                      random = ~ (- 1 + cos(2*pi*Time)) | Mare,
                      data = data, method = 'ML' )
  models <- append(models, list('3' = model))

  model <- lme(fixed = follicles ~ cos(2*pi*Time) + sin(2*pi*Time),
                  random = ~ cos(2*pi*Time) | Mare,
                  data = data, method = 'ML')
  models <- append(models, list('4' = model))


  # GLMM + Autoregressive Process
  model <- lme(fixed = follicles ~ cos(2*pi*Time) + sin(2*pi*Time), random = ~1 | Mare,
                   correlation = corAR1(form = ~1 | Mare),
                   data = data,  method = 'ML')
  models <- append(models, list('5' = model))

  model <- lme(fixed = follicles ~ cos(2*pi*Time) + sin(2*pi*Time), random = ~1 | Mare,
                     correlation = corARMA(form = ~1 | Mare, p = 1, q = 1),
                     data = data,  method = 'ML')
  models <- append(models, list('6' = model))


  # clean-up
  rm(model)


  # metrics
  scores <- function (model, name) {
    values <- data.frame(AIC = stats::AIC(object = model),
               lnlikelihood = stats::logLik(object = model),
               BIC = stats::BIC(object = model), row.names = name)
    return(list(values))
  }

  metrics <- mapply(FUN = scores, model = models, name = '1':'6')
  metrics <- dplyr::bind_rows(metrics)



  return(list(models = models, metrics = metrics))

}