Eingeschränkte maximale Wahrscheinlichkeit mit weniger als dem vollen Spaltenrang von

14

Diese Frage befasst sich mit der eingeschränkten Maximalwahrscheinlichkeitsschätzung (REML) in einer bestimmten Version des linearen Modells, nämlich:

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

Wobei eine durch ; parametrisierte ( ) Matrix ist , wie auch . ist ein unbekannter Vektor von Störparametern; das Interesse liegt in der Schätzung von , und wir haben . Das Modell nach der maximalen Wahrscheinlichkeit abzuschätzen ist kein Problem, aber ich möchte REML verwenden. Es ist bekannt, siehe z. B. LaMotte , dass die Wahrscheinlichkeit , wobei eine beliebige Matrix ist, so dass geschrieben werden kannX(α)n×pαRkΣ(α)βαkpnAYAAX=0

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

wenn volle Spaltenrang ist .X

Mein Problem ist, dass für einige absolut vernünftige und wissenschaftlich interessante die Matrix nicht den vollen Spaltenrang hat. Alle Ableitungen, die ich von der oben genannten eingeschränkten Wahrscheinlichkeit gesehen habe, verwenden Determinantengleichungen, die nicht anwendbar sind, wenn , dh sie nehmen den vollen Spaltenrang von . Dies bedeutet, dass die obige eingeschränkte Wahrscheinlichkeit nur für meine Einstellung auf Teile des Parameterraums richtig ist und daher nicht das ist, was ich optimieren möchte.X ( α ) | X ' X | = 0 XαX(α)|XX|=0X

Frage: Gibt es allgemeinere eingeschränkte Wahrscheinlichkeiten, die in der statistischen Literatur oder anderswo abgeleitet wurden, ohne die Annahme, dass der volle Spaltenrang ist? Wenn ja, wie sehen sie aus?X

Einige Beobachtungen:

  • Das Ableiten des exponentiellen Teils ist für jedes kein Problem, und es kann in Form der Moore-Penrose-Inverse wie oben beschrieben geschrieben werdenX(α)
  • Die Spalten von sind eine (beliebige) orthonormale Basis fürC ( X ) AC(X)
  • Für bekanntes kann die Wahrscheinlichkeit für leicht für jedes Alpha notiert werden , aber natürlich hängt die Anzahl der Basisvektoren, dh der Spalten, in A vom Spaltenrang von X abA ' Y α A XAAYαAX

Wenn jemand, der sich für diese Frage interessiert, der Meinung ist, dass die genaue Parametrisierung von hilfreich ist, lass es mich wissen und ich werde sie aufschreiben. An dieser Stelle bin ich jedoch hauptsächlich an einem REML für ein allgemeines mit den richtigen Abmessungen interessiert .X,Σ X


Eine detailliertere Beschreibung des Modells folgt hier. Sei einryt=μ+Ayt1+vt,t=1,,Tr dimensionale Vektorautoregression erster Ordnung [VAR (1)], wobei . Angenommen, der Prozess wird in einem festen Wert zum Zeitpunkt gestartet .y 0 t = 0vtiidN(0,Ω)y0t=0

Definiere . Das Modell kann in der linearen Modellform Verwendung der folgenden Definitionen und Notation geschrieben werden: Y = X β + εY=[y1,,yT]Y=Xβ+ε

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

wobei bezeichnet1TT dimensionaler Vektor von Einsen und der erste Standard Basisvektor von .e1,TRT

Bezeichne . Beachten Sie, dass, wennα=vec(A)A kein voller Spaltenrang ist kein voller Rang ist. Dies schließt zum Beispiel Fälle ein, in denen eine der Komponenten von nicht von der Vergangenheit abhängt.X(α)yt

Die Idee der Schätzung von VARs unter Verwendung von REML ist beispielsweise in der Literatur zu prädiktiven Regressionen gut bekannt (siehe z. B. Phillips und Chen und die darin enthaltenen Referenzen).

Es kann sich lohnen, klarzustellen, dass die Matrix keine Entwurfsmatrix im üblichen Sinne ist, sondern nur aus dem Modell herausfällt und es keine gibtX a priori - Wissen über ist, soweit ich das beurteilen kann, keine Möglichkeit zu parametrieren es ist voller Rang.A


Ich habe eine Frage zu math.stackexchange gestellt , die sich auf diese in dem Sinne bezieht, dass eine Antwort auf die mathematische Frage hilfreich sein kann, um eine Wahrscheinlichkeit abzuleiten, die diese Frage beantworten würde.

