## Load packages
library(languageR)
library(lattice)
library(Matrix)
library(lme4)
library(lmerTest)
library(car)
library(ggplot2)

## Clear all
rm(list=ls())
objects()

## Read data
data1 <- read.table(file = "c:/.../all_data.csv", header = TRUE, sep = ";", dec = ",") 
data1$freq <- ordered(data1$freq, levels = c("AF", "MF","BF"))
data1$subject <- as.factor(data1$subject)
data1$acc <- as.logical(data1$acc)
data1$rtlog <- log(data1$rt)
data1$rtinv <- 1/data1$rt * -1000
data1$dummy <- 1
str(data1)
head(data1)											
(nrow.data1 <- nrow(data1))

## Pseudowords filter
data2 <- droplevels(subset(data1, type == "word"))
(nrow.data2 <- nrow(data2))
(outlier.data2 <- 100 * (1 - nrow.data2/nrow.data1))
(nrow.data1 - nrow.data2)

## ACC filter
data3 <- droplevels(subset(data2, acc == TRUE))
(nrow.data3 <- nrow(data3))
(outlier.data3 <- 100 * (1 - nrow.data3/nrow.data2))
(nrow.data2 - nrow.data3)

## Subject ACC
data.subject.acc <- xtabs(~subject + acc, data = data2)    		
(data.subject.acc.prop <- 100 * prop.table(data.subject.acc, 1))

data.subject.error.prop <- data.frame(data.subject.acc.prop)
data.subject.error.prop <- droplevels(subset(data.subject.error.prop, acc == FALSE))  			
par(mfrow = c(1,1))
boxplot(data.subject.error.prop$Freq, main = "Subject Error", xlab = "Subjects", ylab = "Error(%)")
text(data.subject.error.prop$acc, data.subject.error.prop$Freq, as.character(data.subject.error.prop$subject), pos = 4)

## Subject RT
boxplot(data3$rt ~ data3$subject, main = "RT vs. Subject", xlab = "Subject", ylab = "RT(ms)")
abline(h = mean(data3$rt))
abline(h = mean(data3$rt) - sd(data3$rt), lty = 3)
abline(h = mean(data3$rt) + sd(data3$rt), lty = 3)

data.subject.rt <- with(data3, aggregate(rt, list(subject, dummy), mean))
colnames(data.subject.rt) <- c("Subject", "Dummy", "MeanRT")
boxplot(data.subject.rt$MeanRT, main = "RT Subject", xlab = "Subject", ylab = "RT(ms)")
text(data.subject.rt$Dummy, data.subject.rt$MeanRT, as.character(data.subject.rt$Subject), pos = 4)

## Reject subjects
data4 <- data3
(nrow.data4 <- nrow(data4))
(outlier.data4 <- 100 * (1 - nrow.data4/nrow.data3))
(nrow.data3 - nrow.data4)
(outlier.general1 <- 100 * (1 - nrow.data4/nrow.data2))

## Item ACC
data.item.acc <- xtabs(~target + acc, data = data2)				
(data.item.acc.prop <- 100 * prop.table(data.item.acc, 1))

data.item.error.prop <- data.frame(data.item.acc.prop)
data.item.error.prop <- droplevels(subset(data.item.error.prop, acc == FALSE))
boxplot(data.item.error.prop$Freq, main = "Item Error", xlab = "Item", ylab = "Error(%)")
text(data.item.error.prop$acc, data.item.error.prop$Freq, as.character(data.item.error.prop$target), pos = 4)
abline(h = 50)

## Item RT
boxplot(data3$rt ~ data3$target, main = "RT vs. Item", xlab = "Item", ylab = "RT(ms)")
abline(h = mean(data3$rt))
abline(h = mean(data3$rt) - sd(data3$rt), lty = 3)
abline(h = mean(data3$rt) + sd(data3$rt), lty = 3)

data.item.rt <- with(data3, aggregate(rtinv, list(target, dummy), mean))  	
colnames(data.item.rt) <- c("Target", "Dummy", "MeanRT")
boxplot(data.item.rt$MeanRT, main = "RT Item", xlab = "Items", ylab = "RT(ms)")
text(data.item.rt$Dummy, data.item.rt$MeanRT, as.character(data.item.rt$Target), pos = 4, cex = 0.7)

## Reject items
data5 <- data4
droplevels(subset(data4, target != "emana"))
(nrow.data5 <- nrow(data5))
(outlier.data5 <- 100 * (1 - nrow.data5/nrow.data4))
(nrow.data4 - nrow.data5)
(outlier.general2 <- 100 * (1 - nrow.data5/nrow.data2))

## RT filter a-priori
data6 <- droplevels(subset(data5, rt > 300 & rt < 1800))
(nrow.data6 <- nrow(data6))
(outlier.data6 <- 100 * (1 - nrow.data6/nrow.data5))
(nrow.data5 - nrow.data6)
(outlier.general4 <- 100 * (1 - nrow.data6/nrow.data2))

