library(tidyverse)
library(ggpubr)
library(rstatix)
library(knitr)
library(kableExtra)
library(nlme)
library(emmeans)
library(standardize)
library(sjPlot)
library(gridExtra)
library(effectsize)

Data Import

raw_data <- read_csv(
  file = "../../Data/TonoRMS.csv",
  # specify data types of colums for later use
  col_types = cols(
    AVRECrms = col_double(),
    AVRECamp = col_double(),
    LIVrms = col_double(),
    LIVamp = col_double())
)

model_data <- raw_data
# new column with combined Measurememnts
model_data$combMeasurement <- "na"

Data Cleaning

# combine trials for post surgery for animals KIV03, KIT09
special_animals <- c("KIV03", "KIT09", "KIV16")

# select max trial from the first period
max_trials <- model_data %>% 
  filter(
    Animal %in% special_animals,
    Measurement == "Pre_6") %>% 
  group_by(Animal) %>% 
  summarise(max_trial = max(Trial))

# add max trial from the first period to the trials in the second period
for(animal in unique(max_trials$Animal)){
  model_data[(model_data$Animal == animal) & (model_data$Measurement == "Pre_7"), "Trial"] = model_data[(model_data$Animal == animal) & (model_data$Measurement == "Pre_7"), "Trial"] + max_trials[max_trials$Animal == animal, "max_trial"] %>% as.integer()
}

# create a new label of combined Measurement
model_data <- model_data %>% 
  mutate(combMeasurement = if_else(
    condition = (Animal %in% special_animals) & (Measurement %in% c("Pre_6", "Pre_7")), 
    true = "postSurgery",
    false = combMeasurement))

# combine trials for post surgery for all other animals
max_trials <- model_data %>% 
  filter(
    !(Animal %in% special_animals),
    Measurement == "Pre_2") %>% 
  group_by(Animal) %>% 
  summarise(max_trial = max(Trial))

for(animal in unique(max_trials$Animal)){
  model_data[(model_data$Animal == animal) & (model_data$Measurement == "Pre_3"), "Trial"] = model_data[(model_data$Animal == animal) & (model_data$Measurement == "Pre_3"), "Trial"] + max_trials[max_trials$Animal == animal, "max_trial"] %>% as.integer()
}

model_data <- model_data %>% 
  mutate(combMeasurement = if_else(
    condition = !(Animal %in% special_animals) & (Measurement %in% c("Pre_2", "Pre_3")), 
    true = "postSurgery",
    false = combMeasurement))

# combine trials for all animals for pre laser
max_trials <- model_data %>% 
  filter(
    Measurement == "preCLtono_3") %>% 
  group_by(Animal) %>% 
  summarise(max_trial = max(Trial))

for(animal in unique(max_trials$Animal)){
  model_data[(model_data$Animal == animal) & (model_data$Measurement == "preCLtono_4"), "Trial"] = model_data[(model_data$Animal == animal) & (model_data$Measurement == "preCLtono_4"), "Trial"] + max_trials[max_trials$Animal == animal, "max_trial"] %>% as.integer()
}

model_data <- model_data %>% 
  mutate(combMeasurement = if_else(
    condition = Measurement %in% c("preCLtono_3", "preCLtono_4"), 
    true = "preLaser",
    false = combMeasurement))

# combine trials for all animals for post laser
max_trials <- model_data %>% 
  filter(
    Measurement == "CLtono_1") %>% 
  group_by(Animal) %>% 
  summarise(max_trial = max(Trial))

for(animal in unique(max_trials$Animal)){
  model_data[(model_data$Animal == animal) & (model_data$Measurement == "AMtono_1"), "Trial"] = model_data[(model_data$Animal == animal) & (model_data$Measurement == "AMtono_1"), "Trial"] + max_trials[max_trials$Animal == animal, "max_trial"] %>% as.integer()
}

model_data <- model_data %>% 
  mutate(combMeasurement = if_else(
    condition = Measurement %in% c("CLtono_1", "AMtono_1"), 
    true = "postLaser",
    false = combMeasurement))

