## ----rt0 --------------------------
load(url("http://pnb.mcmaster.ca/bennett/psy710/datasets/rtData-mw12-1.rda"))
sapply(rt.wide,class)
summary(rt.long)


## ----rt1 ---------------------------
options(width=80,digits=4)
options(contrasts=c("contr.sum","contr.poly"))
# following model is incorrect
rt.aov.00 <- aov(rt ~ 1+Error(subj/(angle*noise)),data=rt.long) # no fixed effects!
# the next 2 models are equivalent:
rt.aov.01 <- aov(rt ~ angle*noise+Error(subj/(angle*noise)),data=rt.long)
# rt.aov.01b <- aov(rt ~ angle*noise+Error(subj/(angle+noise+angle:noise)),data=rt.long)
summary(rt.aov.01)

library(afex)
rt.aov.02 <-aov_car(rt ~ angle*noise+Error(subj/(angle*noise)),data=rt.long)
summary(rt.aov.02)


## ----effectsize -------------------------
omega.sq.angle <- 2*(112711-14460) / (225422 + 260280 + 344423 + 38269)
omega.sq.noise <- 1*(234624-14950) / (234624 + 134551 + 344423 + 38269)
omega.sq.a_x_n <- 2*(93992-12522) / (187983 + 225402 + 344423 + 38269)
c(omega.sq.angle , omega.sq.noise, omega.sq.a_x_n)

cohens.f.angle <- sqrt(omega.sq.angle/(1-omega.sq.angle) )
cohens.f.noise <- sqrt(omega.sq.noise/(1-omega.sq.noise) )
cohens.f.noise <- sqrt(omega.sq.a_x_n/(1-omega.sq.a_x_n) )
c(cohens.f.angle , cohens.f.noise, omega.sq.a_x_n)


## ----effectsize-2 --------------------------
library(effectsize)
omega_squared(rt.aov.02)


## ----effectsize-2b -------------------------
eta_squared(rt.aov.02)
cohens_f(rt.aov.02)


## ----intPlot ------
( tmp <- with(rt.long,tapply(rt,list(noise,angle),mean)) )
absent.means <- tmp["absent",]
present.means <- tmp["present",]
x.angle <- c(0,4,8);
plot(x=c(0,4,8),absent.means,"b",ylim=c(450,800),ylab="RT",xlab="angle")
points(x=c(0,4,8),present.means,"b",pch=19)
legend(x=0,y=750,legend=c("absent","present"),pch=c(21,19))


## ----simpleEffects1 ------------------------------
aov.angle.absent <- aov_car(rt~angle+Error(subj/angle),data=subset(rt.long,noise=="absent"))
aov.angle.present <- aov_car(rt~angle+Error(subj/angle),data=subset(rt.long,noise=="present"))
# Hyunh-Feldt correction; pes == partial eta squared:
anova(aov.angle.absent,correction="HF",es="pes")
anova(aov.angle.present,correction="HF",es="pes")


## ----linearTrend2 -------------------
rt.mat <- as.matrix(rt.wide[,2:7])
rt.mat[1:2,]
lin.C <- c(-1,0,1,-1,0,1)
rt.lin <- rt.mat %*% lin.C
t.test(rt.lin)


## ----linearTrend3 -----------------------
myC <- c(-1,0,1,1,0,-1)
rt.lin.x.noise <- rt.mat %*% myC
t.test(rt.lin.x.noise)


rt.mat[1,] # inspect column names
# note the contrast weights in next 2 lines:
rt.absent.lin <- rt.mat %*% c(-1,0,1,0,0,0)
rt.present.lin <- rt.mat %*% c(0,0,0,-1,0,1)
t.test(rt.absent.lin)
t.test(rt.present.lin)


## ----linearTrend4 ------------------------
N <- length(rt.absent.lin) # number of subjects
linTrend <- c(rt.present.lin,rt.absent.lin)
nz <- as.factor(x=rep(c("present","absent"),each=N))
sid <- factor(x=rep(1:N,times=2),label="s")
linTrend.df <- data.frame(sid,nz,linTrend)
summary(linTrend.df)


## ----linearTrend5 ---------------------------
options(contrasts=c("contr.sum","contr.poly"))
linTrend.aov.01 <- aov_car(linTrend~nz+Error(sid/nz),data=linTrend.df)
summary(linTrend.aov.01)


