gamm_ex.R

# GAM analysis -----------------------------------------------------------------
#
# Resources: 
# - https://arxiv.org/abs/1703.05339
# - https://github.com/soskuthy/gamm_intro


# set modality to ordered variable and relevel so production is reference
prodShadow_gamm <- prodShadow %>% 
  filter(., group == 'midd', picName == 0, 
            labID != "approx" & labID != "error") %>% 
  dplyr::select(., id, session, days, place, voicing, phon, vowelHeight, 
                   word, repN, vot, labID) %>% 
  mutate(., id = as.factor(id), 
            place_ord = as.ordered(place), 
            place_ord = fct_relevel(place_ord, 'bilabial'), 
            phon_ord = as.ordered(phon), 
            phon_ord = fct_relevel(phon_ord, 'b'), 
            voicing_ord = as.ordered(voicing), 
            voicing_ord = fct_relevel(voicing_ord, 'voiceless')) %>% 
  group_by(., voicing) %>% 
  mutate(., votS = (vot - mean(vot)) / sd(vot)) %>% 
  ungroup(.)

# dummy code voicing_ord
contrasts(prodShadow_gamm$voicing_ord) <- "contr.treatment" 

gamm_1 <- readRDS(here('cache', 'gamm', 'gamm_1.rds'))
gamm_0 <- readRDS(here('cache', 'gamm', 'gamm_0.rds'))

if(F){
  # Fit full model with parametric place variable, reference smooth,  
  # difference smooth, and random smooth
  gamm_1 <- gam(
    votS ~ voicing_ord +                                    # parametric predictor
           s(days, bs = "cr", k = 4) +                      # reference smooth
           s(days, by = voicing_ord, bs = "cr", k = 4) +    # difference smooth
           s(days, id, bs = "fs", xt = "cr", m = 1, k = 4), # random smooth
    data = prodShadow_gamm, method = "ML") 

  # Fit null model with no parametric modality predictor nor difference smooth
  # Essentially take out anything related to place
  gamm_0 <- gam(
    votS ~ s(days, bs = "cr", k = 4) +                      # reference smooth
           s(days, id, bs = "fs", xt = "cr", m = 1, k = 4), # random smooth
    data = prodShadow_gamm, method = "ML") 

  saveRDS(gamm_1, "./cache/gamm/gamm_1.rds", compress = 'xz')
  saveRDS(gamm_0, "./cache/gamm/gamm_0.rds", compress = 'xz')
  }


acf_plot(resid(gamm_1), split_by = list(prodShadow_gamm$id))