Imports

library(sjPlot)
library(effects)
library(car)
library(nlme)
library(standardize)
library(tidyverse)
library(gridExtra)
library(ggpubr)
library(knitr)
library(car)
library(kableExtra)
library(lattice)
library(multcomp)

glht.table <- function(x) {

# https://gist.github.com/cheuerde/3acc1879dc397a1adfb0
  pq <- summary(x)$test
  mtests <- cbind(pq$coefficients, pq$sigma, pq$tstat, pq$pvalues)
  error <- attr(pq$pvalues, "error")
  pname <- switch(x$alternativ, less = paste("Pr(<", ifelse(x$df ==0, "z", "t"), ")", sep = ""), 
  greater = paste("Pr(>", ifelse(x$df == 0, "z", "t"), ")", sep = ""), two.sided = paste("Pr(>|",ifelse(x$df == 0, "z", "t"), "|)", sep = ""))
  colnames(mtests) <- c("Estimate", "Std. Error", ifelse(x$df ==0, "z value", "t value"), pname)
  return(mtests)
}

AVRECPeakCL <- read_csv(file = "../../Data_new/AVRECPeakCLST.csv")

Data Preparation

# clean data for modeling
model_df <- AVRECPeakCL %>%
  filter(
    case_when(
      ClickFreq == 2 ~ OrderofClick %in% c(1,2),
      ClickFreq == 5 ~ OrderofClick %in% c(1,3),
      ClickFreq == 10 ~ OrderofClick %in% c(1,5),
      ClickFreq == 20 ~ OrderofClick %in% c(1,10),
      ClickFreq == 40 ~ OrderofClick %in% c(1,20)),
    # Layer != "All",
    Measurement %in% c("preCL_1", "CL_1")) %>%  
  mutate(
    logRMS = log(RMS),
    Measurement = case_when(
      Measurement == "preCL_1" ~ "preLaser",
      TRUE ~ "postLaser"),
    ClickOrder = case_when(
      OrderofClick == 1 ~ "1",
      TRUE ~ "2")
  ) %>% 
  ungroup() %>%   
  drop_na(c(RMS, PeakLat)) %>% 
  dplyr::select(-c(OrderofClick, PeakLat, PeakAmp))

  # filter(logRMS > quantile(logRMS, 0.005))
  # filter(between(logRMS, quantile(logRMS, 0.005), quantile(logRMS, 0.995)))
  

# convert to factors
model_df$Measurement <- factor(
  x = model_df$Measurement, 
  levels = c("preLaser", "postLaser"))
model_df$ClickOrder <- factor(x = model_df$ClickOrder)
model_df$Layer <- factor(
  x = model_df$Layer,
  levels = c("All", "I_II", "IV", "V", "VI"))
model_df$Group <- factor(x = model_df$Group)


# sum-coding
# model_df$Measurement = named_contr_sum(model_df$Measurement, return_contr = FALSE)
# model_df$Group = named_contr_sum(model_df$Group, return_contr = FALSE)
# model_df$Layer = named_contr_sum(model_df$Layer, return_contr = FALSE)

p1 <- ggplot(data = model_df) +
  geom_histogram(
    aes(x = RMS), bins = 50, 
    fill = "lightblue", color = "black") +
  labs(title = "Histogram of RMS")

p2 <- model_df %>% 
  ggplot() +
  geom_histogram(
    aes(x = logRMS), bins = 50, 
    fill = "lightblue", color = "black") +
  labs(title = "Histogram of log(RMS)")

grid.arrange(p1, p2, nrow = 1)

model_df  %>% 
    filter(ClickFreq == 5) %>%  
ggboxplot(
    x = "Group", y = "logRMS", color = "Measurement", 
    palette = "jco", facet.by = c("ClickOrder", "Layer"), main = "5 Hz" #, add = "jitter",
)

head(model_df) %>% 
  kable(escape = F, caption = "Data overview") %>%
  kable_classic_2(full_width = F)

