R.Muralikrishnan, MPI for Empirical Aesthetics. This work is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International Licence.

#Create a variable called behavdata and read into it, the data from the file allresfilter_3words_split.out.

behavdata <- read.table("allres_filtersplit.out")  #Input


#Give names to each field in the data just read.
names(behavdata) <- c("SNo", "Block", "Trial", "Item", "Condition", "Sentence", "WO", "ST", "CT", "RT", "Acceptability", "Answer", "Type")

#The new versions of R do not take a purely numeric input column as a factor,
#unless it is EXPLICITLY DECLARED to be taken as a factor.  Since the 'Item' column in our input files
#is a purely numeric one, this has to be explicity declared to be taken as a factor, 
#as follows: - (Thanks to Petra for pointing this out).
behavdata$Item <- as.factor(behavdata$Item) 


# Example row from allresfilter_3words_split.out
# KU01 1   1   63    DAI   63DAI     DA  I   CQ  1.628 2               2         D
# SNo  Bl  Tr  Item  Cond  Sentence  WO  ST  CT  RT    Acceptability   Accuracy  => Type D: Dative sentence; N for Nominative and F for Filler sentences

#behavdata_Correct <- behavdata[behavdata$F1=="2",]

b_Dt   <- behavdata[behavdata$Type=="D",]  # Only Dative sentences.
b_Dtc  <- b_Dt[b_Dt$Answer=="2",]    # Correctly answered Dative targets;


#Indicate R the file name under which the results are to be saved.
sink("BehavdataDStats.txt")	#Output

print("Total")
print(dim(behavdata))
print("Correctly answered D Targets")
print(dim(b_Dtc))


#Add generous comments to the output file
print("+-------------------------------------------------------------------------+")
print("| KURINJI: Statistics for Behavioural Data                                 |")
print("+-------------------------------------------------------------------------+")
print("")

# For Reaction Time
# print(mean_rt_all <- aggregate(behavdata$RT, by=list(WO=behavdata$WO, AC=behavdata$AC), mean)) 
# print(sd_rt_all <- aggregate(behavdata$RT,, by=list(WO=behavdata$WO, AC=behavdata$AC), sd))
# This is incorrect, because we are taking the means without considering subjects.

#First calculate a mean across conditions for each subject.
rt_per_subj <- aggregate(b_Dtc$RT,list(Subj=b_Dtc$SNo, WO=b_Dtc$WO, ST=b_Dtc$ST, CT=b_Dtc$CT),mean)
# RT in Kurinji in seconds already in the input!
#Then calculate the mean and SD of these means.
print("Mean: Reaction Time (for correct answers) across subjects by condition")
print(rt_m <- aggregate(rt_per_subj$x,list(WO=rt_per_subj$WO, ST=rt_per_subj$ST, CT=rt_per_subj$CT),mean))
print("################################")

print("SD: Reaction Time (for correct answers) across subjects by condition")
print(rt_sd <- aggregate(rt_per_subj$x,list(WO=rt_per_subj$WO, ST=rt_per_subj$ST, CT=rt_per_subj$CT),sd))
print("############## SUBJECTS ANALYSIS ##################")

print(paste("Repeated Measures ANOVA for D targets: Reaction Time -- Factors WO, ST and CT  -- Subjects"))
rt_subj <- aggregate(b_Dtc$RT, by=list(Subj=b_Dtc$SNo, WO=b_Dtc$WO, ST=b_Dtc$ST, CT=b_Dtc$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*ST*CT + Error(Subj/(WO*ST*CT)), data=rt_subj)))
#This is called repeated measures, because it is a within-subjects design, and so each subject would've seen all the conditions...and so we have a measure for each condition (many times...since there are many items in the same condition!) from each participant.  Had it been a design where condition A is seen only by a set of 10 people and condition B by another set, then it will be One-way ANOVA!
#First calculate a mean across conditions for each subject.

#Resolve for Word-order
b_Dtcn <- b_Dtc[b_Dtc$WO=="DA",]
b_Dtcs <- b_Dtc[b_Dtc$WO=="AD",]

print("No. of correctly answered D targets with WO = DA")
print(dim(b_Dtcn))
print("No. of correctly answered D targets with WO = AD")
print(dim(b_Dtcs))


print("############# WO = DA ###################")

print(paste("ANOVA for D targets: Reaction Time -- WO = DA, Factors ST and CT  -- Subjects"))
rt_subj_n <- aggregate(b_Dtcn$RT, by=list(Subj=b_Dtcn$SNo, WO=b_Dtcn$WO, ST=b_Dtcn$ST, CT=b_Dtcn$CT), mean, NA.rm=T)
print(summary (aov(x ~ ST*CT + Error(Subj/(ST*CT)), data=rt_subj_n)))

print("############# WO = AD ###################")

print(paste("ANOVA for D targets: Reaction Time -- WO = AD, Factors ST and CT  -- Subjects"))
rt_subj_s <- aggregate(b_Dtcs$RT, by=list(Subj=b_Dtcs$SNo, WO=b_Dtcs$WO, ST=b_Dtcs$ST, CT=b_Dtcs$CT), mean, NA.rm=T)
print(summary (aov(x ~ ST*CT + Error(Subj/(ST*CT)), data=rt_subj_s)))

