library(dplyr)
library(tibble)
library(tidyr)
library(ppcor) ### install.packages("ppcor")
library(MASS)
library(neuroca)
## load study data
pcstudy <- as_tibble(read.csv("~/Dropbox/New_pc_behavioural_data/study_fixations.csv")) %>%
filter(Image != "." & !(Subject %in% c(28,32, 109))) %>% droplevels()
## load test data to get accuracy values (only taking first fixation row of each trial -- see: 'slice(1)')
pctest <- as_tibble(read.csv("~/Dropbox/New_pc_behavioural_data/delay_fixations.csv")) %>%
filter(Image != "." & !(Subject %in% c(28,32, 109))) %>% droplevels()
pcstudy$ImageNumber <- as.integer(as.character(pcstudy$ImageNumber))
## create table for each study trial (Subject/Image)
study_tab <- eye_table("FixX", "FixY", duration="FixDuration", onset="FixStartTime",
groupvar=c("Image", "Subject", "Block"), data=pcstudy,
clip_bounds=c(112, (112+800), 684, 84),
vars=c("ImageVersion",
"ImageSet", "Block", "Image", "ImageNumber"))
## create table for each study trial (Subject/Image)
delay_tab <- eye_table("FixX", "FixY", duration="FixDuration", onset="FixOffset",
groupvar=c("Image", "Subject"), data=pctest,
clip_bounds=c(112, (112+800), 684, 84),
vars=c("ImageVersion",
"ImageSet", "Image", "ImageNumber"))
fill_mat <- function(imat, Y, levs) {
omat <- matrix(NA, length(levs), ncol(imat))
ind <- match(Y, levs)
omat[ind,] <- imat
row.names(omat) <- levs
omat
}
get_b <- function(sid) {
dfx <- study_tab %>% filter(Subject %in% sid)
dfxdens <- density_by(dfx, groups=c("ImageVersion"),
xbounds=c(0,800),
ybounds=c(0,600),
outdim=c(16,12),
duration_weighted=TRUE, sigma=50)
dfx2 <- delay_tab %>% filter(Subject %in% sid)
dfxdens2 <- density_by(dfx2, groups=c("ImageVersion"),
xbounds=c(0,800),
ybounds=c(0,600),
outdim=c(16,12),
duration_weighted=TRUE, sigma=50)
imat <- do.call(rbind, purrr::map(dfxdens$density, ~ as.vector(.$z)))
Y1 <- droplevels(dfxdens$ImageVersion)
imat2 <- do.call(rbind, purrr::map(dfxdens2$density, ~ as.vector(.$z)))
Y2 <- droplevels(dfxdens2$ImageVersion)
X1 <- fill_mat(imat, Y1, levels(pcstudy$ImageVersion))
X2 <- fill_mat(imat2, Y2, levels(pcstudy$ImageVersion))
list(study=X1, test=X2)
}
get_avgs <- function(Xs) {
res <- lapply(1:nrow(Xs[[1]]), function(i) {
iset <- lapply(1:length(Xs), function(j) {
Xs[[j]][i,]
})
isetavg <- colMeans(do.call(rbind, iset), na.rm = TRUE)
})
Xavg <- do.call(rbind, res)
row.names(Xavg) <- row.names(Xs[[1]])
Xavg
}
get_estimate <- function(Xs, Xavg) {
Xout <- Xs
for (i in 1:length(Xs)) {
for (j in 1:nrow(Xs[[i]])) {
if (any(is.na(Xs[[i]][j,]))) {
Xout[[i]][j,] <- Xavg[j,]
}
}
}
Xout
}
Xl <- lapply(unique(pcstudy$Subject), function(i) {
print(i)
get_block(i)
})
Xstudy <- lapply(Xl, "[[", "study")
Xtest <- lapply(Xl, "[[", "test")
Xavgstudy <- get_avgs(Xstudy)
Xstart <- get_estimate(Xstudy, Xavgstudy)
Xstart <- neuroca::to_block_matrix(do.call(cbind, Xstart),
block_lengths=rep(ncol(Xstart[[1]]), length(Xstart)))
#Xlstudy <- neuroca::to_block_matrix(do.call(cbind, Xstudy),
# block_lengths=rep(ncol(Xstudy[[1]]), length(Xstudy)))
mfa.1 <- neuroca:::mfa(Xstart, 15, center=FALSE, normalization="MFA")
avg_loadings <- function(fit, i) {
lds <- loadings(fit)[,i]
mat <- matrix(lds, length(Xstudy), 16*12, byrow=TRUE)
matrix(colMeans(mat),16,12)
}
res <- imputeMFA(as.data.frame(Xlstudy), group=rep(ncol(Xstudy[[1]]),length(Xstudy)), ncp=15, maxiter=25)
gids <- unlist(lapply(1:length(Xstudy), function(i) rep(i, ncol(Xstudy[[i]]))))
iset <- lapply(1:length(Xstudy), function(i) {
matrix(unlist(res$completeObs[7,gids==i]), 16,12)
})
isetavg <- t(apply(Reduce("+", iset)/length(iset), 1, rev))
image(iset[[1]]^(1/2))
image(iset[[22]]^(1/2))
image(iset[[21]]^(1/2))
Xstudy[[22]][7,]
pres <- neuroca:::mmpca(imat, Y, center=FALSE, knn=10)
refscores <- scores(pres)
sdres <- sda(scores(pres), L=dfx$ImageVersion, diagonal=TRUE)
dfy <- delay_tab %>% filter(Subject %in% c(115))
dfydens <- density_by(dfy, groups=c("ImageVersion"),
xbounds=c(0,800),
ybounds=c(0,600),
outdim=c(40,30),
duration_weighted=TRUE, sigma=80)
dmat <- do.call(rbind, purrr::map(dfydens$density, ~ as.vector(.$z)))
## code match trials
m <- pctest %>% group_by(Subject) %>% rowwise() %>% do( {
sversion <- as.character(filter(pcstudy, Subject == .$Subject & ImageNumber == .$ImageNumber)$ImageSet[1])
Match <- if (sversion == as.character(.$ImageSet)) "match" else "mismatch"
data.frame(Match=Match)
})
pctest$Match <- m$Match
## join encoding sim and test data
pctest2 <- inner_join(pctest, encoding_sim, by=c("Subject", "ImageNumber"))
pctest2$sid <- factor(pctest2$Subject)
## experiment with within-subject lmers... after adding ImageNumber as random effect, sig. effects are gone.
lme.1 <- lmer(Accuracy ~ sim_wavg + Saliency + Duration + (1 | sid), data=subset(pctest2, Match=="mismatch"))
lme.2 <- lmer(Accuracy ~ sim_within + sim_wavg + Saliency + Duration + (1 | sid), data=subset(pctest2, Match=="mismatch"))
lme.3 <- lmer(Accuracy ~ sim_within + sim_wavg + Saliency + Duration + (1 | sid) + (1 | ImageNumber), data=subset(pctest2, Match=="mismatch"))
lme.4 <- lmer(Accuracy ~ sim_wdiff + Saliency + Duration + (1 | sid) + (1 | ImageNumber), data=subset(pctest2, Match=="mismatch"))
## for between subject analysis, average accuracy over subject
pctestacc <- pctest %>% group_by(Subject) %>% summarize(acc=mean(Accuracy))
## compute subject averages for encoding_sim
pcstudysim <- encoding_sim %>% group_by(Subject) %>% summarize(sim_wavg=mean(sim_wavg), sim_within=mean(sim_within), sim_wdiff=mean(sim_wdiff),
sd_fix=mean(sd_fix),
total_fix=mean(total_fix), slope=mean(slope), slope_cor=mean(slope_cor, na.rm=TRUE))
## join accuracy and encoding_sim averages
dfx <- inner_join(pctestacc, pcstudysim, by="Subject")
## linear model
## sim_wavg = the average similarity with the subject's own mean fixation map.
## sim_within = average similarity among repetitions after partialing out the average_map
lm.1 <- lm(acc ~ sim_wavg, data=dfx, subset=Subject < 300 & Subject != 16)
lm.2 <- lm(acc ~ sim_within, data=dfx, subset=Subject < 300 & Subject != 16)
lm.3 <- lm(acc ~ sim_wavg + sim_within, data=dfx, subset=Subject < 300 & Subject != 16)
lm.4 <- lm(acc ~ sim_wdiff, data=dfx, subset=Subject < 300 & Subject != 16)
lm.4 <- lm(acc ~ sd_fix, data=dfx, subset=Subject < 300 & Subject != 16)
lm.5 <- lm(acc ~ total_fix, data=dfx, subset=Subject < 300 & Subject != 16)
lm.6 <- lm(acc ~ slope, data=dfx, subset=Subject < 300 & Subject != 16)
## plot
ggplot(aes(sim_within, acc), data=dfx) + geom_point() + stat_smooth(method=rlm)
ggplot(aes(sim_wavg, acc), data=dfx) + geom_point() + stat_smooth(method=rlm)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.