library(tidyverse)
library(ggpubr)
library(rstatix)
library(knitr)
library(kableExtra)
library(nlme)
library(emmeans)
library(standardize)
library(sjPlot)
library(gridExtra)
library(effectsize)
data_type <- "AMST"
raw_data <- read_csv(
  paste0("../../Data/RELRESrms", data_type,".csv"),
  col_types = cols(RMS = col_double())
)

stimfreq <- c(5, 10)
metrics <- c("RMS", "logRMS")
red_color <- "#ffb3b3"
green_color <- "#d6f5d6"
include_trials <- 
  read_csv("AnimalTrial_AM_10.csv") %>% 
  mutate(ClickFreq = 10) %>% 
  bind_rows(
    read_csv("AnimalTrial_AM_5.csv") %>% 
    mutate(ClickFreq = 5)
  ) %>% 
  unique()
all_data <- raw_data %>% 
  filter(
    OrderofClick == 1, # just 1st click
    grepl('_1', Measurement), # pre/post Laser
    RMS > 0
) %>% 
  inner_join(include_trials) %>% 
  mutate(
    Measurement = case_when(
      grepl('pre', Measurement) ~ "preLaser",
      TRUE ~ "postLaser"),
    Measurement = factor(Measurement, levels = c("postLaser", "preLaser")),
    Animal = factor(Animal),
    TrialNumber = factor(TrialNumber),
    Group = factor(Group, levels = c("KIT", "KIC", "KIV")),
    logRMS = log(RMS)
  )

Results

for (freq in stimfreq){

  cat("##", freq, "Hz Frequency {.tabset .tabset-fade .tabset-pills} \n")
  
  cat('\n\n')
  
  model_data <- all_data %>% 
    filter(ClickFreq == freq)
    
  for(metric in metrics){
    
    cat("###", metric, "{.tabset .tabset-fade .tabset-pills} \n")
    
    cat("#### Data overview\n")
    
    cat('\n\n')
    
    model_data %>% 
      select(-c(OrderofClick, ClickFreq)) %>% 
      sample_n(5) %>% 
      kable(escape = F, digits = 6,
          caption = "Sample of the data") %>%
      kable_classic_2(full_width = F) %>% 
      kable_styling() %>% 
      print()
    
    cat('\n\n')
    
    # plots with data overview
    box_plot <- ggboxplot(
      model_data, x = "Group", y = metric,
      color = "Measurement", palette = "lancet", 
      add = "jitter",
      main = paste0(metric, " ~ Group | Measurement")
    )
    plot(box_plot)
    
    box_plot2 <- ggboxplot(
      model_data, x = "Measurement", y = metric,
      color = "Group", palette = "jco", 
      add = "jitter",
      main = paste0(metric, " ~ Measurement | Group")
    )
    plot(box_plot2)
    
    cat('\n\n')
    
    cat("#### LMM\n")
    
    model_data <- model_data %>% 
      mutate(AnimalTrial = paste(Animal, TrialNumber, sep = "_"))
    
    model_final <- lme(
      fixed = as.formula(paste0(metric, "~ Group*Measurement")),
      random = (~1|Animal/TrialNumber),
      data = model_data)
    
    cat(paste0("`", ". ", capture.output(summary(model_final)), "`  \n"))
    
    cat('\n\n')
    
    cat("#### Residuals check\n")
    
    res <- residuals(model_final)
    
    p1 <- ggplot(data = data.frame(residuals = res)) +
      geom_histogram(aes(x = residuals), bins = 25,
                     fill = "lightblue", color = "black") +
      labs(title = "Histogram of Residuals")
    
    p2 <- ggplot(
      data = data.frame(residuals = res),
      aes(sample = residuals)) +
      stat_qq() +
      stat_qq_line() +
      labs(title = "QQ Plot of Residuals")
    
    grid.arrange(p1, p2, nrow = 1)
    
    p1 <- plot(model_final, resid(.) ~ fitted(.), abline=0, main = "All Data")
    p2 <- plot(model_final, resid(.) ~ fitted(.)|Animal, abline=0, main = "By Animal")
    
    grid.arrange(p1, p2, nrow = 1, widths = c(1.5,2))
    
    cat('\n\n')
    
    cat("#### Effects\n")
    
    p1 <- emmip(model_final, formula = Group ~ Measurement) +
      theme(aspect.ratio = 1)
    
    p2 <- emmip(model_final, formula = Measurement ~ Group) +
      theme(aspect.ratio = 1)
    
    grid.arrange(p1, p2, nrow = 1)
    
    cat('\n\n')
    
    cat("#### Post-hoc tests\n")
    
    
    grp_eff <- pairs(emmeans(model_final, ~ Group | Measurement))
    msrmnt_eff <- pairs(emmeans(model_final, ~ Measurement | Group))
    df1 <- rbind(grp_eff, msrmnt_eff, adjust = "none") %>% 
      data.frame()
    
    df2 <- rbind(grp_eff, msrmnt_eff, adjust = "holm") %>% 
      data.frame() %>% 
      mutate(p.value.adj = p.value) %>% 
      select(-p.value)
    
    df <- df1 %>% 
      left_join(df2)
    
    df <- df %>% bind_cols(
      t_to_d(t = df$t.ratio, df_error = df$df) %>%
        data.frame() %>%
        round(2) %>%
        dplyr::select(d)) %>%
      mutate(eff_size = interpret_d(d)) %>%
      dplyr::select(-c(df))
    
    df %>%
      kable(escape = F, digits = 4, caption = "Multiple Comparisons") %>%
      kable_classic_2(full_width = F) %>%
      footnote(general = "P value adjustment: holm method for 9 tests") %>%
      kable_styling() %>%
      row_spec(
        which(df$p.value.adj < 0.05), bold = T,
        background = "#d6f5d6") %>% 
      print()
        
    cat('\n\n')
  }
  cat('\n\n')
}

