R/ovary/ExploreCyclic.R

In premodelling/mixed: Modelling Multilevel & Longitudinal Problems

# Title     : ExploreCyclic.R
# Objective : Investigating fixed effects (intercept + cyclic effects)
#             a. with random intercept
#             b. without random intercept
# Created by: greyhypotheses
# Created on: 28/03/2022



#' Preliminaries
#'

# the data
data(Ovary)
ovary <- Ovary
rm(Ovary)


# setting-up
model <- NULL
estimates <- NULL


#' III: A fixed cyclic effect, a random intercept
#'
model$a <-glmer(formula = follicles ~ 1 + cos(2*pi*Time) + sin(2*pi*Time) + (1|Mare),
                    family = poisson(link = 'log'), data = ovary)
estimates$a <- summary(object = model$a)

# ... the fixed effects coefficient estimates
estimates$a$coefficients

# ... random effects variance
unlist(estimates$a$varcor['Mare'])
unlist(VarCorr(model$a)['Mare'])

# ... AIC
estimates$a$AICtab['AIC']



#' IV: A fixed cyclic effect, no random intercept
#'
model$b <- glm(follicles ~ 1 + cos(2*pi*Time) + sin(2*pi*Time),
               family = poisson(link = 'log'), data = ovary)
summary(model$b)

# ... cf. standard errors
summary(model$a)$coefficients
summary(model$b)$coefficients

# ... no evidence in favour of the null hypothesis
#     Null Hypothesis: The model fits the data well,
#                      p.value = 1 - pchisqr(residual deviance, residual degrees of freedom, lower.tail = TRUE)
#                              = pchisqr(residual deviance, residual degrees of freedom, lower.tail = FALSE)
#     Alternative Hypothesis: Otherwise
pchisq(q = model$b$deviance, df = model$b$df.residual, lower.tail = FALSE)



#' V: fitted values vs. observations
#'
derivations <- ovary %>%
  dplyr::mutate(fitted.GLMM = as.numeric(fitted.values(object = model$a)),
                fitted.GLM = as.numeric(fitted.values(object = model$b)))

derivations %>%
  gather(key = 'Model', value = 'value', -c('Mare', 'Time')) %>%
  ggplot(mapping = aes(x = Time, y = value, colour = Model)) +
  geom_point(alpha = 0.35) +
  facet_wrap(~Mare) +
  theme_minimal() +
  theme(panel.spacing = unit(x = 2, units = 'lines'),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_line(size = 0.05),
        axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10),
        axis.title.x = element_text(size = 12, face = 'bold'), axis.title.y = element_text(size = 12, face = 'bold'))



#' VI: Diagnostics
#'
diagnostics <- ovary %>%
  dplyr::mutate(fitted.GLMM = as.numeric(fitted.values(object = model$a)),
                deviance.GLMM = residuals(object = model$a, type = 'deviance'))

# ... No discernible pattern, i.e., no evidence of inconstant variance (funnel pattern), no
# evidence of non-linearity (parabolic pattern) [ref.Methodologies, Lecture 5]
plot(x = diagnostics$fitted.GLMM, y = diagnostics$deviance.GLMM, frame.plot = FALSE,
     xlab = 'fitted value', ylab = 'deviance residuals', main = '')

qqnorm(y = diagnostics$deviance.GLMM,frame.plot = FALSE)