Let, \(\mathbf y_i = \mathbf X_i\mathbf\beta + \mathbf Z_i\mathbf b_i+ \mathbf\epsilon_i,\)
where
\(\mathbf y_i\sim N(\mathbf X_i\mathbf\beta, \Sigma_i=\sigma^2\mathbf I_{n_i}+\mathbf Z_i \mathbf G\mathbf Z_i'),\)
\(\mathbf b_i\sim N(\mathbf 0, \mathbf G),\)
\(\mathbf\epsilon_i\sim N(\mathbf 0, \sigma^2\mathbf I_{n_i})\)
\(\mathbf y_i\) is a \(n_i\times 1\) vector of response for \(i^{th}\) individual at \(1,2,\ldots, n_i\) time points, \(\mathbf X_i\) is a \(n_i\times p\) matrix, \(\mathbf \beta\) is a \(p\times 1\) vector of fixed effect parameters, \(\mathbf Z_i\) is a \(n_i\times q\) matrix, \(\mathbf b_i\) is a \(q\times 1\) vector of random effects, \(\mathbf \epsilon_i\) is a \(n_i\times 1\) vector of within errors, \(\mathbf G\) is a \(q\times q\) covariance matrix of between-subject, \(\sigma^2\) is a scalar.
Note that, \(\mathbf X_i\), \(\mathbf Z_i\), and \(\mathbf G\) do NOT involve \(\sigma^2\).
Now I have to find out the Restricted Maximum Likelihood (REML) Estimate of \(\sigma^2\), that is,
\(\hat\sigma^2_R = \frac{1}{N_0-p}\sum_{i=1}^{N}(\mathbf y_i-\mathbf X_i\mathbf\beta)'(\mathbf I_{n_i}+\mathbf Z_i \mathbf G\mathbf Z_i')^{-1}(\mathbf y_i-\mathbf X_i\mathbf\beta),\ldots (1)\)
where \(N_0 = \sum_{i=1}^{N}n_i\).
So first I wrote the Restricted Maximum Log-Likelihood :
\(l_R \propto -\frac{1}{2}\sum_{i=1}^{N}\log\det(\Sigma_i)-\frac{1}{2}\sum_{i=1}^{N}\log\det(\mathbf X_i'\Sigma_i^{-1}\mathbf X_i)-\frac{1}{2}\sum_{i=1}^{N}(\mathbf y_i-\mathbf X_i\mathbf\beta)'\Sigma_i^{-1}(\mathbf y_i-\mathbf X_i\mathbf\beta)\)
Then I have to differentiate log-likelihood, \(l_R\), with respect to \(\sigma^2\) and equate it to zero, i.e.,
\(-\frac{1}{2}\frac{\partial}{\partial\sigma^2}\{\sum_{i=1}^{N}\log\det(\Sigma_i)+\sum_{i=1}^{N}\log\det(\mathbf X_i'\Sigma_i^{-1}\mathbf X_i)\) \(+\sum_{i=1}^{N}(\mathbf y_i-\mathbf X_i\mathbf\beta)'\Sigma_i^{-1}(\mathbf y_i-\mathbf X_i\mathbf\beta)\}|_{\sigma^2=\hat\sigma^2_R}=0\)
But basically I can't differentiate,
\(\frac{\partial}{\partial\sigma^2}\log\det(\Sigma_i)=\frac{\partial}{\partial\sigma^2}\log\det(\sigma^2\mathbf I_{n_i}+\mathbf Z_i \mathbf G\mathbf Z_i')\)
\(\frac{\partial}{\partial\sigma^2}\log\det(\mathbf X_i'\Sigma_i^{-1}\mathbf X_i)=\frac{\partial}{\partial\sigma^2}\log\det(\mathbf X_i'(\sigma^2\mathbf I_{n_i}+\mathbf Z_i \mathbf G\mathbf Z_i')^{-1}\mathbf X_i)\) and
\(\frac{\partial}{\partial\sigma^2}\Sigma_i^{-1}= \frac{\partial}{\partial\sigma^2}(\sigma^2\mathbf I_{n_i}+\mathbf Z_i \mathbf G\mathbf Z_i')^{-1}\).
How can I differentiate the above derivatives and get the REML estimate \(\hat\sigma^2_R\) in equation \((1)\) ?
where
\(\mathbf y_i\sim N(\mathbf X_i\mathbf\beta, \Sigma_i=\sigma^2\mathbf I_{n_i}+\mathbf Z_i \mathbf G\mathbf Z_i'),\)
\(\mathbf b_i\sim N(\mathbf 0, \mathbf G),\)
\(\mathbf\epsilon_i\sim N(\mathbf 0, \sigma^2\mathbf I_{n_i})\)
\(\mathbf y_i\) is a \(n_i\times 1\) vector of response for \(i^{th}\) individual at \(1,2,\ldots, n_i\) time points, \(\mathbf X_i\) is a \(n_i\times p\) matrix, \(\mathbf \beta\) is a \(p\times 1\) vector of fixed effect parameters, \(\mathbf Z_i\) is a \(n_i\times q\) matrix, \(\mathbf b_i\) is a \(q\times 1\) vector of random effects, \(\mathbf \epsilon_i\) is a \(n_i\times 1\) vector of within errors, \(\mathbf G\) is a \(q\times q\) covariance matrix of between-subject, \(\sigma^2\) is a scalar.
Note that, \(\mathbf X_i\), \(\mathbf Z_i\), and \(\mathbf G\) do NOT involve \(\sigma^2\).
Now I have to find out the Restricted Maximum Likelihood (REML) Estimate of \(\sigma^2\), that is,
\(\hat\sigma^2_R = \frac{1}{N_0-p}\sum_{i=1}^{N}(\mathbf y_i-\mathbf X_i\mathbf\beta)'(\mathbf I_{n_i}+\mathbf Z_i \mathbf G\mathbf Z_i')^{-1}(\mathbf y_i-\mathbf X_i\mathbf\beta),\ldots (1)\)
where \(N_0 = \sum_{i=1}^{N}n_i\).
So first I wrote the Restricted Maximum Log-Likelihood :
\(l_R \propto -\frac{1}{2}\sum_{i=1}^{N}\log\det(\Sigma_i)-\frac{1}{2}\sum_{i=1}^{N}\log\det(\mathbf X_i'\Sigma_i^{-1}\mathbf X_i)-\frac{1}{2}\sum_{i=1}^{N}(\mathbf y_i-\mathbf X_i\mathbf\beta)'\Sigma_i^{-1}(\mathbf y_i-\mathbf X_i\mathbf\beta)\)
Then I have to differentiate log-likelihood, \(l_R\), with respect to \(\sigma^2\) and equate it to zero, i.e.,
\(-\frac{1}{2}\frac{\partial}{\partial\sigma^2}\{\sum_{i=1}^{N}\log\det(\Sigma_i)+\sum_{i=1}^{N}\log\det(\mathbf X_i'\Sigma_i^{-1}\mathbf X_i)\) \(+\sum_{i=1}^{N}(\mathbf y_i-\mathbf X_i\mathbf\beta)'\Sigma_i^{-1}(\mathbf y_i-\mathbf X_i\mathbf\beta)\}|_{\sigma^2=\hat\sigma^2_R}=0\)
But basically I can't differentiate,
\(\frac{\partial}{\partial\sigma^2}\log\det(\Sigma_i)=\frac{\partial}{\partial\sigma^2}\log\det(\sigma^2\mathbf I_{n_i}+\mathbf Z_i \mathbf G\mathbf Z_i')\)
\(\frac{\partial}{\partial\sigma^2}\log\det(\mathbf X_i'\Sigma_i^{-1}\mathbf X_i)=\frac{\partial}{\partial\sigma^2}\log\det(\mathbf X_i'(\sigma^2\mathbf I_{n_i}+\mathbf Z_i \mathbf G\mathbf Z_i')^{-1}\mathbf X_i)\) and
\(\frac{\partial}{\partial\sigma^2}\Sigma_i^{-1}= \frac{\partial}{\partial\sigma^2}(\sigma^2\mathbf I_{n_i}+\mathbf Z_i \mathbf G\mathbf Z_i')^{-1}\).
How can I differentiate the above derivatives and get the REML estimate \(\hat\sigma^2_R\) in equation \((1)\) ?