Data overview
Group	Animal	Layer	Measurement	ClickFreq	TrialNumber	RMS	logRMS	ClickOrder
KIC	KIC02	All	preLaser	2	1	0.0010370	-6.871387	1
KIC	KIC02	All	preLaser	2	1	0.0011177	-6.796501	2
KIC	KIC02	All	preLaser	2	2	0.0012743	-6.665325	1
KIC	KIC02	All	preLaser	2	2	0.0014588	-6.530117	2
KIC	KIC02	All	preLaser	2	3	0.0017276	-6.361004	1
KIC	KIC02	All	preLaser	2	3	0.0009644	-6.943973	2

Modelling

Comparing different models

model1 <- lme(
    logRMS ~ Group, 
    random = (~1|Animal), data = model_df,
    method = "ML",
    control = lmeControl(opt = "optim"))
model2 <- lme(
    logRMS ~ Group, 
    random = (~1|Animal/TrialNumber), data = model_df,
    method = "ML",
    control = lmeControl(opt = "optim"))
model3 <- lme(
    logRMS ~ Group + Measurement, 
    random = (~1|Animal/TrialNumber), data = model_df,
    method = "ML",
    control = lmeControl(opt = "optim"))
model4 <- lme(
    logRMS ~ Group*Measurement, 
    random = (~1|Animal/TrialNumber), data = model_df,
    method = "ML",
    control = lmeControl(opt = "optim"))
model5 <- lme(
    logRMS ~ Group*Measurement + Layer, 
    random = (~1|Animal/TrialNumber), data = model_df,
    method = "ML",
    control = lmeControl(opt = "optim"))
model6 <- lme(
    logRMS ~ Group*Measurement + Layer + ClickFreq, 
    random = (~1|Animal/TrialNumber), data = model_df,
    method = "ML",
    control = lmeControl(opt = "optim"))
model7 <- lme(
    logRMS ~ Group*Measurement*Layer + ClickFreq, 
    random = (~1|Animal/TrialNumber), data = model_df,
    method = "ML",
    control = lmeControl(opt = "optim"))
model8 <- lme(
    logRMS ~ Group*Measurement*Layer + ClickFreq + ClickOrder, 
    random = (~1|Animal/TrialNumber), data = model_df,
    method = "ML",
    control = lmeControl(opt = "optim"))

anova(
  model1, model2, model3, model4, 
  model5, model6, model7, model8) %>% 
  dplyr::select(-c(Model, df, L.Ratio)) %>% 
  kable(escape = F, digits = 3, caption = "Models Comparison") %>%
  kable_classic_2()

Models Comparison
	call	AIC	BIC	logLik	Test	p-value
model1	lme.formula(fixed = logRMS ~ Group, data = model_df, random = (~1 \| Animal), method = “ML”, control = lmeControl(opt = “optim”))	165086.3	165133.2	-82538.17		NA
model2	lme.formula(fixed = logRMS ~ Group, data = model_df, random = (~1 \| Animal/TrialNumber), method = “ML”, control = lmeControl(opt = “optim”))	164666.3	164722.5	-82327.13	1 vs 2	0
model3	lme.formula(fixed = logRMS ~ Group + Measurement, data = model_df, random = (~1 \| Animal/TrialNumber), method = “ML”, control = lmeControl(opt = “optim”))	164639.5	164705.1	-82312.76	2 vs 3	0
model4	lme.formula(fixed = logRMS ~ Group * Measurement, data = model_df, random = (~1 \| Animal/TrialNumber), method = “ML”, control = lmeControl(opt = “optim”))	164307.2	164391.6	-82144.60	3 vs 4	0
model5	lme.formula(fixed = logRMS ~ Group * Measurement + Layer, data = model_df, random = (~1 \| Animal/TrialNumber), method = “ML”, control = lmeControl(opt = “optim”))	156407.8	156529.7	-78190.90	4 vs 5	0
model6	lme.formula(fixed = logRMS ~ Group * Measurement + Layer + ClickFreq, data = model_df, random = (~1 \| Animal/TrialNumber), method = “ML”, control = lmeControl(opt = “optim”))	155485.6	155616.8	-77728.79	5 vs 6	0
model7	lme.formula(fixed = logRMS ~ Group * Measurement * Layer + ClickFreq, data = model_df, random = (~1 \| Animal/TrialNumber), method = “ML”, control = lmeControl(opt = “optim”))	154731.0	155049.8	-77331.52	6 vs 7	0
model8	lme.formula(fixed = logRMS ~ Group * Measurement * Layer + ClickFreq + ClickOrder, data = model_df, random = (~1 \| Animal/TrialNumber), method = “ML”, control = lmeControl(opt = “optim”))	152074.6	152402.7	-76002.29	7 vs 8	0