# Resolve for Context pairwise
b_Dtc_c <- b_Dtc[b_Dtc$CT=="CQ",]
b_Dtc_n <- b_Dtc[b_Dtc$CT=="NQ",]
b_Dtc_v <- b_Dtc[b_Dtc$CT=="VQ",]
b_Dtc_w <- b_Dtc[b_Dtc$CT=="WQ",]
# Six possible context pairs possible: cn, cv, cw, nv, nw, vw!
b_Dtc_cn <- merge(b_Dtc_c, b_Dtc_n, all=T)
b_Dtc_cv <- merge(b_Dtc_c, b_Dtc_v, all=T)
b_Dtc_vn <- merge(b_Dtc_v, b_Dtc_n, all=T)
b_Dtc_cw <- merge(b_Dtc_c, b_Dtc_w, all=T)
b_Dtc_wn <- merge(b_Dtc_w, b_Dtc_n, all=T)
b_Dtc_vw <- merge(b_Dtc_v, b_Dtc_w, all=T)

print(paste("########## CT Pairwise Comparison: RT by Subjects #########"))
print(paste("### CT = CQ + NQ: ANOVA for D targets: Reaction Time -- Factors WO, ST and CT  -- Subjects"))
rt_subj_cn <- aggregate(b_Dtc_cn$RT, by=list(Subj=b_Dtc_cn$SNo, WO=b_Dtc_cn$WO, ST=b_Dtc_cn$ST, CT=b_Dtc_cn$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*ST*CT + Error(Subj/(WO*ST*CT)), data=rt_subj_cn)))

print(paste("### CT = CQ + VQ: ANOVA for D targets: Reaction Time -- Factors WO, ST and CT  -- Subjects"))
rt_subj_cv <- aggregate(b_Dtc_cv$RT, by=list(Subj=b_Dtc_cv$SNo, WO=b_Dtc_cv$WO, ST=b_Dtc_cv$ST, CT=b_Dtc_cv$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*ST*CT + Error(Subj/(WO*ST*CT)), data=rt_subj_cv)))

print(paste("### CT = VQ + NQ: ANOVA for D targets: Reaction Time -- Factors WO, ST and CT  -- Subjects"))
rt_subj_vn <- aggregate(b_Dtc_vn$RT, by=list(Subj=b_Dtc_vn$SNo, WO=b_Dtc_vn$WO, ST=b_Dtc_vn$ST, CT=b_Dtc_vn$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*ST*CT + Error(Subj/(WO*ST*CT)), data=rt_subj_vn)))

print(paste("### CT = CQ + WQ: ANOVA for D targets: Reaction Time -- Factors WO, ST and CT  -- Subjects"))
rt_subj_cw <- aggregate(b_Dtc_cw$RT, by=list(Subj=b_Dtc_cw$SNo, WO=b_Dtc_cw$WO, ST=b_Dtc_cw$ST, CT=b_Dtc_cw$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*ST*CT + Error(Subj/(WO*ST*CT)), data=rt_subj_cw)))

print(paste("### CT = WQ + NQ: ANOVA for D targets: Reaction Time -- Factors WO, ST and CT  -- Subjects"))
rt_subj_wn <- aggregate(b_Dtc_wn$RT, by=list(Subj=b_Dtc_wn$SNo, WO=b_Dtc_wn$WO, ST=b_Dtc_wn$ST, CT=b_Dtc_wn$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*ST*CT + Error(Subj/(WO*ST*CT)), data=rt_subj_wn)))

print(paste("### CT = VQ + WQ: ANOVA for D targets: Reaction Time -- Factors WO, ST and CT  -- Subjects"))
rt_subj_vw <- aggregate(b_Dtc_vw$RT, by=list(Subj=b_Dtc_vw$SNo, WO=b_Dtc_vw$WO, ST=b_Dtc_vw$ST, CT=b_Dtc_vw$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*ST*CT + Error(Subj/(WO*ST*CT)), data=rt_subj_vw)))


print("================================================")
print("############# ITEMS ANALYSIS ###################")

print(paste("Repeated Measures ANOVA: Reaction Time -- Factors WO, ST and CT -- Items"))
rt_item <- aggregate(b_Dtc$RT, by=list(Item=b_Dtc$Item, WO=b_Dtc$WO, ST=b_Dtc$ST, CT=b_Dtc$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*ST*CT + Error(Item/(WO*ST*CT)), data=rt_item)))

print("############# WO = DA ###################")

print(paste("Repeated Measures ANOVA: Reaction Time -- WO = DA, Factors ST and CT  -- Items"))
rt_item_n <- aggregate(b_Dtcn$RT, by=list(Item=b_Dtcn$Item, WO=b_Dtcn$WO, ST=b_Dtcn$ST, CT=b_Dtcn$CT), mean, NA.rm=T)
print(summary (aov(x ~ ST*CT + Error(Item/(ST*CT)), data=rt_item_n)))

print("############# WO = AD ###################")

