# March 11, 2024

# install.packages("faux")

library(tidyverse)
library(lmSupport)
library(car)
library(faux) # for creating variables correlated with an existing variable

# fake data
set.seed(2001)
treatment <- round(rnorm(n = 10, mean = 12, sd = 3), 0)
control <- round(rnorm(n = 10, mean = 10, sd = 4), 0)

# group means
mean(treatment); mean(control)

# a t-test to compare the groups
t.test(treatment, control, var.equal = T)

# putting the data into a data frame/tibble
d <- tibble(group = c(rep("tx", 10), rep("ctrl", 10)),
						correct = c(treatment, control))

# adding a contrast code
d <- d %>% 
	mutate(con1 = ifelse(group == "ctrl", -1/2, 1/2))

# modeling
model1 <- lm(correct ~ con1, d)
summary(model1) # hey! it's the t-test again

# creating a fake covariate ("foreign language exposure")
e <- resid(model1)
FLE <- round(rnorm_pre(x = e, mu = 50, sd = 10, r = 0.2), 0)

# add FLE to the tibble, mean-centered
d <- d %>% 
	mutate(FLE = FLE - mean(FLE))

# new model (the ANCOVA)
model2 <- lm(correct ~ con1 + FLE, d)
modelSummary(model2)

# conventionally, ANCOVA is sometimes executed sequentially
step1 <- lm(correct ~ FLE, d)
step2 <- lm(correct ~ con1 + FLE, d)
modelCompare(step1, step2)

# adjusted means via the effects package

# the effects package likes factors to be FACTORs
d <- d %>% 
	mutate(groupF = as.factor(con1))

ANCOVA <- lm(correct ~ groupF + FLE, d)
modelSummary(ANCOVA, t = F) # same as it was!

# install.packages("effects")
library(effects)
effect("groupF", ANCOVA)