Final model

model_full <- lme(
    logRMS ~ Group*Measurement*Layer + ClickFreq + ClickOrder, 
    random = (~1|Animal/TrialNumber), data = model_df,
    control = lmeControl(opt = "optim"))

tab_model(model_full, show.stat = TRUE, show.aic = TRUE, title = "Summary of the LMM")

Summary of the LMM
	logRMS
Predictors	Estimates	CI	Statistic	p
(Intercept)	-6.85	-7.06 – -6.63	-62.01	<0.001
Group [KIT]	-0.15	-0.47 – 0.17	-0.95	0.350
Group [KIV]	-0.03	-0.39 – 0.32	-0.20	0.844
Measurement [postLaser]	0.11	0.08 – 0.14	6.42	<0.001
Layer [I_II]	-0.21	-0.25 – -0.18	-12.20	<0.001
Layer [IV]	-0.01	-0.05 – 0.02	-0.72	0.473
Layer [V]	-0.43	-0.46 – -0.39	-23.96	<0.001
Layer [VI]	-0.54	-0.58 – -0.51	-31.18	<0.001
ClickFreq	-0.01	-0.01 – -0.00	-31.35	<0.001
ClickOrder [2]	-0.20	-0.21 – -0.19	-51.95	<0.001
Group [KIT] * Measurement [postLaser]	-0.26	-0.31 – -0.21	-10.41	<0.001
Group [KIV] * Measurement [postLaser]	-0.01	-0.06 – 0.04	-0.47	0.638
Group [KIT] * Layer [I_II]	-0.05	-0.10 – -0.00	-2.12	0.034
Group [KIV] * Layer [I_II]	0.16	0.10 – 0.21	5.89	<0.001
Group [KIT] * Layer [IV]	-0.00	-0.05 – 0.05	-0.06	0.955
Group [KIV] * Layer [IV]	0.10	0.04 – 0.15	3.65	<0.001
Group [KIT] * Layer [V]	0.07	0.02 – 0.12	2.57	0.010
Group [KIV] * Layer [V]	0.10	0.04 – 0.15	3.59	<0.001
Group [KIT] * Layer [VI]	0.18	0.13 – 0.23	7.24	<0.001
Group [KIV] * Layer [VI]	0.18	0.13 – 0.23	6.74	<0.001
Measurement [postLaser] * Layer [I_II]	-0.04	-0.08 – 0.01	-1.64	0.102
Measurement [postLaser] * Layer [IV]	0.02	-0.02 – 0.07	1.10	0.272
Measurement [postLaser] * Layer [V]	-0.01	-0.05 – 0.04	-0.25	0.805
Measurement [postLaser] * Layer [VI]	-0.03	-0.08 – 0.01	-1.55	0.122
(Group [KIT] Measurement [postLaser]) Layer [I_II]	0.17	0.11 – 0.23	5.29	<0.001
(Group [KIV] Measurement [postLaser]) Layer [I_II]	-0.00	-0.07 – 0.07	-0.01	0.990
(Group [KIT] Measurement [postLaser]) Layer [IV]	-0.07	-0.13 – -0.00	-2.08	0.038
(Group [KIV] Measurement [postLaser]) Layer [IV]	-0.10	-0.16 – -0.03	-2.86	0.004
(Group [KIT] Measurement [postLaser]) Layer [V]	0.10	0.03 – 0.16	2.96	0.003
(Group [KIV] Measurement [postLaser]) Layer [V]	0.03	-0.04 – 0.09	0.84	0.399
(Group [KIT] Measurement [postLaser]) Layer [VI]	0.13	0.07 – 0.20	4.27	<0.001
(Group [KIV] Measurement [postLaser]) Layer [VI]	-0.03	-0.09 – 0.04	-0.77	0.444
Random Effects
σ²	0.33
τ₀₀ _TrialNumber	0.35
τ₀₀ _Animal	0.08
N _TrialNumber	50
N _Animal	27
Observations	87149
Marginal R² / Conditional R²	0.168 / NA
AIC	152293.921

