library(tensorPredictors) library(parallel) library(pROC) #' Mode-Wise PCA preprocessing #' #' @param npc_time Number of Principal Components for time axis #' @param npc_sensor Number of Principal Components for sensor axis #' @param npc_condition Number of Principal Components for stimulus condition axis preprocess <- function(X, npc_time, npc_sensor, npc_condition) { # Mode covariances (for predictor and time point modes) c(Sigma_t, Sigma_s, Sigma_c) %<-% mcov(X) # "predictor" (sensor) and time point principal components V_t <- svd(Sigma_t, npc_time, 0L)$u V_s <- svd(Sigma_s, npc_sensor, 0L)$u V_c <- svd(Sigma_c, npc_condition, 0L)$u # reduce with mode wise PCs mlm(X, list(V_t, V_s, V_c), modes = 1:3, transposed = TRUE) } ### Classification performance measures # acc: Accuracy. P(Yhat = Y). Estimated as: (TP+TN)/(P+N). acc <- function(y.true, y.pred) mean(round(y.pred) == y.true) # err: Error rate. P(Yhat != Y). Estimated as: (FP+FN)/(P+N). err <- function(y.true, y.pred) mean(round(y.pred) != y.true) # fpr: False positive rate. P(Yhat = + | Y = -). aliases: Fallout. fpr <- function(y.true, y.pred) mean((round(y.pred) == 1)[y.true == 0]) # tpr: True positive rate. P(Yhat = + | Y = +). aliases: Sensitivity, Recall. tpr <- function(y.true, y.pred) mean((round(y.pred) == 1)[y.true == 1]) # fnr: False negative rate. P(Yhat = - | Y = +). aliases: Miss. fnr <- function(y.true, y.pred) mean((round(y.pred) == 0)[y.true == 1]) # tnr: True negative rate. P(Yhat = - | Y = -). tnr <- function(y.true, y.pred) mean((round(y.pred) == 0)[y.true == 0]) # auc: Area Under the Curve auc <- function(y.true, y.pred) as.numeric(pROC::roc(y.true, y.pred, quiet = TRUE)$auc) auc.sd <- function(y.true, y.pred) sqrt(pROC::var(pROC::roc(y.true, y.pred, quiet = TRUE))) # Load full EEG dataset of all subjects c(X, y) %<-% readRDS("eeg_data_3d.rds") ##################################### GMLM ##################################### #' Leave-one-out prediction using GMLM #' #' @param X 3D EEG data (preprocessed or not) #' @param y binary responce `y` as a 3D tensor, every obs. is a 1 x 1 matrix loo.predict.gmlm <- function(X, y) { unlist(parallel::mclapply(seq_along(y), function(i) { # Fit with i'th observation removed fit <- gmlm_tensor_normal(X[ , , , -i], as.integer(y[-i]), sample.axis = 4L) # Reduce the entire data set r <- as.vector(mlm(X, fit$betas, modes = 1:3, transpose = TRUE)) # Fit a logit model on reduced data with i'th observation removed logit <- glm(y ~ r, family = binomial(link = "logit"), data = data.frame(y = y[-i], r = r[-i]) ) # predict i'th response given i'th reduced observation y.hat <- predict(logit, newdata = data.frame(r = r[i]), type = "response") # report progress cat(sprintf("dim: (%d, %d, %d) - %3d/%d\n", dim(X)[1], dim(X)[2], dim(X)[3], i, length(y) )) y.hat }, mc.cores = getOption("mc.cores", max(1L, parallel::detectCores() - 1L)))) } # perform preprocessed (reduced) and raw (not reduced) leave-one-out prediction y.hat.3.4 <- loo.predict.gmlm(preprocess(X, 3, 4, 3), y) y.hat.15.15 <- loo.predict.gmlm(preprocess(X, 15, 15, 3), y) y.hat.20.30 <- loo.predict.gmlm(preprocess(X, 20, 30, 3), y) y.hat <- loo.predict.gmlm(X, y) # classification performance measures table by leave-one-out cross-validation (loo.cv <- apply(cbind(y.hat.3.4, y.hat.15.15, y.hat.20.30, y.hat), 2, function(y.pred) { sapply(c("acc", "err", "fpr", "tpr", "fnr", "tnr", "auc", "auc.sd"), function(FUN) { match.fun(FUN)(as.integer(y) - 1L, y.pred) }) })) #> y.hat.3.4 y.hat.15.15 y.hat.20.30 y.hat #> acc 0.83606557 0.80327869 0.80327869 0.80327869 #> err 0.16393443 0.19672131 0.19672131 0.19672131 #> fpr 0.31111111 0.33333333 0.33333333 0.33333333 #> tpr 0.92207792 0.88311688 0.88311688 0.88311688 #> fnr 0.07792208 0.11688312 0.11688312 0.11688312 #> tnr 0.68888889 0.66666667 0.66666667 0.66666667 #> auc 0.88023088 0.87070707 0.87041847 0.86810967 #> auc.sd 0.03124875 0.03244623 0.03248653 0.03295883 ################################## Tensor SIR ################################## #' Leave-one-out prediction using TSIR #' #' @param X 3D EEG data (preprocessed or not) #' @param y binary responce vector loo.predict.tsir <- function(X, y) { unlist(parallel::mclapply(seq_along(y), function(i) { # Fit with i'th observation removed fit <- TSIR(X[ , , , -i], y[-i], c(1L, 1L, 1L), sample.axis = 4L) # Reduce the entire data set r <- as.vector(mlm(X, fit, modes = 1:3, transpose = TRUE)) # Fit a logit model on reduced data with i'th observation removed logit <- glm(y ~ r, family = binomial(link = "logit"), data = data.frame(y = y[-i], r = r[-i]) ) # predict i'th response given i'th reduced observation y.hat <- predict(logit, newdata = data.frame(r = r[i]), type = "response") # report progress cat(sprintf("dim: (%d, %d, %d) - %3d/%d\n", dim(X)[1], dim(X)[2], dim(X)[3], i, length(y) )) y.hat }, mc.cores = getOption("mc.cores", max(1L, parallel::detectCores() - 1L)))) } # perform preprocessed (reduced) and raw (not reduced) leave-one-out prediction y.hat.3.4 <- loo.predict.tsir(preprocess(X, 3, 4, 3), y) y.hat.15.15 <- loo.predict.tsir(preprocess(X, 15, 15, 3), y) y.hat.20.30 <- loo.predict.tsir(preprocess(X, 20, 30, 3), y) y.hat <- loo.predict.tsir(X, y) # classification performance measures table by leave-one-out cross-validation (loo.cv <- apply(cbind(y.hat.3.4, y.hat.15.15, y.hat.20.30, y.hat), 2, function(y.pred) { sapply(c("acc", "err", "fpr", "tpr", "fnr", "tnr", "auc", "auc.sd"), function(FUN) { match.fun(FUN)(as.integer(y) - 1L, y.pred) }) })) #> y.hat.3.4 y.hat.15.15 y.hat.20.30 y.hat #> acc 0.81967213 0.84426230 0.81147541 0.80327869 #> err 0.18032787 0.15573770 0.18852459 0.19672131 #> fpr 0.33333333 0.24444444 0.33333333 0.33333333 #> tpr 0.90909091 0.89610390 0.89610390 0.88311688 #> fnr 0.09090909 0.10389610 0.10389610 0.11688312 #> tnr 0.66666667 0.75555556 0.66666667 0.66666667 #> auc 0.86522367 0.89379509 0.88196248 0.85974026 #> auc.sd 0.03357539 0.03055047 0.02986038 0.03367847