Date <- format(Sys.time(), "%Y%m%d")

source("R-Fonctions.R")

Detacher("DDD")
DDDpro <- read.table("raw-results-pro.txt", quote="", header=TRUE)
	DDDpro$PROBIAS <- 1
DDDanti <- read.table("raw-results-anti.txt", quote="", header=TRUE)
	DDDanti$PROBIAS <- 0
DDD <- rbind(DDDpro, DDDanti)




NB <- gray(c(.9,.5))

DDD$subject <- factor(DDD$subject)

DDD$TYPE <- sub("(.*)==(.*)==(.*)==(.*)", "\\1", DDD$Gname)
DDD$SCENARIO <- sub("(.*)==(.*)==(.*)==(.*)", "\\2", DDD$Gname)
DDD$POS <- sub("(.*)==(.*)==(.*)==(.*)", "\\3", DDD$Gname)
DDD$HI <- sub("(.*)==(.*)==(.*)==(.*)", "\\4", DDD$Gname)

DDD$TYPE[DDD$SCENARIO==5 & DDD$PRO==1] <- "KB+"
DDD$TYPE[DDD$SCENARIO==5 & DDD$PRO==0] <- "KB-"
DDD$HI[DDD$SCENARIO==5 & DDD$Iname==0] <- 1 # BADLY CODED
DDD$HI[DDD$SCENARIO==5 & DDD$Iname==1] <- 0
DDD$POS[DDD$SCENARIO==5] <- 1


# OTHER CODE ERRORS
TEMP <- DDD$POS==DDD$HI
DDD$POS[DDD$SCENARIO==1 & DDD$TYPE=="M" & TEMP] <- 1
DDD$POS[DDD$SCENARIO==1 & DDD$TYPE=="M" & !TEMP] <- 0




DDD$DATE <- sub("(.*)_(.*)", "\\1", DDD$date_time)
DDD$TIME <- sub("(.*)_(.*)", "\\2", DDD$date_time)


DDD$BINARYRESPONSE <- DDD$Ireponse>.5


DD <- subset(DDD, Ireponse >=0 & !(TYPE %in% c("ENT5", "ENT6")))


Mrt <- tapply(DD$Irt, DD$subject, mean)
SDrt <- tapply(DD$Irt, DD$subject, sd)

DD$Irtz <- (DD$Irt - Mrt[DD$subject]) / SDrt[DD$subject]



infos.subject <- data.frame(
	tapply(as.character(DD$DATE), DD$subject, max),
	tapply(as.character(DD$TIME), DD$subject, max),
	tapply(as.character(DD$ip), DD$subject, max),
	tapply(as.character(DD$Qreponse_2_sex), DD$subject, max),
	tapply(as.character(DD$Qreponse_1_age), DD$subject, max),
	tapply(as.character(DD$Qreponse_3_languages), DD$subject, max)
	)
	
colnames(infos.subject) <- c("date", "time", "ip", "sex", "year", "languages")

#infos.subject










DD <- subset(DDD, Ireponse >-1 & !(TYPE %in% c("E")) & Qreponse_3_languages %in% c("American_E", "Other_Engl"))
DDag <- aggregate(DD$Ireponse, by=list(Type=DD$TYPE, Subject=DD$subject, Hi=DD$HI, Pos=DD$POS), FUN=mean)





DDagC <- subset(DDag, Type=="C")

round(tapply(100*DDagC$x, list(paste("POS:", DDagC$Pos, " HI:", DDagC$Hi, sep=""), factor(DDagC$Subject)), mean))
round(tapply(100*DDagC$x, list(factor(DDagC$Subject), paste(DDagC$Pos, "-", DDagC$Hi, sep="")), mean))

BarErr2(100*DDagC$x, paste("POS:", DDagC$Pos, " HI:", DDagC$Hi, sep=""), factor(DDagC$Subject), beside=T, leg=T, ylim=c(0,100), col=c("red", "green", "blue", "black" ))




BarErr2(100*DDag$x, paste("POS:", DDag$Pos, " HI:", DDag$Hi, sep=""), DDag$Type, beside=T, leg=T, ylim=c(0,100))

