rm(list = ls()); dev.off()

setwd("//Volumes/MattSSD Files/Users/matt_a/PHDTHESIS/MA_PHDTHESIS_FILES/DotProbe_Stereo_new/Analyse_Results/")
data = read.csv('allmdls_humandata_varmeandepth.txt',header = T, sep = ",")
data = data[data$expt == 2,]
data = data[data$RT > 500,]
data = data[data$RT < 12000,]

# most of the preprocessing has already been done.
data = data[c('Ppt','Prestime','gt','human_resp','human_respratio','gtratio','ddiff','ID','scene','view','stereo','meandepth','elev','depthmu')]
data = data[data$gt == data$human_resp,]

head(data); sapply(data, class)

data$Ppt = as.factor(data$Ppt)
data$ID = as.factor(data$ID)
data$scene = as.factor(data$scene)
data$view = as.factor(data$view)
data$elev = as.factor(data$elev)
data$meandepth = as.factor(data$meandepth)
data$stereo = as.factor(data$stereo)

# convert to numeric? 
mapping = c("meandist_10m" = 10, "meandist_20m" = 20, "meandist_40m" = 40)
data$meandepth = mapping[as.character(data$meandepth)]
data$meandepth = as.numeric(data$meandepth)
data$depthmu = as.numeric(data$depthmu)

data$ddiff = as.numeric(data$ddiff)
data$Prestime = as.numeric(data$Prestime)
data$gtratio = as.numeric(data$gtratio)
data$human_resp = as.numeric(data$human_resp)

sapply(data, class)

# remove any outlier subjects
submeans = tapply(data$human_respratio, data$Ppt, mean)
plot(submeans)

# now our dv could be mean-squared error or just error
data$abserror = abs(data$human_respratio-data$gtratio)
submeans = tapply(data$abserror, data$Ppt, mean)
plot(submeans)

# let's test whether log odds varies linearly with presentation time
logodds_prestime = tapply(data$abserror, data$Prestime,mean)
plot(c(1,2,4,8,16),logodds_prestime, type="l")
plot(log(c(1,2,4,8,16)), logodds_prestime, type="l")

# log odds varies more linearly with the presentation time than log prestime
data$Prestime = log(data$Prestime)

# and let's test whether log odds varies linearly with gtratio
gtratbin = round(data$gtratio*10)
logodds_gtrat = tapply(data$abserror, gtratbin, mean)
plot(0:10,logodds_gtrat, type="l") # no monotonic relationship

# and finally the depth difference
gtratbin = round((data$ddiff/max(data$ddiff))*10)
logodds_gtrat = tapply(data$abserror, gtratbin, mean)
plot(0:10,logodds_gtrat, type="l") # no monotonic relationship

# check whether variance is homogenous across groups. No obvious outliers can be detected
plot(tapply(data$abserror, data$Ppt, var))
plot(tapply(data$abserror, data$ID, var))
plot(tapply(data$abserror, data$scene, var))


# ---------------------- now we set up the contrasts ---------------------- #

# For post-hocs & main effects, Set up the appropriate contrasts. There will be 3*2 = 6 unique conditions
# generalized inverse function:
library(MASS)
ginv2 = function(x) fractions(provideDimnames(ginv(x), base = dimnames(x)[2:1]))

data$elev = factor(data$elev, labels = c("same","var"))
data$stereo = factor(data$stereo, labels = c("mono","stereo"))

h_matrix1 = rbind(elevationH1 = c(same = -1, var = 1))
h_matrix2 = rbind(stereoH1 = c(mono = -1, stereo = 1))
c_matrix_elev = ginv2(h_matrix1)
c_matrix_stereo = ginv2(h_matrix2)


# ---------------------- now we run the LMM ---------------------- #

library(lme4); library(optimx)
library(lmerTest)

data$Prestime = (data$Prestime-mean(data$Prestime))/sd(data$Prestime)
data$gtratio = (data$gtratio-mean(data$gtratio))/sd(data$gtratio)
data$ddiff = (data$ddiff-mean(data$ddiff))/sd(data$ddiff)
data$meandepth = (data$meandepth-mean(data$meandepth))/sd(data$meandepth)
data$depthmu = (data$depthmu-mean(data$depthmu))/sd(data$depthmu)
data$human_respratio = (data$human_respratio-mean(data$human_respratio))/sd(data$human_respratio)

logitlmertest2 = lmer(human_respratio ~ meandepth * elev * stereo * Prestime * ddiff + (1 + Prestime|Ppt) + (1 + Prestime + gtratio|scene),
                      data = data, contrasts = list(elev = c_matrix_elev, stereo = c_matrix_stereo))

logitlmertest = lmer(human_respratio ~ meandepth * elev * stereo * Prestime * gtratio + (1 + Prestime|Ppt) + (1 + Prestime + gtratio|scene),
                     data = data, contrasts = list(elev = c_matrix_elev, stereo = c_matrix_stereo))
summary(logitlmertest)
AIC(logitlmertest)
AIC(logitlmertest2)

# export results to csv? 
write.csv(summary(logitlmertest)$coefficients, "intervaltaskFixedEstimates.txt")

## note that it matters very little whether the depth difference is coded as a ratio, or as a raw metric difference
print("------------------------ LMM OUTPUT ---------------------------------")
summary(logitlmertest)
x = confint.merMod(object = logitlmertest, level = .95, method = "Wald")
x


library(arm)
designmatrix = as.data.frame(model.matrix(logitlmertest))
fixedestimates = fixef(logitlmertest)

sstot = sum((data$human_respratio - mean(data$human_respratio))^2)
r2 = 0

for (i in 1:length(fixedestimates)){
  singpred = as.vector(fixedestimates[i] %*% t(designmatrix[,i]))
  ssresid = sum((data$human_respratio - singpred)^2)
  r2[i] = 1-(ssresid/sstot)
}


x = data.frame(r2)
x

# at what point does performance exceed chance? 
data$Pres_Time = as.factor(data$Prestime)

data$Pres_Time = factor(data$Pres_Time, labels = c("P1", "P2", "P3", "P4","P5"))

# let's try a different coding
library(MASS)
ginv2 = function(x) fractions(provideDimnames(ginv(x), base = dimnames(x)[2:1]))
data2 = data[data$Pres_Time == "P1",]
logitlmertest_prestime = lmer(human_respratio ~ gtratio + (1|Ppt) + (1|scene), data = data2)
summ1 = summary(logitlmertest_prestime)
print(summ1)