print(paste("Repeated Measures ANOVA: Reaction Time -- WO = AD, Factors ST and CT  -- Items"))
rt_item_s <- aggregate(b_Dtcs$RT, by=list(Item=b_Dtcs$Item, WO=b_Dtcs$WO, ST=b_Dtcs$ST, CT=b_Dtcs$CT), mean, NA.rm=T)
print(summary (aov(x ~ ST*CT + Error(Item/(ST*CT)), data=rt_item_s)))

print(paste("########## CT Pairwise Comparison: RT by Items #########"))
print(paste("### CT = CQ + NQ: ANOVA for D targets: Reaction Time -- Factors WO, ST and CT  -- Items"))
rt_item_cn <- aggregate(b_Dtc_cn$RT, by=list(Item=b_Dtc_cn$Item, WO=b_Dtc_cn$WO, ST=b_Dtc_cn$ST, CT=b_Dtc_cn$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*ST*CT + Error(Item/(WO*ST*CT)), data=rt_item_cn)))

print(paste("### CT = CQ + VQ: ANOVA for D targets: Reaction Time -- Factors WO, ST and CT  -- Items"))
rt_item_cv <- aggregate(b_Dtc_cv$RT, by=list(Item=b_Dtc_cv$Item, WO=b_Dtc_cv$WO, ST=b_Dtc_cv$ST, CT=b_Dtc_cv$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*ST*CT + Error(Item/(WO*ST*CT)), data=rt_item_cv)))

print(paste("### CT = VQ + NQ: ANOVA for D targets: Reaction Time -- Factors WO, ST and CT  -- Items"))
rt_item_vn <- aggregate(b_Dtc_vn$RT, by=list(Item=b_Dtc_vn$Item, WO=b_Dtc_vn$WO, ST=b_Dtc_vn$ST, CT=b_Dtc_vn$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*ST*CT + Error(Item/(WO*ST*CT)), data=rt_item_vn)))

print(paste("### CT = CQ + WQ: ANOVA for D targets: Reaction Time -- Factors WO, ST and CT  -- Items"))
rt_item_cw <- aggregate(b_Dtc_cw$RT, by=list(Item=b_Dtc_cw$Item, WO=b_Dtc_cw$WO, ST=b_Dtc_cw$ST, CT=b_Dtc_cw$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*ST*CT + Error(Item/(WO*ST*CT)), data=rt_item_cw)))

print(paste("### CT = WQ + NQ: ANOVA for D targets: Reaction Time -- Factors WO, ST and CT  -- Items"))
rt_item_wn <- aggregate(b_Dtc_wn$RT, by=list(Item=b_Dtc_wn$Item, WO=b_Dtc_wn$WO, ST=b_Dtc_wn$ST, CT=b_Dtc_wn$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*ST*CT + Error(Item/(WO*ST*CT)), data=rt_item_wn)))

print(paste("### CT = VQ + WQ: ANOVA for D targets: Reaction Time -- Factors WO, ST and CT  -- Items"))
rt_item_vw <- aggregate(b_Dtc_vw$RT, by=list(Item=b_Dtc_vw$Item, WO=b_Dtc_vw$WO, ST=b_Dtc_vw$ST, CT=b_Dtc_vw$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*ST*CT + Error(Item/(WO*ST*CT)), data=rt_item_vw)))

print("=================================================================================================")

# Acceptability: For rating the Context+Stimulus combination acceptable or otherwise
# Only those trials where the comprehensions task was correctly performed are considered here!!!
# (1 was Unacceptable; 2 was Acceptable, but to avoid confusion, 
# (F1-1)x100 or F1*100 - 100 will give the percentage accuracy,
# and this is what will be reported.  
# Eg:- Acceptability = 1 => (1-1)x100 or 1*100 -100 = 0%; Acceptability = 2 => (2-1)x100 or 2*100 - 100 = 100%.
# Doing this already here DOES NOT change the anova results.
# What changes is just the mean and SD (both multiplied by 100!), which then 
# can be directly reported!
#First calculate a mean across conditions for each subject.
f1_per_subj <- aggregate(b_Dtc$Acceptability*100-100,list(Subj=b_Dtc$SNo, WO=b_Dtc$WO, ST=b_Dtc$ST, CT=b_Dtc$CT),mean)

#Then calculate the mean and SD of these means.
print("Mean: Acceptability (for correct answers) across subjects by condition")
print(f1_m <- aggregate(f1_per_subj$x,list(WO=f1_per_subj$WO, ST=f1_per_subj$ST, CT=f1_per_subj$CT),mean))
print("################################")

print("SD: Acceptability (for correct answers) across subjects by condition")
print(f1_sd <- aggregate(f1_per_subj$x,list(WO=f1_per_subj$WO, ST=f1_per_subj$ST, CT=f1_per_subj$CT),sd))
print("################################")

print("########### SUBJECTS ANALYSIS, Acceptability ##########")
print(paste("Repeated Measures ANOVA for D targets: Acceptability -- Factors WO, ST and CT  -- Subjects"))
f1_subj <- aggregate(b_Dtc$Acceptability*100-100, by=list(Subj=b_Dtc$SNo, WO=b_Dtc$WO, ST=b_Dtc$ST, CT=b_Dtc$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*ST*CT + Error(Subj/(WO*ST*CT)), data=f1_subj)))