par(mai=c(.8,.65,.2,0))
BarErr2(100*DDag$x,
	factor(DDag$Hi, levels=c(0,1), labels=c("LOW", "HIGH")),
	factor(DDag$Pos, levels=c(0,1), labels=c("NEG", "POS")),
	DDag$Type=="C",
	beside=T, leg=F, ylim=c(0,100), axes=F, col=rep(c("red", "green"), each=2))
BarErr2(100*DDag$x,
	factor(DDag$Hi, levels=c(0,1), labels=c("LOW", "HIGH")),
	factor(DDag$Pos, levels=c(0,1), labels=c("NEG", "POS")),
	DDag$Type=="C",
	beside=T, leg=F, ylim=c(0,100), axes=F, density=15, add=T, col=c("black", "red", "black", "green"))
axis(2, at=c(0:4)*25, labels=paste(c(0:4)*25, "%", sep=""), las=1)
legend("topright", density=25, fill=c("black", "white"), legend=c("LOW", "HIGH"), pt.lwd=2, inset=.05)






par(mai=c(.8,.65,.2,0))
SCENARIO <- "K"
BarErr2(100*DDag$x,
	factor(DDag$Hi, levels=c(0,1), labels=c("LOW", "HIGH")),
	factor(DDag$Pos, levels=c(0,1), labels=c("Negative", "Positive")),
	DDag$Type==SCENARIO,
	beside=T, leg=F, ylim=c(0,100), axes=F, col=rep(c("green", "red"), 2))
BarErr2(100*DDag$x,
	factor(DDag$Hi, levels=c(0,1), labels=c("", "")),
	factor(DDag$Pos, levels=c(0,1), labels=c("", "")),
	DDag$Type==SCENARIO,
	beside=T, leg=F, ylim=c(0,100), axes=F, density=15, add=T, col=c("black", "black", "green", "red"))
axis(2, at=c(0:4)*25, labels=paste(c(0:4)*25, "%", sep=""), las=1)
legend("topright", fill=c("green", "red"), legend=c("LOW", "HIGH"), pt.lwd=2, inset=.05)





DDagHIDE <- subset(DDag, paste(DDag$Pos,DDag$Hi,sep="") %in% c("10", "01"))
par(mai=c(.8,.65,.2,0))
SCENARIO <- "K"
BarErr2(100*DDagHIDE$x,
	factor(DDagHIDE$Hi, levels=c(0,1), labels=c("LOW", "HIGH")),
	factor(DDagHIDE$Pos, levels=c(0,1), labels=c("Negative", "Positive")),
	DDagHIDE$Type==SCENARIO,
	beside=T, leg=F, ylim=c(0,100), axes=F, col=rep(c("green", "red"), 2))
BarErr2(100*DDagHIDE$x,
	factor(DDagHIDE$Hi, levels=c(0,1), labels=c("", "")),
	factor(DDagHIDE$Pos, levels=c(0,1), labels=c("", "")),
	DDagHIDE$Type==SCENARIO,
	beside=T, leg=F, ylim=c(0,100), axes=F, density=15, add=T, col=c("black", "black", "green", "red"))
axis(2, at=c(0:4)*25, labels=paste(c(0:4)*25, "%", sep=""), las=1)
legend("topright", fill=c("green", "red"), legend=c("LOW", "HIGH"), pt.lwd=2, inset=.05)









BarErr2(100*DDag$x, paste("POS:", DDag$Pos, " HI:", DDag$Hi, sep=""), factor(DDag$Type, levels=c("CONT", "C", "M", "K", "KB+", "KB-"), labels=c("Controls", "Color", "Misc.", "Knowl", "BiasPro", "BiasAnti")), beside=T, leg=T, ylim=c(0,100), col=c("red", "green", "blue", "black" ), xlim=c(0,40))



BarErr2(100*DD$Ireponse, paste("POS:", DD$POS, " HI:", DD$HI, sep=""), paste(DD$TYPE, DD$SCENARIO), beside=T, leg=T, ylim=c(0,100), col=c("red", "green", "blue", "black" ), xlim=c(0,70))


