Analysis.R

library(tidyverse)

# file path management.
library(here)

# predictions from regression models
library(marginaleffects)

# t test summaries
library(rempsyc)

# plotting
library(patchwork)
theme_set(theme_bw())

#load data
all_data <- read.csv(here("Data","Data.csv"))

#rename data
all_data <- all_data %>% 
  rename(
    language_status = `TD.0.LT.1`,
    atypical = `T3Atypical.Error`,
    omission = `T3Omission.Error`,
    total_atypical = Atypical,
    total_omissions = Omissions,
    elicited = `T3_Consonants_Elicited`,
    t4pcc = `T4_PCC`,
    t4ss = 'T4_SS',
    t4pipa_ss_r = 'T4_pipa_ss_r',
    t4pipa_ra_r = 'T4_pipa_ra_r',
    t4pipa_al_r = 'T4_pipa_al_r',
    t4pipa_pi_r = 'T4_pipa_pi_r',
    t4pipa_ps_r = 'T4_pipa_ps_r',
    t4pipa_lk_r = 'T4_pipa_lk_r',
    t4pipa_ss_s = 'T4_pipa_ss_s',
    t4pipa_ra_s = 'T4_pipa_ra_s',
    t4pipa_al_s = 'T4_pipa_al_s',
    t4pipa_pi_s = 'T4_pipa_pi_s',
    t4pipa_ps_s = 'T4_pipa_ps_s'
  ) %>% 
  select(
    Participant, language_status, atypical, omission, total_atypical, 
    total_omissions, elicited, t4pcc, t4ss,
    t4pipa_ss_r, t4pipa_ra_r, t4pipa_al_r, t4pipa_pi_r, t4pipa_ps_r,
    t4pipa_lk_r, t4pipa_ss_s, t4pipa_ra_s, t4pipa_al_s, t4pipa_pi_s,
    t4pipa_ps_s, av_pipa_s, av_alpips_s
  )

# Remove participants with identified hearing loss at T4
all_data <- all_data |> 
  filter(
    !Participant %in% c(4, 83, 93, 109, 113, 166)
  )

#exploratory visualisation
summary(all_data)

all_data <- all_data %>% 
  mutate(
    omission = as.numeric(omission),
    atypical = as.numeric(atypical),
    t4pcc = as.numeric(t4pcc),
    t4ss = as.numeric(t4ss),
    av_pipa_s = as.numeric(av_pipa_s),
    av_alpips_s = as.numeric(av_alpips_s),
    language_status = factor(language_status)
  ) %>% 
  filter(
    !is.na(omission),
    !is.na(atypical),
    !is.na(language_status),#check participant 91
    t4ss != 0
  )

all_data %>%
  ggplot(
    aes(
      x = factor(language_status),
      y = as.numeric(omission)
    )
  ) +
  geom_boxplot() +
  labs(
    x = "Group",
    y = "Omission Patterns"
  ) +
  scale_x_discrete(labels = c("TD", "LT")
  )

all_data %>%
  ggplot(
    aes(
      x = factor(language_status),
      y = as.numeric(atypical)
    )
  ) +
  geom_boxplot() +
  labs(
    x = "Group",
    y = "Atypical Patterns"
  ) +
  scale_x_discrete(labels = c("TD", "LT")
  )

# Scale predictor variables
all_data <- all_data %>% 
  filter(complete.cases(all_data)) %>%
  mutate(
    omission_s = scale(omission)[,1],
    atypical_s = scale(atypical)[,1]
  )

## t-tests for RQ1

t.test.results <- nice_t_test(
  data = all_data,
  response = c("atypical", "omission", "t4ss", "av_pipa_s"),
  group = "language_status",
  var.equal = TRUE
)
t.test.results
# Dependent Variable         t df            p          d   CI_lower   CI_upper
# 1           atypical -3.506832 99 6.833078e-04 -0.7500256 -1.1802174 -0.3162610
# 2           omission -6.400145 99 5.205450e-09 -1.3688346 -1.8263564 -0.9055343
# 3               t4ss  3.513664 99 6.678851e-04  0.7514868  0.3176700  1.1817247
# 4          av_pipa_s  2.528972 99 1.301892e-02  0.5408853  0.1136739  0.9654496