# Resolve for CT alone!
print(paste("########## CT Pairwise Comparison: Acceptability by Subjects #########"))
print(paste("### CT = CQ + NQ: ANOVA for D targets: Acceptability -- Factors WO, ST and CT  -- Subjects"))
f1_subj_cn <- aggregate(b_Dtc_cn$Acceptability*100-100, by=list(Subj=b_Dtc_cn$SNo, WO=b_Dtc_cn$WO, ST=b_Dtc_cn$ST, CT=b_Dtc_cn$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*ST*CT + Error(Subj/(WO*ST*CT)), data=f1_subj_cn)))

print(paste("### CT = CQ + VQ: ANOVA for D targets: Acceptability -- Factors WO, ST and CT  -- Subjects"))
f1_subj_cv <- aggregate(b_Dtc_cv$Acceptability*100-100, by=list(Subj=b_Dtc_cv$SNo, WO=b_Dtc_cv$WO, ST=b_Dtc_cv$ST, CT=b_Dtc_cv$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*ST*CT + Error(Subj/(WO*ST*CT)), data=f1_subj_cv)))

print(paste("### CT = VQ + NQ: ANOVA for D targets: Acceptability -- Factors WO, ST and CT  -- Subjects"))
f1_subj_vn <- aggregate(b_Dtc_vn$Acceptability*100-100, by=list(Subj=b_Dtc_vn$SNo, WO=b_Dtc_vn$WO, ST=b_Dtc_vn$ST, CT=b_Dtc_vn$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*ST*CT + Error(Subj/(WO*ST*CT)), data=f1_subj_vn)))

print(paste("### CT = CQ + WQ: ANOVA for D targets: Acceptability -- Factors WO, ST and CT  -- Subjects"))
f1_subj_cw <- aggregate(b_Dtc_cw$Acceptability*100-100, by=list(Subj=b_Dtc_cw$SNo, WO=b_Dtc_cw$WO, ST=b_Dtc_cw$ST, CT=b_Dtc_cw$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*ST*CT + Error(Subj/(WO*ST*CT)), data=f1_subj_cw)))

print(paste("### CT = WQ + NQ: ANOVA for D targets: Acceptability -- Factors WO, ST and CT  -- Subjects"))
f1_subj_wn <- aggregate(b_Dtc_wn$Acceptability*100-100, by=list(Subj=b_Dtc_wn$SNo, WO=b_Dtc_wn$WO, ST=b_Dtc_wn$ST, CT=b_Dtc_wn$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*ST*CT + Error(Subj/(WO*ST*CT)), data=f1_subj_wn)))

print(paste("### CT = VQ + WQ: ANOVA for D targets: Acceptability -- Factors WO, ST and CT  -- Subjects"))
f1_subj_vw <- aggregate(b_Dtc_vw$Acceptability*100-100, by=list(Subj=b_Dtc_vw$SNo, WO=b_Dtc_vw$WO, ST=b_Dtc_vw$ST, CT=b_Dtc_vw$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*ST*CT + Error(Subj/(WO*ST*CT)), data=f1_subj_vw)))


#Resolve for ST, and then, for CT pairwise!  Because the interaction ST x CT is significant!
b_DtcDS <- b_Dtc[b_Dtc$ST=="S",]
b_DtcDI <- b_Dtc[b_Dtc$ST=="I",]

# Six possible context pairs possible: cn, cv, cw, nv, nw, vw!
# For DS
b_DtcDS_1 <- b_DtcDS[b_DtcDS$CT!="VQ",]
b_DtcDS_cn <- b_DtcDS_1[b_DtcDS_1$CT!="WQ",]
b_DtcDS_cw <- b_DtcDS_1[b_DtcDS_1$CT!="NQ",]
b_DtcDS_wn <- b_DtcDS_1[b_DtcDS_1$CT!="CQ",]

b_DtcDS_2 <- b_DtcDS[b_DtcDS$CT!="WQ",]
b_DtcDS_vn <- b_DtcDS_2[b_DtcDS_2$CT!="CQ",]
b_DtcDS_cv <- b_DtcDS_2[b_DtcDS_2$CT!="NQ",]

b_DtcDS_3 <- b_DtcDS[b_DtcDS$CT!="CQ",]
b_DtcDS_vw <- b_DtcDS_3[b_DtcDS_3$CT!="NQ",]

# For DI
b_DtcDI_1 <- b_DtcDI[b_DtcDI$CT!="VQ",]
b_DtcDI_cn <- b_DtcDI_1[b_DtcDI_1$CT!="WQ",]
b_DtcDI_cw <- b_DtcDI_1[b_DtcDI_1$CT!="NQ",]
b_DtcDI_wn <- b_DtcDI_1[b_DtcDI_1$CT!="CQ",]

b_DtcDI_2 <- b_DtcDI[b_DtcDI$CT!="WQ",]
b_DtcDI_vn <- b_DtcDI_2[b_DtcDI_2$CT!="CQ",]
b_DtcDI_cv <- b_DtcDI_2[b_DtcDI_2$CT!="NQ",]

b_DtcDI_3 <- b_DtcDI[b_DtcDI$CT!="CQ",]
b_DtcDI_vw <- b_DtcDI_3[b_DtcDI_3$CT!="NQ",]