BarErr2(100*DD$Ireponse,
	factor(paste("POS:", DD$POS, " HI:", DD$HI, sep=""), 
		levels=c("POS:0 HI:0", "POS:0 HI:1", "POS:1 HI:0", "POS:1 HI:1"),
		labels=c("Neg-Low", "Neg-Hi", "Pos-Low", "Pos-Hi")),
	factor(paste(DD$TYPE, DD$SCENARIO), 
		levels=c("K 1", "K 2", "K 3", "K 4", "C 1", "C 2" , "C 3", "C 4", "M 1", "M 2" ),
		labels=c("Knowl1", "Knowl2", "Knowl3", "Knowl4", "Col 1", "Col 2", "Col 3", "Col 4", "Misc 1", "Misc 2"  )), 
	beside=T, leg=F, ylim=c(0,100), col=c("green", "red" ), xlim=c(0,60), axes=F)

BarErr2(100*DD$Ireponse,
	factor(paste("POS:", DD$POS, " HI:", DD$HI, sep=""), 
		levels=c("POS:0 HI:0", "POS:0 HI:1", "POS:1 HI:0", "POS:1 HI:1"),
		labels=c("Neg-Low", "Neg-Hi", "Pos-Low", "Pos-Hi")),
	factor(paste(DD$TYPE, DD$SCENARIO), 
		levels=c("K 1", "K 2", "K 3", "K 4", "C 1", "C 2" , "C 3", "C 4", "M 1", "M 2" ),
		labels=c("Knowl1", "Knowl2", "Knowl3", "Knowl4", "Col 1", "Col 2", "Col 3", "Col 4", "Misc 1", "Misc 2"  )), 
	beside=T, leg=F, ylim=c(0,100), density=25, col=c("black", "black", "green", "red" ), xlim=c(0,60), add=T, axes=F)
axis(2, at=c(0:4)*25, labels=paste(c(0:4)*25, "%", sep=""), las=1)
legend("topright", density=25, fill=c("black", "white"), legend=c("Negative", "Positive"), pt.lwd=2, inset=.05)
legend("right", fill=c("green", "red"), legend=c("Low", "High"), pt.lwd=2, inset=.05)










op <- options(warn=(-1)) # suppress warnings

summary(aov(x~Pos*Hi + Error(Subject/(Pos*Hi)), data=subset(DDag, Type %in% c("C", "K", "M"))))
summary(aov(x~Hi+Type + Error(Subject/(Hi*Type)), data=subset(DDag, Type %in% c("K", "C", "M") & Pos==1)))
summary(aov(x~Hi*Type + Error(Subject/(Hi)), data=subset(DDag, Type %in% c("KB+", "KB-") & Pos==1)))
summary(aov(x~Pos*Hi*Type + Error(Subject/(Pos*Hi*Type)), data=subset(DDag, Type %in% c("C", "K", "M"))))


summary(aov(x~Pos*Hi + Error((Subject:Type)/(Pos*Hi)), data=subset(DDag, Type %in% c("C", "K", "M"))))


SCENARIO <- "K"
summary(aov(x~Hi + Error(Subject/Hi), data=subset(DDag, Type == SCENARIO & Pos==1)))
summary(aov(x~Hi + Error(Subject/Hi), data=subset(DDag, Type == SCENARIO & Pos==0)))
summary(aov(x~Pos*Hi + Error(Subject/(Pos*Hi)), data=subset(DDag, Type == SCENARIO)))


options(op)








# BLOCK ANALYSES

DDagB <- aggregate(DD$Ireponse, by=list(Type=DD$TYPE, Subject=DD$subject, Hi=DD$HI, Pos=DD$POS, Block=DD$Bplace), FUN=mean)

DDagBT <- subset(DDagB, Type %in% c("K", "KB+", "KB-"))
BarErr2(100*DDagBT$x, 
        factor(paste("POS:", DDagBT$Pos, " HI:", DDagBT$Hi, sep="")), 
      	factor(DDagBT$Block), beside=T, leg=T, ylim=c(0,100), col=c("red", "green", "blue", "black" )
)

DDagBT <- subset(DDagB, Type %in% c("K"))
BarErr2(100*DDagBT$x, 
	factor(paste("POS:", DDagBT$Pos, " HI:", DDagBT$Hi, sep=""), 
		levels=c("POS:0 HI:0", "POS:0 HI:1", "POS:1 HI:0", "POS:1 HI:1"),
		labels=c("POS:0 HI:0", "POS:0 HI:1", "POS:1 HI:0", "POS:1 HI:1"),		), 
	factor(DDagBT$Block, levels=2:5, labels=1:4), beside=T, leg=F, ylim=c(0,100), 
	col=c("red", "green", "red", "green" ), axes=F)