Anova(model_full) %>% 
  kable(
    escape = F, digits = 3,
    caption = "Analysis of Deviance Table (Type II tests)") %>% 
  kable_classic_2(full_width = F)

Analysis of Deviance Table (Type II tests)
	Chisq	Df	Pr(>Chisq)
Group	3.413	2	0.182
Measurement	31.969	1	0.000
Layer	7899.725	4	0.000
ClickFreq	982.768	1	0.000
ClickOrder	2699.239	1	0.000
Group:Measurement	389.508	2	0.000
Group:Layer	686.173	8	0.000
Measurement:Layer	28.346	4	0.000
Group:Measurement:Layer	101.881	8	0.000

Post-hoc tests

glht(model_full, linfct = mcp(Group = "Tukey"), test = adjusted("holm")) %>% 
  glht.table() %>% 
  kable(escape = F, digits = 3, caption = "Groups Comparison") %>%
  kable_classic_2(full_width = F)

Groups Comparison
	Estimate	Std. Error	z value	Pr(>\|z\|)
KIT - KIC	-0.149	0.156	-0.953	0.606
KIV - KIC	-0.034	0.172	-0.199	0.978
KIV - KIT	0.115	0.172	0.666	0.783

glht(model_full, linfct = mcp(Layer = "Tukey"), test = adjusted("holm")) %>% 
  glht.table() %>% 
  kable(escape = F, digits = 3, caption = "Layers Comparison") %>%
  kable_classic_2(full_width = F)

Layers Comparison
	Estimate	Std. Error	z value	Pr(>\|z\|)
I_II - All	-0.215	0.018	-12.195	0.000
IV - All	-0.013	0.018	-0.718	0.952
V - All	-0.425	0.018	-23.958	0.000
VI - All	-0.542	0.017	-31.176	0.000
IV - I_II	0.202	0.016	12.452	0.000
V - I_II	-0.211	0.016	-12.996	0.000
VI - I_II	-0.327	0.016	-20.746	0.000
V - IV	-0.413	0.016	-25.172	0.000
VI - IV	-0.529	0.016	-33.141	0.000
VI - V	-0.116	0.016	-7.305	0.000

model_df$IntTerm <- interaction(model_df$Group, model_df$Measurement)
lme(
  logRMS ~ IntTerm*Layer + ClickFreq + ClickOrder, 
  random = (~1|Animal/TrialNumber), data = model_df,
  control = lmeControl(opt = "optim")) %>% 
  glht(linfct = mcp(IntTerm = "Tukey"), test = adjusted("holm")) %>% 
  glht.table() %>% 
  kable(escape = F, digits = 3, caption = "Group*Measurement Comparison") %>%
  kable_classic_2(full_width = F)

Group*Measurement Comparison
	Estimate	Std. Error	z value	Pr(>\|z\|)
KIT.preLaser - KIC.preLaser	-0.149	0.156	-0.953	0.895
KIV.preLaser - KIC.preLaser	-0.034	0.172	-0.199	1.000
KIC.postLaser - KIC.preLaser	0.110	0.017	6.420	0.000
KIT.postLaser - KIC.preLaser	-0.297	0.156	-1.906	0.293
KIV.postLaser - KIC.preLaser	0.064	0.172	0.373	0.998
KIV.preLaser - KIT.preLaser	0.115	0.172	0.666	0.976
KIC.postLaser - KIT.preLaser	0.259	0.156	1.661	0.445
KIT.postLaser - KIT.preLaser	-0.148	0.018	-8.270	0.000
KIV.postLaser - KIT.preLaser	0.213	0.172	1.240	0.740
KIC.postLaser - KIV.preLaser	0.145	0.172	0.842	0.936
KIT.postLaser - KIV.preLaser	-0.263	0.172	-1.531	0.537
KIV.postLaser - KIV.preLaser	0.098	0.019	5.245	0.000
KIT.postLaser - KIC.postLaser	-0.407	0.156	-2.617	0.055
KIV.postLaser - KIC.postLaser	-0.046	0.171	-0.270	1.000
KIV.postLaser - KIT.postLaser	0.361	0.171	2.107	0.196