5 Hz Frequency

RMS

Data overview

Sample of the data
Group Animal Measurement TrialNumber RMS logRMS
KIC KIC14 preLaser 27 0.353820 -1.038968
KIT KIT11 preLaser 7 0.235442 -1.446292
KIT KIT03 postLaser 3 0.128098 -2.054957
KIV KIV02 preLaser 12 0.802764 -0.219695
KIV KIV16 postLaser 36 0.466241 -0.763053

LMM

. Linear mixed-effects model fit by REML
. Data: model_data
. AIC BIC logLik
. -1276.671 -1226.5 647.3354
.
. Random effects:
. Formula: ~1 | Animal
. (Intercept)
. StdDev: 0.09003444
.
. Formula: ~1 | TrialNumber %in% Animal
. (Intercept) Residual
. StdDev: 1.325261e-05 0.1688084
.
. Fixed effects: as.formula(paste0(metric, "~ Group*Measurement"))
. Value Std.Error DF t-value p-value
. (Intercept) 0.3739030 0.02960485 1201 12.629788 0.0000
. GroupKIC -0.0335800 0.04183226 24 -0.802730 0.4300
. GroupKIV -0.0473761 0.04612000 24 -1.027236 0.3146
. MeasurementpreLaser -0.0617815 0.01304107 723 -4.737457 0.0000
. GroupKIC:MeasurementpreLaser 0.0032164 0.01834245 723 0.175353 0.8609
. GroupKIV:MeasurementpreLaser 0.0890216 0.02011984 723 4.424571 0.0000
. Correlation:
. (Intr) GrpKIC GrpKIV MsrmnL GKIC:M
. GroupKIC -0.708
. GroupKIV -0.642 0.454
. MeasurementpreLaser -0.170 0.120 0.109
. GroupKIC:MeasurementpreLaser 0.121 -0.167 -0.078 -0.711
. GroupKIV:MeasurementpreLaser 0.110 -0.078 -0.171 -0.648 0.461
.
. Standardized Within-Group Residuals:
. Min Q1 Med Q3 Max
. -2.62694992 -0.71074448 -0.04270875 0.64632602 3.59133767
.
. Number of Observations: 1954
. Number of Groups:
. Animal TrialNumber %in% Animal
. 27 1228

