From bb6e10a8096fe40a275a7776b2a4f63ba13089a1 Mon Sep 17 00:00:00 2001
From: Paul McCarthy <pauld.mccarthy@gmail.com>
Date: Fri, 4 Nov 2016 14:13:48 +0000
Subject: [PATCH] Further tweaks to FEAT model fit - trying to get visually
sensible scaling on partial model fits. Also added new method
isFirstLevelAnalysis.
---
fsl/data/featanalysis.py | 10 ++++++
fsl/data/featimage.py | 70 +++++++++++++++++++++++++--------------
fsl/data/melodiclabels.py | 3 +-
3 files changed, 58 insertions(+), 25 deletions(-)
diff --git a/fsl/data/featanalysis.py b/fsl/data/featanalysis.py
index 63cbcc9d0..8aa24c4a7 100644
--- a/fsl/data/featanalysis.py
+++ b/fsl/data/featanalysis.py
@@ -19,6 +19,7 @@ following functions are provided:
hasMelodicDir
getAnalysisDir
getTopLevelAnalysisDir
+ isFirstLevelAnalysis
loadDesign
loadContrasts
loadSettings
@@ -279,6 +280,15 @@ def getThresholds(settings):
}
+def isFirstLevelAnalysis(settings):
+ """Returns ``True`` if the FEAT analysis described by ``settings``
+ is a first level analysis, ``False`` otherwise.
+
+ :arg settings: A FEAT settings dictionary (see :func:`loadSettings`).
+ """
+ return settings['level'] == '1'
+
+
def loadClusterResults(featdir, settings, contrast):
"""If cluster thresholding was used in the FEAT analysis, this function
will load and return the cluster results for the specified (0-indexed)
diff --git a/fsl/data/featimage.py b/fsl/data/featimage.py
index d9fa82cf4..d5857d50a 100644
--- a/fsl/data/featimage.py
+++ b/fsl/data/featimage.py
@@ -45,7 +45,7 @@ class FEATImage(fslimage.Image):
# [23, 30, 42] (in this example, we
# have 4 EVs - the first argument
# is a contrast vector).
- img.fit([1, 1, 1, 1], [23, 30, 42], fullmodel=True)
+ img.fit([1, 1, 1, 1], [23, 30, 42], fullModel=True)
"""
@@ -111,6 +111,13 @@ class FEATImage(fslimage.Image):
return self.__analysisName
+ def isFirstLevelAnalysis(self):
+ """Returns ``True`` if the FEAT analysis described by ``settings``
+ is a first level analysis, ``False`` otherwise.
+ """
+ return featanalysis.isFirstLevelAnalysis(self.__settings)
+
+
def getTopLevelAnalysisDir(self):
"""Returns the path to the higher level analysis directory of
which this FEAT analysis is a part, or ``None`` if this analysis
@@ -287,12 +294,12 @@ class FEATImage(fslimage.Image):
return self.__clustMasks[con]
- def fit(self, contrast, xyz, fullmodel=False):
+ def fit(self, contrast, xyz, fullModel=False):
"""Calculates the model fit for the given contrast vector
at the given voxel.
- Passing in a contrast of all 1s, and ``fullmodel=True`` will
- get you the full model fit. Pass in ``fullmodel=False`` for
+ Passing in a contrast of all 1s, and ``fullModel=True`` will
+ get you the full model fit. Pass in ``fullModel=False`` for
all other contrasts, otherwise the model fit values will not
be scaled correctly.
@@ -302,13 +309,20 @@ class FEATImage(fslimage.Image):
:arg xyz: Coordinates of the voxel to calculate the model fit
for.
- :arg fullmodel: Set to ``True`` for a full model fit, ``False``
+ :arg fullModel: Set to ``True`` for a full model fit, ``False``
otherwise.
"""
if self.__design is None:
raise RuntimeError('No design')
+ x, y, z = xyz
+ firstLevel = self.isFirstLevelAnalysis()
+ numEVs = self.numEVs()
+
+ if len(contrast) != numEVs:
+ raise ValueError('Contrast is wrong length')
+
# Here we are basically trying to
# replicate the behaviour of tsplot.
# There are some differences though -
@@ -316,19 +330,21 @@ class FEATImage(fslimage.Image):
# data by Z statistics. We're not
# doing that here.
- # If we are not plotting a full model fit,
- # normalise the contrast vector to unit
- # length. Why? Because that's what tsplot
- # does.
- if not fullmodel:
- contrast = np.array(contrast)
- contrast = contrast / np.sqrt((contrast ** 2).sum())
-
- x, y, z = xyz
- numEVs = self.numEVs()
-
- if len(contrast) != numEVs:
- raise ValueError('Contrast is wrong length')
+ # Normalise the contrast vector.
+ # The scaling factor is arbitrary,
+ # but should result in a visually
+ # sensible scaling of the model fit.
+ # For a vector of all 1's (i.e. a
+ # full model fit) this is a no-op.
+ #
+ # We also take the absolute value
+ # of the values in the contrast
+ # vector, as the parameter estimates
+ # should be appropriately signed.
+ contrast = np.array(contrast)
+ nonzero = sum(~np.isclose(contrast, 0))
+ contrast = contrast / np.sqrt((contrast ** 2).sum())
+ contrast = np.abs(contrast * np.sqrt(nonzero))
X = self.__design.getDesign(xyz)
data = self[x, y, z, :]
@@ -340,13 +356,19 @@ class FEATImage(fslimage.Image):
pe = self.getPE(i)[x, y, z]
modelfit += ev * pe * contrast[i]
- # Make sure the full model fit
- # has the same mean as the data
- if fullmodel: return modelfit - modelfit.mean() + data.mean()
- else: return modelfit
+ # Make sure the model fit has an appropriate mean.
+ #
+ # Full model fits should, but won't necessarily,
+ # have the same mean as the data.
+ #
+ # The data in first level analyses is demeaned before
+ # model fitting, so we need to add it back in.
+ if fullModel: return modelfit - modelfit.mean() + data.mean()
+ elif firstLevel: return modelfit - modelfit.mean() + data.mean()
+ else: return modelfit
- def partialFit(self, contrast, xyz, fullmodel=False):
+ def partialFit(self, contrast, xyz, fullModel=False):
"""Calculates and returns the partial model fit for the specified
contrast vector at the specified voxel.
@@ -355,6 +377,6 @@ class FEATImage(fslimage.Image):
x, y, z = xyz
residuals = self.getResiduals()[x, y, z, :]
- modelfit = self.fit(contrast, xyz, fullmodel)
+ modelfit = self.fit(contrast, xyz, fullModel=fullModel)
return residuals + modelfit
diff --git a/fsl/data/melodiclabels.py b/fsl/data/melodiclabels.py
index 241827914..1b09858a8 100644
--- a/fsl/data/melodiclabels.py
+++ b/fsl/data/melodiclabels.py
@@ -320,7 +320,8 @@ def loadLabelFile(filename, includeLabel=None, excludeLabel=None):
[2, 5, 6, 7]
The square brackets may or may not be present, i.e. the
- following format is also accepted::
+ following format is also accepted (this format is generated
+ by ICA-AROMA)::
2, 5, 6, 7
--
GitLab