Estimación de Máxima Verosimilitud Restringida (REML) del Componente de Varianza

Question

Deja, $$\mathbf y_i = \mathbf X_i\mathbf\beta + \mathbf Z_i\mathbf b_i+ \mathbf\epsilon_i,$$

donde

$\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$ es un vector de $n_i\times 1$ de respuesta para el individuo $i$ en $1,2,\ldots, n_i$ puntos temporales, $\mathbf X_i$ es una matriz de $n_i\times p$, $\mathbf \beta$ es un vector de $p\times 1$ de parámetros de efecto fijo, $\mathbf Z_i$ es una matriz de $n_i\times q$, $\mathbf b_i$ es un vector de $q\times 1$ de efectos aleatorios, $\mathbf \epsilon_i$ es un vector de $n_i\times 1$ de errores dentro, $\mathbf G$ es una matriz de covarianza de $q\times q$ entre-sujetos, $\sigma^2$ es un escalar.

Nótese que, $\mathbf X_i$, $\mathbf Z_i$, y $\mathbf G$ NO involucran $\sigma^2$.

Ahora tengo que encontrar la Estimación del Máximo Verosimilitud Restringido (REML) de $\sigma^2$, es decir,

$$\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)$$

donde $N_0 = \sum_{i=1}^{N}n_i$.

Entonces primero escribí el Logaritmo de Verosimilitud Máximo Restringido:

$$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).$$

Luego tengo que diferenciar el logaritmo de la verosimilitud, $l_R$, con respecto a $\sigma^2$ y igualarlo a cero, es decir,

$-\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.$

Pero básicamente no puedo diferenciar,

$\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)$ y

$\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}$.

Aquí hay una respuesta para diferenciar esa derivada pero no puedo combinarla para obtener la estimación REML $\hat\sigma^2_R$ en la ecuación $(1)$.

¿Cómo puedo obtener la estimación REML $\hat\sigma^2_R$ en la ecuación $(1)$?

Answer 1

NB. Simplifico la notación un poco y no uso formato en negrita.

---

Las siguientes reglas para diferenciales de matrices son útiles:

\begin{align} d\log \vert A\vert &= \mathrm{tr}(A^{-1}dA) \\ dA^{-1} & = -A^{-1}(dA)A^{-1}. \end{align}

Una buena fuente para tales reglas y cómo derivarlas es Magnus Neudecker. Ese texto también explica cómo ir entre diferenciales y derivadas. También puedes consultar este documento de Minka. The Matrix Cookbook es otra referencia estándar para tales reglas. Además de derivadas, esta respuesta también hace un uso repetido de la invariancia cíclica y la linealidad del operador de traza.

Basándonos en estas reglas, obtenemos las siguientes expresiones para los diferenciales en cuestión: $$ d\Sigma_i^{-1} = -\Sigma_i^{-1}d\Sigma_i \Sigma_i^{-1} = -\Sigma^{-2}d\sigma^2, $$ $$ d\log \vert \Sigma_i \vert = \mathrm{tr}\Sigma_i^{-1}d\Sigma_i = \mathrm{tr}\Sigma_i^{-1}Id\sigma^2, $$

y

\begin{align} d\log\vert X_i'\Sigma_i^{-1}X_i\vert &= \mathrm{tr}[(X_i'\Sigma_iX)^{-1}X_i'd\Sigma_i^{-1}X_i] \\ &= -\mathrm{tr}[(X_i'\Sigma_iX)^{-1}X_i'(\Sigma_i^{-1}d\Sigma_i \Sigma_i^{-1})X_i]\\ &= -\mathrm{tr}[\Sigma_i^{-1}X_i(X_i'\Sigma_iX)^{-1}X_i'\Sigma_i^{-1}Id\sigma^2] \end{align}

Así, el gradiente del log-likelihood restringido hasta términos escalares y aditivos que no involucran $\sigma^2$ es:

\begin{align} \sum_i \frac{\partial}{\partial \sigma^2}\left[-\log \vert \Sigma_i \vert - \log \vert X_i'\Sigma_i^{-1}X_i \vert - (y_i - X_i\beta)'\Sigma_i^{-1}(y_i - X_i\beta) \right] = \\ \sum_i \mathrm{tr}\left[-\Sigma_i^{-1} + [\Sigma_i^{-1}X_i(X_i'\Sigma_i^{-1}X)^{-1}X_i'\Sigma_i^{-1}] + (y_i - X_i \beta)'\Sigma_i^{-2}(y_i - X_i\beta)\right]. \end{align}

Dado que las $\Sigma_i$ son (asumo) definitivas positivas, y en particular de rango completo, poner el gradiente en cero es equivalente a: $$ \sum_i \mathrm{tr}\left[I - [X_i(X_i'\Sigma_i^{-1}X)^{-1}X_i'\Sigma_i^{-1}] - \Sigma_i^{-1}(y_i - X_i\beta)(y_i - X_i \beta)'\right] = 0. $$

Ahora, nota que $Q_i := I - [X_i(X_i'\Sigma_i^{-1}X)^{-1}X_i'\Sigma_i^{-1}]$ es una matriz de proyección en el complemento ortogonal del espacio columna de $X_i$. Por lo tanto, su traza es igual a la dimensión del espacio en el que se proyecta: $n_i - p$. Factorizando $\sigma^2$, se obtiene la estimación REML de $\sigma^2$ como la solución a $$\sigma^2 = \sum_i (y_i - X\beta)'[I + Z_iGZ_i' / \sigma^2]^{-1}(y_i - X_i\beta) / (N_0 - p),$$

suponiendo que el problema de optimización es convexo.

Esto es diferente de tu expresión. Es posible que haya cometido un error en algún lugar y, en ese caso, espero que tú u otro lector lo pueda identificar. Después de leer un poco más, tengo la impresión de que los modelos mixtos comúnmente están parametrizados en términos de $\tilde{G} = G/\sigma^2$. Parametrizar la condición de primer orden derivada aquí en términos de $\tilde{G}$ haría que coincidiera con tu condición de primer orden.

Answer 2

Este es el código del Modelo del Constructor de Anuncios para el modelo. Esta versión estimará los efectos fijos.

DATA_SECTION
  init_number N   // número de individuos
  init_int p
  init_int q
  init_vector n(1,N) // número de puntos temporales
  init_3darray X(1,N,1,n,1,p)
  init_3darray Z(1,N,1,n,1,q)
  init_matrix obsy(1,N,1,n)
PARAMETER_SECTION
  init_number log_sigma
 !! int q1=q*(q-1)/2;
  init_vector chcoff(1,q1)
  init_vector logs(1,q)
  init_vector beta(1,p)
  random_effects_matrix eps(1,N,1,q)
  objective_function_value f
PROCEDURE_SECTION
  for (int i=1;i<=N;i++)
  {
    f+=f1(i,eps(i),beta,log_sigma,logs,chcoff);
  }
SEPARABLE_FUNCTION   dvariable f1(int i,const dvar_vector& epsi,const dvar_vector beta,const prevariable & log_sigma,const dvar_vector& logs,const dvar_vector& chcoff)
  dvariable sigma=exp(log_sigma);
  dvar_vector s=exp(logs);
  dvariable var=square(sigma);
  dvar_matrix ch(1,q,1,q);
  ch.initialize();
  // ch es la descomposición choleski de G
  ch(1,1)=s(1);
  int offset=0;
  for (int j=2;j<=q;j++)
  {
    ch(j)(1,j-1)=chcoff(1+offset,j-1+offset).shift(1);
    ch(j,j)=1.0;
    ch(j)/=norm(ch(j));
    offset+=j-1;
    ch(j)*=s(j);
  }
  dvar_vector predy=X(i)*beta+Z(i)*(ch*epsi);
  f+= 0.5*n(i)*log(var)+0.5*norm2(obsy(i)-predy)/var;
  f+= 0.5*norm2(epsi);

La versión REML se calcula integrando los efectos fijos. Esto implica cambiar una línea en el código anterior.

 init_vector beta(1,p)

se cambia a

 random_effects_vector beta(1,p)

Ahora supongamos que quiero modificar este código para usar la distribución t de Student. La línea

f+= 0.5*n(i)*log(var)+0.5*norm2(obsy(i)-predy)/var;

se cambia a

f-=ld_student(obsy(i),predy,var,dof);

y las líneas

init_number log_dof_coff

dvariable dof=1.0+exp(log_dof_coff);

se agregan al código en los lugares apropiados para incluir el parámetro de grados de libertad.
Ten en cuenta que mientras que hay una solución de forma cerrada para integrar sobre los parámetros para la distribución normal, esto no es cierto para la distribución t de Student. La integración se realiza a través de la aproximación de Laplace.