library(grt)
library(tidyverse)

#load the data to be fit
#This is a restricted set, only one subject and one block
setwd('/home/tzhu/Desktop/pilot_data/CR_task')
d<- read.csv("data3_cr.csv")
cols <- c(4:6, 25:33)
d<-d[,cols]
d <- d[complete.cases(d), ]
d$block<-rep(1:17, each=30)
#relabel the category as 1 and 2
d<-(d %>% mutate(corrAns = case_when(
  corrAns == "a" ~ 1,
  corrAns == "b" ~ 2))
)

#relabel the key_resp_3.keys as 1 and 2
d<-(d %>% mutate(key_resp_3.keys = case_when(
  key_resp_3.keys == "a" ~ 1,
  key_resp_3.keys == "b" ~ 2))
)

Test code with block 1

d1<-d[d$block==1, ]

fit.ori <- glc(key_resp_3.keys ~ ori_transform, data=d1, category=d1$corrAns, zlimit=7)
fit.len <- glc(key_resp_3.keys ~ length, data=d1, category=d1$corrAns, zlimit=7)
fit.2dl <- glc(key_resp_3.keys ~ ori_transform + length, data=d1, category=d1$corrAns, zlimit=7)
fit.conj <- gcjc(key_resp_3.keys ~ ori_transform + length, data=d1, category=d1$corrAns, zlimit=7)
fit.guess <- grg(d$key_resp_3.keys, fixed=TRUE, k=2)

#This code compares the three by calculating an AIC. Lower is better
AIC(fit.ori,fit.len,fit.2dl, fit.conj)
##          df      AIC
## fit.ori   2 35.25931
## fit.len   2 46.85566
## fit.2dl   3 34.43228
## fit.conj  3 29.45281
extractAIC(fit.guess) #AIC argument does not work for guess model
## [1]   0.0000 707.0101
#Now just plot the boundary.
#Circles are category 1, triangles are category2
#Open are category 1 responses, filled are category 2 responses
plot.glc(fit.ori)

plot.glc(fit.len)

plot.glc(fit.2dl)

plot.gcjc(fit.conj)

AIC table for 4 models across 17 blocks

AIC_table <- lapply(1:17, function(i){
  k <- d$block == i
  fit.ori <- glc(key_resp_3.keys ~ ori_transform, data = d[k, ], category = d$corrAns[k], zlimit = 7)
  fit.len <- glc(key_resp_3.keys ~ length, data = d[k, ], category = d$corrAns[k], zlimit = 7)
  fit.2dl <- glc(key_resp_3.keys ~ ori_transform + length, data = d[k, ], category = d$corrAns[k], zlimit = 7)
  fit.conj <- gcjc(key_resp_3.keys ~ ori_transform + length, data= d[k,], category=d$corrAns[k], zlimit=7)
  accuracy<-mean(d[k,]$key_resp_3.corr)
  data.frame(Accuracy= accuracy, AIC_ori = AIC(fit.ori), AIC_len = AIC(fit.len), AIC_2dl = AIC(fit.2dl), AIC_conj = AIC(fit.conj))
})

AIC_table <- do.call(rbind, AIC_table)

# find the minimus AIC value in each row, then return column name as block strategy
AIC_table$strategy<-apply(AIC_table[, 2:5], 1, function(x)which(x==min(x)))

AIC_table
##     Accuracy  AIC_ori  AIC_len   AIC_2dl AIC_conj strategy
## 1  0.5333333 35.25931 46.85566 34.432282 29.45281        4
## 2  0.7000000 31.13986 39.66358 31.520196 33.76653        1
## 3  0.8333333 25.82446 42.33881 26.558980 28.21877        1
## 4  0.8000000 19.69266 46.11335 19.950742 21.69045        1
## 5  0.5666667 25.12632 44.22716 27.104466 25.81603        1
## 6  0.7666667 24.38431 45.04774 25.877140 26.15590        1
## 7  0.5666667 24.81727 43.54058 26.781696 26.71505        1
## 8  0.5666667 27.41906 45.40984 28.467868 30.07352        1
## 9  0.7000000 23.75327 46.20208 20.514740 25.91689        3
## 10 0.7000000 27.50840 44.45289 29.248128 29.73962        1
## 11 0.8000000 10.03191 43.35437  9.905092 12.03309        3
## 12 0.7666667 19.34971 45.49179 21.349699 21.34958        1
## 13 0.7000000 19.61067 41.37415 20.788665 22.01121        1
## 14 0.7666667 31.56453 43.24881 32.586302 34.13098        1
## 15 0.6666667 29.96392 45.40342 31.772704 32.65796        1
## 16 0.7666667 37.79191 41.72087 36.458396 40.63903        3
## 17 0.7000000 32.68212 42.31376 32.718841 34.60942        1