## RT distribution
data <- data6
par(mfrow = c(3,3))
plot(density(data$rt), main = "RT Density", xlab = "RT(ms)", ylab = "Density")
plot(sort(data$rt), main = "RT Ordered", ylab = "RT(ms)")
boxplot(data$rt, main = "RT Boxplot", ylab = "RT(ms)")
plot(density(data$rtlog), main = "RTlog Density", xlab = "RTlog", ylab = "Density")
plot(sort(data$rtlog), main = "RTlog Ordered", ylab = "RTlog")
boxplot(data$rtlog, main = "RTlog Boxplot", ylab = "RTlog")
plot(density(data$rtinv), main = "RTinv Density", xlab = "RTinv", ylab = "Density")
plot(sort(data$rtinv), main = "RTinv Ordered", ylab = "RTinv")
boxplot(data$rtinv, main = "RTinv Boxplot", ylab = "RTinv")

par(mfrow = c(2,2))
qqnorm(rnorm(length(data$rt), 10, 3))
qqline(rnorm(length(data$rt), 10, 3))
qqnorm(data$rt, main = "RT Q-Q Plot")
qqline(data$rt, main = "RT Q-Q Plot")
qqnorm(data$rtlog, main = "RTlog Q-Q Plot")
qqline(data$rtlog, main = "RTlog Q-Q Plot")
qqnorm(data$rtinv, main = "RTinv Q-Q Plot")
qqline(data$rtinv, main = "RTinv Q-Q Plot")

## Kolmogorov test
ks.test(jitter(data$rt), "pnorm", mean(data$rt), sd(data$rt))
ks.test(jitter(data$rtlog), "pnorm", mean(data$rtlog), sd(data$rtlog))
ks.test(jitter(data$rtinv), "pnorm", mean(data$rtinv), sd(data$rtinv))

## RT vs. type
round(tapply(data1$rt, data1$type, mean))    
round(tapply(data1$rt, data1$type, sd))

## RT vs. freq
round(with(data, tapply(rt, list(freq), mean)))    
round(with(data, tapply(rt, list(freq), sd)))    

## RT vs. neigh
round(with(data, tapply(rt, list(neigh), mean)))    
round(with(data, tapply(rt, list(neigh), sd))) 

## RT vs. cat
round(with(data, tapply(rt, list(cat), mean)))    
round(with(data, tapply(rt, list(cat), sd))) 

## RT vs. letters
round(with(data, tapply(rt, list(letters), mean)))    
round(with(data, tapply(rt, list(cat), sd))) 

## RTs and SDs
round(with(data, tapply(rt, list(freq, neigh, cat), mean)))    
round(with(data, tapply(rt, list(freq, neigh, cat), sd)))    

## RT T-test
rt.ttest <- with(data1, aggregate(rtinv, list(subject, type), mean))
colnames(rt.ttest) <- c("subject", "type", "rtinv")
rt.ttest <- droplevels(subset(rt.ttest, subject != 20 & subject != 35))
t.test(rtinv ~ type, rt.ttest)

## RT ANOVA
rt.sub <- with(data, aggregate(rtinv, list(subject, cat, freq, neigh), mean))
colnames(rt.sub) <- c("subject", "cat", "freq", "neigh", "rtinv")
summary(aov(rtinv ~ cat * freq * neigh + Error(subject), rt.sub))

freq.ttest <- with(data, aggregate(rtinv, list(subject, freq), mean))
colnames(freq.ttest) <- c("subject", "freq", "rtinv")
freq.ttest <- droplevels(subset(freq.ttest, freq != "HF"))
t.test(rtinv ~ freq, freq.ttest)

cat.ttest <- with(data, aggregate(rtinv, list(subject, freq, cat), mean))
colnames(cat.ttest) <- c("subject", "freq", "cat", "rtinv")
cat.ttest <- droplevels(subset(cat.ttest, freq == "HF"))
t.test(rtinv ~ cat, cat.ttest)

## RT linear regression
rt.lm <- with(data, aggregate(rtinv, list(target, letters, syllables), mean))
colnames(rt.lm) <- c("target", "letters", "syllables", "rtinv")
anova(lm(rtinv ~ letters + syllables, rt.lm))
summary(lm(rtinv ~ letters + syllables, rt.lm))

## ACC type
acc.ttest <- data.frame(xtabs(~subject + type + acc, data = data1))      	
round(data.acc <- 100 * data.acc/sum(data.acc), 2)

## ACCs
data.acc <- xtabs(~neigh + freq + cat + acc, data = data2)      	
round(data.acc <- 100 * data.acc/sum(data.acc), 2)

## ACC T-test
acc.ttest <- data.frame(xtabs(~subject + type + acc, data = data1))      	
acc.ttest <- droplevels(subset(acc.ttest, acc == FALSE))
t.test(Freq ~ type, acc.ttest)

bwplot(data$rt ~ data$neigh | data$freq, main = "A: RT x Vizinhança x Frequência", xlab = "Vizinhança", ylab = "RT(ms)")
bwplot(data$rt ~ data$cat | data$freq, main = "B: RT x Categoria Gram. x Frequência", xlab = "Categoria Gram.", ylab = "RT(ms)")
boxplot(data$rt ~ data$letters, main = "A: RT x Letras", xlab = "Letras", ylab = "RT(ms)")
boxplot(data$rt ~ data$syllables, main = "B: RT x Sílabas", xlab = "Sílabas", ylab = "RT(ms)")