print("############## Acceptability Subjects Analysis, ST = DS ##################")
f1_subj_DS <- aggregate(b_DtcDS$Acceptability*100-100, by=list(Subj=b_DtcDS$SNo, WO=b_DtcDS$WO, ST=b_DtcDS$ST, CT=b_DtcDS$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*CT + Error(Subj/(WO*CT)), data=f1_subj_DS)))

print(paste("########## Acceptability by Subjects, ST = DS, CT Pairwise Comparison #########"))
print(paste("### ST = DS; CT = CQ + NQ: ANOVA for D targets: Acceptability -- Factors WO, ST and CT  -- Subjects"))
f1_subjDS_cn <- aggregate(b_DtcDS_cn$Acceptability*100-100, by=list(Subj=b_DtcDS_cn$SNo, WO=b_DtcDS_cn$WO, CT=b_DtcDS_cn$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*CT + Error(Subj/(WO*CT)), data=f1_subjDS_cn)))

print(paste("### ST = DS; CT = CQ + VQ: ANOVA for D targets: Acceptability -- Factors WO, ST and CT  -- Subjects"))
f1_subjDS_cv <- aggregate(b_DtcDS_cv$Acceptability*100-100, by=list(Subj=b_DtcDS_cv$SNo, WO=b_DtcDS_cv$WO, CT=b_DtcDS_cv$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*CT + Error(Subj/(WO*CT)), data=f1_subjDS_cv)))

print(paste("### ST = DS; CT = VQ + NQ: ANOVA for D targets: Acceptability -- Factors WO, ST and CT  -- Subjects"))
f1_subjDS_vn <- aggregate(b_DtcDS_vn$Acceptability*100-100, by=list(Subj=b_DtcDS_vn$SNo, WO=b_DtcDS_vn$WO, CT=b_DtcDS_vn$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*CT + Error(Subj/(WO*CT)), data=f1_subjDS_vn)))

print(paste("### ST = DS; CT = CQ + WQ: ANOVA for D targets: Acceptability -- Factors WO, ST and CT  -- Subjects"))
f1_subjDS_cw <- aggregate(b_DtcDS_cw$Acceptability*100-100, by=list(Subj=b_DtcDS_cw$SNo, WO=b_DtcDS_cw$WO, CT=b_DtcDS_cw$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*CT + Error(Subj/(WO*CT)), data=f1_subjDS_cw)))

print(paste("### ST = DS; CT = WQ + NQ: ANOVA for D targets: Acceptability -- Factors WO, ST and CT  -- Subjects"))
f1_subjDS_wn <- aggregate(b_DtcDS_wn$Acceptability*100-100, by=list(Subj=b_DtcDS_wn$SNo, WO=b_DtcDS_wn$WO, CT=b_DtcDS_wn$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*CT + Error(Subj/(WO*CT)), data=f1_subjDS_wn)))

print(paste("### ST = DS; CT = VQ + WQ: ANOVA for D targets: Acceptability -- Factors WO, ST and CT  -- Subjects"))
f1_subjDS_vw <- aggregate(b_DtcDS_vw$Acceptability*100-100, by=list(Subj=b_DtcDS_vw$SNo, WO=b_DtcDS_vw$WO, CT=b_DtcDS_vw$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*CT + Error(Subj/(WO*CT)), data=f1_subjDS_vw)))

print("############## Acceptability Subjects Analysis, ST = DI ##################")
f1_subj_DI <- aggregate(b_DtcDI$Acceptability*100-100, by=list(Subj=b_DtcDI$SNo, WO=b_DtcDI$WO, CT=b_DtcDI$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*CT + Error(Subj/(WO*CT)), data=f1_subj_DI)))

print(paste("########## Acceptability by Subjects, ST = DI, CT Pairwise Comparison #########"))
print(paste("### ST = DI; CT = CQ + NQ: ANOVA for D targets: Acceptability -- Factors WO, ST and CT  -- Subjects"))
f1_subjDI_cn <- aggregate(b_DtcDI_cn$Acceptability*100-100, by=list(Subj=b_DtcDI_cn$SNo, WO=b_DtcDI_cn$WO, CT=b_DtcDI_cn$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*CT + Error(Subj/(WO*CT)), data=f1_subjDI_cn)))

print(paste("### ST = DI; CT = CQ + VQ: ANOVA for D targets: Acceptability -- Factors WO, ST and CT  -- Subjects"))
f1_subjDI_cv <- aggregate(b_DtcDI_cv$Acceptability*100-100, by=list(Subj=b_DtcDI_cv$SNo, WO=b_DtcDI_cv$WO,  CT=b_DtcDI_cv$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*CT + Error(Subj/(WO*CT)), data=f1_subjDI_cv)))

print(paste("### ST = DI; CT = VQ + NQ: ANOVA for D targets: Acceptability -- Factors WO, ST and CT  -- Subjects"))
f1_subjDI_vn <- aggregate(b_DtcDI_vn$Acceptability*100-100, by=list(Subj=b_DtcDI_vn$SNo, WO=b_DtcDI_vn$WO, CT=b_DtcDI_vn$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*CT + Error(Subj/(WO*CT)), data=f1_subjDI_vn)))