# combine trials for baseline except special animals
max_trials <- model_data %>% 
  filter(
    !(Animal %in% special_animals),
    Measurement == "Pre_7") %>% 
  group_by(Animal) %>% 
  summarise(max_trial = max(Trial))

for(animal in unique(max_trials$Animal)){
  model_data[(model_data$Animal == animal) & (model_data$Measurement == "Pre_8"), "Trial"] = model_data[(model_data$Animal == animal) & (model_data$Measurement == "Pre_8"), "Trial"] + max_trials[max_trials$Animal == animal, "max_trial"] %>% as.integer()
}

model_data <- model_data %>% 
  mutate(combMeasurement = if_else(
    condition = !(Animal %in% special_animals) & (Measurement %in% c("Pre_7", "Pre_8")), 
    true = "Baseline",
    false = combMeasurement))

#  baseline for special animals
model_data <- model_data %>% 
  mutate(combMeasurement = if_else(
    condition = (Animal %in% special_animals) & (Measurement  == "Pre_8"), 
    true = "Baseline",
    false = combMeasurement))

# select only relevant columns
model_data <- model_data %>% 
  filter(combMeasurement %in% c("postSurgery", "Baseline", "preLaser", "postLaser")) %>%
  select(Group, Animal, Trial, combMeasurement, AVRECrms, AVRECamp, LIVrms, LIVamp)

Modelling

AVREC logRMS

avrec_data <- model_data %>% 
  filter(!is.na(AVRECamp)) %>% # only those trials where Amplitude was detected
  mutate(
    combMeasurement = factor(
      combMeasurement, 
      levels = c("postSurgery", "Baseline", "preLaser", "postLaser")),
    Group = factor(
      Group, 
      levels = c("KIT", "KIC", "KIV")),
    # Trial = factor(Trial),
    Animal = factor(Animal),
    logAVRECrms = log(AVRECrms)
  )

Overview of the Data

ggboxplot(
  avrec_data, x = "combMeasurement", y = "logAVRECrms",
  color = "Group", palette = "lancet", 
  add = "jitter",
  main = "logAVRECrms ~ Group | Measurement"
)

ggboxplot(
  avrec_data, x = "combMeasurement", y = "AVRECrms",
  color = "Group", palette = "lancet", 
  add = "jitter",
  main = "AVRECrms ~ Group | Measurement"
)

LMM

model_final <- lme(
    fixed = logAVRECrms ~ Group*combMeasurement,
    random = (~1|Animal/Trial),
    data = avrec_data)

tab_model(model_final)

	logAVRECrms
Predictors	Estimates	CI	p
(Intercept)	-7.18	-7.42 – -6.94	<0.001
Group [KIC]	0.09	-0.28 – 0.45	0.632
Group [KIV]	0.04	-0.35 – 0.43	0.825
combMeasurement [Baseline]	0.14	0.11 – 0.17	<0.001
combMeasurement [preLaser]	0.04	0.01 – 0.06	0.017
combMeasurement [postLaser]	0.01	-0.02 – 0.04	0.362
Group [KIC] * combMeasurement [Baseline]	0.09	0.05 – 0.14	<0.001
Group [KIV] * combMeasurement [Baseline]	0.13	0.08 – 0.18	<0.001
Group [KIC] * combMeasurement [preLaser]	0.18	0.14 – 0.22	<0.001
Group [KIV] * combMeasurement [preLaser]	0.22	0.18 – 0.27	<0.001
Group [KIC] * combMeasurement [postLaser]	0.29	0.24 – 0.33	<0.001
Group [KIV] * combMeasurement [postLaser]	0.18	0.14 – 0.23	<0.001
Random Effects
σ²	0.10
τ₀₀ _Trial	0.38
τ₀₀ _Animal	0.00
N _Trial	100
N _Animal	26
Observations	9342
Marginal R² / Conditional R²	0.155 / NA

Residuals Check

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))

Effects

emmip(model_final, formula = Group ~ combMeasurement) +
  theme(aspect.ratio = 0.5)