ekvall
quelle
1
Eine Möglichkeit, die Frage zu beantworten, besteht darin, zu fragen, was in linearen gemischten Modellen geschieht, wenn die Modellmatrix nicht den vollen Spaltenrang hat.
Greenparker
Danke für das Kopfgeld @Greenparker. Und ja, wenn eine eingeschränkte Wahrscheinlichkeit für ein lineares gemischtes Modell mit weniger als der vollen Spaltenrang-Entwurfsmatrix für feste Effekte aufgeschrieben werden könnte, würde dies helfen.
Freitag,

Antworten:

2

Das Ableiten des exponentiellen Teils ist für jedes X (α) X (α) kein Problem, und es kann in Bezug auf die Moore-Penrose-Inverse wie oben geschrieben werden

Ich habe Zweifel, dass diese Beobachtung richtig ist. Die verallgemeinerte Inverse schränkt Ihre Schätzer zusätzlich linear ein [Rao & Mitra], daher sollten wir die gemeinsame Wahrscheinlichkeit als Ganzes betrachten, anstatt zu raten, dass die Moore-Penrose-Inverse für den exponentiellen Teil funktioniert. Dies scheint formal korrekt zu sein, aber Sie verstehen das gemischte Modell wahrscheinlich nicht richtig.

(1) Wie man richtig denkt?

Sie müssen das Mixed-Effekt-Modell auf eine andere Weise überlegen, bevor Sie versuchen, die g-Inverse (ODER Moore-Penrose-Inverse, eine spezielle Art der reflexiven g-Inverse [Rao & Mitra]) mechanisch in die von RMLE (Restricted) gegebene Formel einzufügen Maximum Likelihood Estimator (siehe unten).

X=(fixedeffectrandomeffect)

Eine übliche Denkweise für einen gemischten Effekt ist, dass der zufällige Effektteil in der Entwurfsmatrix durch einen Messfehler eingefügt wird, der einen anderen Namen "stochastischer Prädiktor" trägt, wenn wir uns mehr für die Vorhersage als für die Schätzung interessieren. Dies ist auch eine historische Motivation für das Studium der stochastischen Matrix in der Statistik.

Mein Problem ist, dass für einige absolut vernünftige und wissenschaftlich interessante αα die Matrix X (α) X (α) nicht den vollen Spaltenrang hat.

Unter Berücksichtigung der Wahrscheinlichkeit ist die Wahrscheinlichkeit, dass nicht den vollen Rang hat, Null. Dies liegt daran, dass die Determinantenfunktion in Matrixeinträgen stetig ist und die Normalverteilung eine stetige Verteilung ist, die einem einzelnen Punkt eine Wahrscheinlichkeit von Null zuweist. Die Wahrscheinlichkeit eines fehlerhaften Rangs ist positiv, wenn Sie ihn auf pathologische Weise wie folgt parametrisieren: .X(α)X(α)(ααααrandomeffect)

Die Lösung Ihrer Frage ist also auch ziemlich einfach: Sie stören einfach Ihre Entwurfsmatrix (stören Sie nur den Part mit festem Effekt) und verwenden Sie die gestörte Matrix (die den vollen Rang hat), um alle Ableitungen durchzuführen. Wenn Ihr Modell keine komplizierten Hierarchien hat oder selbst nahezu singulär ist, sehe ich kein ernstes Problem, wenn Sie im Endergebnis nehmen, da die Determinantenfunktion stetig ist und wir die Grenze innerhalb der Determinantenfunktion nehmen können. . Und in Störungsform die Umkehrung vonXϵ(α)=X(α)+ϵ(I000)Xϵ0limϵ0|Xϵ|=|limϵ0Xϵ|Xϵkann durch Sherman-Morrision-Woodbury Theorem erhalten werden. Die Determinante von Matrix ist in einem Standardbuch zur linearen Algebra wie [Horn & Johnson] angegeben. Natürlich können wir die Determinante in Bezug auf jeden Eintrag in der Matrix schreiben, aber eine Störung wird immer bevorzugt [Horn & Johnson].I+X

(2) Wie sollen wir mit Störparametern in einem Modell umgehen?

Wie Sie sehen, sollten wir den Zufallseffektteil im Modell als eine Art "Störparameter" betrachten. Das Problem ist: Ist RMLE der geeignetste Weg, um einen Störparameter zu eliminieren? Auch bei GLM- und Mixed-Effect-Modellen ist RMLE bei weitem nicht die einzige Wahl. [Basu] wies darauf hin, dass viele andere Möglichkeiten zur Beseitigung von Parametern bei der Festlegung der Schätzung. Heutzutage tendieren die Menschen dazu, zwischen RMLE- und Bayes-Modellierung zu wählen, da sie zwei gängigen computergestützten Lösungen entsprechen: EM und MCMC.

