This document forms part of the data and code deposited at:
https://github.com/acp29/Elmasri_GRIN2A

Load package requirements

if (!require(package="tidyverse")) utils::install.packages("tidyverse")
library(tidyverse) 
if (!require(package="lme4")) utils::install.packages("lme4")
library(lme4)  
if (!require(package="lmerTest")) utils::install.packages("lmerTest")
library(lmerTest)  
if (!require(package="parameters")) utils::install.packages("parameters")
library(parameters) 
if (!require(package="car")) utils::install.packages("car")
library(car)  
if (!require(package="stats")) utils::install.packages("stats")
library(stats)
if (!require(package="pCalibrate")) utils::install.packages("pCalibrate")
library(pCalibrate)
if (!require(package="afex")) utils::install.packages("afex")
library(afex)
if (!require(package="emmeans")) utils::install.packages("emmeans")
library(emmeans)
if (!require(package="knitr")) utils::install.packages("knitr")
library(knitr)
if (!require(package="kableExtra")) utils::install.packages("kableExtra")
library(kableExtra)

Read text in from file

Data <- read.delim("../data/n2a_dose_rescue_nmdar.dat", header = TRUE)

Factor encoding

Data$dose <- as.factor(Data$dose)

Set genotype WT and cre - as reference levels

Data$dose<- factor(Data$dose, levels=c("Unt","0","1.1","3.3","3.8","4.4","5.6","10.9"))
#Data$doseinterpret_bf(1)factor(Data$dose, ordered = TRUE, levels=c("Unt","0","1.1","3.3","3.8","4.4","5.6","10.9"))

Fit a linear model

# Initialize
variates <- c("peak","decay")
l <- length(variates)

for (i in 1:l) {
  
variates[i] -> resp
  
cat('\n\n\n# Analysis of',resp,'\n\n')
    
# Model fitting  
formula <- sprintf("log(%s) ~ dose", resp)
model <- lm(formula, data=Data) 

# ANOVA table
aov <- car::Anova(model, type = 3, test.statistic = "F")
aov %>%
  rename(p.value = "Pr(>F)") %>%
  mutate(p.value = afex::round_ps_apa(p.value)) %>%
  knitr::kable(caption = sprintf("**ANOVA table: %s**",resp), digits = 2) %>% 
  kableExtra::kable_styling(bootstrap_options = c("striped", "hover")) %>%
  print()

# Checking model assumptions
# Plot of standardized residuals to assess for outliers and homo/heteroskedasticity of residuals
resid = residuals(model)
plot(fitted(model),resid/sd(resid), 
     main = sprintf("Residuals Plot: log%s",resp), 
     xlab = "Fitted values", ylab = 
       "Standardized Residuals")
abline(h = 0)
# Histogram to asses Normality of residuals
hist(resid, col="darkgray", main = sprintf("Histogram of residuals: log%s",resp))
# Normal Q-Q plot of the residuals
qqnorm(resid, main = sprintf("Normal Q-Q Plot: log%s",resp))
qqline(resid)  # Do the points fall nicely onto the line?

# Extract model parameters and return them on the response scale 
# After backtransforming to the response scale, the highlighted coefficients correspond to the fold change 
# between mutant and WT 
model_parameters(model, exponentiate = TRUE) %>%   
  dplyr::select(-c(p,SE)) %>%
  knitr::kable(caption = sprintf("**Model parameters (on the response scale) with 95%% confidence intervals: %s**",resp), digits = 2) %>% 
  kableExtra::kable_styling(bootstrap_options = c("striped", "hover")) %>%
  print()

# Use emmmeans to obtain estimated marginal means from our model
emm <- emmeans(model, ~ dose, tran = 'log', type = 'response')
emm %>% 
  as.data.frame() %>%
  dplyr::select(-SE) %>%
  knitr::kable(caption = sprintf("**Standardized effect sizes (Cohen's *d*) with marginal 95%% confidence intervals: %s**",resp), digits = 2) %>% 
  kableExtra::kable_styling(bootstrap_options = c("striped", "hover")) %>%
  print()

}

Analysis of peak

ANOVA table: peak
Sum Sq Df F value p.value
(Intercept) 0.39 1 1.71 .193
dose 6.96 7 4.31 <.001
Residuals 33.88 147 NA NA
Model parameters (on the response scale) with 95% confidence intervals: peak
Parameter Coefficient CI CI_low CI_high t df_error
(Intercept) 0.90 0.95 0.76 1.06 -1.31 147
dose0 0.65 0.95 0.52 0.82 -3.69 147
dose1.1 0.80 0.95 0.56 1.14 -1.24 147
dose3.3 1.03 0.95 0.76 1.39 0.17 147
dose3.8 1.13 0.95 0.84 1.52 0.82 147
dose4.4 1.06 0.95 0.79 1.42 0.36 147
dose5.6 0.67 0.95 0.50 0.89 -2.78 147
dose10.9 0.96 0.95 0.72 1.28 -0.30 147
Standardized effect sizes (Cohen’s d) with marginal 95% confidence intervals: peak
dose response df lower.CL upper.CL
Unt 0.90 147 0.76 1.06
0 0.59 147 0.50 0.69
1.1 0.72 147 0.52 0.98
3.3 0.92 147 0.71 1.18
3.8 1.01 147 0.79 1.29
4.4 0.95 147 0.74 1.21
5.6 0.60 147 0.47 0.76
10.9 0.86 147 0.68 1.09

Analysis of decay

ANOVA table: decay
Sum Sq Df F value p.value
(Intercept) 0.00 1 0.06 .801
dose 19.01 7 37.70 <.001
Residuals 10.59 147 NA NA
Model parameters (on the response scale) with 95% confidence intervals: decay
Parameter Coefficient CI CI_low CI_high t df_error
(Intercept) 0.99 0.95 0.90 1.08 -0.25 147
dose0 1.94 0.95 1.70 2.20 10.27 147
dose1.1 1.65 0.95 1.35 2.01 4.95 147
dose3.3 1.21 0.95 1.03 1.44 2.27 147
dose3.8 1.03 0.95 0.87 1.21 0.32 147
dose4.4 0.78 0.95 0.66 0.92 -2.91 147
dose5.6 0.77 0.95 0.66 0.91 -3.12 147
dose10.9 0.79 0.95 0.67 0.93 -2.89 147
Standardized effect sizes (Cohen’s d) with marginal 95% confidence intervals: decay
dose response df lower.CL upper.CL
Unt 0.99 147 0.90 1.08
0 1.91 147 1.75 2.09
1.1 1.63 147 1.36 1.94
3.3 1.20 147 1.04 1.38
3.8 1.01 147 0.88 1.16
4.4 0.77 147 0.68 0.89
5.6 0.77 147 0.67 0.87
10.9 0.78 147 0.68 0.89