print(paste("### ST = DI; CT = CQ + WQ: ANOVA for D targets: Acceptability -- Factors WO, ST and CT  -- Subjects"))
f1_subjDI_cw <- aggregate(b_DtcDI_cw$Acceptability*100-100, by=list(Subj=b_DtcDI_cw$SNo, WO=b_DtcDI_cw$WO, CT=b_DtcDI_cw$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*CT + Error(Subj/(WO*CT)), data=f1_subjDI_cw)))

print(paste("### ST = DI; CT = WQ + NQ: ANOVA for D targets: Acceptability -- Factors WO, ST and CT  -- Subjects"))
f1_subjDI_wn <- aggregate(b_DtcDI_wn$Acceptability*100-100, by=list(Subj=b_DtcDI_wn$SNo, WO=b_DtcDI_wn$WO, CT=b_DtcDI_wn$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*CT + Error(Subj/(WO*CT)), data=f1_subjDI_wn)))

print(paste("### ST = DI; CT = VQ + WQ: ANOVA for D targets: Acceptability -- Factors WO, ST and CT  -- Subjects"))
f1_subjDI_vw <- aggregate(b_DtcDI_vw$Acceptability*100-100, by=list(Subj=b_DtcDI_vw$SNo, WO=b_DtcDI_vw$WO, ST=b_DtcDI_vw$ST, CT=b_DtcDI_vw$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*CT + Error(Subj/(WO*CT)), data=f1_subjDI_vw)))

print("===============================================")
print("########### ITEMS ANALYSIS, Acceptability ##########")
print(paste("Repeated Measures ANOVA for D targets: Acceptability -- Factors WO, ST and CT -- Items"))
f1_item <- aggregate(b_Dtc$Acceptability*100-100, by=list(Item=b_Dtc$Item, WO=b_Dtc$WO, ST=b_Dtc$ST, CT=b_Dtc$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*ST*CT + Error(Item/(WO*ST*CT)), data=f1_item)))

# Resolve for CT alone!
print(paste("########## CT Pairwise Comparison: Acceptability by Items #########"))
print(paste("### CT = CQ + NQ: ANOVA for D targets: Acceptability -- Factors WO, ST and CT  -- Items"))
f1_item_cn <- aggregate(b_Dtc_cn$Acceptability*100-100, by=list(Item=b_Dtc_cn$Item, WO=b_Dtc_cn$WO, ST=b_Dtc_cn$ST, CT=b_Dtc_cn$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*ST*CT + Error(Item/(WO*ST*CT)), data=f1_item_cn)))

print(paste("### CT = CQ + VQ: ANOVA for D targets: Acceptability -- Factors WO, ST and CT  -- Items"))
f1_item_cv <- aggregate(b_Dtc_cv$Acceptability*100-100, by=list(Item=b_Dtc_cv$Item, WO=b_Dtc_cv$WO, ST=b_Dtc_cv$ST, CT=b_Dtc_cv$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*ST*CT + Error(Item/(WO*ST*CT)), data=f1_item_cv)))

print(paste("### CT = VQ + NQ: ANOVA for D targets: Acceptability -- Factors WO, ST and CT  -- Items"))
f1_item_vn <- aggregate(b_Dtc_vn$Acceptability*100-100, by=list(Item=b_Dtc_vn$Item, WO=b_Dtc_vn$WO, ST=b_Dtc_vn$ST, CT=b_Dtc_vn$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*ST*CT + Error(Item/(WO*ST*CT)), data=f1_item_vn)))

print(paste("### CT = CQ + WQ: ANOVA for D targets: Acceptability -- Factors WO, ST and CT  -- Items"))
f1_item_cw <- aggregate(b_Dtc_cw$Acceptability*100-100, by=list(Item=b_Dtc_cw$Item, WO=b_Dtc_cw$WO, ST=b_Dtc_cw$ST, CT=b_Dtc_cw$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*ST*CT + Error(Item/(WO*ST*CT)), data=f1_item_cw)))

print(paste("### CT = WQ + NQ: ANOVA for D targets: Acceptability -- Factors WO, ST and CT  -- Items"))
f1_item_wn <- aggregate(b_Dtc_wn$Acceptability*100-100, by=list(Item=b_Dtc_wn$Item, WO=b_Dtc_wn$WO, ST=b_Dtc_wn$ST, CT=b_Dtc_wn$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*ST*CT + Error(Item/(WO*ST*CT)), data=f1_item_wn)))

print(paste("### CT = VQ + WQ: ANOVA for D targets: Acceptability -- Factors WO, ST and CT  -- Items"))
f1_item_vw <- aggregate(b_Dtc_vw$Acceptability*100-100, by=list(Item=b_Dtc_vw$Item, WO=b_Dtc_vw$WO, ST=b_Dtc_vw$ST, CT=b_Dtc_vw$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*ST*CT + Error(Item/(WO*ST*CT)), data=f1_item_vw)))

print("############## Acceptability Items Analysis, ST = DS ##################")
f1_item_DS <- aggregate(b_DtcDS$Acceptability*100-100, by=list(Item=b_DtcDS$Item, WO=b_DtcDS$WO, CT=b_DtcDS$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*CT + Error(Item/(WO*CT)), data=f1_item_DS)))


