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','fdist','pdist','stereo','meandepth','elev','depthmu')]

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$correct = as.numeric(data$gt == data$human_resp)

# 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)

submeans = tapply(data$correct, data$Ppt, mean)
plot(submeans)

# let's test whether log odds varies linearly with presentation time
logodds_prestime = tapply(data$correct, 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$correct, 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$correct, gtratbin, mean)
plot(0:10,logodds_gtrat, type="l") # no monotonic relationship

# ---------------------- 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)

logitlmertest = glmer(correct ~ meandepth * elev * stereo * Prestime * gtratio + (1 + Prestime|Ppt) + (1 + Prestime + gtratio|scene),
                      data = data, family = binomial, contrasts = list(elev = c_matrix_elev, stereo = c_matrix_stereo), 
                      control=glmerControl(optimizer = "optimx", optCtrl=list(method = 'bobyqa', maxfun=1000000)))
# logitlmertest2 = glmer(correct ~ meandepth * elev * stereo * Prestime * ddiff + (1 + Prestime|Ppt) + (1 + Prestime + gtratio|scene),
#                       data = data, family = binomial, contrasts = list(elev = c_matrix_elev, stereo = c_matrix_stereo), 
#                       control=glmerControl(optimizer = "optimx", optCtrl=list(method = 'bobyqa', maxfun=1000000)))
summary(logitlmertest)
AIC(logitlmertest)
AIC(logitlmertest2)
write.csv(summary(logitlmertest)$coefficients, "ordinaltaskFixedEstimates.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")
data.frame(x)

library(arm)
designmatrix = as.data.frame(model.matrix(logitlmertest))
fixedestimates = fixef(logitlmertest)
fixedpreds = var(as.vector(fixedestimates %*% t(designmatrix)))

# Overall model fit:
randomeffect_variance = as.data.frame(VarCorr(logitlmertest))
randomeffect_variance = sum(randomeffect_variance$vcov)
residual = pi^2/3

rsquared = fixedpreds / (fixedpreds + randomeffect_variance + residual)
print(rsquared)

# effect sizes for all our tested effects
sumrsquared = 0
for (f in 1:length(fixedestimates)){
  single_pred = var(as.vector(fixedestimates[f] %*% t(designmatrix[,f])))
  rsquared[f] = single_pred / (fixedpreds + randomeffect_variance + residual)
  sumrsquared = sumrsquared + rsquared[f]
  print(fixedestimates[f])
  print(rsquared[f])
}

x = data.frame(rsquared)
x

# let's do two post-hoc interactions: mean depth x depth offset for the two elevation conditions
data2 = data[data$elev == "var",]
logitlmertest = glmer(correct ~ meandepth * gtratio + (1|Ppt) + (1 + Prestime + gtratio|scene),
                      data = data2, family = binomial)
summary(logitlmertest)
confint.merMod(object = logitlmertest, level = .95, method = "Wald")

data2 = data[data$elev == "same",]
logitlmertest = glmer(correct ~ meandepth * gtratio + (1|Ppt) + (1 + Prestime + gtratio|scene),
                      data = data2, family = binomial)
summary(logitlmertest)
confint.merMod(object = logitlmertest, level = .95, method = "Wald")


# test when performance exceeds 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
data2 = data[data$Pres_Time == "P1",]
logitlmertest_prestime = glmer(correct ~ 1 + (1|Ppt) + (1|scene), data = data2, family = binomial)
summ1 = summary(logitlmertest_prestime)
print(summ1)

# logitlmer automatically uses a baseline of 0.5. We want chance to be 1/C
chance = -log(-1 + (1/(1/2)))
oddsratios = summ1$coefficients[,1]-chance
zscores = (summ1$coefficients[,1]-chance)/(summ1$coefficients[,2])
pvalues = (1 - pnorm(abs(zscores))) * 2
print(oddsratios)
print(zscores)
print(pvalues)
