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