print(paste("########## Acceptability by Items, ST = DS, CT Pairwise Comparison #########"))
print(paste("### ST = DS; CT = CQ + NQ: ANOVA for D targets: Acceptability -- Factors WO, ST and CT  -- Items"))
f1_itemDS_cn <- aggregate(b_DtcDS_cn$Acceptability*100-100, by=list(Item=b_DtcDS_cn$Item, WO=b_DtcDS_cn$WO, CT=b_DtcDS_cn$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*CT + Error(Item/(WO*CT)), data=f1_itemDS_cn)))

print(paste("### ST = DS; CT = CQ + VQ: ANOVA for D targets: Acceptability -- Factors WO, ST and CT  -- Items"))
f1_itemDS_cv <- aggregate(b_DtcDS_cv$Acceptability*100-100, by=list(Item=b_DtcDS_cv$Item, WO=b_DtcDS_cv$WO, CT=b_DtcDS_cv$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*CT + Error(Item/(WO*CT)), data=f1_itemDS_cv)))

print(paste("### ST = DS; CT = VQ + NQ: ANOVA for D targets: Acceptability -- Factors WO, ST and CT  -- Items"))
f1_itemDS_vn <- aggregate(b_DtcDS_vn$Acceptability*100-100, by=list(Item=b_DtcDS_vn$Item, WO=b_DtcDS_vn$WO, CT=b_DtcDS_vn$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*CT + Error(Item/(WO*CT)), data=f1_itemDS_vn)))

print(paste("### ST = DS; CT = CQ + WQ: ANOVA for D targets: Acceptability -- Factors WO, ST and CT  -- Items"))
f1_itemDS_cw <- aggregate(b_DtcDS_cw$Acceptability*100-100, by=list(Item=b_DtcDS_cw$Item, WO=b_DtcDS_cw$WO, CT=b_DtcDS_cw$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*CT + Error(Item/(WO*CT)), data=f1_itemDS_cw)))

print(paste("### ST = DS; CT = WQ + NQ: ANOVA for D targets: Acceptability -- Factors WO, ST and CT  -- Items"))
f1_itemDS_wn <- aggregate(b_DtcDS_wn$Acceptability*100-100, by=list(Item=b_DtcDS_wn$Item, WO=b_DtcDS_wn$WO, CT=b_DtcDS_wn$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*CT + Error(Item/(WO*CT)), data=f1_itemDS_wn)))

print(paste("### ST = DS; CT = VQ + WQ: ANOVA for D targets: Acceptability -- Factors WO, ST and CT  -- Items"))
f1_itemDS_vw <- aggregate(b_DtcDS_vw$Acceptability*100-100, by=list(Item=b_DtcDS_vw$Item, WO=b_DtcDS_vw$WO, CT=b_DtcDS_vw$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*CT + Error(Item/(WO*CT)), data=f1_itemDS_vw)))

print("############## Acceptability Items Analysis, ST = DI ##################")
f1_item_DI <- aggregate(b_DtcDI$Acceptability*100-100, by=list(Item=b_DtcDI$Item, WO=b_DtcDI$WO, CT=b_DtcDI$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*CT + Error(Item/(WO*CT)), data=f1_item_DI)))


print(paste("########## Acceptability by Items, ST = DI, CT Pairwise Comparison #########"))
print(paste("### ST = DI; CT = CQ + NQ: ANOVA for D targets: Acceptability -- Factors WO, ST and CT  -- Items"))
f1_itemDI_cn <- aggregate(b_DtcDI_cn$Acceptability*100-100, by=list(Item=b_DtcDI_cn$Item, WO=b_DtcDI_cn$WO, CT=b_DtcDI_cn$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*CT + Error(Item/(WO*CT)), data=f1_itemDI_cn)))

print(paste("### ST = DI; CT = CQ + VQ: ANOVA for D targets: Acceptability -- Factors WO, ST and CT  -- Items"))
f1_itemDI_cv <- aggregate(b_DtcDI_cv$Acceptability*100-100, by=list(Item=b_DtcDI_cv$Item, WO=b_DtcDI_cv$WO, CT=b_DtcDI_cv$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*CT + Error(Item/(WO*CT)), data=f1_itemDI_cv)))

print(paste("### ST = DI; CT = VQ + NQ: ANOVA for D targets: Acceptability -- Factors WO, ST and CT  -- Items"))
f1_itemDI_vn <- aggregate(b_DtcDI_vn$Acceptability*100-100, by=list(Item=b_DtcDI_vn$Item, WO=b_DtcDI_vn$WO, CT=b_DtcDI_vn$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*CT + Error(Item/(WO*CT)), data=f1_itemDI_vn)))

print(paste("### ST = DI; CT = CQ + WQ: ANOVA for D targets: Acceptability -- Factors WO, ST and CT  -- Items"))
f1_itemDI_cw <- aggregate(b_DtcDI_cw$Acceptability*100-100, by=list(Item=b_DtcDI_cw$Item, WO=b_DtcDI_cw$WO, CT=b_DtcDI_cw$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*CT + Error(Item/(WO*CT)), data=f1_itemDI_cw)))

