Probabilité maximale restreinte avec un rang de colonne inférieur à la totalité de

14

Cette question traite de l'estimation du maximum de vraisemblance restreint (REML) dans une version particulière du modèle linéaire, à savoir:

Y=X(α)β+ϵ,ϵNn(0,Σ(α)),

X(α) est une matrice ( n×p ) paramétrée par αRk , tout comme Σ(α) . β est un vecteur inconnu de paramètres de nuisance; l'intérêt est d'estimer α , et on a kpn . L'estimation du modèle par maximum de vraisemblance n'est pas un problème, mais je veux utiliser REML. Il est bien connu, voir par exemple LaMotte , que la vraisemblance AY , où A est une matrice semi-orthogonale telle queAX=0 peut s'écrire

LREML(αY)|XX|1/2|Σ|1/2|XΣ1X|1/2exp{12rΣ1r},r=(IX(XΣ1X)+XΣ1)Y,

lorsque est le rang de colonne completX .

Mon problème est que pour certains parfaitement raisonnables et scientifiquement intéressants, la matrice X ( α ) n'est pas de plein rang de colonne. Toutes les dérivations que j'ai vues de la probabilité restreinte ci-dessus utilisent des égalités déterminantes qui ne sont pas applicables lorsque | X X | = 0 , à savoir qu'ils assument rang de la colonne complète de X . Cela signifie que la probabilité restreinte ci-dessus n'est correcte que pour mon réglage sur des parties de l'espace des paramètres, et n'est donc pas ce que je souhaite optimiser.αX(α)|XX|=0X

Question: Y a-t-il des probabilités restreintes plus générales, dérivées, dans la littérature statistique ou ailleurs, sans l'hypothèse que soit le rang complet des colonnes? Si oui, à quoi ressemblent-ils?X

Quelques observations:

  • Dériver la partie exponentielle n'est pas un problème pour tout et elle peut être écrite en termes de l'inverse de Moore-Penrose comme ci-dessusX(α)
  • Les colonnes de sont une base orthonormée (quelconque) pour C ( X ) AC(X)
  • Pour connu , la probabilité pour A Y peut facilement être écrite pour chaque α , mais bien sûr le nombre de vecteurs de base, c'est-à-dire de colonnes, dans A dépend du rang de colonne de XAAYαAX

Si quelqu'un intéressé par cette question croit que le paramétrage exact de aiderait, faites-le moi savoir et je les noterai. À ce stade, je suis surtout intéressé par un REML pour un X général des dimensions correctes.X,Σ X


Une description plus détaillée du modèle suit ici. Soit une autorégression vectorielle de premier ordre r- dimensionnelle [VAR (1)] où v t i i d N ( 0 , Ω ) . Supposons que le processus démarre à une valeur fixe y 0 au temps t = 0 .yt=μ+Ayt1+vt,t=1,,TrvtiidN(0,Ω)y0t=0

Définissez . Le modèle peut être écrit sous la forme du modèle linéaire Y = X β + ε en utilisant les définitions et la notation suivantes:Y=[y1,,yT]Y=Xβ+ε

X=[1TIr,C1B]β=[μ,y0μ]var(ε)1=C(ITΩ1)CC=[Ir00AIr00AIr]B=e1,TA,

désigne un T - vectoriel de dimension de uns et de e 1 , T le premier vecteur de base de type R T .1TTe1,TRT

Notons . Notez que si A n'est pas un rang complet, alors X ( α ) n'est pas un rang complet de colonne. Cela inclut, par exemple, les cas où l'une des composantes de y t ne dépend pas du passé.α=vec(A)AX(α)yt

L'idée d'estimer les VAR en utilisant REML est bien connue, par exemple, dans la littérature sur les régressions prédictives (voir par exemple Phillips et Chen et les références qui y sont contenues).

Il peut être utile de préciser que la matrice n'est pas une matrice de conception au sens habituel, elle tombe juste hors du modèle et à moins qu'il n'y ait une connaissance a priori de A, il n'y a, pour autant que je sache , aucun moyen de reparamétrer que ce soit complet.XA


J'ai posté une question sur math.stackexchange qui est liée à celle-ci dans le sens où une réponse à la question mathématique peut aider à dériver une probabilité qui répondrait à cette question.

ekvall
la source
1
Peut-être qu'une façon de répondre à la question est de demander, que se passe-t-il dans les modèles mixtes linéaires lorsque la matrice du modèle n'est pas de rang de colonne complet?
Greenparker
Merci pour la prime @Greenparker. Et, oui, si une probabilité restreinte pouvait être notée pour un modèle mixte linéaire, avec une matrice de conception d'effets fixes de rang inférieur à la pleine colonne, cela aiderait.
ekvall

Réponses:

2

Dériver la partie exponentielle n'est pas un problème pour tout X (α) X (α) et elle peut être écrite en termes de l'inverse de Moore-Penrose comme ci-dessus

Je doute que cette observation soit correcte. L'inverse généralisé impose en fait une restriction linéaire supplémentaire à vos estimateurs [Rao & Mitra], donc nous devrions considérer la vraisemblance conjointe dans son ensemble au lieu de deviner "l'inverse de Moore-Penrose fonctionnera pour la partie exponentielle". Cela semble formellement correct mais vous ne comprenez probablement pas correctement le modèle mixte.

(1) Comment penser correctement les modèles à effets mixtes?