model_df$IntTerm <- interaction(model_df$Group, model_df$Layer)
lme(
  logRMS ~ IntTerm*Measurement + ClickFreq + ClickOrder, 
  random = (~1|Animal/TrialNumber), data = model_df,
  control = lmeControl(opt = "optim")) %>% 
  glht(linfct = mcp(IntTerm = "Tukey"), test = adjusted("holm")) %>% 
  glht.table() %>% 
  kable(escape = F, digits = 3, caption = "Group*Layer Comparison") %>%
  kable_classic_2(full_width = F)

Group*Layer Comparison
	Estimate	Std. Error	z value	Pr(>\|z\|)
KIT.All - KIC.All	-0.149	0.156	-0.953	0.999
KIV.All - KIC.All	-0.034	0.172	-0.199	1.000
KIC.I_II - KIC.All	-0.215	0.018	-12.195	0.000
KIT.I_II - KIC.All	-0.418	0.156	-2.679	0.187
KIV.I_II - KIC.All	-0.094	0.172	-0.545	1.000
KIC.IV - KIC.All	-0.013	0.018	-0.718	1.000
KIT.IV - KIC.All	-0.163	0.156	-1.045	0.998
KIV.IV - KIC.All	0.050	0.172	0.291	1.000
KIC.V - KIC.All	-0.425	0.018	-23.958	0.000
KIT.V - KIC.All	-0.509	0.156	-3.263	0.034
KIV.V - KIC.All	-0.364	0.172	-2.121	0.567
KIC.VI - KIC.All	-0.542	0.017	-31.176	0.000
KIT.VI - KIC.All	-0.510	0.156	-3.271	0.033
KIV.VI - KIC.All	-0.400	0.172	-2.326	0.405
KIV.All - KIT.All	0.115	0.172	0.666	1.000
KIC.I_II - KIT.All	-0.066	0.156	-0.424	1.000
KIT.I_II - KIT.All	-0.269	0.019	-14.510	0.000
KIV.I_II - KIT.All	0.055	0.172	0.321	1.000
KIC.IV - KIT.All	0.136	0.156	0.872	1.000
KIT.IV - KIT.All	-0.014	0.018	-0.777	1.000
KIV.IV - KIT.All	0.199	0.172	1.157	0.995
KIC.V - KIT.All	-0.277	0.156	-1.773	0.823
KIT.V - KIT.All	-0.360	0.018	-19.696	0.000
KIV.V - KIT.All	-0.216	0.172	-1.254	0.988
KIC.VI - KIT.All	-0.393	0.156	-2.519	0.272
KIT.VI - KIT.All	-0.361	0.018	-20.145	0.000
KIV.VI - KIT.All	-0.251	0.172	-1.460	0.956
KIC.I_II - KIV.All	-0.181	0.172	-1.052	0.998
KIT.I_II - KIV.All	-0.384	0.172	-2.233	0.475
KIV.I_II - KIV.All	-0.059	0.020	-3.018	0.074
KIC.IV - KIV.All	0.021	0.172	0.125	1.000
KIT.IV - KIV.All	-0.129	0.172	-0.750	1.000
KIV.IV - KIV.All	0.084	0.020	4.263	0.001
KIC.V - KIV.All	-0.391	0.172	-2.278	0.442
KIT.V - KIV.All	-0.475	0.172	-2.763	0.150
KIV.V - KIV.All	-0.330	0.020	-16.677	0.000
KIC.VI - KIV.All	-0.508	0.172	-2.956	0.088
KIT.VI - KIV.All	-0.476	0.172	-2.771	0.149
KIV.VI - KIV.All	-0.365	0.020	-18.659	0.000
KIT.I_II - KIC.I_II	-0.203	0.156	-1.302	0.984
KIV.I_II - KIC.I_II	0.121	0.172	0.707	1.000
KIC.IV - KIC.I_II	0.202	0.016	12.452	0.000
KIT.IV - KIC.I_II	0.052	0.156	0.334	1.000
KIV.IV - KIC.I_II	0.265	0.172	1.544	0.931
KIC.V - KIC.I_II	-0.211	0.016	-12.996	0.000
KIT.V - KIC.I_II	-0.294	0.156	-1.886	0.748
KIV.V - KIC.I_II	-0.149	0.172	-0.870	1.000
KIC.VI - KIC.I_II	-0.327	0.016	-20.746	0.000
KIT.VI - KIC.I_II	-0.295	0.156	-1.894	0.743
KIV.VI - KIC.I_II	-0.185	0.172	-1.076	0.998
KIV.I_II - KIT.I_II	0.324	0.172	1.888	0.747
KIC.IV - KIT.I_II	0.405	0.156	2.599	0.228
KIT.IV - KIT.I_II	0.255	0.017	15.278	0.000
KIV.IV - KIT.I_II	0.468	0.172	2.725	0.167
KIC.V - KIT.I_II	-0.008	0.156	-0.049	1.000
KIT.V - KIT.I_II	-0.091	0.017	-5.454	0.000
KIV.V - KIT.I_II	0.054	0.172	0.312	1.000
KIC.VI - KIT.I_II	-0.124	0.156	-0.795	1.000
KIT.VI - KIT.I_II	-0.092	0.016	-5.669	0.000
KIV.VI - KIT.I_II	0.018	0.172	0.106	1.000
KIC.IV - KIV.I_II	0.081	0.172	0.471	1.000
KIT.IV - KIV.I_II	-0.069	0.172	-0.404	1.000
KIV.IV - KIV.I_II	0.144	0.019	7.657	0.000
KIC.V - KIV.I_II	-0.332	0.172	-1.933	0.715
KIT.V - KIV.I_II	-0.415	0.172	-2.419	0.337
KIV.V - KIV.I_II	-0.271	0.019	-14.400	0.000
KIC.VI - KIV.I_II	-0.448	0.172	-2.611	0.222
KIT.VI - KIV.I_II	-0.416	0.172	-2.426	0.332
KIV.VI - KIV.I_II	-0.306	0.019	-16.480	0.000
KIT.IV - KIC.IV	-0.150	0.156	-0.964	0.999
KIV.IV - KIC.IV	0.063	0.172	0.366	1.000
KIC.V - KIC.IV	-0.413	0.016	-25.172	0.000
KIT.V - KIC.IV	-0.496	0.156	-3.184	0.045
KIV.V - KIC.IV	-0.352	0.172	-2.048	0.625
KIC.VI - KIC.IV	-0.529	0.016	-33.141	0.000
KIT.VI - KIC.IV	-0.497	0.156	-3.192	0.043
KIV.VI - KIC.IV	-0.387	0.172	-2.254	0.459
KIV.IV - KIT.IV	0.213	0.172	1.241	0.990
KIC.V - KIT.IV	-0.262	0.156	-1.684	0.872
KIT.V - KIT.IV	-0.346	0.016	-21.012	0.000
KIV.V - KIT.IV	-0.201	0.172	-1.173	0.994
KIC.VI - KIT.IV	-0.379	0.156	-2.432	0.331
KIT.VI - KIT.IV	-0.347	0.016	-21.634	0.000
KIV.VI - KIT.IV	-0.237	0.172	-1.378	0.973
KIC.V - KIV.IV	-0.476	0.172	-2.770	0.148
KIT.V - KIV.IV	-0.559	0.172	-3.256	0.034
KIV.V - KIV.IV	-0.414	0.019	-21.938	0.000
KIC.VI - KIV.IV	-0.592	0.172	-3.449	0.018
KIT.VI - KIV.IV	-0.560	0.172	-3.263	0.034
KIV.VI - KIV.IV	-0.450	0.019	-24.098	0.000
KIT.V - KIC.V	-0.083	0.156	-0.535	1.000
KIV.V - KIC.V	0.061	0.172	0.356	1.000
KIC.VI - KIC.V	-0.116	0.016	-7.305	0.000
KIT.VI - KIC.V	-0.085	0.156	-0.543	1.000
KIV.VI - KIC.V	0.026	0.172	0.151	1.000
KIV.V - KIT.V	0.145	0.172	0.842	1.000
KIC.VI - KIT.V	-0.033	0.156	-0.211	1.000
KIT.VI - KIT.V	-0.001	0.016	-0.072	1.000
KIV.VI - KIT.V	0.109	0.172	0.637	1.000
KIC.VI - KIV.V	-0.177	0.172	-1.034	0.998
KIT.VI - KIV.V	-0.146	0.172	-0.849	1.000
KIV.VI - KIV.V	-0.035	0.019	-1.888	0.747
KIT.VI - KIC.VI	0.032	0.156	0.204	1.000
KIV.VI - KIC.VI	0.142	0.172	0.829	1.000
KIV.VI - KIT.VI	0.110	0.172	0.643	1.000

