No parallel computing under Windows

[Ede1994]

I've tried the following lien of code:

 stats <- lme_mass_fit_vw(X, Zcols, Y, ni, numcore=8)

and get this error:

Error in parallel::mclapply(c(1:nv), build_the_model, Y, X, Xcols, Zcols,  : 
  'mc.cores' > 1 is not supported under Windows

If I omit numcores, i.e. only work with one core (without parallel computing), then the code runs, but takes a very long time. Is there a way to use parallel computing under Windows as well? Or what other alternatives are there?

PS: I am not that experienced in R, only in Python.

[kdiers]

Hello, thanks for your report. I am afraid that we cannot support Windows systems currently; one option may be to try the original Matlab implementation: https://surfer.nmr.mgh.harvard.edu/fswiki/LinearMixedEffectsModels if that is an alternative. Sorry to not have any better new at the moment.

[Ede1994]

Okay, thanks for the feedback. I thought maybe you just need to install settings or another package to be able to use the function.

If the whole thing depends on the operating system, would it run under Linux (Ubuntu)?

[kdiers]

Yes, linux / ubuntu should work.

[Ede1994]

I've installed it under WSL2 in RStudio Server and now it runs with more cores. Thx.
I now have a different problem or an error message.

This is my code so far:

# load the package
library(fslmer)

# read thickness file
lh.thickness <- lme_openmgh("/mnt/d/FastSurfer/lh.thickness_sm10.mgh")

# read template surface
lh.sphere <- lme_readsurf("/mnt/d/FastSurfer/raw/fsaverage6/surf/lh.sphere")

# read cortical label file
lh.cortex <- lme_readlabel("/mnt/d/FastSurfer/raw/fsaverage6/label/lh.cortex.label")[,1]

# read qdec table
dat <- read.csv("/mnt/d/FastSurfer/qdec.long.csv", header=TRUE, sep = ";")

# order qdec table based on fsid.base and time
dat <- dat[order(dat$fsid.base, dat$days), ]

# create vector of number of observations per subject
ni <- matrix(unname(table(dat$fsid.base)), ncol=1)

# create nSubjects-times-nVertices matrix for thickness data
Y <- t(drop(lh.thickness$x))

# create (1-based) vector of cortical label indices (to exclude non-cortical vertices)
maskvtx <- sort(lh.cortex)+1

# create new column 'dX' in 'dat'
dat$DX <- factor(dat$days, 
                 levels = c(0,1,4), 
                 labels = c("day0","day1","day4"))


# Creating the design matrix
X <- model.matrix(~DX, dat)

# create contrasts
# C <- matrix(c(0, 0, 0, 1), nrow = 1)
C <- rbind("day1-day0" = c(0,  1,  0),
           "day4-day0" = c(0,  0,  1),
           "day4-day1" = c(0, -1,  1))

# determine the type of mixed-effects model (assuming the intercept is column 1 
# in the design matrix and the time variable is column 2).
Zcols <- 1 # random-intercept 
# Zcols <- c(1, 2) # random-slope, random-intercept

# Estimate the model
stats <- lme_mass_fit_vw(X, Zcols, Y, ni, numcore=6)

# obtain initial estimates
FitInit <- lme_mass_fit_init(X=X, Zcols=Zcols, Y, ni, maskvtx=maskvtx, numcore=6)

# run algorithm to identify spatially homogeneous regions
RgGrow <- lme_mass_RgGrow(lh.sphere, FitInit$Re0, FitInit$Theta0, maskvtx=maskvtx, nst=2, prc=95)

# fit model
FitRgw <- lme_mass_fit_Rgw(X, Zcols, Y, ni, FitInit$Theta0, RgGrow$Regions, lh.sphere, prs=6)

# optional: random-effects coefficients
# rfx <- lme_mass_rfx(FitRgw$stats, X, Zcols, Y, ni, maskvtx, prs=6)

# inference
F_C <- lme_mass_F(FitRgw$stats, C)

And I get this error:

> F_C <- lme_mass_F(FitRgw$stats, C)
Error in solve.default((C %*% CBhat %*% t(C))) : 
  system is computationally singular: reciprocal condition number = 4.50134e-17

My data (dat) looks like this:

fsid  | fsid.base    | days 
xx_00 | xx           | 0 
xx_01 | xx           | 1 
xx_02 | xx           | 4 
yy_00 | yy           | 0 
yy_01 | yy           | 1
yy_02 | yy           | 4

The data are balanced and there are no missing data or so.
For me, the differences from day=0 to day=1, day=1 to day=4 and day=0 to day=4 are particularly interesting. Is my contrast possibly defined incorrectly?

Maybe I didn't understand something about the example in the README.md correctly...

[kdiers]

Yes, I think your suspicion about the contrast is correct, because contrast 4-1 can be represented as a linear combination of contrasts 1-0 and 4-0, and therefore is redundant and cannot be combined into a single F-contrast. What should work is to have three separate contrasts, i.e. 1-0, 4-0, and 4-1, and run the corresponding code 3 times.

[Ede1994]

Thx for your help. Now my code runs without errors. So I just compared 0-1:

library(fslmer)

# read thickness file
lh.thickness <- lme_openmgh("/mnt/d/FastSurfer/lh.thickness_sm10.mgh")

# read template surface
lh.sphere <- lme_readsurf("/mnt/d/FastSurfer/raw/fsaverage6/surf/lh.sphere")

# read cortical label file
lh.cortex <- lme_readlabel("/mnt/d/FastSurfer/raw/fsaverage6/label/lh.cortex.label")[,1]

# read qdec table
dat <- read.csv("/mnt/d/SleepyQSM/FastSurfer/qdec.long.csv",
                header=TRUE, sep = ";")
#dat <- dat[order(dat$fsid.base, dat$days), ]

# create nSubjects x nVertices thickness matrix
Y <- t(drop(lh.thickness$x))
# 1-based cortical label index
maskvtx <- sort(lh.cortex) + 1

# keep days
days_keep=c(0,1)

# create data subset
dat_sub <- subset(dat, days %in% days_keep)
# dat_sub <- dat_sub[order(dat_sub$fsid.base, dat_sub$days), ]
  
# create nSubjects-times-nVertices matrix for thickness data
Y_sub <- Y[dat$days %in% days_keep, ]
  
# create vector of number of observations per subject
ni_sub <- matrix(unname(table(dat_sub$fsid.base)), ncol=1)
  
# create new column 'dX' in 'dat' (labels)
day_labels <- paste0("day", days_keep)
dat_sub$DX <- factor(dat_sub$days,
                     levels = days_keep,
                     labels = day_labels)
  
# Creating the design matrix
X <- model.matrix(~DX, dat_sub)
  
# create contrasts "dayX - dayY" => (0 1) test effects of "DX != baseline"
C <- matrix(c(0, 1), nrow = 1)
  
# random-intercept 
Zcols <- 1
  
# Estimate the model
stats <- lme_mass_fit_vw(X, Zcols, Y_sub, ni_sub, numcore=8)
  
# obtain initial estimates
init <- lme_mass_fit_init(X, Zcols, Y_sub, ni_sub,
                          maskvtx=maskvtx, numcore=8)
  
# run algorithm to identify spatially homogeneous regions
RgGrow <- lme_mass_RgGrow(lh.sphere, init$Re0, init$Theta0,
                          maskvtx=maskvtx, nst=2, prc=95)
  
# fit model
FitRgw <- lme_mass_fit_Rgw(X, Zcols, Y_sub, ni_sub,
                           init$Theta0, RgGrow$Regions,
                           lh.sphere, prs=8)
  
# inference
F_C <- lme_mass_F(FitRgw$stats, C)

# multiple comparison correction using classical FDR
thr_pFDR  <- lme_mass_FDR(F_C$pval, 0.05)

# multiple comparison correction using two-stage approach
FDR2_C <- lme_mass_FDR2(F_C$pval, F_C$sgn)

# save F-values
vol <- NULL
vol$ndim1 <- length(F_C$F)
vol$ndim2 <- 1
vol$ndim3 <- 1
vol$nframes <- 1
vol$x <- array(data=F_C$F, dim=c(length(F_C$F), 1, 1, 1))
lme_savemgh(vol=vol, fname="/mnt/d/FastSurfer/lh.thickness_sm10.f_val.mgh")

# save uncorrected p-values
vol <- NULL
vol$ndim1 <- length(F_C$pval)
vol$ndim2 <- 1
vol$ndim3 <- 1
vol$nframes <- 1
vol$x <- array(data=F_C$pval, dim=c(length(F_C$pval), 1, 1, 1))
lme_savemgh(vol=vol, fname="/mnt/d/FastSurfer/lh.thickness_sm10.uncorr_p_val.mgh")

# save FDR2-corrected p-values
vol <- NULL
vol$ndim1 <- length(FDR2_C$sided_pval)
vol$ndim2 <- 1
vol$ndim3 <- 1
vol$nframes <- 1
vol$x <- array(data=FDR2_C$sided_pval, dim=c(length(FDR2_C$sided_pval), 1, 1, 1))
lme_savemgh(vol=vol, fname="/mnt/d/FastSurfer/lh.thickness_sm10.corr_p_val.mgh")

Hopefully one last problem/question 😄
I now have the .mgh file and tried the following:

freeview \
  -f /mnt/d/FastSurfer/raw/fsaverage6/surf/lh.inflated \
       :overlay=/mnt/d/FastSurfer/lh.thickness_sm10.f_val.mgh \
       :overlay_threshold=2,5

But then I get the error:

Did not find any volume info
error: No such file or directory
error: MRISread(): could not open file
MRISread failed

I've checked the path twice, so there is no error. Freeview shows just the fsaverage6 file.
In general, I lack the idea of how to visualise the data and ‘recognise’ the results.

[kdiers]

Can you open the surface file and the mgh overlay file in the Freeview GUI directly, i.e. without using the command line? I could also imagine (but have not checked yet) that there is a problem with whitespace characters in the command line arguments; maybe try to put everything on one line and no whitespace characters before or after the colons.

[Ede1994]

Unfortunately, opening via the Freeview GUI does not work either. I first loaded lh.inflated from fsaverage6 (worked) and when loading one of the created files Freeview crashed with the following message in the shell:

Did not find any volume info
Segmentation fault (core dumped)

I've checked with mri_info one of the results:

>> mri_info /mnt/d/FastSurfer/lh.thickness_sm10.f_val.mgh
Volume information for /mnt/d/FastSurfer/lh.thickness_sm10.f_val.mgh
          type: MGH
    dimensions: 40962 x 1 x 1
   voxel sizes: 1.000000, 1.000000, 1.000000
          type: FLOAT (3)
           fov: 40962.000
           dof: 0
        xstart: -20481.0, xend: 20481.0
        ystart: -0.5, yend: 0.5
        zstart: -0.5, zend: 0.5
            TR: 0.00 msec, TE: 0.00 msec, TI: 0.00 msec, flip angle: 0.00 degrees
       nframes: 1
       PhEncDir: UNKNOWN
       FieldStrength: 0.000000
ras xform not present
    xform info: x_r =  -1.0000, y_r =   0.0000, z_r =   0.0000, c_r =     0.0000
              : x_a =   0.0000, y_a =   0.0000, z_a =   1.0000, c_a =     0.0000
              : x_s =   0.0000, y_s =  -1.0000, z_s =   0.0000, c_s =     0.0000

talairach xfm :
Orientation   : ???
Primary Slice Direction: unknown

voxel to ras transform:
               -1.0000   0.0000   0.0000 20481.0000
                0.0000   0.0000   1.0000    -0.5000
                0.0000  -1.0000   0.0000     0.5000
                0.0000   0.0000   0.0000     1.0000

voxel-to-ras determinant -1

ras to voxel transform:
               -1.0000  -0.0000  -0.0000 20481.0000
               -0.0000  -0.0000  -1.0000     0.5000
               -0.0000   1.0000  -0.0000     0.5000
               -0.0000  -0.0000  -0.0000     1.0000

So there are no informations about orientation and direction in the header.

EDIT:
Okay, I've found the problem. Since my operating system is Windows and I installed FreeSurfer etc. in WSL2, Freeview could not work properly with the /mnt/d... path and this caused the error. When I copy the files into my WSL environment and then open them with Freeview it loads both files and now ‘only’ the error occurs:

Did not find any volume info
NaN found at voxel (39682, 0, 0, 0)
NaN found at voxel (28162, 0, 0, 0)
NaN found at voxel (3855, 0, 0, 0)
NaN found at voxel (29450, 0, 0, 0)
NaN found at voxel (32067, 0, 0, 0)
NaN found at voxel (9022, 0, 0, 0)
NaN found at voxel (5354, 0, 0, 0)
NaN found at voxel (24441, 0, 0, 0)
NaN found at voxel (26977, 0, 0, 0)
NaN found at voxel (26179, 0, 0, 0)
NaN found at voxel (35699, 0, 0, 0)
NaN found at voxel (14385, 0, 0, 0)
NaN found at voxel (1304, 0, 0, 0)
NaN found at voxel (2597, 0, 0, 0)
NaN found at voxel (2597, 0, 0, 0)

error: WARNING: 3486 NaNs found in volume /home/einspaen/Downloads/lh.thickness_sm10.corr_p_val.mgh

But I can see the overlay in Freeview.