## ----linearTrend-emmeans -----------------------
library(emmeans)
rt.emm <- emmeans(rt.aov.02,specs="angle")
lin.C <- c(-1,0,1) # trend weights

 # linear trend ignoring noise
contrast(rt.emm,method=list(linear=lin.C))

# linear trend x noise interaction:
rt.emm.2 <- emmeans(rt.aov.02,specs="angle",by="noise")
contrast(rt.emm.2,interaction=list(angle=list(lin.C),noise=list(c(-1,1))),by=NULL)


# linear trend separately for each noise
contrast(rt.emm.2,method=list(linear=lin.C)) 


## ----r-alt-1 ----------------------------
# using lmer:
library(lmerTest)
rt.lmer.01 <- lmer(rt~angle*noise+(1|subj),data=rt.long)
anova(rt.lmer.01) # assumes sphericity


## ----r-alt-2 --------------------------------
rt.lmer.00 <- lmer(rt~angle+noise+(1|subj),data=rt.long)
anova(rt.lmer.00,rt.lmer.01) 


## ----r-alt-3 -----------------------------------
# conservative F tests:
N <- 10 # 10 subjects
F.angle <- 8.18
F.angle.x.noise <- 6.82
# conservative F tests for angle and angle x noise:
1-pf(F.angle,1,N) # p value for angle
1-pf(F.angle.x.noise,1,N) # p value for angle x noise


## ----r-alt-4 -------------------------------
library(effectsize)
omega_squared(rt.lmer.01)
cohens_f(rt.lmer.01)


## ----r-alt-5 ----------------------------
# random components
ranova(rt.lmer.01) # anova-like table
# variance components
print(VarCorr(rt.lmer.01),comp=c("Variance","Std.Dev."))
# association strength
library(performance)
icc(rt.lmer.01)


## ----r-alt-6 ---------------------------------
# using lme in nlme
library(nlme)
# assumes independence:
rt.lme.01 <- lme(rt~angle*noise,random=~1|subj,
                 data=rt.long)
# assumes sphericity:
rt.lme.02 <- lme(rt~angle*noise,random=~1|subj,
                 data=rt.long,
                 correlation=corCompSymm(value=0.3,form=~1|subj))
# does not assume sphericity:
rt.lme.03 <- lme(rt~angle*noise,random=~1|subj,
                 data=rt.long,
                 correlation=corSymm(form=~1|subj))
# no significant difference in fit:
anova(rt.lme.01,rt.lme.02,rt.lme.03) 


## ----r-alt-7 -------------------------------
# fixed effects for model 1
anova(rt.lme.01)
omega_squared(rt.lme.01)
cohens_f(rt.lme.01)


## ----r-alt-8 ------------------------------------
# variance components (independence)
print(VarCorr(rt.lme.01),comp=c("Variance","Std.Dev."))
( icc.subj <- 4081/(4081+13783) )

# variance components (compound symmetry)
print(VarCorr(rt.lme.02),comp=c("Variance","Std.Dev."))
( icc.subj <- 1142/(1142+16722) )


## ----split2 -----------------------------
load(url("http://pnb.mcmaster.ca/bennett/psy710/datasets/rtVisual.rda"))
sapply(rtVisual.wide,class)
summary(rtVisual)


## ----split3 ---------------------------------
rtAge.aov.01 <- aov(rt ~ group*angle + Error(subj/angle),data=rtVisual)
summary(rtAge.aov.01)


## ----split3b -------------------------
library(afex)
options(contrasts=c("contr.sum","contr.poly"))
rtAge.aov.02 <- aov_car(rt ~ group*angle + Error(subj/angle),data=rtVisual)
summary(rtAge.aov.02)


## ----split5 ------------------------------
rt.mat <- rtVisual.wide[,3:5] # extract rt measures
rtVisual.wide$rtAvg <- rowMeans(rt.mat) # calculate average
summary(aov(rtAvg~group,data=rtVisual.wide))


## ----split6 ------------------------------
rt.young <- aov(rt~angle+Error(subj/angle),data=subset(rtVisual,group=="young"))
rt.old <- aov(rt~angle+Error(subj/angle),data=subset(rtVisual,group=="old"))


## ----split7 --------------------------------
summary(rt.young)[[2]]
summary(rt.old)[[2]]


## ----split8 ---------------------------------
library(effectsize)
# eta_squared(rtAge.aov.02)
omega_squared(rtAge.aov.02)
cohens_f(rtAge.aov.02)


