library(tidyverse)
library(ggpubr)
library(rstatix)
library(knitr)
library(kableExtra)
library(nlme)
library(emmeans)
library(standardize)
library(sjPlot)
library(gridExtra)
library(effectsize)
data_type <- "CLST"
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_CL_10.csv") %>%
  mutate(ClickFreq = 10) %>%
  bind_rows(
    read_csv("AnimalTrial_CL_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
KIT KIT01 postLaser 22 0.405662 -0.902235
KIC KIC03 postLaser 39 0.263341 -1.334305
KIC KIC11 postLaser 3 0.233298 -1.455437
KIC KIC12 postLaser 14 0.220786 -1.510559
KIV KIV12 postLaser 24 0.470421 -0.754128

LMM

. Linear mixed-effects model fit by REML
. Data: model_data
. AIC BIC logLik
. -1282.051 -1231.657 650.0257
.
. Random effects:
. Formula: ~1 | Animal
. (Intercept)
. StdDev: 0.09573827
.
. Formula: ~1 | TrialNumber %in% Animal
. (Intercept) Residual
. StdDev: 4.087625e-06 0.1699367
.
. Fixed effects: as.formula(paste0(metric, "~ Group*Measurement"))
. Value Std.Error DF t-value p-value
. (Intercept) 0.3883310 0.03131665 1232 12.400145 0.0000
. GroupKIC -0.0385617 0.04428551 24 -0.870751 0.3925
. GroupKIV -0.0347429 0.04875097 24 -0.712660 0.4829
. MeasurementpreLaser -0.0188273 0.01300505 741 -1.447693 0.1481
. GroupKIC:MeasurementpreLaser 0.0325509 0.01825582 741 1.783042 0.0750
. GroupKIV:MeasurementpreLaser 0.0464327 0.01998205 741 2.323720 0.0204
. Correlation:
. (Intr) GrpKIC GrpKIV MsrmnL GKIC:M
. GroupKIC -0.707
. GroupKIV -0.642 0.454
. MeasurementpreLaser -0.157 0.111 0.101
. GroupKIC:MeasurementpreLaser 0.112 -0.158 -0.072 -0.712
. GroupKIV:MeasurementpreLaser 0.102 -0.072 -0.155 -0.651 0.464
.
. Standardized Within-Group Residuals:
. Min Q1 Med Q3 Max
. -3.00284346 -0.69779741 -0.02727847 0.70848387 2.99575165
.
. Number of Observations: 2003
. Number of Groups:
. Animal TrialNumber %in% Animal
. 27 1259

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.0386 0.0443 0.8708 0.3925 1.0000 0.36 small
postLaser . KIT - KIV 0.0347 0.0488 0.7127 0.4829 1.0000 0.29 small
postLaser . KIC - KIV -0.0038 0.0487 -0.0783 0.9382 1.0000 -0.03 very small
preLaser . KIT - KIC 0.0060 0.0452 0.1331 0.8952 1.0000 0.05 very small
preLaser . KIT - KIV -0.0117 0.0497 -0.2350 0.8162 1.0000 -0.10 very small
preLaser . KIC - KIV -0.0177 0.0497 -0.3562 0.7248 1.0000 -0.15 very small
. KIT postLaser - preLaser 0.0188 0.0130 1.4477 0.1481 1.0000 0.11 very small
. KIC postLaser - preLaser -0.0137 0.0128 -1.0712 0.2844 1.0000 -0.08 very small
. KIV postLaser - preLaser -0.0276 0.0152 -1.8196 0.0692 0.6229 -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.205316 -1.583206 KIC04_4
KIC KIC14 preLaser 10 0.140318 -1.963842 KIC14_10
KIC KIC04 postLaser 20 0.330247 -1.107914 KIC04_20
KIT KIT03 postLaser 13 0.201392 -1.602504 KIT03_13
KIV KIV12 postLaser 42 0.180145 -1.713994 KIV12_42

LMM

. Linear mixed-effects model fit by REML
. Data: model_data
. AIC BIC logLik
. 3804.091 3854.485 -1893.045
.
. Random effects:
. Formula: ~1 | Animal
. (Intercept)
. StdDev: 0.2875882
.
. Formula: ~1 | TrialNumber %in% Animal
. (Intercept) Residual
. StdDev: 6.617417e-05 0.6084282
.
. Fixed effects: as.formula(paste0(metric, "~ Group*Measurement"))
. Value Std.Error DF t-value p-value
. (Intercept) -1.1080811 0.09535708 1232 -11.620333 0.0000
. GroupKIC -0.1183684 0.13484309 24 -0.877823 0.3887
. GroupKIV -0.1638829 0.14838244 24 -1.104463 0.2803
. MeasurementpreLaser -0.0624081 0.04656171 741 -1.340332 0.1805
. GroupKIC:MeasurementpreLaser 0.1434290 0.06536112 741 2.194408 0.0285
. GroupKIV:MeasurementpreLaser 0.1715574 0.07154180 741 2.398002 0.0167
. Correlation:
. (Intr) GrpKIC GrpKIV MsrmnL GKIC:M
. GroupKIC -0.707
. GroupKIV -0.643 0.454
. MeasurementpreLaser -0.185 0.131 0.119
. GroupKIC:MeasurementpreLaser 0.132 -0.186 -0.085 -0.712
. GroupKIV:MeasurementpreLaser 0.120 -0.085 -0.182 -0.651 0.464
.
. Standardized Within-Group Residuals:
. Min Q1 Med Q3 Max
. -6.5421985 -0.4721583 0.2071221 0.7170518 2.1572014
.
. Number of Observations: 2003
. Number of Groups:
. Animal TrialNumber %in% Animal
. 27 1259

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.1184 0.1348 0.8778 0.3887 1.0000 0.36 small
postLaser . KIT - KIV 0.1639 0.1484 1.1045 0.2803 1.0000 0.45 small
postLaser . KIC - KIV 0.0455 0.1484 0.3068 0.7617 1.0000 0.13 very small
preLaser . KIT - KIC -0.0251 0.1385 -0.1810 0.8579 1.0000 -0.07 very small
preLaser . KIT - KIV -0.0077 0.1526 -0.0503 0.9603 1.0000 -0.02 very small
preLaser . KIC - KIV 0.0174 0.1524 0.1141 0.9101 1.0000 0.05 very small
. KIT postLaser - preLaser 0.0624 0.0466 1.3403 0.1805 1.0000 0.10 very small
. KIC postLaser - preLaser -0.0810 0.0459 -1.7663 0.0778 0.6221 -0.13 very small
. KIV postLaser - preLaser -0.1091 0.0543 -2.0095 0.0448 0.4036 -0.15 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 KIC04 postLaser 35 0.218636 -1.520345
KIC KIC15 postLaser 7 0.230981 -1.465421
KIC KIC10 preLaser 5 0.432388 -0.838433
KIT KIT05 postLaser 23 0.597950 -0.514247
KIT KIT12 preLaser 14 0.077900 -2.552333

LMM

. Linear mixed-effects model fit by REML
. Data: model_data
. AIC BIC logLik
. -992.1759 -943.4502 505.088
.
. Random effects:
. Formula: ~1 | Animal
. (Intercept)
. StdDev: 0.09098412
.
. Formula: ~1 | TrialNumber %in% Animal
. (Intercept) Residual
. StdDev: 1.574144e-05 0.1730518
.
. Fixed effects: as.formula(paste0(metric, "~ Group*Measurement"))
. Value Std.Error DF t-value p-value
. (Intercept) 0.3166270 0.03040683 1083 10.413022 0.0000
. GroupKIC -0.0003945 0.04280569 24 -0.009216 0.9927
. GroupKIV -0.0113823 0.04698500 24 -0.242255 0.8106
. MeasurementpreLaser -0.0280404 0.01508349 552 -1.859013 0.0636
. GroupKIC:MeasurementpreLaser -0.0074434 0.02067698 552 -0.359983 0.7190
. GroupKIV:MeasurementpreLaser 0.0378916 0.02226736 552 1.701668 0.0894
. Correlation:
. (Intr) GrpKIC GrpKIV MsrmnL GKIC:M
. GroupKIC -0.710
. GroupKIV -0.647 0.460
. MeasurementpreLaser -0.192 0.137 0.124
. GroupKIC:MeasurementpreLaser 0.140 -0.188 -0.091 -0.729
. GroupKIV:MeasurementpreLaser 0.130 -0.093 -0.180 -0.677 0.494
.
. Standardized Within-Group Residuals:
. Min Q1 Med Q3 Max
. -2.7479994 -0.6970997 -0.0600754 0.6278760 3.0529226
.
. Number of Observations: 1665
. Number of Groups:
. Animal TrialNumber %in% Animal
. 27 1110

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.0004 0.0428 0.0092 0.9927 1.0000 0.00 very small
postLaser . KIT - KIV 0.0114 0.0470 0.2423 0.8106 1.0000 0.10 very small
postLaser . KIC - KIV 0.0110 0.0468 0.2348 0.8164 1.0000 0.10 very small
preLaser . KIT - KIC 0.0078 0.0439 0.1786 0.8598 1.0000 0.07 very small
preLaser . KIT - KIV -0.0265 0.0482 -0.5495 0.5877 1.0000 -0.22 small
preLaser . KIC - KIV -0.0343 0.0480 -0.7158 0.4810 1.0000 -0.29 small
. KIT postLaser - preLaser 0.0280 0.0151 1.8590 0.0636 0.5085 0.16 very small
. KIC postLaser - preLaser 0.0355 0.0141 2.5089 0.0124 0.1116 0.21 small
. KIV postLaser - preLaser -0.0099 0.0164 -0.6014 0.5478 1.0000 -0.05 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 35 0.218636 -1.520345 KIC04_35
KIC KIC15 postLaser 7 0.230981 -1.465421 KIC15_7
KIC KIC10 preLaser 5 0.432388 -0.838433 KIC10_5
KIT KIT05 postLaser 23 0.597950 -0.514247 KIT05_23
KIT KIT12 preLaser 14 0.077900 -2.552333 KIT12_14

LMM

. Linear mixed-effects model fit by REML
. Data: model_data
. AIC BIC logLik
. 4205.466 4254.192 -2093.733
.
. Random effects:
. Formula: ~1 | Animal
. (Intercept)
. StdDev: 0.2875863
.
. Formula: ~1 | TrialNumber %in% Animal
. (Intercept) Residual
. StdDev: 0.0001571872 0.8334988
.
. Fixed effects: as.formula(paste0(metric, "~ Group*Measurement"))
. Value Std.Error DF t-value p-value
. (Intercept) -1.4013678 0.10242139 1083 -13.682374 0.0000
. GroupKIC -0.0343380 0.14356827 24 -0.239175 0.8130
. GroupKIV -0.1057731 0.15695262 24 -0.673917 0.5068
. MeasurementpreLaser -0.1543809 0.07261651 552 -2.125975 0.0339
. GroupKIC:MeasurementpreLaser 0.0552460 0.09956311 552 0.554884 0.5792
. GroupKIV:MeasurementpreLaser 0.2457530 0.10722609 552 2.291914 0.0223
. Correlation:
. (Intr) GrpKIC GrpKIV MsrmnL GKIC:M
. GroupKIC -0.713
. GroupKIV -0.653 0.466
. MeasurementpreLaser -0.276 0.197 0.180
. GroupKIC:MeasurementpreLaser 0.201 -0.271 -0.131 -0.729
. GroupKIV:MeasurementpreLaser 0.187 -0.133 -0.259 -0.677 0.494
.
. Standardized Within-Group Residuals:
. Min Q1 Med Q3 Max
. -6.5963456 -0.4255285 0.2192718 0.6704474 1.7367839
.
. Number of Observations: 1665
. Number of Groups:
. Animal TrialNumber %in% Animal
. 27 1110

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.0343 0.1436 0.2392 0.8130 1.0000 0.10 very small
postLaser . KIT - KIV 0.1058 0.1570 0.6739 0.5068 1.0000 0.28 small
postLaser . KIC - KIV 0.0714 0.1558 0.4586 0.6507 1.0000 0.19 very small
preLaser . KIT - KIC -0.0209 0.1509 -0.1385 0.8910 1.0000 -0.06 very small
preLaser . KIT - KIV -0.1400 0.1655 -0.8456 0.4061 1.0000 -0.35 small
preLaser . KIC - KIV -0.1191 0.1639 -0.7266 0.4745 1.0000 -0.30 small
. KIT postLaser - preLaser 0.1544 0.0726 2.1260 0.0339 0.3055 0.18 very small
. KIC postLaser - preLaser 0.0991 0.0681 1.4554 0.1461 1.0000 0.12 very small
. KIV postLaser - preLaser -0.0914 0.0789 -1.1582 0.2473 1.0000 -0.10 very small
Note:
P value adjustment: holm method for 9 tests