model_df$IntTerm <- interaction(model_df$Measurement, model_df$Layer)
lme(
  logRMS ~ IntTerm*Group + ClickFreq + ClickOrder, 
  random = (~1|Animal/TrialNumber), data = model_df,
  control = lmeControl(opt = "optim")) %>% 
  glht(linfct = mcp(IntTerm = "Tukey"), test = adjusted("holm")) %>% 
  glht.table() %>% 
  kable(escape = F, digits = 3, caption = "Measurement*Layer Comparison") %>%
  kable_classic_2(full_width = F)

Measurement*Layer Comparison
	Estimate	Std. Error	z value	Pr(>\|z\|)
postLaser.All - preLaser.All	0.110	0.017	6.420	0.000
preLaser.I_II - preLaser.All	-0.215	0.018	-12.195	0.000
postLaser.I_II - preLaser.All	-0.141	0.016	-8.669	0.000
preLaser.IV - preLaser.All	-0.013	0.018	-0.718	0.999
postLaser.IV - preLaser.All	0.122	0.016	7.462	0.000
preLaser.V - preLaser.All	-0.425	0.018	-23.958	0.000
postLaser.V - preLaser.All	-0.321	0.016	-19.616	0.000
preLaser.VI - preLaser.All	-0.542	0.017	-31.176	0.000
postLaser.VI - preLaser.All	-0.466	0.016	-28.894	0.000
preLaser.I_II - postLaser.All	-0.325	0.016	-20.854	0.000
postLaser.I_II - postLaser.All	-0.251	0.014	-18.344	0.000
preLaser.IV - postLaser.All	-0.123	0.016	-7.792	0.000
postLaser.IV - postLaser.All	0.012	0.014	0.864	0.997
preLaser.V - postLaser.All	-0.536	0.016	-34.023	0.000
postLaser.V - postLaser.All	-0.431	0.014	-31.271	0.000
preLaser.VI - postLaser.All	-0.652	0.015	-42.583	0.000
postLaser.VI - postLaser.All	-0.576	0.014	-42.648	0.000
postLaser.I_II - preLaser.I_II	0.074	0.015	5.076	0.000
preLaser.IV - preLaser.I_II	0.202	0.016	12.452	0.000
postLaser.IV - preLaser.I_II	0.337	0.015	22.995	0.000
preLaser.V - preLaser.I_II	-0.211	0.016	-12.996	0.000
postLaser.V - preLaser.I_II	-0.106	0.015	-7.229	0.000
preLaser.VI - preLaser.I_II	-0.327	0.016	-20.746	0.000
postLaser.VI - preLaser.I_II	-0.251	0.014	-17.469	0.000
preLaser.IV - postLaser.I_II	0.128	0.015	8.696	0.000
postLaser.IV - postLaser.I_II	0.263	0.013	20.832	0.000
preLaser.V - postLaser.I_II	-0.284	0.015	-19.312	0.000
postLaser.V - postLaser.I_II	-0.180	0.013	-14.246	0.000
preLaser.VI - postLaser.I_II	-0.401	0.014	-28.152	0.000
postLaser.VI - postLaser.I_II	-0.324	0.012	-26.439	0.000
postLaser.IV - preLaser.IV	0.135	0.015	9.068	0.000
preLaser.V - preLaser.IV	-0.413	0.016	-25.172	0.000
postLaser.V - preLaser.IV	-0.308	0.015	-20.733	0.000
preLaser.VI - preLaser.IV	-0.529	0.016	-33.141	0.000
postLaser.VI - preLaser.IV	-0.453	0.015	-31.067	0.000
preLaser.V - postLaser.IV	-0.548	0.015	-36.917	0.000
postLaser.V - postLaser.IV	-0.443	0.013	-34.759	0.000
preLaser.VI - postLaser.IV	-0.664	0.014	-46.249	0.000
postLaser.VI - postLaser.IV	-0.588	0.012	-47.337	0.000
postLaser.V - preLaser.V	0.105	0.015	7.070	0.000
preLaser.VI - preLaser.V	-0.116	0.016	-7.305	0.000
postLaser.VI - preLaser.V	-0.040	0.015	-2.755	0.149
preLaser.VI - postLaser.V	-0.221	0.014	-15.424	0.000
postLaser.VI - postLaser.V	-0.145	0.012	-11.689	0.000
postLaser.VI - preLaser.VI	0.076	0.014	5.438	0.000

