library(tidyverse)
library(nlme)
library(lme4)

# ------------------------------------------------------------------------------------------------------------
# Experiment 1

exp1.pre <- read.csv("experiment1.csv")

# Model building and final model selection
M1.Full <- lme(rt ~ hand * response, random = ~1 | subject, method="ML", data = exp1.pre)
M1.A <- update(M1.Full, .~. - hand:response)
anova(M1.Full, M1.A)
M1.Final <- lme(rt ~ hand * response, random = ~1 | subject, method="REML", data = exp1.pre)
summary(M1.Final)

# Compute effect size
lmerfit1 <- lmer(rt ~ hand * response + (1|subject) + (1|trial), data = exp1.pre)
randomvar <- VarCorr(lmerfit1)
print(randomvar, comp=c("Variance"))
randomvar <- as.data.frame(randomvar)
trialvar <- randomvar$vcov[1]
subjectvar <- randomvar$vcov[2]
resid <- randomvar$vcov[3]
fixed <- getME(lmerfit1, "beta")
diff <- fixed[4]
effectsize <- diff/sqrt(trialvar + subjectvar + resid)
effectsize

# ------------------------------------------------------------------------------------------------------------
# Experiment 2

exp2.pre <- read.csv("experiment2.csv")

# Model building and model testing
M1.Full <- lme(rt ~ hand * response, random = ~1 | subject, method="ML", data = exp2.pre)
M1.A <- update(M1.Full, .~. - hand:response)
anova(M1.Full, M1.A)
M1.Final <- lme(rt ~ hand * response, random = ~1 | subject, method="REML", data = exp2.pre)
summary(M1.Final)

# Compute effect size
lmerfit1 <- lmer(rt ~ hand * response + (1|subject) + (1|trial), data = exp2.pre)
randomvar <- VarCorr(lmerfit1)
print(randomvar, comp=c("Variance"))
randomvar <- as.data.frame(randomvar)
trialvar <- randomvar$vcov[1]
subjectvar <- randomvar$vcov[2]
resid <- randomvar$vcov[3]
fixed <- getME(lmerfit1, "beta")
diff <- fixed[4]
effectsize <- diff/sqrt(trialvar + subjectvar + resid)
effectsize

# ------------------------------------------------------------------------------------------------------------
# Experiment 3

exp3.pre <- read.csv("experiment3.csv")

# Current linear model and assumption testing

M1.1 <- lme(rt ~ hand * response * space, random = ~1 | subject, method="ML", data = exp3.pre)
summary(M1.1)
M1.2 <- update(M1.1, .~. - hand:space:response)
summary(M1.2)
anova(M1.1, M1.2)
M1.3 <- update(M1.2, .~. - space:response)
anova(M1.2, M1.3)
M1.4 <- update(M1.3, .~. - hand:response)
anova(M1.3, M1.4)
M1.5 <- update(M1.3, .~. - hand:space)
anova(M1.3, M1.5)

M1.Full <- lme(rt ~ hand * response + space * hand, random = ~1 | subject, method="REML", data = exp3.pre)
summary(M1.Full)

#calculate effect size
lmerfit1 <- lmer(rt ~ hand * response + space * hand + (1|subject) + (1|trial), data = exp3.pre)
randomvar <- VarCorr(lmerfit1)
print(randomvar, comp=c("Variance"))
randomvar <- as.data.frame(randomvar)
trialvar <- randomvar$vcov[1]
subjectvar <- randomvar$vcov[2]
resid <- randomvar$vcov[3]
fixed <- getME(lmerfit1, "beta")
diff1 <- fixed[6]
diff2 <- fixed[7]
effectsize_response_interaction <- diff1/sqrt(trialvar + subjectvar + resid)
effectsize_space_interaction <- diff2/sqrt(trialvar + subjectvar + resid)

op <- par(mfrow = c(2,2))
#check assumptions
plot(M1.Final, add.smooth = FALSE)