# Group means
all_data |> 
  group_by(language_status) |> 
  summarise(
    atypical_mean = mean(atypical),
    atypical_sd = sd(atypical),
    omission_mean = mean(omission),
    omission_sd = sd(omission),
    t4ss_mean = mean(t4ss),
    t4ss_sd = sd(t4ss),
    av_pipa_s_mean = mean(av_pipa_s),
    av_pipa_s_sd = sd(av_pipa_s)
  )
# language_status atypical_mean atypical_sd omission_mean omission_sd t4ss_mean t4ss_sd av_pipa_s_mean av_pipa_s_sd
# <fct>                   <dbl>       <dbl>         <dbl>       <dbl>     <dbl>   <dbl>          <dbl>        <dbl>
#   1 0                      0.0817      0.0918         0.106       0.104      9.51    3.54           11.4         2.20
# 2 1                      0.157       0.116          0.269       0.147      6.91    3.29           10.2         2.36


## Regression models for RQ2

# Check vif for collinearity in simple model (no interactions).
simple_model <- lm(t4ss ~ language_status + omission_s + atypical_s, data = all_data)

car::vif(simple_model)
# language_status      omission_s      atypical_s
# 1.454953        2.181255        1.664774
# OK (FWIW).

# Fit complex model with interactions
complex_model <- lm(t4ss ~ language_status * omission_s * atypical_s, data = all_data)

summary(complex_model)
# Call:
# lm(formula = t4ss ~ language_status * omission_s * atypical_s, 
#     data = all_data)
# 
# Residuals:
#     Min      1Q  Median      3Q     Max 
# -7.9420 -2.1454 -0.0128  2.3560  6.7910 
# 
# Coefficients:
#                                        Estimate Std. Error t value Pr(>|t|)    
# (Intercept)                              8.5776     0.4870  17.615  < 2e-16 ***
# language_status1                        -1.0853     0.8905  -1.219 0.226013    
# omission_s                              -2.3928     0.6182  -3.871 0.000201 ***
# atypical_s                               0.4964     0.5226   0.950 0.344569    
# language_status1:omission_s              1.6978     1.0900   1.558 0.122723    
# language_status1:atypical_s             -1.1357     1.0103  -1.124 0.263863    
# omission_s:atypical_s                    0.4211     0.4869   0.865 0.389324    
# language_status1:omission_s:atypical_s  -0.1852     0.5984  -0.309 0.757693    
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 3.263 on 93 degrees of freedom
# Multiple R-squared:  0.2576,	Adjusted R-squared:  0.2018 
# F-statistic: 4.611 on 7 and 93 DF,  p-value: 0.0001839

## Diagnostic plots
complex_resid <- tibble(
  Fitted = predict(complex_model),
  Residuals = resid(complex_model)
)

residual_plot <- complex_resid |> 
  ggplot(
    aes(
      x = Fitted,
      y = Residuals
    )
  ) +
  geom_point() +
  geom_smooth() +
  labs(
    title = "Residual plot for `complex_model`"
  )

residual_plot

ggsave(plot = residual_plot, filename = here('Plots', 'Model1_Residuals.png'))

qq_plot <- complex_resid |> 
  ggplot(
    aes(
      sample = Residuals
    )
  ) +
  stat_qq() +
  stat_qq_line() +
  labs(
    title = "QQ plot for `complex_model`"
  )

qq_plot
  
ggsave(plot = qq_plot, filename = here('Plots', 'Model1_QQ.png'))


# Condition number
X <- model.matrix(complex_model)
kappa(X, exact = TRUE)
# 11.41 - weak collinearity

complex_model_pipa <- lm(
  av_pipa_s ~ language_status * omission_s * atypical_s, 
  data = all_data
)