# model_df$IntTerm <- interaction(model_df$Group, model_df$Measurement, model_df$Layer)
# lme(
#   logRMS ~ IntTerm + ClickFreq + ClickOrder, 
#   random = (~1|Animal/TrialNumber), data = model_df,
#   control = lmeControl(opt = "optim")) %>% 
#   glht(linfct = mcp(IntTerm = "Tukey"), test = adjusted("holm")) %>% 
#   glht.table() %>% 
#   kable(escape = F, digits = 3, caption = "Group*Measurement*Layer Comparison") %>%
#   kable_classic_2(full_width = F)

Model Diagnostics

Random Intercepts

p1 <- ggplot(
  data = data.frame(intercept = model_full$coefficients$random$Animal)) +
  geom_histogram(
    aes(x = X.Intercept.), bins = 25, 
    fill = "lightblue", color = "black") +
  labs(title = "Intercept by Animal")

p2 <- ggplot(
  data = data.frame(intercept = model_full$coefficients$random$Animal), 
  aes(sample = X.Intercept.)) +
  stat_qq() + 
  stat_qq_line() +
  labs(title = "QQ Plot (Animal)")

p3 <- ggplot(
  data = data.frame(intercept = model_full$coefficients$random$TrialNumber)) +
  geom_histogram(
    aes(x = X.Intercept.), bins = 25, 
    fill = "lightblue", color = "black") +
  labs(title = "Intercept by Animal/TrialNumber")