BarErr2(100*DDagBT$x, 
	factor(paste("POS:", DDagBT$Pos, " HI:", DDagBT$Hi, sep=""), 
		levels=c("POS:0 HI:0", "POS:0 HI:1", "POS:1 HI:0", "POS:1 HI:1"),
		labels=c("POS:0 HI:0", "POS:0 HI:1", "POS:1 HI:0", "POS:1 HI:1"),		), 
	factor(DDagBT$Block, levels=2:5, labels=1:4), beside=T, leg=F, ylim=c(0,100), 
	density=15, add=T,
	col=c("black", "black", "red", "green" ), axes=F)
axis(2, at=c(0:4)*25, labels=paste(c(0:4)*25, "%", sep=""), las=1)
axis(1, at=c(1:4)*25, labels=1:4)
legend("topright", density=25, fill=c("black", "white"), legend=c("NEGATIVE", "POSITIVE"), pt.lwd=2, inset=.05)
legend("top", fill=c("red", "green"), legend=c("LOW", "HIGH"), pt.lwd=2, inset=.05)






DDagBTA <- subset(DDagB, Type %in% c("M") & Block<3)
par(mai=c(.8,.65,.2,0))
BarErr2(100*DDagBTA$x,
	factor(DDagBTA$Hi, levels=c(0,1), labels=c("LOW", "HIGH")),
	factor(DDagBTA$Pos, levels=c(0,1), labels=c("NEG", "POS")),
	beside=T, leg=F, ylim=c(0,100), axes=F, col=rep(c("red", "green"), times=2))
BarErr2(100*DDagBTA$x,
	factor(DDagBTA$Hi, levels=c(0,1), labels=c("LOW", "HIGH")),
	factor(DDagBTA$Pos, levels=c(0,1), labels=c("NEG", "POS")),
	beside=T, leg=F, ylim=c(0,100), axes=F, density=15, add=T, col=c("black", "black", "red", "green"))
axis(2, at=c(0:4)*25, labels=paste(c(0:4)*25, "%", sep=""), las=1)
legend("topright", fill=c("red", "green"), legend=c("LOW", "HIGH"), pt.lwd=2, inset=.05)

















DDagBT <- subset(DDagB, Type %in% c("K", "KB+", "KB-"))
BarErr2(100*DDagBT$x, 
	factor(paste("POS:", DDagBT$Pos, " HI:", DDagBT$Hi, sep=""), 
		levels=c("POS:0 HI:0", "POS:0 HI:1", "POS:1 HI:0", "POS:1 HI:1"),
		labels=c("POS:0 HI:0", "POS:0 HI:1", "POS:1 HI:0", "POS:1 HI:1"),		), 
	factor(DDagBT$Block, levels=2:5, labels=1:4), beside=T, leg=T, ylim=c(0,100), 
	col=c("red", "green", "blue", "black" ))






SCENARIO <- "K"
summary(aov(x~Hi + Error(Subject/Hi), data=subset(DDagB, Type == SCENARIO & Pos==1 & Block==2)))
summary(aov(x~Hi + Error(Subject/Hi), data=subset(DDagB, Type == SCENARIO & Pos==0 & Block==2)))
summary(aov(x~Pos*Hi + Error(Subject/(Pos*Hi)), data=subset(DDagB, Type == SCENARIO & Block==2)))











## PREDICTION GRAPHS








par(mai=c(.8,.65,.2,0))
BarErr2(c(0,0,90,10),
	factor(c(0,1,0,1), levels=c(0,1), labels=c("LOW", "HIGH")),
	factor(c(0,0,1,1), levels=c(0,1), labels=c("", "POS")),
	beside=T, leg=F, ylim=c(0,100), axes=F, col=rep(c("red", "green"), times=2), xpd=F)
BarErr2(c(0,0,90,10),
	factor(c(0,1,0,1), levels=c(0,1), labels=c("LOW", "HIGH")),
	factor(c(0,0,1,1), levels=c(0,1), labels=c("", "POS")),
	beside=T, leg=F, ylim=c(0,100), axes=F, density=15, add=T, col=c("black", "black", "red", "green"), xpd=F)
