Add GLM read write support

[27-apizzuti]

Hi Faruk, it would be great to have a read/write .glm function.
Thank you in advance,
Alessandra

[ofgulban]

Initial implementation is done, "read_GLM" function is working. See an example script (read_glm_export_nifti.py) at wip folder.

However, there are 2 remaining issues:

1. There are 9 unknown bytes associated with the predictors. I can not find documentation on what these bytes are stand for. Therefore, I am currently reading and skipping them. See

   bvbabel/bvbabel/glm.py

   Lines 211 to 212 in 6fc13b8

   	# TODO: Unknown bytes, ask to senior dev. for documentation
   	f.read(9)

2. The docstring descriptions needs to be double checked. See

   bvbabel/bvbabel/glm.py

   Lines 18 to 67 in 6fc13b8

   	Returns
   	-------
   	header : dictionary
   	Pre-data headers.
   	data_R2 : 3D (if volume) OR 1D (if vertices) numpy.array
   	Multiple correlation coefficient R indicating the goodness-of-fit for
   	the respective voxel's time course and to allow to calculate the
   	proportion of explained (R^2) and unexplained (1 - R^2) variance.
   	data_SS : 3D (if volume) OR 1D (if vertices) numpy.array
   	Sum-of-squares term (SS_total) that can be used together with the R
   	value to calculate the variance of the residuals.
   	data_beta : 4D (if volume) OR 2D (if vertices) numpy.array
   	Estimated beta values. One value for each predictor of the design
   	matrix ("Nr all predictors" values).
   	data_XY = 4D (if volume) OR 2D (if vertices) numpy.array
   	Sum-of-squares indicating the covariation of each predictor with the
   	time course data (SS_XiY). These values are stored to allow easy
   	calculation of explained variance terms for restricted models (i.e. to
   	allow application of the extra-sum-of-squares principle); these values
   	may be ignored (not stored) for custom processing.
   	data_ARlag : 3D or 4D (if volume) OR 1D or 2D (if vertices) numpy.array
   	Auto-regression lag value. If "serial correlation" is 1, this will be a
   	3D numpy.array (if volume) or 1D numpy.array (if vertices). If "serial
   	correlation" is 2, this will be a 4D numpy.array (if volume) or 2D
   	numpy.array (if vertices). If "serial correlation" is zero, this will
   	be a 3D numpy.array (if volume) or 1D numpy.array (if vertices)
   	containing all zeros (can be ignored)
   	Notes
   	-----
   	- NOTE[Faruk]: "Developer Guide - The Format Of GLM Files (v4)" has
   	interesting details about calculating standard errors for beta and contrast
   	values. In order to retain this extra information, I am including these
   	notes below.
   	- The (non-RFX) GLM file stores enough values to allow calculation of
   	standard errors for beta and contrast values for each voxel. The stored
   	multiple correlation coefficient R (data_R2) together with the overall
   	sum-of-squares term (data_SS) can be used to calculate the variance of
   	the residuals as follows:
   	`VAR_residuals = data_SS * (1 - data_R2) / (header["Nr time points"] - header["Nr all predictors"])`
   	Together with the stored inverted X'X matrix, this allows calculating the
   	standard error for any beta or contrast t value using the usual equation
   	(c is the contrast vector and b is the voxel's vector of stored beta
   	values):
   	`t = c'b / sqrt(VAR_residuals * c' * header["Inverted X'X matrix"] * c)`
   	However, note that in case of performed serial correlation correction, the
   	inverted X'X matrix needs to be recalculated for each voxel from the stored
   	design matrix X using the voxel-specific autocorrelation function term(s);
   	furthermore the number of time points (NTimePoints) needs to be corrected
   	[subtraction of 1 for AR(1) model, subtraction of 2 for AR(2) model].

   
   @nausikaa8 , would you have time to help me with the second item?

[nausikaa8]

Sure, @ofgulban, I'll get back within a few days.

[nausikaa8]

Hi @ofgulban, I've emailed you with some initial comments.