# Computes by-subject and by-item ANOVA,
# minF',
# repeated-measures regression across subjects,
# linear mixed-effect model across subjects
# and across both subjects and items
# for the data in "ventoneR.txt"

############# Load data #############

vt = read.table("ventoneR.txt",T)
attach(vt)

############# By-subject ANOVA #############

# Create by-subject data frame
SToneX = NULL
SBegin = NULL
SF0 = NULL
SSubj = NULL
for (s in 1:10) {
 for (t in c(-1,1)) {
  for (b in c(-1,1)) {
   SToneX = c(SToneX,t)
   SBegin = c(SBegin,b)
   SF0 = c(SF0,mean(F0[Subj == s & ToneX == t & Begin == b]))
   SSubj = c(SSubj, s)
  }
 }
}
Svt = data.frame(SSubj,SToneX,SBegin,SF0)

# Run by-subject analysis
SSubj.f = as.factor(SSubj)
SToneX.f = as.factor(SToneX)
SBegin.f = as.factor(SBegin)
summary(aov(SF0~SToneX.f*SBegin.f+Error(SSubj.f/(SToneX.f*SBegin.f)),data=Svt))

############# By-item ANOVA #############

# Create by-item data frame
WToneX = NULL
WBegin = NULL
WF0 = NULL
WWord = NULL
for (w in 1:20) {
 for (b in c(-1,1)) {
  WBegin = c(WBegin,b)
  if (w < 11) {
   WToneX = c(WToneX,1)
  }
  if (w > 10) {
   WToneX = c(WToneX,-1)
  }
  WF0 = c(WF0,mean(F0[Word == w & Begin == b]))
  WWord = c(WWord, w)
 }
}
Wvt = data.frame(WWord,WToneX,WBegin,WF0)

# Run by-item analysis
WWord.f = as.factor(WWord)
WToneX.f = as.factor(WToneX)
WBegin.f = as.factor(WBegin)
summary(aov(WF0~WToneX.f*WBegin.f+Error(WWord.f/(WBegin.f)),data=Svt))

############# minF' by hand #############

# I can't figure out how to extract F and df automatically from aov outputs
# so I'll copy them by hand from R's displayed ANOVA tables

F1ToneX = 1.3298
F2ToneX = 3.3129
dfToneX = 1
df1ToneX = 9
df2ToneX = 18

F1Begin = 6.2678
F2Begin = 26.826
dfBegin = 1
df1Begin = 9
df2Begin = 18

F1TxB = 21.995
F2TxB = 81.927
dfTxB = 1
df1TxB = 9
df2TxB = 18

minFToneX = (F1ToneX*F2ToneX)/(F1ToneX+F2ToneX)
minFBegin = (F1Begin*F2Begin)/(F1Begin+F2Begin)
minFTxB = (F1TxB*F2TxB)/(F1TxB+F2TxB)

dfnumToneX = dfToneX
dfdenToneX = (F1ToneX+F2ToneX)^2/(F1ToneX^2/df2ToneX+F2ToneX^2/df1ToneX)
dfnumBegin =dfBegin
dfdenBegin =(F1Begin+F2Begin)^2/(F1Begin^2/df2Begin+F2Begin^2/df1Begin)
dfnumTxB =dfTxB
dfdenTxB =(F1TxB+F2TxB)^2/(F1TxB^2/df2TxB+F2TxB^2/df1TxB)

p.ToneX = pf(minFToneX,dfnumToneX,dfdenToneX, lower.tail = F)
p.Begin = pf(minFBegin,dfnumBegin,dfdenBegin, lower.tail = F)
p.TxB = pf(minFTxB,dfnumTxB,dfdenTxB, lower.tail = F)

minFToneX
minFBegin
minFTxB

p.ToneX
p.Begin
p.TxB

############# By-subject repeated-measures regression #############

# Prepare by-subject repeated-measures regression
int.coef = NULL
ToneX.coef = NULL
Begin.coef = NULL
TxB.coef = NULL
for (i in 1:10) {
 F0.i = F0[Subj == i]
 ToneX.i = ToneX[Subj == i]
 Begin.i = Begin[Subj == i]
 dat.i = data.frame(F0.i,ToneX.i,Begin.i)
 lm.i = summary(lm(F0.i~ToneX.i*Begin.i, data = dat.i))
 int.coef = c(int.coef,lm.i$coef[1,1])
 ToneX.coef = c(ToneX.coef,lm.i$coef[2,1])
 Begin.coef = c(Begin.coef,lm.i$coef[3,1])
 TxB.coef = c(TxB.coef,lm.i$coef[4,1])
}

# Compute one-sample t tests on coefficients
M.int = mean(int.coef)
SE.int = sd(int.coef)/sqrt(10)
t.int = M.int/SE.int
p.int = 2*pt(-abs(t.int),9)
cat("int.coef = ", M.int, "\n")
cat("int SE = ", SE.int, "\n")
cat("int t = ", t.int, "\n")
cat("int p = ", p.int, "\n")

M.ToneX = mean(ToneX.coef)
SE.ToneX = sd(ToneX.coef)/sqrt(10)
t.ToneX = M.ToneX/SE.ToneX
p.ToneX = 2*pt(-abs(t.ToneX),9)
cat("ToneX.coef = ", M.ToneX, "\n")
cat("ToneX SE = ", SE.ToneX, "\n")
cat("ToneX t = ", t.ToneX, "\n")
cat("ToneX p = ", p.ToneX, "\n")

M.Begin = mean(Begin.coef)
SE.Begin = sd(Begin.coef)/sqrt(10)
t.Begin = M.Begin/SE.Begin
p.Begin = 2*pt(-abs(t.Begin),9)
cat("Begin.coef = ", M.Begin, "\n")
cat("Begin SE = ", SE.Begin, "\n")
cat("Begin t = ", t.Begin, "\n")
cat("Begin p = ", p.Begin, "\n")

M.TxB = mean(TxB.coef)
SE.TxB = sd(TxB.coef)/sqrt(10)
t.TxB = M.TxB/SE.TxB
p.TxB = 2*pt(-abs(t.TxB),9)
cat("TxB.coef = ", M.TxB, "\n")
cat("TxB SE = ", SE.TxB, "\n")
cat("TxB t = ", t.TxB, "\n")
cat("TxB p = ", p.TxB, "\n")

############# By-subject LME using "nlme" package #############

library("nlme") # Make sure you install it first.

vt.g = groupedData(F0~1|Subj, data=vt)
vt.g.s = lme(F0~ToneX*Begin, random=~1|Subj, data=vt.g)
summary(vt.g.s)

############# LME using "lme4" package #############

library("lme4")  # Make sure you install it first; loading takes about 30 seconds

vt.s = lmer(F0 ~ ToneX*Begin + (1|Subj), data = vt)
vt.s
vt.sw = lmer(F0 ~ ToneX*Begin + (1|Subj) + (1|Word), data = vt)
vt.sw
anova(vt.s,vt.sw) # If significant, use vt.sw; otherwise vt.s is sufficient.