summary(complex_model_pipa)
# Call:
# lm(formula = av_pipa_s ~ language_status * omission_s * atypical_s, 
#     data = all_data)
# 
# Residuals:
#     Min      1Q  Median      3Q     Max 
# -6.3925 -1.5367  0.1291  1.5824  3.9178 
# 
# Coefficients:
#                                        Estimate Std. Error t value Pr(>|t|)    
# (Intercept)                            10.99302    0.32846  33.468   <2e-16 ***
# language_status1                       -0.62486    0.60068  -1.040    0.301    
# omission_s                             -0.94311    0.41698  -2.262    0.026 *  
# atypical_s                             -0.24765    0.35248  -0.703    0.484    
# language_status1:omission_s             0.50406    0.73525   0.686    0.495    
# language_status1:atypical_s             0.63670    0.68146   0.934    0.353    
# omission_s:atypical_s                  -0.02314    0.32841  -0.070    0.944    
# language_status1:omission_s:atypical_s -0.01248    0.40363  -0.031    0.975    
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 2.201 on 93 degrees of freedom
# Multiple R-squared:  0.1553,	Adjusted R-squared:  0.09174 
# F-statistic: 2.443 on 7 and 93 DF,  p-value: 0.02413

## Diagnostic plots
complex_resid <- tibble(
  Fitted = predict(complex_model_pipa),
  Residuals = resid(complex_model_pipa)
)

residual_plot <- complex_resid |> 
  ggplot(
    aes(
      x = Fitted,
      y = Residuals
    )
  ) +
  geom_point() +
  geom_smooth() +
  labs(
    title = "Residual plot for `complex_model_pipa`"
  )

residual_plot

ggsave(plot = residual_plot, filename = here('Plots', 'Model2_Residuals.png'))

qq_plot <- complex_resid |> 
  ggplot(
    aes(
      sample = Residuals
    )
  ) +
  stat_qq() +
  stat_qq_line() +
  labs(
    title = "QQ plot for `complex_model_pipa`"
  )

qq_plot

ggsave(plot = qq_plot, filename = here('Plots', 'Model2_QQ.png'))

# kappa (i.e. condition number) is the same, as model predictors are the same!

# Post-hoc test test of pipa subsets.
complex_model_t4pipa_al_s <- lm(t4pipa_al_s ~ language_status * omission_s * atypical_s, data = all_data)

summary(complex_model_t4pipa_al_s)
# Call:
# lm(formula = av_pipa_s ~ language_status * omission_s * atypical_s, 
#     data = all_data)
# 
# Residuals:
#     Min      1Q  Median      3Q     Max 
# -6.3925 -1.5367  0.1291  1.5824  3.9178 
# 
# Coefficients:
#                                        Estimate Std. Error t value Pr(>|t|)    
# (Intercept)                            10.99302    0.32846  33.468   <2e-16 ***
# language_status1                       -0.62486    0.60068  -1.040    0.301    
# omission_s                             -0.94311    0.41698  -2.262    0.026 *  
# atypical_s                             -0.24765    0.35248  -0.703    0.484    
# language_status1:omission_s             0.50406    0.73525   0.686    0.495    
# language_status1:atypical_s             0.63670    0.68146   0.934    0.353    
# omission_s:atypical_s                  -0.02314    0.32841  -0.070    0.944    
# language_status1:omission_s:atypical_s -0.01248    0.40363  -0.031    0.975    
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 2.201 on 93 degrees of freedom
# Multiple R-squared:  0.1553,	Adjusted R-squared:  0.09174 
# F-statistic: 2.443 on 7 and 93 DF,  p-value: 0.02413

complex_model_t4pipa_pi_s <- lm(t4pipa_pi_s ~ language_status * omission_s * atypical_s, data = all_data)

summary(complex_model_t4pipa_pi_s)
# Call:
# lm(formula = t4pipa_pi_s ~ language_status * omission_s * atypical_s, 
#     data = all_data)
# 
# Residuals:
#     Min      1Q  Median      3Q     Max 
# -8.3822 -1.2246  0.7859  1.8372  3.9053 
# 
# Coefficients:
#                                        Estimate Std. Error t value Pr(>|t|)    
# (Intercept)                             12.8730     0.4387  29.342   <2e-16 ***
# language_status1                        -0.5189     0.8023  -0.647   0.5194    
# omission_s                              -1.1453     0.5570  -2.056   0.0426 *  
# atypical_s                              -0.5429     0.4708  -1.153   0.2518    
# language_status1:omission_s              1.0292     0.9821   1.048   0.2974    
# language_status1:atypical_s              0.5628     0.9102   0.618   0.5379    
# omission_s:atypical_s                    0.1989     0.4387   0.453   0.6514    
# language_status1:omission_s:atypical_s  -0.2179     0.5391  -0.404   0.6871    
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 2.94 on 93 degrees of freedom
# Multiple R-squared:  0.1348,	Adjusted R-squared:  0.0697 
# F-statistic:  2.07 on 7 and 93 DF,  p-value: 0.05443

