Merge branch 'master' of https://git.art-ist.cc/daniel/tensor_predictors
This commit is contained in:
commit
6cbb5029d4
2
.gitattributes
vendored
Normal file
2
.gitattributes
vendored
Normal file
@ -0,0 +1,2 @@
|
||||
*.rds filter=lfs diff=lfs merge=lfs -text
|
||||
*.Rds filter=lfs diff=lfs merge=lfs -text
|
||||
4
.gitignore
vendored
4
.gitignore
vendored
@ -45,8 +45,8 @@
|
||||
|
||||
## R environment, data and package build intermediate files/folders
|
||||
# R Data Object files
|
||||
*.Rds
|
||||
*.rds
|
||||
# *.Rds # git-lfs
|
||||
# *.rds # git-lfs
|
||||
*.Rdata
|
||||
|
||||
# Example code in package build process
|
||||
|
||||
@ -25,17 +25,6 @@ c(X, F, y) %<-% local({
|
||||
list(X, array(y, c(n, 1L, 1L)), y)
|
||||
})
|
||||
|
||||
# fit a tensor normal model to the data sample axis 1 indexes persons)
|
||||
fit.gmlm <- gmlm_tensor_normal(X, F, sample.axis = 1L)
|
||||
|
||||
# plot the fitted mode wise reductions (for time and sensor axis)
|
||||
with(fit.gmlm, {
|
||||
par.reset <- par(mfrow = c(2, 1))
|
||||
plot(seq(0, 1, len = 256), betas[[1]], main = "Time", xlab = "Time [s]", ylab = expression(beta[1]))
|
||||
plot(betas[[2]], main = "Sensors", xlab = "Sensor Index", ylab = expression(beta[2]))
|
||||
par(par.reset)
|
||||
})
|
||||
|
||||
|
||||
#' (2D)^2 PCA preprocessing
|
||||
#'
|
||||
@ -54,31 +43,6 @@ preprocess <- function(X, tpc, ppc) {
|
||||
}
|
||||
|
||||
|
||||
#' Leave-one-out prediction
|
||||
#'
|
||||
#' @param X 3D EEG data (preprocessed or not)
|
||||
#' @param F binary responce `y` as a 3D tensor, every obs. is a 1 x 1 matrix
|
||||
loo.predict <- function(X, F) {
|
||||
sapply(seq_len(dim(X)[1L]), function(i) {
|
||||
# Fit with i'th observation removes
|
||||
fit <- gmlm_tensor_normal(X[-i, , ], F[-i, , , drop = FALSE], sample.axis = 1L)
|
||||
|
||||
# Reduce the entire data set
|
||||
r <- as.vector(mlm(X, fit$betas, modes = 2: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) - %3d/%d\n", dim(X)[2L], dim(X)[3L], i, dim(X)[1L]))
|
||||
|
||||
y.hat
|
||||
})
|
||||
}
|
||||
|
||||
|
||||
### 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)
|
||||
@ -97,11 +61,48 @@ auc <- function(y.true, y.pred) as.numeric(pROC::roc(y.true, y.pred, quiet = TRU
|
||||
auc.sd <- function(y.true, y.pred) sqrt(pROC::var(pROC::roc(y.true, y.pred, quiet = TRUE)))
|
||||
|
||||
|
||||
##################################### GMLM #####################################
|
||||
|
||||
# fit a tensor normal model to the data sample axis 1 indexes persons)
|
||||
fit.gmlm <- gmlm_tensor_normal(X, F, sample.axis = 1L)
|
||||
|
||||
# plot the fitted mode wise reductions (for time and sensor axis)
|
||||
with(fit.gmlm, {
|
||||
par.reset <- par(mfrow = c(2, 1))
|
||||
plot(seq(0, 1, len = 256), betas[[1]], main = "Time", xlab = "Time [s]", ylab = expression(beta[1]))
|
||||
plot(betas[[2]], main = "Sensors", xlab = "Sensor Index", ylab = expression(beta[2]))
|
||||
par(par.reset)
|
||||
})
|
||||
|
||||
#' Leave-one-out prediction using GMLM
|
||||
#'
|
||||
#' @param X 3D EEG data (preprocessed or not)
|
||||
#' @param F binary responce `y` as a 3D tensor, every obs. is a 1 x 1 matrix
|
||||
loo.predict.gmlm <- function(X, F) {
|
||||
unlist(parallel::mclapply(seq_len(dim(X)[1L]), function(i) {
|
||||
# Fit with i'th observation removed
|
||||
fit <- gmlm_tensor_normal(X[-i, , ], F[-i, , , drop = FALSE], sample.axis = 1L)
|
||||
|
||||
# Reduce the entire data set
|
||||
r <- as.vector(mlm(X, fit$betas, modes = 2: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) - %3d/%d\n", dim(X)[2L], dim(X)[3L], i, dim(X)[1L]))
|
||||
|
||||
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(preprocess(X, 3, 4), F)
|
||||
y.hat.15.15 <- loo.predict(preprocess(X, 15, 15), F)
|
||||
y.hat.20.30 <- loo.predict(preprocess(X, 20, 30), F)
|
||||
y.hat <- loo.predict(X, F)
|
||||
y.hat.3.4 <- loo.predict.gmlm(preprocess(X, 3, 4), F)
|
||||
y.hat.15.15 <- loo.predict.gmlm(preprocess(X, 15, 15), F)
|
||||
y.hat.20.30 <- loo.predict.gmlm(preprocess(X, 20, 30), F)
|
||||
y.hat <- loo.predict.gmlm(X, F)
|
||||
|
||||
# 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) {
|
||||
@ -117,3 +118,51 @@ y.hat <- loo.predict(X, F)
|
||||
#> tnr 0.64444444 0.60000000 0.60000000 0.60000000
|
||||
#> auc 0.85108225 0.83838384 0.83924964 0.83896104
|
||||
#> auc.sd 0.03584791 0.03760531 0.03751307 0.03754553
|
||||
|
||||
|
||||
################################## 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_len(dim(X)[1L]), function(i) {
|
||||
# Fit with i'th observation removed
|
||||
fit <- TSIR(X[-i, , ], y[-i], c(1L, 1L), sample.axis = 1L)
|
||||
|
||||
# Reduce the entire data set
|
||||
r <- as.vector(mlm(X, fit, modes = 2: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) - %3d/%d\n", dim(X)[2L], dim(X)[3L], i, dim(X)[1L]))
|
||||
|
||||
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), y)
|
||||
y.hat.15.15 <- loo.predict.tsir(preprocess(X, 15, 15), y)
|
||||
y.hat.20.30 <- loo.predict.tsir(preprocess(X, 20, 30), 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)(y, y.pred) })
|
||||
}))
|
||||
#> y.hat.3.4 y.hat.15.15 y.hat.20.30 y.hat
|
||||
#> acc 0.79508197 0.78688525 0.7540984 0.67213115
|
||||
#> err 0.20491803 0.21311475 0.2459016 0.32786885
|
||||
#> fpr 0.33333333 0.37777778 0.3777778 0.53333333
|
||||
#> tpr 0.87012987 0.88311688 0.8311688 0.79220779
|
||||
#> fnr 0.12987013 0.11688312 0.1688312 0.20779221
|
||||
#> tnr 0.66666667 0.62222222 0.6222222 0.46666667
|
||||
#> auc 0.84646465 0.83376623 0.8040404 0.68946609
|
||||
#> auc.sd 0.03596227 0.04092069 0.0446129 0.05196611
|
||||
|
||||
BIN
dataAnalysis/eeg/eeg_data.rds
(Stored with Git LFS)
Normal file
BIN
dataAnalysis/eeg/eeg_data.rds
(Stored with Git LFS)
Normal file
Binary file not shown.
Loading…
x
Reference in New Issue
Block a user