Residuals check

Effects

Post-hoc tests

Multiple Comparisons
Measurement Group contrast estimate SE t.ratio p.value p.value.adj d eff_size
postLaser . KIT - KIC 0.0336 0.0418 0.8027 0.4300 1.0000 0.33 small
postLaser . KIT - KIV 0.0474 0.0461 1.0272 0.3146 1.0000 0.42 small
postLaser . KIC - KIV 0.0138 0.0461 0.2993 0.7673 1.0000 0.12 very small
preLaser . KIT - KIC 0.0304 0.0428 0.7098 0.4846 1.0000 0.29 small
preLaser . KIT - KIV -0.0416 0.0471 -0.8848 0.3850 1.0000 -0.36 small
preLaser . KIC - KIV -0.0720 0.0471 -1.5303 0.1390 0.8341 -0.62 medium
. KIT postLaser - preLaser 0.0618 0.0130 4.7375 0.0000 0.0000 0.35 small
. KIC postLaser - preLaser 0.0586 0.0129 4.5404 0.0000 0.0001 0.34 small
. KIV postLaser - preLaser -0.0272 0.0153 -1.7779 0.0758 0.5308 -0.13 very small
Note:
P value adjustment: holm method for 9 tests

logRMS

Data overview

Sample of the data
Group Animal Measurement TrialNumber RMS logRMS AnimalTrial
KIC KIC04 postLaser 4 0.357351 -1.029036 KIC04_4
KIC KIC14 preLaser 21 0.325659 -1.121903 KIC14_21
KIC KIC04 postLaser 21 0.574317 -0.554573 KIC04_21
KIT KIT03 postLaser 23 0.085666 -2.457301 KIT03_23
KIV KIV16 postLaser 14 0.304169 -1.190172 KIV16_14

LMM

. Linear mixed-effects model fit by REML
. Data: model_data
. AIC BIC logLik
. 4116.343 4166.514 -2049.171
.
. Random effects:
. Formula: ~1 | Animal
. (Intercept)
. StdDev: 0.2844165
.
. Formula: ~1 | TrialNumber %in% Animal
. (Intercept) Residual
. StdDev: 6.688009e-05 0.6757027
.
. Fixed effects: as.formula(paste0(metric, "~ Group*Measurement"))
. Value Std.Error DF t-value p-value
. (Intercept) -1.1981302 0.09562255 1201 -12.529786 0.0000
. GroupKIC -0.0594068 0.13505613 24 -0.439867 0.6640
. GroupKIV -0.1292120 0.14893721 24 -0.867560 0.3942
. MeasurementpreLaser -0.1788854 0.05219831 723 -3.427034 0.0006
. GroupKIC:MeasurementpreLaser -0.0397170 0.07341882 723 -0.540965 0.5887
. GroupKIV:MeasurementpreLaser 0.2421202 0.08052760 723 3.006673 0.0027
. Correlation:
. (Intr) GrpKIC GrpKIV MsrmnL GKIC:M
. GroupKIC -0.708
. GroupKIV -0.642 0.455
. MeasurementpreLaser -0.210 0.149 0.135
. GroupKIC:MeasurementpreLaser 0.150 -0.207 -0.096 -0.711
. GroupKIV:MeasurementpreLaser 0.136 -0.097 -0.211 -0.648 0.461
.
. Standardized Within-Group Residuals:
. Min Q1 Med Q3 Max
. -6.2939578 -0.4763477 0.2100404 0.6855370 1.9957147
.
. Number of Observations: 1954
. Number of Groups:
. Animal TrialNumber %in% Animal
. 27 1228

Residuals check

Effects

Post-hoc tests