complex_model_t4pipa_ps_s <- lm(t4pipa_ps_s ~ language_status * omission_s * atypical_s, data = all_data)

summary(complex_model_t4pipa_ps_s)
# Call:
# lm(formula = t4pipa_ps_s ~ language_status * omission_s * atypical_s, 
#     data = all_data)
# 
# Residuals:
#     Min      1Q  Median      3Q     Max 
# -5.3670 -2.3412 -0.5521  2.1293  5.8575 
# 
# Coefficients:
#                                        Estimate Std. Error t value Pr(>|t|)    
# (Intercept)                             12.0047     0.4722  25.423   <2e-16 ***
# language_status1                        -2.1762     0.8635  -2.520   0.0134 *  
# omission_s                               0.0923     0.5995   0.154   0.8780    
# atypical_s                              -0.5006     0.5067  -0.988   0.3258    
# language_status1:omission_s             -1.4402     1.0570  -1.363   0.1763    
# language_status1:atypical_s              1.6244     0.9797   1.658   0.1007    
# omission_s:atypical_s                   -0.4744     0.4721  -1.005   0.3176    
# language_status1:omission_s:atypical_s   0.5583     0.5803   0.962   0.3385    
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 3.164 on 93 degrees of freedom
# Multiple R-squared:  0.1685,	Adjusted R-squared:  0.1059 
# F-statistic: 2.692 on 7 and 93 DF,  p-value: 0.01385

complex_model_av_alpips_s <- lm(av_alpips_s ~ language_status * omission_s * atypical_s, data = all_data)

summary(complex_model_av_alpips_s)
# Call:
# lm(formula = av_alpips_s ~ language_status * omission_s * atypical_s, 
#     data = all_data)
# 
# Residuals:
#     Min      1Q  Median      3Q     Max 
# -6.3397 -1.3465  0.0557  1.5662  4.3908 
# 
# Coefficients:
#                                        Estimate Std. Error t value Pr(>|t|)    
# (Intercept)                            12.10936    0.37363  32.410   <2e-16 ***
# language_status1                       -1.11506    0.68327  -1.632    0.106    
# omission_s                             -0.73569    0.47432  -1.551    0.124    
# atypical_s                             -0.45300    0.40094  -1.130    0.261    
# language_status1:omission_s             0.28095    0.83634   0.336    0.738    
# language_status1:atypical_s             0.89256    0.77516   1.151    0.252    
# omission_s:atypical_s                   0.04229    0.37356   0.113    0.910    
# language_status1:omission_s:atypical_s -0.11805    0.45913  -0.257    0.798    
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 2.504 on 93 degrees of freedom
# Multiple R-squared:  0.1664,	Adjusted R-squared:  0.1036 
# F-statistic: 2.652 on 7 and 93 DF,  p-value: 0.01517

# Group means by language status.
aggregate(all_data$t4pipa_pi_s, list(all_data$language_status), FUN=mean)

aggregate(all_data$t4pipa_al_s, list(all_data$language_status), FUN=mean)

aggregate(all_data$t4pipa_ps_s, list(all_data$language_status), FUN=mean)

aggregate(all_data$av_alpips_s, list(all_data$language_status), FUN=mean)

# Plots by status
omission_by_status <- all_data %>%
  ggplot(
    aes(
      x = factor(language_status),
      y = as.numeric(omission)
    )
  ) +
  geom_boxplot() +
  labs(
    x = "Group",
    y = "Omission Patterns"
  ) +
  scale_x_discrete(labels = c("TD", "LT")
  )

atypical_by_status <- all_data %>%
  ggplot(
    aes(
      x = factor(language_status),
      y = as.numeric(atypical)
    )
  ) +
  geom_boxplot() +
  labs(
    x = "Group",
    y = "Atypical Patterns"
  ) +
  scale_x_discrete(labels = c("TD", "LT")
  )

ss_by_status <- all_data %>%
  ggplot(
    aes(
      x = factor(language_status),
      y = t4ss
    )
  ) +
  geom_boxplot() +
  labs(
    x = "Group",
    y = "DEAP Standard Score"
  ) +
  scale_x_discrete(labels = c("TD", "LT")
  )

