Random subject effect

Firstly, notice that the random subject effect model could use this same model with a covariance structure (block-diagonal) which would need to be estimated in the first place. One would get a EGLS but this model cannot be used to estimate the covariance matrix (this is because writing this model as a mixed model makes $\beta$ random and fixed at the same time). One has to rewrite the model, for example with the form:

$$\left ( \begin{array}{c } y_1 \\ y_2 \\ \vdots\\ ... ...silon'_2 \\ \vdots\\ \epsilon'_n \end{array} \right )$$ (35)

or

$$y=(1_n \otimes X) \beta +(I\!d_n \otimes 1_T) \gamma +\epsilon'$$

where $Z=(I\!d_n \otimes 1_T)$ identifying the subjects (dummy variables)²⁶ , with also $E(\gamma)=0$ and $var(\gamma)= G$ as random effects. One can easily see that OLS gives $\hat{\beta}=\overline{\hat{\beta_{{(i)}}}}_{[1..n]}$. Depending on the structure of covariance chosen an estimate will be found directly ($e.g.$ with compound symmetry) or algorithmically (for general REML estimation). For the random analysis with simple compound symmetry (i.e. with $var(y)=var(\epsilon')+ZG{\;}^tZ=\sigma_\epsilon^2I\!d_{nT} + \sigma_s^2I\!d_n\otimes J_T=I\!d_n\otimes(\sigma_\epsilon^2I\!d_T+\sigma_s^2J_T)$, where $\sigma_\epsilon^2$ is the scan-to-scan variation, supposed to be the same for all the subjects, and $\sigma_s^2$ is the random subject variation²⁷ introducing correlations into the errors) then with estimation for $\beta$ as OLS (here equivalent to GLS) ``ANOVA" estimators (equivalent to REML) can be given by:

$$\hat{\sigma}_\epsilon^2=MSE={}^tyP_{[1_n \otimes X,I\!d_n \otimes 1_T ]^\bot}y/[n(T-1) -rank(X)+1]$$ (36)

and

$$\hat{\sigma}_s^2=\frac{MSs-MSE}{T} \mbox{ where }MSs={}^ty... ... \otimes X,I\!d_n \otimes 1_T ]}-P_{[1_n \otimes X]})y/(n-1)$$ (37)

Our interest here is only on the subject error $\hat{\sigma}_\epsilon^2+(T\hat{\sigma}_s^2)=MSs$ which is the two-stage approach estimation seen before (obviously dividing by $T$ to go on second level).
If all parameters in $X$ are considered as random, making $Z=(I\!d_n \otimes X)$, the estimation can be done in the same way with: $var(y)=var(\epsilon')+ZG{\;}^tZ=I\!d_n\otimes(\sigma_\epsilon^2I\!d_T+\sum_{j=1}^p\sigma_j^2X_j{\;}^tX_j)$ and considering the random parameters independent (non-correlated). Using Hendenson's Method [12,3] 36 and 37 (for each random parameters) become:

$$\hat{\sigma}_\epsilon^2=MSE={}^tyP_{[1_n \otimes X,I\!d_n \otimes X]^\bot}y/[n(T-rank(X))]$$ (38)

and

$$\hat{\sigma}_j^2=\frac{SSs -MSE {\;}(n-1)}{trace((P_{[1_n ... ...otimes X,I\!d_n \otimes X_{(-j)}]}) Id_n \otimes X_j{}^tX_j)}$$ (39)

where $X_{(-j)}$ is $X$ without the column $X_j$ and

$$\nonumber SSs={}^ty(P_{[1_n \otimes X,I\!d_n \otimes X]}-P_{[1_n \otimes X,I\!d_n \otimes X_{(-j)}]})y \nonumber$$

The error considered for a parameter is similarly $\hat{\sigma}_\epsilon^2+(\theta \hat{\sigma}_j^2)=SSs/(n-1)=MSs$, with $\theta = trace((P_{[1_n \otimes X,I\!d_n \otimes X]}-P_{[1_n \otimes X,I\!d_n \otimes X_{(-j)}]}) Id_n \otimes X_j{}^tX_j)/(n-1)$ equal to $T$ under orthogonality of $X_j$'s. For general correlation structure, or one derived from time series analysis, instead of using an REML estimation algorithm, a ``two-stage strategy'' could be to estimate the autocorrelation matrix $C_i$ from every subject's time series, then pool these to give an estimate of the common auto-correlation $C$, to be able to define the covariance structure of the form $V=I\!d_n \otimes \sigma^2_s J_T +I\!d_n \otimes \sigma_w^2 \hat{C}$, then solve REML (or Henderson's method III [12]) for $\sigma^2_w$ and $\sigma^2_s$. Notice a redundancy in the previous formula as the covariance structure describes at the same time a constant covariance for the same subject ($\sigma^2_s$) and the auto-covariance coming from time series autocorrelation; this suggests a better covariance model of the form $V=\sigma^2I\!d_n \otimes \hat{C}$ which then take subject the variance component off the model (no need of $Z$ as well). At this point, time series methods are used to estimate the $C_i$ in the first place, and robust estimation such as in Woolrich et al.[16] including the non-linear spatial smoothing would improve the result.
Remark: The model (35) without $Z$ could be considered as an intermediary model (Pme) between the fixed model (Fix) ( the model (35) but $var(\epsilon')=\sigma^2I\!d_{nT}$), i.e. $var(\gamma)=0$) and the random subject model (Ran) as it can be shown that the natural estimate of $\sigma^2$ is the pooled error variance of the fixed version ($MSE$ given in (36)) of and random approaches (subject error $MSs$ (37), but because of the $df$ differences between the two approaches this will underestimate the error variance. This method called here``pooling model errors" or Pme may be interesting when too few subjects are in the study : the error variance is going to be larger than for a fixed analysis and the increase of degrees of freedom ($(n-1)$) not too large:

$$\mbox{model Pme: } y = (1_n \otimes X)\beta + \epsilon_{\mbox{\tiny {Pme}}}$$

$$= \mathcal{X} \beta +\epsilon_{\mbox{\tiny {Pme}}}$$

$$\mbox{model Fix: } y = (\mathcal{X} \quad Z )B + \epsilon_{\mbox{\tiny {Fix}}}$$

$$= F B + \epsilon_{\mbox{\tiny {Fix}}}$$

$$\mbox{model Ran: } y = \mathcal{X}\beta + Z \gamma +\epsilon_{\mbox{\tiny {Ran}}} \mbox{ and }var(\gamma)\neq 0$$

then from model Fix ,

$$df_{\mbox{\tiny {Fix}}}MSE_{\mbox{\tiny {Fix}}} = {\;}^ty(Id -P_{F} y$$

$$= {\;}^ty(Id -P_\mathcal{X} +P_\mathcal{X} -P_{F} y$$

$$= df_{\mbox{\tiny {Pme}}}MSE_{\mbox{\tiny {Pme}}} -(n-1)MSs$$

thereby

$$\widehat{\sigma}_{\mbox{\tiny {Pme}}}^2= MSE_{\mbox{\tiny... ...MSE_{\mbox{\tiny {Fix}}}+(n-1)MSs }{df_{\mbox{\tiny {Pme}}}}$$ (40)

and one can easily check $df_{\mbox{\tiny {Pme}}}=df_{\mbox{\tiny {Fix}}} +(n-1)=nT-rank(X)$.