Meiner Meinung nach ist es definitiv besser, einen Prior in die Situation des fehlerhaften Ranges in den Fixeffektteil einzuführen. Oder Sie können Ihr Modell neu parametrisieren, um es zu einem vollständigen Modell zu machen.

Falls Ihr fester Effekt nicht den vollen Rang hat, können Sie sich über der falsch spezifizierten Kovarianzstruktur Sorgen machen, da die Freiheitsgrade in festen Effekten in den Fehlerteil eingehen sollten. Um diesen Punkt klarer zu sehen, sollten Sie das MLE (auch LSE) für das GLS (General least squre) berücksichtigen. wobei die Kovarianzstruktur des Fehlerausdrucks ist, für den Fall, dass nicht den vollen Rang hat.β^=(XΣ1X)1Σ1yΣX(α)

(3) Weitere Kommentare

Das Problem besteht nicht darin, wie Sie das RMLE ändern, damit es funktioniert, wenn ein Teil der Matrix mit festen Effekten nicht den vollen Rang hat. Das Problem ist, dass in diesem Fall Ihr Modell selbst problematisch sein kann, wenn ein nicht vollständiger Fall eine positive Wahrscheinlichkeit hat.

Ein relevanter Fall, auf den ich gestoßen bin, ist, dass die Leute im räumlichen Fall den Rang eines Teils mit festem Effekt aufgrund von rechnerischen Überlegungen reduzieren möchten [Wikle].

Ich habe in einer solchen Situation keinen "wissenschaftlich interessanten" Fall gesehen. Können Sie auf Literatur verweisen, in der der nicht vollständige Fall von größter Bedeutung ist? Ich würde gerne weiter wissen und diskutieren, danke.

Referenz

[Rao & Mitra] Rao, Calyampudi Radhakrishna und Sujit Kumar Mitra. Verallgemeinerte Inverse von Matrizen und ihren Anwendungen. Vol. 7. New York: Wiley, 1971.

[Basu] Basu, Debabrata. "Über die Beseitigung von Störparametern." Journal of the American Statistical Association 72.358 (1977): 355 & ndash; 366.

[Horn & Johnson] Horn, Roger A. und Charles R. Johnson. Matrixanalyse. Cambridge University Press, 2012.

[Wikle] Wikle, Christopher K. "Geringwertige Darstellungen für räumliche Prozesse." Handbook of Spatial Statistics (2010): 107-118.

Henry.L
quelle
Vielen Dank für Ihr Interesse und Ihre durchdachte Antwort. + 1 für Ihre Mühe. Ich werde es ausführlicher lesen und mit einigen Erläuterungen zurückkommen. Ich denke, als Erstes muss ich klarstellen, dass dieses Modell keine zufälligen Effekte enthält und die Matrix keine Entwurfsmatrix ist, außer vielleicht beim Namen, weil es kein besseres Wort gibt. es ist eine stark nichtlineare Funktion (deterministisch) des Parameters α, die aus der Vektorisierung der Koeffizientenmatrix in einem vektorautoregressiven Prozess besteht, daher ist das Konzept der Wahrscheinlichkeit eines niedrigen Ranges nicht aussagekräftig. Xα
ekvall
@ Student001 Ja, zögern Sie nicht, dies zu klären, da ich es auch eher als GLM empfinde als als gemischtes Modell. Ich werde versuchen, noch einmal zu antworten, wenn ich kann :)
Henry.L
@ Student001 Wenn du kannst, schreibe das ganze Modell und ich würde gerne einen solchen Fall studieren, möglicherweise AR (1) in räumlicher Umgebung, denke ich.
Henry.L
"Angesichts dieser Betrachtungsweise der Wahrscheinlichkeit ist die Wahrscheinlichkeit, dass nicht den vollen Rang hat, Null." Richtige Antwort, falsches Problem. Die Wahrscheinlichkeit, dass es in endlicher Genauigkeit numerisch nicht den vollen Rang hat, ist ungleich Null. X(α)
Mark L. Stone
@ MarkL.Stone Als Lösung für das sorgfältige Lesen von Zeilen habe ich bereits eine Störung bereitgestellt. Dies ist eine Standardlösung für die numerische Singularität. Und das OP sagte, dass er die Beschreibung aktualisieren wird, also denke ich, dass wir eine Einigung über das richtig formulierte Problem erzielen werden.
Henry.L