pipa_by_status <- all_data %>%
  ggplot(
    aes(
      x = factor(language_status),
      y = av_pipa_s
    )
  ) +
  geom_boxplot() +
  labs(
    x = "Group",
    y = "Average PIPA Standard Score"
  ) +
  scale_x_discrete(labels = c("TD", "LT")
  )

combined_plot <- (atypical_by_status + omission_by_status) / 
  (ss_by_status + pipa_by_status) + 
  plot_annotation(tag_levels = "A")

ggsave(
  filename = here("Plots", "combined_plot.png"),
  plot = combined_plot,
  dpi = 300,
  unit = "mm",
  width = 150,
  height = 150
)

# Generate predictions from models and plot.
omission_ss_preds <- avg_predictions(
  complex_model,
  newdata = datagrid(
    "omission_s" = seq(min(all_data$omission_s), max(all_data$omission_s), 0.1),
    "atypical_s" = seq(min(all_data$atypical_s), max(all_data$atypical_s), 0.1),
    "language_status" = c(0, 1)
  ),
  by = "omission_s"
)

omission_ss_preds <- omission_ss_preds |> 
  mutate(
    omission = omission_s * sd(all_data$omission) + mean(all_data$omission)
  )

ss_preds <- all_data |> 
  ggplot(
    aes(
      x = omission, 
      y = t4ss
    )
  ) +
  geom_point() +
  geom_ribbon(
    aes(
      x = omission,
      ymin = conf.low,
      ymax = conf.high
    ),
    data = omission_ss_preds,
    inherit.aes = FALSE,
    fill = "grey",
    alpha = 0.3,
    linetype = 0
  ) +
  geom_line(
    aes(
      x = omission,
      y = estimate
    ),
    inherit.aes = FALSE,
    data = omission_ss_preds,
    linewidth = 1,
    colour = "black"
  ) +
  labs(
    x = "Omission Patterns",
    y = "DEAP Standard Score"
  )

ggsave(
  filename = here("Plots", "t4ss.png"),
  plot = ss_preds,
  dpi = 300,
  unit = "mm",
  width = 100,
  height = 100
)

omission_pipa_preds <- avg_predictions(
  complex_model_pipa,
  newdata = datagrid(
    "omission_s" = seq(min(all_data$omission_s), max(all_data$omission_s), 0.1),
    "atypical_s" = seq(min(all_data$atypical_s), max(all_data$atypical_s), 0.1),
    "language_status" = c(0, 1)
  ),
  by = "omission_s"
)

omission_pipa_preds <- omission_pipa_preds |> 
  mutate(
    omission = omission_s * sd(all_data$omission) + mean(all_data$omission)
  )

pipa_preds <- all_data |> 
  ggplot(
    aes(
      x = omission, 
      y = av_pipa_s
    )
  ) +
  geom_point() +
  geom_ribbon(
    aes(
      x = omission,
      ymin = conf.low,
      ymax = conf.high
    ),
    data = omission_pipa_preds,
    inherit.aes = FALSE,
    fill = "grey",
    alpha = 0.3,
    linetype = 0
  ) +
  geom_line(
    aes(
      x = omission,
      y = estimate
    ),
    inherit.aes = FALSE,
    data = omission_pipa_preds,
    linewidth = 1,
    colour = "black"
  ) +
  labs(
    x = "Omission Patterns",
    y = "Average PIPA Standard Score"
  )

ggsave(
  filename = here("Plots", "pipa.png"),
  plot = pipa_preds,
  dpi = 300,
  unit = "mm",
  width = 100,
  height = 100
)

# Additional analysis using omissions/elicited and atypical/elicited.
all_data <- all_data |> 
  mutate(
    omissions_alt = total_omissions/elicited,
    atypical_alt = total_atypical/elicited,
    omissions_alt_s = scale(omissions_alt)[,1],
    atypical_alt_s = scale(atypical_alt)[,1]
  )

# How much do they correlate with our previous measures?
cor.test(all_data$omission, all_data$omissions_alt)
# Pearson's product-moment correlation
# 
# data:  all_data$omission and all_data$omissions_alt
# t = 3.8207, df = 82, p-value = 0.0002579
# alternative hypothesis: true correlation is not equal to 0
# 95 percent confidence interval:
#  0.1901979 0.5567357
# sample estimates:
#       cor 
# 0.3887412 