Multiple Comparisons
Measurement Group contrast estimate SE t.ratio p.value p.value.adj d eff_size
postLaser . KIT - KIC 0.0594 0.1351 0.4399 0.6640 1.0000 0.18 very small
postLaser . KIT - KIV 0.1292 0.1489 0.8676 0.3942 1.0000 0.35 small
postLaser . KIC - KIV 0.0698 0.1488 0.4692 0.6432 1.0000 0.19 very small
preLaser . KIT - KIC 0.0991 0.1397 0.7095 0.4849 1.0000 0.29 small
preLaser . KIT - KIV -0.1129 0.1536 -0.7350 0.4694 1.0000 -0.30 small
preLaser . KIC - KIV -0.2120 0.1536 -1.3808 0.1801 1.0000 -0.56 medium
. KIT postLaser - preLaser 0.1789 0.0522 3.4270 0.0006 0.0052 0.25 small
. KIC postLaser - preLaser 0.2186 0.0516 4.2340 0.0000 0.0002 0.31 small
. KIV postLaser - preLaser -0.0632 0.0613 -1.0312 0.3028 1.0000 -0.08 very small
Note:
P value adjustment: holm method for 9 tests

10 Hz Frequency

RMS

Data overview

Sample of the data
Group Animal Measurement TrialNumber RMS logRMS
KIC KIC11 preLaser 8 0.002103 -6.164472
KIC KIC16 postLaser 42 0.282593 -1.263749
KIC KIC11 preLaser 29 0.202880 -1.595139
KIT KIT09 postLaser 19 0.007780 -4.856205
KIV KIV03 postLaser 22 0.197421 -1.622418

LMM

. Linear mixed-effects model fit by REML
. Data: model_data
. AIC BIC logLik
. -768.7005 -721.0145 393.3503
.
. Random effects:
. Formula: ~1 | Animal
. (Intercept)
. StdDev: 0.0770658
.
. Formula: ~1 | TrialNumber %in% Animal
. (Intercept) Residual
. StdDev: 6.653761e-06 0.1799049
.
. Fixed effects: as.formula(paste0(metric, "~ Group*Measurement"))
. Value Std.Error DF t-value p-value
. (Intercept) 0.3444814 0.02681245 966 12.847815 0.0000
. GroupKIC -0.0763192 0.03767623 24 -2.025658 0.0541
. GroupKIV -0.0359151 0.04104865 24 -0.874939 0.3903
. MeasurementpreLaser -0.0684064 0.01657861 488 -4.126182 0.0000
. GroupKIC:MeasurementpreLaser 0.0353278 0.02306283 488 1.531807 0.1262
. GroupKIV:MeasurementpreLaser 0.0576645 0.02385459 488 2.417334 0.0160
. Correlation:
. (Intr) GrpKIC GrpKIV MsrmnL GKIC:M
. GroupKIC -0.712
. GroupKIV -0.653 0.465
. MeasurementpreLaser -0.256 0.182 0.167
. GroupKIC:MeasurementpreLaser 0.184 -0.249 -0.120 -0.719
. GroupKIV:MeasurementpreLaser 0.178 -0.127 -0.236 -0.695 0.500
.
. Standardized Within-Group Residuals:
. Min Q1 Med Q3 Max
. -2.7098611 -0.7248810 -0.1225964 0.6500029 3.4105718
.
. Number of Observations: 1484
. Number of Groups:
. Animal TrialNumber %in% Animal
. 27 993

Residuals check

Effects

Post-hoc tests

Multiple Comparisons
Measurement Group contrast estimate SE t.ratio p.value p.value.adj d eff_size
postLaser . KIT - KIC 0.0763 0.0377 2.0257 0.0541 0.3784 0.83 large
postLaser . KIT - KIV 0.0359 0.0410 0.8749 0.3903 1.0000 0.36 small
postLaser . KIC - KIV -0.0404 0.0408 -0.9897 0.3322 1.0000 -0.40 small
preLaser . KIT - KIC 0.0410 0.0390 1.0517 0.3034 1.0000 0.43 small
preLaser . KIT - KIV -0.0217 0.0423 -0.5139 0.6120 1.0000 -0.21 small
preLaser . KIC - KIV -0.0627 0.0422 -1.4880 0.1498 0.8986 -0.61 medium
. KIT postLaser - preLaser 0.0684 0.0166 4.1262 0.0000 0.0004 0.37 small
. KIC postLaser - preLaser 0.0331 0.0160 2.0632 0.0396 0.3170 0.19 very small
. KIV postLaser - preLaser 0.0107 0.0172 0.6263 0.5314 1.0000 0.06 very small
Note:
P value adjustment: holm method for 9 tests

