Towards Multivariable Multivariate GLM

Repeated-measures GLM is also available to describe fMRI studies; the so-called growth curve model could be the best way of conceptually modelling a multi-subject experiment. It has the great advantage of being very close to the two-stage understanding of the random subject analysis, and of being clearly explicit for a g groups analysis. A multivariate GLM model is usually written like the general form we have seen before but $y=Y$, $\beta$ and $\epsilon$ are matrices. A multivariate linear model can be in fact rewritten as a classical GLM and for example resolved this way as well. The growth curve model is an extension of the multivariate GLM to allow a model of the curve response. It is in fact also equivalent to a classical GLM considering a tensor product of the designs:

$$\begin{array}{c } Y\\ N \times T \\ \end{array} = \... ...\begin{array}{c } \epsilon\\ N \times T \\ \end{array}$$ (43)

For a 1-group analysis $N=n$ and $X_N$ reduces to the ``constant" $1_n$ (or $I\!d_n$ for the fixed model), otherwise $N=\sum_k n_k$ and $X_N$ identifies the groups. Notice that $X_T$ is in fact the design we had before (then denoted $X$) to identify the paradigm and covariates (within the time series). The model (46) is in fact equivalent to the univariate model (where underlined means all of the columns stacked, $i.e.$ $NT \times 1$ or $pm \times 1$ vectors)

$$\underline{Y}=(X_N \otimes X_T) \underline{\beta}+ \underline{\epsilon}$$ (44)

from which one recognises the model (35) if one has only one group. Some algebra gives the result for the least squares estimate:

$$(X_N\otimes X_T)\underline{\hat{\beta}}=(P_{X_N} \otimes P_{X_T})\underline{Y} \equiv P_{X_N} Y {\;}^tP_{X_T}$$ (45)

(where $\equiv$ means algebraically equivalent, $i.e.$ a vector $NT \times 1$ or a matrix $N \times T$ of the same thing) ²⁸.
The normal distributional assumption comes with the restriction of ``separability'', $i.e.$ of the form $var(\underline{\epsilon})=V_N \otimes V_T$, with usually the form $var(\underline{\epsilon})=(I\!d_N \otimes \sigma^2 \Sigma)$ where $\Sigma$ is a $T \times T$ matrix expressing the common (to all rows) correlation of the columns (random variables) of $Y$. Historically presented by Potthof and Roy(1964) (see reference in [14]), the model (46) was named GMANOVA (Growth curve Multivariate ANalysis Of VAriance) and was ``solved'' with a two-level analysis: projecting first according to $X_T$ and then solving the second level model. Let the projector onto $X_T$ expressed with a metric $Q_0$, $P_{X_T}=X_T({\;}^tX_TQ_0X_T)^{-1}{\;}^tX_TQ_0$, then writing $Y{\;}^tP_{X_T}=Y_0{\;}^tX_T$ one obtains $E(Y_0)=X_N\beta$ as MANOVA model. A known statistical result is that, under normality, the least squares estimate of $\beta$ in the MANOVA model, which is in fact $\hat{\beta}$ (48), is also the maximum likelihood of $\beta$ in the MANOVA model for the choice of $Q_0=\widehat{\sigma^2\Sigma}^{-1}=({\;}^tYP_{X_N\bot}Y/(rank(P_{X_N\bot}))^{-1}$. In fMRI this estimation is not sensible as $N$ is small comparatively to $T$, one would prefer time series methods (as mentioned before) to estimate the correlation structure $\Sigma$. Estimating $\sigma^2$ brings back to the choices between random, fixed subject analysis or the ``pooling model errors" (43). In the random analysis as in the usual Growth curve analysis $\sigma^2$ relates to common variance (across subject population) for every time measurement, the estimate $1/T [trace(Q_0)]$ could be used i.e. the pooling (across time) of residual variances of the subject model ( $y_t= X_N\beta_t +\epsilon_t$ at t fixed), or could be better estimated from the MANOVA model so with $1/m [trace(S_0)]$ (where $S_0={\;}^tY_0P_{X_N\bot}Y_0/ (N-rank(X_N))$). Notice (if $X_N=1_N$) the diagonal of $S_0$ contains the variances of the parameters ²⁹ used for the random subject analysis (9). When $X_N=I\!d_N$ as in the fixed subject analysis, the previous estimate does not make any sense (it is zero) and the natural estimate is used:

$$\hat{\sigma}^2=\frac{{\;}^t\underline{Y}[ {\;}^t(I\!d_{NT}-P... ...\!d_{NT}-P_{X_N} \otimes P_{X_T})I\!d_N\otimes \hat{\Sigma})}$$ (46)

with $X_N=I\!d_N$ and is the pooled subject errors with autocorrelation taken into account. This last estimate is also the ``pooling model errors" estimate of $\sigma^2$ when $X_N\neq I\!d_N$ (basically $rank(X_N)<N$). Model (47) could come with assumption $var(\underline{\epsilon})=(I\!d_N \otimes D_{\sigma^2} \Sigma)$ and an estimation related to the later comment would be $D_{\sigma^2}=diag(Q_0)$ or $D_{\sigma^2}=diag(X_TS_0{\;}^tX_T)=diag(P_{X_T}Q_0{\;}^tP_{X_T})$. To retrieve completely the mixed model, one would have to write the covariance with the form $I\!d_N \otimes \sigma^2\Sigma=I\!d_N \otimes (\sigma_\epsilon^2C + \sigma_s^2J_T)$, but after all it is not necessary in fMRI studies to split or separate the covariance structure as the interest is only on the fixed effects and not actually on the variance components. The Lawley-Hotelling statistic already described (21) can be used and a general hypothesis takes the form $L \beta M=0$ or $L \otimes M \underline{\beta}=0$ where $L$ will contrast the groups and $M$ contrasts the paradigm. Mathematically it is possible to extend further this representation to any ``dimension'' in the data, for example incorporating the spatial dimension:

$$\underline{Y}=(X_B \otimes X_N \otimes X_T) \underline{\beta}+ \underline{\epsilon}$$ (47)

where $X_B$ would relate a model on the voxels (or $X_B=I\!d_v$ if spatial modelling is not carried out, where $v$ is the number of voxels), $Y$ can be interpreted as a tensor of order 3 with a $(vNT) \times 1$ vector representation. To be useful this model would have to take into account the spatial covariance structure. Remark: Note that if the groups constitute a repeated experiment ($e.g.$ with different doses of a drug), $V_N$ can be chosen to be for the form $V_N=I\!d_n \otimes \sigma_g^2 J_g$ as in compound symmetry, and similar techniques described for mixed models could be derived. This can be also implemented as a three-stage model: seeing the random model as a GMANOVA model then resolved as a two-stage model allow you to carry on levels in presence of a repeated design on the subjects giving a third model. The estimation in multi-stage for multilevel models is valid as long as the designs are well balanced which is the case in fMRI. Iterative GLS ³⁰ may be needed to improve estimation. Note a dense literature on multilevel modelling allows considerations of more complex situations [9].