Vous devez penser le modèle à effets mixtes d'une manière différente avant d'essayer de brancher mécaniquement l'inverse g (OU Moore-Penrose inverse, qui est un type spécial d'inverse g réflexif [Rao & Mitra]) dans la formule donnée par RMLE (Restreint Estimateur du maximum de vraisemblance, même ci-dessous.).

X=(fixedeffectrandomeffect)

Une façon courante de penser l'effet mixte est que la partie de l'effet aléatoire dans la matrice de conception est introduite par une erreur de mesure, qui porte un autre nom de "prédicteur stochastique" si nous nous soucions davantage de la prédiction plutôt que de l'estimation. Il s'agit également d'une motivation historique de l'étude de la matrice stochastique dans l'établissement des statistiques.

Mon problème est que pour certains parfaitement raisonnables et scientifiquement intéressants, αα la matrice X (α) X (α) n'est pas de rang de colonne complet.

Compte tenu de cette façon de penser la vraisemblance, la probabilité que ne soit pas de rang complet est nulle. En effet, la fonction déterminante est continue dans les entrées de la matrice et la distribution normale est une distribution continue qui attribue une probabilité nulle à un seul point. La probabilité de défaut de rang X ( α ) est positive si vous ne la paramétrez de manière pathologique comme ( α α α αX(α)X(α).(ααααrandomeffect)

Donc, la solution à votre question est également assez simple, vous perturbez simplement votre matrice de conception (perturbe la partie à effet fixe uniquement), et utilisez la matrice perturbée (qui est de rang complet) pour effectuer toutes les dérivations. À moins que votre modèle n'ait des hiérarchies compliquées ou que X lui-même soit presque singulier, je ne vois pas de problème sérieux lorsque vous prenez ϵ 0 dans le résultat final car la fonction déterminante est continue et nous pouvons prendre la limite à l'intérieur de la fonction déterminante. l iXϵ(α)=X(α)+ϵ(I000)Xϵ0. Et sous forme de perturbation, l'inverse de X ϵ peut être obtenu par le théorème de Sherman-Morrision-Woodbury. Et le déterminant de la matrice I + X est donné dans un livre d'algèbre linéaire standard comme [Horn & Johnson]. Bien sûr, nous pouvons écrire le déterminant en termes de chaque entrée de la matrice, mais la perturbation est toujours préférée [Horn & Johnson].limϵ0|Xϵ|=|limϵ0Xϵ|XϵI+X

(2) Comment traiter les paramètres de nuisance dans un modèle?

Comme vous le voyez, pour traiter la partie à effet aléatoire du modèle, nous devons la considérer comme une sorte de "paramètre de nuisance". Le problème est le suivant: RMLE est-il le moyen le plus approprié d'éliminer un paramètre de nuisance? Même dans les modèles GLM et à effets mixtes, RMLE est loin d'être le seul choix. [Basu] a souligné que de nombreuses autres façons d'éliminer les paramètres dans le cadre de l'estimation. Aujourd'hui, les gens ont tendance à choisir entre RMLE et la modélisation bayésienne car ils correspondent à deux solutions informatiques populaires: EM et MCMC respectivement.

À mon avis, il est certainement plus approprié d'introduire un a priori dans la situation de rang défectueux dans la partie à effet fixe. Ou vous pouvez re-paramétrer votre modèle afin d'en faire un rang complet.

β^=(XΣ1X)1Σ1yΣX(α)

Le problème n'est pas de savoir comment modifier le RMLE pour le faire fonctionner dans le cas où une partie à effet fixe de la matrice n'est pas de plein rang; le problème est que dans ce cas, votre modèle lui-même peut être problématique si le cas non complet a une probabilité positive.

Un cas pertinent que j'ai rencontré est que dans le cas spatial, les gens peuvent vouloir réduire le rang de la partie à effet fixe en raison de considérations de calcul [Wikle].

Je n'ai pas vu de cas "scientifiquement intéressant" dans une telle situation, pouvez-vous signaler certains documents où le cas non complet est un problème majeur? Je voudrais savoir et discuter davantage, merci.

[Rao et Mitra] Rao, Calyampudi Radhakrishna et Sujit Kumar Mitra. Inverse généralisé des matrices et de ses applications. Vol. 7. New York: Wiley, 1971.

[Basu] Basu, Debabrata. "Sur l'élimination des paramètres de nuisance." Journal de l'American Statistical Association 72.358 (1977): 355-366.

[Horn & Johnson] Horn, Roger A. et Charles R. Johnson. Analyse matricielle. Cambridge University Press, 2012.

[Wikle] Wikle, Christopher K. "Représentations de bas rang pour les processus spatiaux." Handbook of Spatial Statistics (2010): 107-118.

Henry.L
la source
Xα
@ Student001 Oui, n'hésitez pas à apporter des précisions, car je le ressens également plus comme un GLM plutôt que comme un modèle mixte. Je vais essayer de répondre à nouveau si je peux :)
Henry.L
@ Student001 Si vous le pouvez, écrivez tout le modèle et j'aimerais étudier un tel cas, peut-être AR (1) dans un cadre spatial, je suppose.
Henry.L
X(α)
@ MarkL.Stone J'ai déjà fourni la perturbation comme solution si vous lisez attentivement les lignes, ce qui est une solution standard à la singularité numérique. Et l'OP a dit qu'il mettra à jour la description, donc je suppose que nous parviendrons à un consensus sur le problème correctement formulé.
Henry.L