logRMS

Data overview

Sample of the data
Group Animal Measurement TrialNumber RMS logRMS AnimalTrial
KIC KIC11 preLaser 8 0.002103 -6.164472 KIC11_8
KIC KIC16 postLaser 42 0.282593 -1.263749 KIC16_42
KIC KIC11 preLaser 29 0.202880 -1.595139 KIC11_29
KIT KIT09 postLaser 19 0.007780 -4.856205 KIT09_19
KIV KIV03 postLaser 22 0.197421 -1.622418 KIV03_22

LMM

. Linear mixed-effects model fit by REML
. Data: model_data
. AIC BIC logLik
. 3896.056 3943.743 -1939.028
.
. Random effects:
. Formula: ~1 | Animal
. (Intercept)
. StdDev: 0.2580884
.
. Formula: ~1 | TrialNumber %in% Animal
. (Intercept) Residual
. StdDev: 0.0001546594 0.8763101
.
. Fixed effects: as.formula(paste0(metric, "~ Group*Measurement"))
. Value Std.Error DF t-value p-value
. (Intercept) -1.3821308 0.09787888 966 -14.120828 0.0000
. GroupKIC -0.2467508 0.13692513 24 -1.802085 0.0841
. GroupKIV -0.0801055 0.14793448 24 -0.541493 0.5932
. MeasurementpreLaser -0.1951635 0.08072828 488 -2.417536 0.0160
. GroupKIC:MeasurementpreLaser 0.0750735 0.11229999 488 0.668509 0.5041
. GroupKIV:MeasurementpreLaser 0.1422671 0.11616079 488 1.224742 0.2213
. Correlation:
. (Intr) GrpKIC GrpKIV MsrmnL GKIC:M
. GroupKIC -0.715
. GroupKIV -0.662 0.473
. MeasurementpreLaser -0.341 0.244 0.225
. GroupKIC:MeasurementpreLaser 0.245 -0.333 -0.162 -0.719
. GroupKIV:MeasurementpreLaser 0.237 -0.169 -0.319 -0.695 0.500
.
. Standardized Within-Group Residuals:
. Min Q1 Med Q3 Max
. -6.7521150 -0.4950631 0.2075273 0.6957252 1.9374197
.
. Number of Observations: 1484
. Number of Groups:
. Animal TrialNumber %in% Animal
. 27 993

Residuals check

Effects

Post-hoc tests

Multiple Comparisons
Measurement Group contrast estimate SE t.ratio p.value p.value.adj d eff_size
postLaser . KIT - KIC 0.2468 0.1369 1.8021 0.0841 0.6729 0.74 medium
postLaser . KIT - KIV 0.0801 0.1479 0.5415 0.5932 1.0000 0.22 small
postLaser . KIC - KIV -0.1666 0.1465 -1.1372 0.2667 1.0000 -0.46 small
preLaser . KIT - KIC 0.1717 0.1454 1.1810 0.2492 1.0000 0.48 small
preLaser . KIT - KIV -0.0622 0.1562 -0.3978 0.6943 1.0000 -0.16 very small
preLaser . KIC - KIV -0.2338 0.1553 -1.5057 0.1452 0.8723 -0.61 medium
. KIT postLaser - preLaser 0.1952 0.0807 2.4175 0.0160 0.1439 0.22 small
. KIC postLaser - preLaser 0.1201 0.0781 1.5383 0.1246 0.8723 0.14 very small
. KIV postLaser - preLaser 0.0529 0.0835 0.6333 0.5268 1.0000 0.06 very small
Note:
P value adjustment: holm method for 9 tests