cor.test(all_data$atypical, all_data$atypical_alt)
# Pearson's product-moment correlation
# 
# data:  all_data$atypical and all_data$atypical_alt
# t = 2.0971, df = 82, p-value = 0.03907
# alternative hypothesis: true correlation is not equal to 0
# 95 percent confidence interval:
#  0.01179026 0.41970868
# sample estimates:
#       cor 
# 0.2256151 

# RQ1
t.test.results <- nice_t_test(
  data = all_data,
  response = c(
    "atypical", "atypical_alt", "omission", "omissions_alt", 
    "t4ss", "av_pipa_s"
  ),
  group = "language_status",
  var.equal = TRUE
)
t.test.results
# Dependent Variable         t df            p          d    CI_lower   CI_upper
# 1           atypical -2.518128 82 1.374306e-02 -0.6838418 -1.22425765 -0.1394101
# 2       atypical_alt -4.611165 82 1.456955e-05 -1.2522424 -1.81446259 -0.6832314
# 3           omission -4.884897 82 5.045934e-06 -1.3265793 -1.89258834 -0.7534704
# 4      omissions_alt -6.116032 82 3.114443e-08 -1.6609155 -2.24643189 -1.0670756
# 5               t4ss  2.861096 82 5.354819e-03  0.7769806  0.22937786  1.3200671
# 6          av_pipa_s  1.895638 82 6.153098e-02  0.5147936 -0.02478601  1.0513030

# Group means (supplemented)
all_data |> 
  group_by(language_status) |> 
  summarise(
    atypical_mean = mean(atypical),
    atypical_sd = sd(atypical),
    alt_atypical_mean = mean(atypical_alt),
    alt_atypical_sd = sd(atypical_alt),
    omission_mean = mean(omission),
    omission_sd = sd(omission),
    alt_omission_mean = mean(omissions_alt),
    alt_omission_sd = sd(omissions_alt),
    t4ss_mean = mean(t4ss),
    t4ss_sd = sd(t4ss),
    av_pipa_s_mean = mean(av_pipa_s),
    av_pipa_s_sd = sd(av_pipa_s)
  ) |> 
  t()
# [,1]         [,2]        
# language_status   "0"          "1"         
# atypical_mean     "0.08186873" "0.15014833"
# atypical_sd       "0.09315422" "0.12368434"
# alt_atypical_mean "0.03339244" "0.08760951"
# alt_atypical_sd   "0.03732553" "0.06212970"
# omission_mean     "0.1032513"  "0.2429937" 
# omission_sd       "0.1026642"  "0.1157275" 
# alt_omission_mean "0.08783669" "0.18606592"
# alt_omission_sd   "0.05845176" "0.06190607"
# t4ss_mean         "9.507463"   "6.764706"  
# t4ss_sd           "3.539274"   "3.491587"  
# av_pipa_s_mean    "11.55522"   "10.51765"  
# av_pipa_s_sd      "1.988397"   "2.123745"  

# RQ2 - alternative measures

# Check vif for collinearity in simple model (no interactions).
simple_model <- lm(
  t4ss ~ language_status + omissions_alt_s + atypical_alt_s, 
  data = all_data
)

car::vif(simple_model)
# language_status omissions_alt_s  atypical_alt_s 
# 1.466368        2.460665        2.127997 
# OK (FWIW).

# Fit complex model with interactions
complex_model_alt <- lm(
  t4ss ~ language_status * omissions_alt_s * atypical_alt_s, 
  data = all_data
)