## ----splitSimple1 -------------------------------
levels(rtVisual$angle)
summary(aov(rt~1+group,data=subset(rtVisual,angle=="a1")))
summary(aov(rt~1+group,data=subset(rtVisual,angle=="a2")))
summary(aov(rt~1+group,data=subset(rtVisual,angle=="a3")))


## ----splitSimple2 -------------------------------
anova(aov_car(rt~angle+Error(subj/angle),
              correction="HF",
              es="pes",
              data=subset(rtVisual,group=="young")))

anova(aov_car(rt~angle+Error(subj/angle),
              correction="HF",
              es="pes",
              data=subset(rtVisual,group=="old")))


## ----lin-con-between -----------------------------
library(emmeans)
rtAge.emm <- emmeans(rtAge.aov.02,specs="group")
contrast(rtAge.emm,method="pairwise") # ignores angle


# compare age at each angle:
rtAge.emm.2 <- emmeans(rtAge.aov.02,specs="group",by="angle")
contrast(rtAge.emm.2,method="pairwise")


## ----linC1 -----------------------
sapply(rtVisual.wide,class)
rt.mat<-as.matrix(rtVisual.wide[,3:5] )
lin.C <- c(-1,0,1)
rtVisual.wide$linTrend <- rt.mat %*% lin.C
rtVisual.wide[1:3,]


## ----linC2 ----------------------
t.test(linTrend~group,data=rtVisual.wide)


## ----linC3 -----------------------------------------
t.test(rtVisual.wide$linTrend)


## ----linC4 ------------------------------------------
options(contrasts=c("contr.sum","contr.poly"))
lin.scores.aov <- aov(linTrend~group,data=rtVisual.wide)
summary(lin.scores.aov,intercept=T)


## ----linC5-emmeans -----------------------------
library(emmeans)
# create aov_car object:
rtAge.aov.10 <- aov_car(rt ~ group*angle + Error(subj/angle),data=rtVisual)
linTrendWeights <- c(-1,0,1)

# linear trend ignoring group
angle.emm <- emmeans(rtAge.aov.10,specs="angle")
contrast(angle.emm,method=list(linTrendWeights))

# linear trend for each group
angle.emm.2 <- emmeans(rtAge.aov.10,specs="angle",by="group")
contrast(angle.emm.2,method=list(linTrendWeights))

# linear trend x group interaction
contrast(angle.emm.2,interaction=c("poly","consec"),by=NULL)


## ----sp-lmer-1 -----------------------------
# fixed effects:
library(lmerTest)
library(effectsize)
rtAge.lmer.01 <- lmer(rt ~ group*angle + (1|subj),data=rtVisual)
anova(rtAge.lmer.01) # assumes sphericity
omega_squared(rtAge.lmer.01)
cohens_f(rtAge.lmer.01)


## ----sp-lmer-2 -------------------------------------
ranova(rtAge.lmer.01)
print(VarCorr(rtAge.lmer.01),comp=c("Variance","Std.Dev."))
library(performance)
icc(rtAge.lmer.01) # intraclass correlation


## ----sp-lmer-3 -----------------------------
library(nlme)
rtAge.lme.01 <- lme(rt ~ group*angle,
                    random=~1|subj,
                    data=rtVisual,
                    correlation=corCompSymm(value=0.3,form=~1|subj))
# fixed effects:
anova(rtAge.lme.01)
omega_squared(rtAge.lme.01)
cohens_f(rtAge.lme.01)


## ----sp-lmer-4 ------------------------
print(VarCorr(rtAge.lme.01),comp=c("Variance","Std.Dev."))
554/(554+6597) # small intraclass correlation


# assume independence:
rtAge.lme.00 <- lme(rt ~ group*angle,
                    random=~1|subj,
                    data=rtVisual)
print(VarCorr(rtAge.lme.00),comp=c("Variance","Std.Dev."))


## ----sp-lmer-5 -----------------------
shapiro.test(residuals(rtAge.lmer.01))
shapiro.test(residuals(rtAge.lme.01))


## ----C12qqLMER ------
op <- par(no.readonly = T)
par(mfrow=c(1,2),cex.axis=0.5,cex.main=0.8)
qqnorm(residuals(rtAge.lmer.01),main="LMER model");qqline(residuals(rtAge.lmer.01))
qqnorm(residuals(rtAge.lme.01),main="LME model 01");qqline(residuals(rtAge.lme.01))