LogisticTaskRelIrr

Ref: Furuki & Takiyama, 2020, Sci Rep

https://www.nature.com/articles/s41598-020-59257-z

Copyright: Ken Takiyama

This is a R Markdown file as a referenece of our paper, Furuki & Takiyama, 2020, Sci Rep.

Although we analyzed our data in MATLAB using glmnet for MATLAB, the software does not work in new versins of MATLAB. We thus open our code as R markdown.

If you want to analyze your own data, plese arrange motion data \(\boldsymbol{X}\) to be a matrix whose size is \(T \times D\) and categorical data \(\boldsymbol{y}\) to be a vector whose size is \(T \times 1\). \(T\) and \(D\) mean the number of trials and the dimension of motion data, respectively. \(\boldsymbol{y}\) should include wither 0 or 1. Further, the number of 0s and 1s should be almost the same.

# for simulation
NumTri = 50
NumDim = 150
TotNum = NumTri*NumDim
W = matrix(rnorm(NumDim, 0, 1/NumDim), NumDim, 1)
X = matrix(rnorm(TotNum, 0, 1), nr = NumTri)
y = matrix(0.5*(sign(X%*%W)+1))
# for simulation

## NOTE 1: if you want to analyze your own data, replace the above X and y with your data. 
## NOTE 2: y should be defined as either 0 or 1. 


# find trial numbers with 0 and 1
fzero = which(y==0)
fone  = which(y==1)

# Calculate number of trials and dimensions
Xdim   = dim(X)
numtri = Xdim[1]
numdim = Xdim[2]
# Normalization (to make all the weight values be comparable)
Xn = (X-matrix(1,numtri,1)%*%colMeans(X))/(matrix(1,numtri,1)%*%apply(X, 2, sd))

# Don't use R 3.5.x. Pls use 3.6.2. for MacOSX Catalina. 
# install.packages("glmnet") # if necessary
# find the best regularization parameter
library(glmnet)

## Loading required package: Matrix

## Loaded glmnet 3.0-2

cvfit = cv.glmnet(Xn, y, family = "binomial", standarize = TRUE, alpha = 0)

# plot cross validation accuracy
plot(cvfit)

# fit again with the optimal regularization paramters
fit = glmnet(Xn, y, family = "binomial", lambda = cvfit$lambda.min, alpha = 0)

# weight values for classification
weight = fit$beta

# Task-relevant & task-irrelevant components
Xrel = Xn%*%weight%*%t(weight)/norm(weight)/norm(weight)
Xirr = X - Xrel

# Task-relevant differences
DistRel = sqrt(sum((apply(Xrel[fzero,], 2, mean) - apply(Xrel[fone,], 2, mean))**2))

# Task-relevant variabilities
VarRel  = mean(0.5*(apply(Xrel[fzero,], 2, sd) + apply(Xrel[fone,], 2, sd)))