summary(complex_model_alt)
# Call:
# lm(formula = t4ss ~ language_status * omissions_alt_s * atypical_alt_s, 
#     data = all_data)
# 
# Residuals:
#     Min      1Q  Median      3Q     Max 
# -6.7251 -2.6863  0.0555  2.3231  7.2171 
# 
# Coefficients:
#                                                 Estimate Std. Error t value Pr(>|t|)    
# (Intercept)                                       8.8241     0.5692  15.503   <2e-16 ***
# language_status1                                 -2.4238     1.5047  -1.611   0.1114    
# omissions_alt_s                                  -0.1419     0.7439  -0.191   0.8492    
# atypical_alt_s                                   -1.4963     0.8154  -1.835   0.0704 .  
# language_status1:omissions_alt_s                  1.3567     1.6091   0.843   0.4018    
# language_status1:atypical_alt_s                   3.2051     1.7732   1.808   0.0746 .  
# omissions_alt_s:atypical_alt_s                    0.6428     0.7524   0.854   0.3956    
# language_status1:omissions_alt_s:atypical_alt_s  -2.1334     1.2117  -1.761   0.0823 .  
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 3.432 on 76 degrees of freedom
# Multiple R-squared:  0.2036,	Adjusted R-squared:  0.1303 
# F-statistic: 2.776 on 7 and 76 DF,  p-value: 0.01262

## Diagnostic plots
complex_resid <- tibble(
  Fitted = predict(complex_model_alt),
  Residuals = resid(complex_model_alt)
)

residual_plot <- complex_resid |> 
  ggplot(
    aes(
      x = Fitted,
      y = Residuals
    )
  ) +
  geom_point() +
  geom_smooth() +
  labs(
    title = "Residual plot for `complex_model` (alt measures)"
  )

residual_plot

ggsave(plot = residual_plot, filename = here('Plots', 'Model1_Alt_Residuals.png'))

qq_plot <- complex_resid |> 
  ggplot(
    aes(
      sample = Residuals
    )
  ) +
  stat_qq() +
  stat_qq_line() +
  labs(
    title = "QQ plot for `complex_model` (alt measures)"
  )

qq_plot

ggsave(plot = qq_plot, filename = here('Plots', 'Model1_Alt_QQ.png'))


# Condition number
X <- model.matrix(complex_model_alt)
kappa(X, exact = TRUE)
# 13.12 - weak collinearity

complex_model_pipa_alt <- lm(
  av_pipa_s ~ language_status * omissions_alt_s * atypical_alt_s, 
  data = all_data
)

summary(complex_model_pipa_alt)
# Call:
# lm(formula = av_pipa_s ~ language_status * omissions_alt_s * 
#     atypical_alt_s, data = all_data)
# 
# Residuals:
#     Min      1Q  Median      3Q     Max 
# -4.1010 -1.2566  0.1278  1.4077  3.5168 
# 
# Coefficients:
#                                                 Estimate Std. Error t value Pr(>|t|)    
# (Intercept)                                     11.02929    0.32913  33.510   <2e-16 ***
# language_status1                                -0.83721    0.87008  -0.962    0.339    
# omissions_alt_s                                 -0.48305    0.43017  -1.123    0.265    
# atypical_alt_s                                  -0.06299    0.47151  -0.134    0.894    
# language_status1:omissions_alt_s                 0.81457    0.93047   0.875    0.384    
# language_status1:atypical_alt_s                 -0.26489    1.02533  -0.258    0.797    
# omissions_alt_s:atypical_alt_s                   0.79666    0.43507   1.831    0.071 .  
# language_status1:omissions_alt_s:atypical_alt_s -0.64626    0.70067  -0.922    0.359    
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 1.984 on 76 degrees of freedom
# Multiple R-squared:  0.1394,	Adjusted R-squared:  0.06013 
# F-statistic: 1.759 on 7 and 76 DF,  p-value: 0.1082

## Diagnostic plots
complex_resid <- tibble(
  Fitted = predict(complex_model_pipa_alt),
  Residuals = resid(complex_model_pipa_alt)
)

residual_plot <- complex_resid |> 
  ggplot(
    aes(
      x = Fitted,
      y = Residuals
    )
  ) +
  geom_point() +
  geom_smooth() +
  labs(
    title = "Residual plot for `complex_model_pipa` (alt measures)"
  )

residual_plot

ggsave(plot = residual_plot, filename = here('Plots', 'Model2_Alt_Residuals.png'))

qq_plot <- complex_resid |> 
  ggplot(
    aes(
      sample = Residuals
    )
  ) +
  stat_qq() +
  stat_qq_line() +
  labs(
    title = "QQ plot for `complex_model_pipa` (alt measures)"
  )

qq_plot

ggsave(plot = qq_plot, filename = here('Plots', 'Model2_Alt_QQ.png'))

# kappa (i.e. condition number) is the same, as model predictors are the same!