p4 <- ggplot(
  data = data.frame(intercept = model_full$coefficients$random$TrialNumber), 
  aes(sample = X.Intercept.)) +
  stat_qq() + 
  stat_qq_line() +
  labs(title = "QQ Plot (Animal/TrialNumber)")

grid.arrange(p1, p2, p3, p4, nrow = 2)

Residuals Plots

res <- residuals(model_full)

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_full, resid(.) ~ fitted(.), abline=0, main = "All Data")
p2 <- plot(model_full, resid(.) ~ fitted(.)|Animal, abline=0, main = "By Animal")

grid.arrange(p1, p2, nrow = 1, widths = c(1.5,2))

Effects

p1 <- plot(Effect(c("Group"), model_full), multiline = TRUE)
p2 <- plot(Effect(c("Measurement"), model_full), multiline = TRUE)
p3 <- plot(Effect(c("Layer"), model_full), multiline = TRUE)
p4 <- plot(Effect(c("ClickFreq"), model_full), multiline = TRUE)
p5 <- plot(Effect(c("ClickOrder"), model_full), multiline = TRUE)

grid.arrange(p1, p2, p3, p4, p5, nrow = 2)

p1 <- plot(Effect(c("Group", "Measurement"), model_full), multiline = TRUE)
p2 <- plot(Effect(c("Measurement", "Layer"), model_full), multiline = TRUE)
p3 <- plot(Effect(c("Group", "Layer"), model_full), multiline = TRUE)
p4 <- plot(Effect(c("Layer", "Group", "Measurement"), model_full), multiline = TRUE)

grid.arrange(p1, p2, p3, p4, nrow = 2)

Linear Mixed Models

Ruslan Klymentiev

6/9/2021

Imports

Data Preparation

Modelling

Comparing different models

Final model

Post-hoc tests

Model Diagnostics

Random Intercepts

Residuals Plots

Effects