Warum muss man REML (anstelle von ML) verwenden, um unter verschachtelten Var-Covar-Modellen zu wählen?

Verschiedene Beschreibungen zur Modellauswahl für zufällige Effekte von linearen gemischten Modellen weisen an, REML zu verwenden. Ich kenne den Unterschied zwischen REML und ML auf einer bestimmten Ebene, aber ich verstehe nicht, warum REML verwendet werden sollte, weil ML voreingenommen ist. Ist es beispielsweise falsch, mit ML eine LRT für einen Varianzparameter eines Normalverteilungsmodells durchzuführen (siehe folgenden Code)? Ich verstehe nicht, warum es bei der Modellauswahl wichtiger ist, unvoreingenommen zu sein als ML. Ich denke, die ultimative Antwort muss lauten: "Weil die Modellauswahl mit REML besser funktioniert als mit ML", aber ich möchte ein bisschen mehr darüber wissen. Ich habe die Ableitungen von LRT und AIC nicht gelesen (ich bin nicht gut genug, um sie gründlich zu verstehen), aber wenn REML explizit in den Ableitungen verwendet wird, muss ich nur wissen, dass dies tatsächlich ausreicht (z. B.

n <- 100
a <- 10
b <- 1
alpha <- 5
beta <- 1
x <- runif(n,0,10)
y <- rnorm(n,a+b*x,alpha+beta*x)

loglik1 <- function(p,x,y){
   a <- p[1]
   b <- p[2]
   alpha <- p[3]
  -sum(dnorm(y,a+b*x,alpha,log=T))
}

loglik2 <- function(p,x,y){
   a <- p[1]
   b <- p[2]
   alpha <- p[3]
   beta <- p[4]
  -sum(dnorm(y,a+b*x,alpha+beta*x,log=T))
}

m1 <- optim(c(a,b,alpha),loglik1,x=x,y=y)$value
m2 <- optim(c(a,b,alpha,beta),loglik2,x=x,y=y)$value
D <- 2*(m1-m2)
1-pchisq(D,df=1) # p-value

mixed-model maximum-likelihood unbiased-estimator likelihood-ratio reml streiten
quelle

Über REML und AIC sollten Sie sich diese Frage ansehen .

Elvis,

Antworten:

Eine sehr kurze Antwort: Das REML ist ein ML, der auf REML basierende Test ist also trotzdem korrekt. Da die Schätzung der Varianzparameter mit REML besser ist, ist es natürlich, diese zu verwenden.

Warum ist REML eine ML? Betrachten wir zum Beispiel ein Modell mit , und der Vektor der feststehenden Effekte, ist der Vektor von Zufallseffekten und

Y = X β + Z u + e

$Y = X\beta + Zu + e \def\R{\mathbb{R}}$

X \in R^{n \times p}

$X\in\R^{n\times p}$

Z \in R^{n \times q}

$Z\in\R^{n\times q}$

β \in R^{p}

$\beta \in \R^p$

u \sim N (0, τ I_{q})

$u \sim \mathcal N(0, \tau I_q)$

e \sim N (0, σ^{2} I_{n})

$e \sim \mathcal N(0, \sigma^2 I_n)$ . Die eingeschränkte Wahrscheinlichkeit kann erhalten werden, indem

Kontraste berücksichtigt werden , um die festen Effekte "zu entfernen". Genauer gesagt sei

, so dass

und

(d. H. Die Spalten von

sind eine orthonormale Basis des Vektorraums orthognal zu Raum, der durch die Spalten von

) erzeugt wird; dann ist

n - p

$n-p$

C \in R^{(n - p) \times n}

$C \in \R^{(n-p)\times n}$

C X = 0

$CX = 0$

C C^{'} = I_{n - p}

$CC' = I_{n-p}$

C^{'}

$C'$

X

$X$

C Y = C Z u + ϵ

$CY = CZ u + \epsilon$ mit

und der Wahrscheinlichkeit für

bei

die eingeschränkte Wahrscheinlichkeit .

ϵ \sim N (0, σ^{2} I_{n - p})

$\epsilon \sim \mathcal N(0, \sigma^2 I_{n-p})$

τ, σ^{2}

$\tau, \sigma^2$

C Y

$CY$

Elvis
quelle

Schöne Antwort (+1), kann ich zu Recht sagen, dass die Matrix

für den Durchschnitt vom Modell abhängt? Sie können also nur REML-Schätzungen für dieselbe

Matrix vergleichen?

C

$C$

C

$C$

Ja,

hängt von

(ich bearbeite die Antwort in einer Minute, um dies zu verdeutlichen), sodass Ihre verschachtelten Modelle dieselben Variablen mit festen Effekten haben müssen.

C

$C$

X

$X$

Elvis,

REML ist nicht ein ML! Die ML ist für ein gegebenes Wahrscheinlichkeitsmodell eindeutig definiert, die REML ist jedoch von der Parametrisierung mit festen Effekten abhängig. Siehe zB diesen Kommentar von Doug Bates (sowie viele historische zu R-SIG-Mixed-Modellen).

Livius

@Livius Ich denke, meine Antwort gibt deutlich genug an, wie die eingeschränkte Wahrscheinlichkeit aufgebaut ist. Es ist eine Wahrscheinlichkeit, es ist nur nicht die Wahrscheinlichkeit, wenn das beobachtete

in dem in der ersten angezeigten Gleichung geschriebenen Modell gegeben ist, sondern wenn der projizierte Vektor

in dem in der zweiten angezeigten Gleichung geschriebenen Modell gegeben ist. Die REML ist die aus dieser Wahrscheinlichkeit erhaltene ML.

Y

$Y$

C Y

$CY$

Elvis

Ich denke, das ist irgendwie der Grund für die Proteste von DBates zu diesem Thema: Es ist ein anderes Modell, und es ist ein Modell, für das Vergleiche schwierig sind, weil das Modell und die Parametrisierung miteinander verflochten sind. Sie berechnen also nicht die ML für Ihr Originalmodell, sondern die ML für ein anderes Modell, das sich aus einer bestimmten Parametrisierung Ihres Originalmodells ergibt. Daher sind REML-angepasste Modelle mit verschachtelten Strukturen mit festen Effekten keine verschachtelten Modelle mehr (wie oben erwähnt). ML-angepasste Modelle sind jedoch weiterhin verschachtelt, da Sie die Wahrscheinlichkeit für das angegebene Modell maximieren.

Livius

Likelihood-Ratio-Tests sind statistische Hypothesentests, die auf einem Verhältnis von zwei Wahrscheinlichkeiten basieren. Ihre Eigenschaften sind mit der Maximum Likelihood Estimation (MLE) verknüpft. (siehe zB Maximum Likelihood Estimation (MLE) in Laienbegriffen ).

In Ihrem Fall (siehe Frage) möchten Sie zwischen zwei verschachtelten var-covar-Modellen wählen. Nehmen wir an, Sie möchten zwischen einem Modell mit dem var-covar-Wert und einem Modell mit dem var-covar-Wert wählen wobei das zweite (einfaches Modell) ein Sonderfall des ersten (allgemeines Modell) ist. $\Sigma_g$ $\Sigma_s$

Der Test basiert auf der Likelihood - Verhältnis . Wo und $LR=-2 (log(\mathcal{L}_s(\hat{\Sigma}_s)) - log(\mathcal{L}_g(\hat{\Sigma}_g) )$ $\hat{\Sigma}_s$ $\hat{\Sigma}_g$ sind die Maximum-Likelihood-Schätzer.

Die Statistik ist asymptotisch (!) $LR$ $\chi^2$ .

Maximum-Likelihood-Schätzer sind bekanntermaßen konsistent, in vielen Fällen jedoch voreingenommen. Dies ist der Fall für den MLE Schätzer für die und , kann es zeigen, dass sie vorgespannt sind. Dies liegt daran, dass sie mit einem Mittelwert berechnet werden, der aus den Daten abgeleitet wurde, sodass die Streuung um diesen 'geschätzten Durchschnitt' kleiner ist als die Streuung um den wahren Mittelwert (siehe z. B. intuitive Erklärung zur Division durch bei der Berechnung der Standardabweichung) ? ) $\hat{\Sigma}_s$ $\hat{\Sigma}_g$ $n-1$

$LR$ $\chi^2$ $\hat{\Sigma}_s$ $\hat{\Sigma}_g$

$\hat{\Sigma}_s$ $\hat{\Sigma}_g$ $LR$ $\chi^2$ $\Sigma_s$ $\Sigma_g$ $LR$ $\chi^2$

Note that REML should only be used to choose among nested var-covar structures of models with the same mean, for models with different means, the REML is not appropriate, for models with different means one should use ML.

quelle

The statement "The statistic LR is , asymptotically (!) χ2" is not true in this case. This is because if

Σ_{s}

$\Sigma_s$ is nested in

Σ_{g}

$\Sigma_g$ , then

Σ_{s}

$\Sigma_s$ is on the boundary of

Σ_{g}

$\Sigma_g$ . In this case, the

χ^{2}

$\chi^2$ distribution does not hold. For example, see here

Cliff AB

@Cliff AB, this is what is explained below that statement and it is the reason you have to use REML.

-4

I have an answer that has more to do with common sense than with Statistics. If you take a look at PROC MIXED in SAS, the estimation can be performed with six methods:

http://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_mixed_sect008.htm

but REML is the default. Why? Apparently, the practical experience showed it has the best performance (e.g., the smallest chance of convergence problems). Therefore, if your goal is achievable with REML, then it makes sense to use REML as opposed to the other five methods.

James
quelle

It has to to with 'large sample theory' and the biasedness of the MLE estimates, see my answer.

"It's the default in SAS" is not an acceptable answer to a "why" question on this site.

Paul

p-values for mixed models provided by SAS by default are not available by design in lme4 library for R because being untrustworthy (stat.ethz.ch/pipermail/r-help/2006-May/094765.html). So "default SAS" can be even wrong.

Tim