axis(2, at=c(0:4)*25, labels=paste(c(0:4)*25, "%", sep=""), las=1)
legend("topright", fill=c("red", "green"), legend=c("LOW", "HIGH"), pt.lwd=2, inset=.05)




par(mai=c(.8,.65,.2,0))
BarErr2(c(0,90,90,0),
	factor(c(0,1,0,1), levels=c(0,1), labels=c("LOW", "HIGH")),
	factor(c(0,0,1,1), levels=c(0,1), labels=c("NEG", "POS")),
	beside=T, leg=F, ylim=c(0,100), axes=F, col=rep(c("red", "green"), times=2), xpd=F)
BarErr2(c(0,90,90,0),
	factor(c(0,1,0,1), levels=c(0,1), labels=c("LOW", "HIGH")),
	factor(c(0,0,1,1), levels=c(0,1), labels=c("NEG", "POS")),
	beside=T, leg=F, ylim=c(0,100), axes=F, density=15, add=T, col=c("black", "black", "red", "green"), xpd=F)
axis(2, at=c(0:4)*25, labels=paste(c(0:4)*25, "%", sep=""), las=1)
legend("topright", fill=c("red", "green"), legend=c("LOW", "HIGH"), pt.lwd=2, inset=.05)


par(mai=c(.8,.65,.2,0))
BarErr2(c(0,90,90,90),
	factor(c(0,1,0,1), levels=c(0,1), labels=c("LOW", "HIGH")),
	factor(c(0,0,1,1), levels=c(0,1), labels=c("NEG", "POS")),
	beside=T, leg=F, ylim=c(0,100), axes=F, col=rep(c("red", "green"), times=2), xpd=F)
BarErr2(c(0,90,90,90),
	factor(c(0,1,0,1), levels=c(0,1), labels=c("LOW", "HIGH")),
	factor(c(0,0,1,1), levels=c(0,1), labels=c("NEG", "POS")),
	beside=T, leg=F, ylim=c(0,100), axes=F, density=15, add=T, col=c("black", "black", "red", "green"), xpd=F)
axis(2, at=c(0:4)*25, labels=paste(c(0:4)*25, "%", sep=""), las=1)
legend("top", fill=c("red", "green"), legend=c("LOW", "HIGH"), pt.lwd=2, inset=.05)







par(mai=c(.8,.65,.2,0))
BarErr2(c(10,90,90,90),
	factor(c(0,1,0,1), levels=c(0,1), labels=c("LOW", "HIGH")),
	factor(c(0,0,1,1), levels=c(0,1), labels=c("NEG", "POS")),
	beside=T, leg=F, ylim=c(0,100), axes=F, col=rep(c("red", "green"), times=2), xpd=F)
BarErr2(c(10,90,90,90),
	factor(c(0,1,0,1), levels=c(0,1), labels=c("LOW", "HIGH")),
	factor(c(0,0,1,1), levels=c(0,1), labels=c("NEG", "POS")),
	beside=T, leg=F, ylim=c(0,100), axes=F, density=15, add=T, col=c("black", "black", "red", "green"), xpd=F)
axis(2, at=c(0:4)*25, labels=paste(c(0:4)*25, "%", sep=""), las=1)
legend("top", fill=c("red", "green"), legend=c("LOW", "HIGH"), pt.lwd=2, inset=.05)













#NUMBERS FOR LATEX

meansderr <- function(x) {return (paste("\\drawbar{", round(mean(x),1), "}{", round(StdErr(x),1), "}", sep=""));    }

tapply(100*DDag$x, 
	list(
		factor(paste("POS:", DDag$Pos, " HI:", DDag$Hi, sep="")),
		factor(DDag$Type, levels=c("K", "C", "M"))
		),
	meansderr)

tapply(100*DD$Ireponse,	
	list(
		paste("POS:", DD$POS, " HI:", DD$HI, sep=""),
		factor(paste(DD$TYPE, DD$SCENARIO), levels=paste("C",1:4))
		), 
	meansderr)


DDagBTA <- subset(DDagB, Type %in% c("K", "C", "M") & Block<3)
tapply(100*DDagBTA$x, 
	list(
		factor(paste("POS:", DDagBTA$Pos, " HI:", DDagBTA$Hi, sep="")),
		factor(DDagBTA$Type, levels=c("K", "C", "M"))
		),
	meansderr)