print(paste("### ST = DI; CT = WQ + NQ: ANOVA for D targets: Acceptability -- Factors WO, ST and CT  -- Items"))
f1_itemDI_wn <- aggregate(b_DtcDI_wn$Acceptability*100-100, by=list(Item=b_DtcDI_wn$Item, WO=b_DtcDI_wn$WO, CT=b_DtcDI_wn$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*CT + Error(Item/(WO*CT)), data=f1_itemDI_wn)))

print(paste("### ST = DI; CT = VQ + WQ: ANOVA for D targets: Acceptability -- Factors WO, ST and CT  -- Items"))
f1_itemDI_vw <- aggregate(b_DtcDI_vw$Acceptability*100-100, by=list(Item=b_DtcDI_vw$Item, WO=b_DtcDI_vw$WO, CT=b_DtcDI_vw$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*CT + Error(Item/(WO*CT)), data=f1_itemDI_vw)))

print("=================================================================================================")


# Answering Accuracy: For answering correctly or otherwise
# (1 was Incorrect; 2 was Correct, but to avoid confusion, 
# (F1-1)x100 or F1*100 - 100 will give the percentage accuracy,
# and this is what will be reported.  
# Eg:- Accuracy = 1 => (1-1)x100 or 1*100 -100 = 0%; Accuracy = 2 => (2-1)x100 or 2*100 - 100 = 100%.
# Doing this already here DOES NOT change the anova results.
# What changes is just the mean and SD (both multiplied by 100!), which then 
# can be directly reported!

#First calculate a mean across conditions for each subject.
f2_per_subj <- aggregate(b_Dt$Answer*100-100,list(Subj=b_Dt$SNo, WO=b_Dt$WO, ST=b_Dt$ST, CT=b_Dt$CT),mean)

#Then calculate the mean and SD of these means.
print("Mean: Answer accuracy(correct/incorrect/TO) across subjects by condition")
print(f2_m <- aggregate(f2_per_subj$x,list(WO=f2_per_subj$WO, ST=f2_per_subj$ST, CT=f2_per_subj$CT),mean))
print("################################")

print("SD: Answer accuracy(correct/incorrect/TO) across subjects by condition")
print(f2_sd <- aggregate(f2_per_subj$x,list(WO=f2_per_subj$WO, ST=f2_per_subj$ST, CT=f2_per_subj$CT),sd))
print("################################")

print("########### SUBJECTS ANALYSIS, ACCURACY ##########")
print(paste("Repeated Measures ANOVA for D targets: Accuracy -- Factors WO, ST and CT  -- Subjects"))
f2_subj <- aggregate(b_Dt$Answer*100-100, by=list(Subj=b_Dt$SNo, WO=b_Dt$WO, ST=b_Dt$ST, CT=b_Dt$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*ST*CT + Error(Subj/(WO*ST*CT)), data=f2_subj)))


#Resolve for Word-order
b_Dtn <- b_Dt[b_Dt$WO=="DA",]
b_Dts <- b_Dt[b_Dt$WO=="AD",]

print("############## Accuracy Subjects Analysis, WO = DA ##################")
f2_subj_n <- aggregate(b_Dtn$Answer*100-100, by=list(Subj=b_Dtn$SNo, WO=b_Dtn$WO, ST=b_Dtn$ST, CT=b_Dtn$CT), mean, NA.rm=T)
print(summary (aov(x ~ ST*CT + Error(Subj/(ST*CT)), data=f2_subj_n)))

print("############## Accuracy Subjects Analysis, WO = AD ##################")
f2_subj_s <- aggregate(b_Dts$Answer*100-100, by=list(Subj=b_Dts$SNo, WO=b_Dts$WO, ST=b_Dts$ST, CT=b_Dts$CT), mean, NA.rm=T)
print(summary (aov(x ~ ST*CT + Error(Subj/(ST*CT)), data=f2_subj_s)))


print("===============================================")
print("########### ITEMS ANALYSIS, ACCURACY ##########")
print(paste("Repeated Measures ANOVA for D targets: Accuracy -- Factors WO, ST and CT -- Items"))
f2_item <- aggregate(b_Dt$Answer*100-100, by=list(Item=b_Dt$Item, WO=b_Dt$WO, ST=b_Dt$ST, CT=b_Dt$CT), mean, NA.rm=T)
print(summary (aov(x ~ WO*ST*CT + Error(Item/(WO*ST*CT)), data=f2_item)))

print("############## Accuracy Items Analysis, WO = DA ##################")
f2_item_n <- aggregate(b_Dtn$Answer*100-100, by=list(Item=b_Dtn$Item, WO=b_Dtn$WO, ST=b_Dtn$ST, CT=b_Dtn$CT), mean, NA.rm=T)
print(summary (aov(x ~ ST*CT + Error(Item/(ST*CT)), data=f2_item_n)))

print("############## Accuracy Items Analysis, WO = AD ##################")
f2_item_s <- aggregate(b_Dts$Answer*100-100, by=list(Item=b_Dts$Item, WO=b_Dts$WO, ST=b_Dts$ST, CT=b_Dts$CT), mean, NA.rm=T)
print(summary (aov(x ~ ST*CT + Error(Item/(ST*CT)), data=f2_item_s)))

sink()