Máxima verosimilitud restringida con un rango de columna inferior al completo de $X$

Question

Esta pregunta trata de la estimación por máxima verosimilitud restringida (REML) en una versión particular del modelo lineal, a saber:

$$Y = X(\alpha)\beta + \epsilon, \\ \epsilon\sim N_n(0, \Sigma(\alpha)),$$

donde $X(\alpha)$ es un ( $n \times p$ ) parametrizada por $\alpha \in \mathbb R^k$ , tal y como está $\Sigma(\alpha)$ . $\beta$ es un vector desconocido de parámetros molestos; lo que interesa es estimar $\alpha$ y tenemos $k\leq p\ll n$ . Estimar el modelo por máxima verosimilitud no es un problema, pero quiero utilizar REML. Es bien conocido, véase por ejemplo LaMotte que la probabilidad $A'Y$ , donde $A$ es cualquier matriz semiortogonal tal que $A'X=0$ se puede escribir

$$L_{\text{REML}}(\alpha\mid Y) \propto\vert X'X\vert^{1/2} \vert \Sigma\vert^{-1/2}\vert X'\Sigma^{-1}X\vert^{-1/2}\exp\left\{-\frac{1}{2} r'\Sigma^{-1}r \right\}, \\ r = (I - X(X'\Sigma^{-1}X)^+X'\Sigma^{-1})Y,$$

cuando $X$ es el rango de la columna completa .

Mi problema es que para algunos perfectamente razonable, y científicamente interesante, $\alpha$ la matriz $X(\alpha)$ no es de rango de columna completo. Todas las derivaciones que he visto de la probabilidad restringida anterior hacen uso de igualdades determinantes que no son aplicables cuando $\vert X'X\vert=0$ es decir, asumen el rango de columna completo de $X$ . Esto significa que la probabilidad restringida anterior sólo es correcta para mi ajuste en partes del espacio de parámetros, y por lo tanto no es lo que quiero optimizar.

Pregunta: ¿Existen probabilidades restringidas más generales, derivadas, en la literatura estadística o en otros lugares, sin la suposición de que $X$ ¿se trata de un rango de columna completo? Si es así, ¿qué aspecto tienen?

Algunas observaciones:

• Derivar la parte exponencial no es un problema para cualquier $X(\alpha)$ y se puede escribir en términos de la inversa de Moore-Penrose como se ha indicado anteriormente
• Las columnas de $A$ son una base ortonormal (cualquiera) para $C(X)^\bot$
• Para los conocidos $A$ la probabilidad de $A'Y$ puede escribirse fácilmente para cada $\alpha$ pero, por supuesto, el número de vectores base, es decir, de columnas, en $A$ depende del rango de la columna de $X$

Si alguien interesado en esta cuestión cree que la parametrización exacta de $X,\Sigma$ ayudaría, hágamelo saber y los anotaré. En este momento, sin embargo, estoy principalmente interesado en un REML para un general $X$ de las dimensiones correctas.

---

A continuación se ofrece una descripción más detallada del modelo. Dejemos que $y_t = \mu + Ay_{t - 1} + v_t, t = 1, \dots, T$ ser un $r$ -Autoregresión vectorial de primer orden [VAR(1)] donde $v_t \overset{iid}{\sim}N(0, \Omega)$ . Supongamos que el proceso se inicia en algún valor fijo $y_0$ en el momento $t = 0$ .

Definir $Y = [y_1', \dots, y_T']'$ . El modelo puede escribirse en la forma de modelo lineal $Y = X\beta + \varepsilon$ utilizando las siguientes definiciones y notación:

\begin{align} X &= [1_T \otimes I_r, C^{-1}B] \\ \beta &= [\mu', y_0' - \mu']' \\ \mathrm{var}(\varepsilon)^{-1} &= C'(I_T \otimes \Omega^{-1})C \\ C &= $$\begin{bmatrix} I_r & 0 & 0 & \cdots \\ -A & I_r & 0 & \cdots \\ 0 & -A & I_r & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{bmatrix}$$ \\ B &= e_{1, T} \N a veces A, \fin

donde $1_T$ denota una $T-$ vector dimensional de unos y $e_{1,T}$ el primer vector base estándar de $\mathbb R^T$ .

Denote $\alpha = \mathrm{vec}(A)$ . Tenga en cuenta que si $A$ no es de rango completo entonces $X(\alpha)$ no es un rango de columna completo. Esto incluye, por ejemplo, los casos en los que uno de los componentes de $y_t$ no depende del pasado.

La idea de estimar VARs usando REML es bien conocida, por ejemplo, en la literatura de regresiones predictivas (ver por ejemplo Phillips y Chen y sus referencias).

Quizá convenga aclarar que la matriz $X$ no es una matriz de diseño en el sentido habitual, simplemente se sale del modelo y a menos que haya a priori conocimiento sobre $A$ no hay, hasta donde yo sé, ninguna forma de reparametrizarlo para que sea de rango completo.

---

He publicado una pregunta en math.stackexchange que está relacionada con ésta en el sentido de que una respuesta a la pregunta matemática puede ayudar a derivar una probabilidad que responda a esta pregunta.

Answer 1

Derivar la parte exponencial no es un problema para cualquier X()X() y puede escribirse en términos de la inversa de Moore-Penrose como se ha indicado anteriormente

Tengo dudas de que esta observación sea correcta. La inversa generalizada en realidad pone una restricción lineal adicional en sus estimadores [Rao&Mitra], por lo tanto deberíamos considerar la probabilidad conjunta como un todo en lugar de adivinar "la inversa de Moore-Penrose funcionará para la parte exponencial". Esto parece formalmente correcto, pero probablemente no entienda el modelo mixto correctamente.

$\blacksquare$ (1)¿Cómo pensar correctamente los modelos de efectos mixtos?

Hay que pensar en el modelo de efectos mixtos de una manera diferente antes de tratar de enchufar el g-inverso (O el inverso de Moore-Penrose, que es un tipo especial de g-inverso reflexivo [Rao&Mitra]) mecánicamente en la fórmula dada por RMLE (Estimador de Máxima Verosimilitud Restringido, lo mismo de abajo).

$$\boldsymbol{X}=\left(\begin{array}{cc} fixed\quad effect\\ & random\quad effect \end{array}\right)$$

Una forma común de pensar en el efecto mixto es que la parte del efecto aleatorio en la matriz de diseño es introducida por el error de medición, que lleva otro nombre de "predictor estocástico" si nos importa más la predicción que la estimación. Esta es también una motivación histórica del estudio de la matriz estocástica en el ámbito de la estadística.

Mi problema es que para algunos perfectamente razonable, y científicamente interesante, la matriz X()X() no es de rango de columna completo.

Teniendo en cuenta esta forma de pensar la probabilidad, la probabilidad de que $X(\alpha)$ no es de rango completo es cero. Esto se debe a que la función determinante es continua en las entradas de la matriz y la distribución normal es una distribución continua que asigna probabilidad cero a un solo punto. La probabilidad de rango defectuoso $X(\alpha)$ es positivo si se parametriza de forma patológica como $\left(\begin{array}{ccc} \alpha & \alpha\\ \alpha & \alpha\\ & & random\quad effect \end{array}\right)$ .

Así que la solución a tu pregunta también es bastante sencilla, simplemente perturbas tu matriz de diseño $X_\epsilon(\alpha)=X(\alpha)+\epsilon\left(\begin{array}{cc} I & 0\\ 0 & 0 \end{array}\right)$ (perturbar sólo la parte del efecto fijo), y utilizar la matriz perturbada (que es de rango completo) para llevar a cabo todas las derivaciones. A menos que su modelo tenga jerarquías complicadas o $X$ sí es casi singular, no veo que haya un problema serio cuando se toma $\epsilon\rightarrow 0$ en el resultado final ya que la función determinante es continua y podemos tomar el límite dentro de la función determinante. $lim_{\epsilon\rightarrow 0}|X_\epsilon|=|lim_{\epsilon\rightarrow 0}X_\epsilon|$ . Y en forma de perturbación la inversa de $X_\epsilon$ se puede obtener mediante el teorema de Sherman-Morrision-Woodbury. Y el determinante de la matriz $I+X$ se da en un libro estándar de álgebra lineal como [Horn&Johnson]. Por supuesto, podemos escribir el determinante en términos de cada entrada de la matriz, pero siempre se prefiere la perturbación [Horn&Johnson].

$\blacksquare$ (2)¿Cómo debemos tratar los parámetros molestos en un modelo?

Como ve, para tratar la parte del efecto aleatorio en el modelo, debemos considerarla como una especie de "parámetro molesto". El problema es: ¿es la RMLE la forma más adecuada de eliminar un parámetro molesto? Incluso en los modelos GLM y de efectos mixtos, el RMLE está lejos de ser la única opción. [Basu] señaló que hay muchas otras formas de eliminar parámetros en el marco de la estimación. Hoy en día, la gente tiende a elegir entre el RMLE y la modelización bayesiana porque corresponden a dos soluciones populares basadas en el ordenador: EM y MCMC respectivamente.

En mi opinión, es definitivamente más adecuado introducir una prioridad en la situación de rango defectuoso en la parte de efectos fijos. O puede reparametrizar su modelo para convertirlo en uno de rango completo.

Además, en caso de que su efecto fijo no sea de rango completo, podría preocuparse por una estructura de covarianza mal especificada, ya que los grados de libertad de los efectos fijos deberían ir a la parte del error. Para ver este punto más claramente, puede considerar el MLE (también LSE) para el GLS (General least squre) $\hat{\beta}=(X\Sigma^{-1} X')^{-1}\Sigma^{-1}y$ donde $\Sigma$ es la estructura de covarianza del término de error, para el caso en que $X(\alpha)$ no es de rango completo.

$\blacksquare$ (3)Otros comentarios

El problema no es cómo modificar el RMLE para que funcione en el caso de que la parte de efectos fijos de la matriz no sea de rango completo; el problema es que en ese caso su modelo en sí mismo puede ser problemático si el caso de rango no completo tiene una probabilidad positiva.

Un caso relevante que he encontrado es que en el caso espacial la gente puede querer reducir el rango de la parte de efectos fijos debido a consideraciones computacionales[Wikle].

No he visto ningún caso "científicamente interesante" en tal situación, ¿puede señalar alguna literatura en la que el caso sin rango completo sea de mayor preocupación? Me gustaría saberlo y debatirlo más a fondo, gracias.

$\blacksquare$ Referencia

[Rao&Mitra]Rao, Calyampudi Radhakrishna, y Sujit Kumar Mitra. Generalized inverse of matrices and its applications. Vol. 7. Nueva York: Wiley, 1971.

[Basu]Basu, Debabrata. "Sobre la eliminación de parámetros molestos". Journal of the American Statistical Association 72.358 (1977): 355-366.

[Horn&Johnson]Horn, Roger A., y Charles R. Johnson. Matrix analysis. Cambridge university press, 2012.

[Wikle]Wikle, Christopher K. "Low-rank representations for spatial processes". Handbook of Spatial Statistics (2010): 107-118.