emmip(model_final, formula = combMeasurement ~ Group) +
  theme(aspect.ratio = 0.5)

# grid.arrange(p1, p2, nrow = 2)

Post-hoc tests

grp_eff <- pairs(emmeans(model_final, ~ Group | combMeasurement))
msrmnt_eff <- pairs(emmeans(model_final, ~ combMeasurement | 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 30 tests") %>%
  kable_styling() %>%
  row_spec(
    which(df$p.value.adj < 0.05), bold = T,
    background = "#d6f5d6")

Multiple Comparisons
combMeasurement	Group	contrast	estimate	SE	t.ratio	p.value	p.value.adj	d	eff_size
postSurgery	.	KIT - KIC	-0.0856	0.1765	-0.4853	0.6321	1.0000	-0.20	small
postSurgery	.	KIT - KIV	-0.0424	0.1893	-0.2243	0.8245	1.0000	-0.09	very small
postSurgery	.	KIC - KIV	0.0432	0.1935	0.2231	0.8254	1.0000	0.09	very small
Baseline	.	KIT - KIC	-0.1785	0.1765	-1.0111	0.3225	1.0000	-0.42	small
Baseline	.	KIT - KIV	-0.1705	0.1895	-0.8996	0.3777	1.0000	-0.38	small
Baseline	.	KIC - KIV	0.0080	0.1937	0.0412	0.9675	1.0000	0.02	very small
preLaser	.	KIT - KIC	-0.2654	0.1765	-1.5039	0.1462	1.0000	-0.63	medium
preLaser	.	KIT - KIV	-0.2674	0.1893	-1.4130	0.1711	1.0000	-0.59	medium
preLaser	.	KIC - KIV	-0.0020	0.1935	-0.0104	0.9918	1.0000	0.00	very small
postLaser	.	KIT - KIC	-0.3738	0.1766	-2.1169	0.0453	0.7248	-0.88	large
postLaser	.	KIT - KIV	-0.2264	0.1894	-1.1955	0.2441	1.0000	-0.50	medium
postLaser	.	KIC - KIV	0.1474	0.1936	0.7613	0.4542	1.0000	0.32	small
.	KIT	postSurgery - Baseline	-0.1417	0.0153	-9.2848	0.0000	0.0000	-0.23	small
.	KIT	postSurgery - preLaser	-0.0354	0.0149	-2.3822	0.0172	0.2930	-0.06	very small
.	KIT	postSurgery - postLaser	-0.0142	0.0156	-0.9115	0.3621	1.0000	-0.02	very small
.	KIT	Baseline - preLaser	0.1063	0.0153	6.9312	0.0000	0.0000	0.17	very small
.	KIT	Baseline - postLaser	0.1275	0.0160	7.9487	0.0000	0.0000	0.19	very small
.	KIT	preLaser - postLaser	0.0212	0.0156	1.3602	0.1738	1.0000	0.03	very small
.	KIC	postSurgery - Baseline	-0.2346	0.0155	-15.1104	0.0000	0.0000	-0.37	small
.	KIC	postSurgery - preLaser	-0.2152	0.0155	-13.8662	0.0000	0.0000	-0.34	small
.	KIC	postSurgery - postLaser	-0.3023	0.0160	-18.9185	0.0000	0.0000	-0.46	small
.	KIC	Baseline - preLaser	0.0194	0.0156	1.2459	0.2129	1.0000	0.03	very small
.	KIC	Baseline - postLaser	-0.0678	0.0160	-4.2255	0.0000	0.0005	-0.10	very small
.	KIC	preLaser - postLaser	-0.0871	0.0160	-5.4383	0.0000	0.0000	-0.13	very small
.	KIV	postSurgery - Baseline	-0.2698	0.0195	-13.8171	0.0000	0.0000	-0.34	small
.	KIV	postSurgery - preLaser	-0.2604	0.0173	-15.0427	0.0000	0.0000	-0.37	small
.	KIV	postSurgery - postLaser	-0.1981	0.0177	-11.1747	0.0000	0.0000	-0.27	small
.	KIV	Baseline - preLaser	0.0094	0.0195	0.4804	0.6310	1.0000	0.01	very small
.	KIV	Baseline - postLaser	0.0717	0.0200	3.5892	0.0003	0.0063	0.09	very small
.	KIV	preLaser - postLaser	0.0623	0.0177	3.5119	0.0004	0.0081	0.09	very small
Note:
P value adjustment: holm method for 30 tests

Layer IV logRMS

liv_data <- model_data %>% 
  filter(!is.na(LIVamp)) %>% # only those trials where Amplitude was detected
  mutate(
    combMeasurement = factor(
      combMeasurement, 
      levels = c("postSurgery", "Baseline", "preLaser", "postLaser")),
    Group = factor(
      Group, 
      levels = c("KIT", "KIC", "KIV")),
    Trial = factor(Trial),
    Animal = factor(Animal),
    logLIVrms = log(LIVrms)
  )

Overview of the Data

ggboxplot(
  liv_data, x = "combMeasurement", y = "logLIVrms",
  color = "Group", palette = "lancet", 
  add = "jitter",
  main = "logLIVrms ~ Group | Measurement"
)

ggboxplot(
  liv_data, x = "combMeasurement", y = "LIVrms",
  color = "Group", palette = "lancet", 
  add = "jitter",
  main = "LIVrms ~ Group | Measurement"
)

LMM

model_final <- lme(
    fixed = logLIVrms ~ Group*combMeasurement,
    random = (~1|Animal/Trial),
    data = liv_data)

tab_model(model_final)

	logLIVrms
Predictors	Estimates	CI	p
(Intercept)	-7.16	-7.49 – -6.83	<0.001
Group [KIC]	0.07	-0.44 – 0.58	0.792
Group [KIV]	0.08	-0.47 – 0.63	0.767
combMeasurement [Baseline]	0.26	0.21 – 0.31	<0.001
combMeasurement [preLaser]	0.05	-0.00 – 0.10	0.060
combMeasurement [postLaser]	0.01	-0.04 – 0.06	0.738
Group [KIC] * combMeasurement [Baseline]	0.11	0.04 – 0.18	0.003
Group [KIV] * combMeasurement [Baseline]	0.12	0.03 – 0.20	0.007
Group [KIC] * combMeasurement [preLaser]	0.20	0.13 – 0.27	<0.001
Group [KIV] * combMeasurement [preLaser]	0.33	0.25 – 0.41	<0.001
Group [KIC] * combMeasurement [postLaser]	0.47	0.39 – 0.54	<0.001
Group [KIV] * combMeasurement [postLaser]	0.17	0.09 – 0.25	<0.001
Random Effects
σ²	0.32
τ₀₀ _Trial	0.53
τ₀₀ _Animal	0.00
N _Trial	100
N _Animal	26
Observations	9611
Marginal R² / Conditional R²	0.102 / NA

Residuals Check

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))

