R/cor_select.R

In fdrennan/test: What the Package Does (One Line, Title Case)

# This corresponds to the fourth box
cor_select <- function(transformed_data, variable) {
  num_times <- transformed_data %>%
    distinct(Time) %>%
    nrow()
  # Try AR1, CS, corSymm (unstructured), corExp (SP) correlations structures
  # and let user know which one is 'best' based on smallest AIC

  AR1 <- gls(
    data = transformed_data %>% filter(basic_model),
    model = as.formula(paste(variable, "~ TreatmentNew * Time")),
    correlation = corAR1(form = ~ 1 | SubjectID)
  )
  aic_AR1 <- c("AR1" = AIC(AR1))
  ARH1 <- gls(
    data = transformed_data %>% filter(basic_model),
    model = as.formula(paste(variable, "~ TreatmentNew * Time")),
    correlation = corAR1(form = ~ 1 | SubjectID),
    weights = varIdent(form = ~ 1 | Time)
  )
  aic_ARH1 <- c("ARH1" = AIC(ARH1))
  CS <- gls(
    data = transformed_data %>% filter(basic_model),
    model = as.formula(paste(variable, "~ TreatmentNew * Time")),
    correlation = corCompSymm(form = ~ 1 | SubjectID)
  )
  aic_CS <- c("CS" = AIC(CS))
  CSH <- gls(
    data = transformed_data %>% filter(basic_model),
    model = as.formula(paste(variable, "~ TreatmentNew * Time")),
    correlation = corCompSymm(form = ~ 1 | SubjectID),
    weights = varIdent(form = ~ 1 | Time)
  )
  aic_CSH <- c("CSH" = AIC(CSH))
  TOEP <- gls(
    data = transformed_data %>% filter(basic_model),
    model = as.formula(paste(variable, "~ TreatmentNew * Time")),
    correlation = corARMA(p = num_times - 1, q = 0, form = ~ 1 | SubjectID),
    weights = varIdent(form = ~ 1 | Time)
  )
  aic_TOEP <- c("TOEP" = AIC(TOEP))
  UN <- gls(
    data = transformed_data %>% filter(basic_model),
    model = as.formula(paste(variable, "~ TreatmentNew * Time")),
    correlation = corSymm(form = ~ 1 | SubjectID),
    weights = varIdent(form = ~ 1 | Time)
  )
  aic_UN <- c("UN" = AIC(UN))

  aic <- c(aic_AR1, aic_ARH1, aic_CS, aic_CSH, aic_TOEP, aic_UN)
  print(aic)
  message(paste("The best correlation structure is", names(aic[which.min(aic)])))
  return(names(aic[which.min(aic)]))
}