Effects

emmip(model_final, formula = Group ~ combMeasurement) +
  theme(aspect.ratio = 0.5)

emmip(model_final, formula = combMeasurement ~ Group) +
  theme(aspect.ratio = 0.5)

# grid.arrange(p1, p2, nrow = 2)

Post-hoc tests

grp_eff <- pairs(emmeans(model_final, ~ Group | combMeasurement))
msrmnt_eff <- pairs(emmeans(model_final, ~ combMeasurement | 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 30 tests") %>%
  kable_styling() %>%
  row_spec(
    which(df$p.value.adj < 0.05), bold = T,
    background = "#d6f5d6")

Multiple Comparisons
combMeasurement	Group	contrast	estimate	SE	t.ratio	p.value	p.value.adj	d	eff_size
postSurgery	.	KIT - KIC	-0.0658	0.2466	-0.2670	0.7919	1.0000	-0.11	very small
postSurgery	.	KIT - KIV	-0.0795	0.2645	-0.3004	0.7666	1.0000	-0.13	very small
postSurgery	.	KIC - KIV	-0.0136	0.2705	-0.0503	0.9603	1.0000	-0.02	very small
Baseline	.	KIT - KIC	-0.1764	0.2468	-0.7150	0.4818	1.0000	-0.30	small
Baseline	.	KIT - KIV	-0.1969	0.2651	-0.7427	0.4652	1.0000	-0.31	small
Baseline	.	KIC - KIV	-0.0205	0.2710	-0.0755	0.9404	1.0000	-0.03	very small
preLaser	.	KIT - KIC	-0.2681	0.2466	-1.0871	0.2883	1.0000	-0.45	small
preLaser	.	KIT - KIV	-0.4095	0.2645	-1.5480	0.1353	1.0000	-0.65	medium
preLaser	.	KIC - KIV	-0.1414	0.2705	-0.5226	0.6062	1.0000	-0.22	small
postLaser	.	KIT - KIC	-0.5331	0.2469	-2.1596	0.0415	0.6636	-0.90	large
postLaser	.	KIT - KIV	-0.2524	0.2647	-0.9533	0.3503	1.0000	-0.40	small
postLaser	.	KIC - KIV	0.2807	0.2707	1.0372	0.3104	1.0000	0.43	small
.	KIT	postSurgery - Baseline	-0.2622	0.0264	-9.9199	0.0000	0.0000	-0.24	small
.	KIT	postSurgery - preLaser	-0.0479	0.0254	-1.8824	0.0598	0.8973	-0.04	very small
.	KIT	postSurgery - postLaser	-0.0089	0.0265	-0.3339	0.7385	1.0000	-0.01	very small
.	KIT	Baseline - preLaser	0.2143	0.0265	8.0751	0.0000	0.0000	0.19	very small
.	KIT	Baseline - postLaser	0.2534	0.0277	9.1610	0.0000	0.0000	0.22	small
.	KIT	preLaser - postLaser	0.0390	0.0266	1.4647	0.1431	1.0000	0.03	very small
.	KIC	postSurgery - Baseline	-0.3728	0.0269	-13.8424	0.0000	0.0000	-0.33	small
.	KIC	postSurgery - preLaser	-0.2502	0.0269	-9.3018	0.0000	0.0000	-0.22	small
.	KIC	postSurgery - postLaser	-0.4761	0.0278	-17.1079	0.0000	0.0000	-0.41	small
.	KIC	Baseline - preLaser	0.1227	0.0269	4.5615	0.0000	0.0001	0.11	very small
.	KIC	Baseline - postLaser	-0.1033	0.0278	-3.7125	0.0002	0.0035	-0.09	very small
.	KIC	preLaser - postLaser	-0.2260	0.0278	-8.1411	0.0000	0.0000	-0.19	very small
.	KIV	postSurgery - Baseline	-0.3797	0.0344	-11.0242	0.0000	0.0000	-0.26	small
.	KIV	postSurgery - preLaser	-0.3779	0.0305	-12.3897	0.0000	0.0000	-0.30	small
.	KIV	postSurgery - postLaser	-0.1818	0.0311	-5.8523	0.0000	0.0000	-0.14	very small
.	KIV	Baseline - preLaser	0.0017	0.0345	0.0506	0.9596	1.0000	0.00	very small
.	KIV	Baseline - postLaser	0.1979	0.0351	5.6345	0.0000	0.0000	0.13	very small
.	KIV	preLaser - postLaser	0.1962	0.0311	6.2980	0.0000	0.0000	0.15	very small
Note:
P value adjustment: holm method for 30 tests

Clustering voltage-gated calcium channels in the auditory cortex suppresses impulse response

Linear Mixed Models. Tonotopy Analysis

Katrina Deane, Ruslan Klymentiev

06 November 2021

Data Import

Data Cleaning

Modelling

AVREC logRMS

Overview of the Data

LMM

Residuals Check

Effects

Post-hoc tests

Layer IV logRMS

Overview of the Data

LMM

Residuals Check

Effects

Post-hoc tests