großer Unterschied zwischen den Schätzungen der binomialen Regression bei Einbeziehung des Zufallseffekts und der Nicht-Regression

Ich versuche, die durchschnittliche Punktzahl für zwei Gruppen von Studenten zu schätzen. Ich benutze ein binomiales Regressionsmodell. Dies total_ansist die Gesamtfrage, die sie beantwortet haben und die für verschiedene Schüler unterschiedlich sein kann.

Modell 1 schätzt direkt

model <- glm(cbind(total_correct, total_ans-total_correct) ~ student_type,family= binomial, data = df)

Call:  glm(formula = cbind(total_correct, total_ans - total_correct) ~ student_type, family = binomial, data = df)

Coefficients:
                  (Intercept)  student_group_2  
                      -1.9684           0.2139  

Degrees of Freedom: 1552 Total (i.e. Null);  1551 Residual Null
Deviance:       1480  Residual Deviance: 1477   AIC: 1764


lsmeans(model,~ student_type, type="response")

 student_type           prob         SE df asymp.LCL asymp.UCL
 student_group_1   0.1225627 0.00654160 NA 0.1103074 0.1359715
 student_group_2   0.1474774 0.01275231 NA 0.1241918 0.1742602

In Modell 2 verwende ich einen zufälligen Effekt, um einzelne Abweichungen besser zu berücksichtigen.

model <- glmer(cbind(total_correct, total_ans-total_correct) ~ (1|student)  + student_type, family= binomial, data = sub_df, control=glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE))

Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: binomial  ( logit )
Formula: cbind(total_correct, total_ans - total_correct) ~ (1 | student) +  
    student_type
   Data: sub_df
      AIC       BIC    logLik  deviance  df.resid 
1653.9049 1669.9488 -823.9525 1647.9049      1550 
Random effects:
 Groups Name        Std.Dev.
 student (Intercept) 1.881   
Number of obs: 1553, groups:  author, 1553
Fixed Effects:
                  (Intercept)  student_group_2  
                      -3.0571   0.3915  


lsmeans(model,~ student_type, type="response")

 student_type            prob          SE df  asymp.LCL asymp.UCL
 student_group_1   0.04491007 0.004626728 NA 0.03666574 0.0549025
 student_group_2   0.06503249 0.015117905 NA 0.04097546 0.1017156

Ich bin überrascht, dass es einen so großen Unterschied zwischen den Ergebnissen in beiden Gruppen gibt. Was könnte der Grund dafür sein?

Weitere Informationen: Gruppe 1 hat 1434 Schüler, Gruppe 2 hat 119 Schüler. Dies sind natürlich vorkommende Gruppen

generalized-linear-model lme4-nlme random-effects-model lsmeans user926321
quelle

@PeterFlom, das ist nur eine Möglichkeit, die Anzahl der Erfolge und die Anzahl der Fehler im glm anzugeben. Die Summe der beiden ist die Gesamtzahl der Versuche

user926321

Der Hauptgrund, warum sich die Ergebnisse so stark unterscheiden, ist, dass das zweite Modell zeigt, dass das erste grob falsch ist. Die erste geht von einer Binomialantwort mit konstanter Wahrscheinlichkeit innerhalb jeder Gruppe aus, während die große Standardabweichung der Zufallskomponente (1,8) im zweiten Modell zeigt, dass die Daten nicht mit dieser Annahme einer konstanten Wahrscheinlichkeit übereinstimmen. Histogramme der nach Gruppen korrekten Proportionen sollten dies deutlich machen. Ein paar gut ausgewählte Diagramme der Daten würden helfen, dieses Problem zu lösen.

whuber

Interessanter Thread, obwohl eine wichtige Information zu fehlen scheint: Haben die Schüler in jeder Gruppe den Test mehrmals beantwortet, sodass Sie jedes Mal die Anzahl der richtigen Antworten im Verhältnis zur Anzahl der beantworteten Fragen aufgezeichnet haben? Wenn nicht, verstehe ich nicht, warum Sie einen zufälligen Schülereffekt in Ihr Modell aufnehmen würden.

Isabella Ghement

Das andere Problem ist, dass es sich um natürlich vorkommende Gruppen handelt, sodass sie sich wahrscheinlich in vielerlei Hinsicht unterscheiden (was die Unterschiede zwischen den Gruppen erklären kann). Wenn Sie nach verwirrenden Variablen suchen können, die Sie möglicherweise gemessen haben, kann dies bei der Interpretation Ihrer Ergebnisse hilfreich sein.

Isabella Ghement

Kurzer Kommentar: lsmeans ist veraltet. Sie sollten zu wechseln emmeans. Die Syntax wird höchstwahrscheinlich nahezu identisch sein.

COOLSerdash

Antworten:

Die Frage enthält genügend Informationen, um dies zu klären. lsmeans verwendet einfach die Koeffizienten, um die vorhergesagten Gruppenwahrscheinlichkeiten zu erhalten.

Für das GLM lautet das implizite Modell von OP also:

\hat{π} = \frac{1}{1 + e^{- - (- - 1,9684 + 0,2139 \times d)}}

$\begin{equation} \hat\pi=\frac{1}{1+e^{-(-1.9684+0.2139\times d)}} \end{equation}$

wo $d$ ist ein Indikator für die Mitgliedschaft in Gruppe 2. Die vorhergesagten Wahrscheinlichkeiten sind also: $(1+e^{-(-1.9684)})^{-1}$ und $(1+e^{-(-1.9684+0.2139)})^{-1}$ für die Gruppen 1 bzw. 2. Diese führen zu vorhergesagten Wahrscheinlichkeiten von etwa $12.25\%$ und $14.75\%$ beziehungsweise.

Für das Mehrebenenmodell (oder GLMM) lautet das implizite Modell von OP:

\hat{π} = \frac{1}{1 + e^{- - (- - 3,0571 + 0,3915 \times d + 1,881 \times \hat{u})}}

$\begin{equation} \hat\pi=\frac{1}{1+e^{-(-3.0571+0.3915\times d+1.881\times \hat{u})}} \end{equation}$

wo $\hat{u}$ ist der zufällige Schnittpunkt, der als Standardnormal angenommen wird. Die vorhergesagten Wahrscheinlichkeiten von lsmeans nehmen einen zufälligen Schnittwert von Null an $(\hat{u}=0)$ ergebend: $(1+e^{-(-3.0571)})^{-1}$ und $(1+e^{-(-3.0571+0.3915)})^{-1}$ für die Gruppen 1 bzw. 2. Diese führen zu vorhergesagten Wahrscheinlichkeiten von etwa $4.49\%$ und $6.50\%$ beziehungsweise. Dies sind die lsmeans GLMM-Ergebnisse.

Ein Problem ist, dass der Abschnitt in GLMM die erwarteten logarithmischen Erfolgsaussichten für jemanden ist, der im Vergleich zu anderen "durchschnittlich" ist. Daher ist es ein Problem, dies als Grundlage für die Berichterstattung über das gesamte Modell zu verwenden. In Bezug auf den anderen Koeffizienten ist ein Vorschlag, warum der Gruppendifferenzkoeffizient zunimmt, dass die Modellqualität besser ist, so dass der Koeffizient zunimmt, auf dieser Website nach zusammenklappbar zu suchen oder zu sehen. Ist eine Metaanalyse der Quotenverhältnisse im Wesentlichen hoffnungslos? .

Um die lsmeans GLMM-Ergebnisse mit den lsmeans GLM-Ergebnissen vergleichbar zu machen. Wir müssen die beobachteten Werte des zufälligen Abschnitts verwenden, $\hat{u}$ . Man kann zufällige Abschnitte simulieren, um dem spezifischen Modell von OP zu entsprechen:

set.seed(12345)
u_hat <- rnorm(1553, 0, 1.881)
d <- c(rep(0, 1434), rep(1, 119))
# predicted log odds:
pred_log_odds <- -3.0571 + 0.3915 * d + u_hat
pred_prob <- plogis(pred_log_odds)
coef(lm(pred_prob ~ 0 + factor(d)))

factor(d)0 factor(d)1 
 0.1189244  0.1490395

In diesem simulierten Beispiel liegen diese Werte viel näher an den GLM-Ergebnissen von lsmeans. Wenn Sie für Ihre Daten so etwas wie die folgende Syntax ausführen:

lm(fitted(model) ~ 0 + df$student_type)

Wo modelsich das GLMM befindet, sollten Sie Werte erhalten, die den lsmeans-Werten sehr nahe kommen. Ich gehe davon aus, dass , wenn Sie rufen fitted()auf glmer(), sondern schließt auch die Zufalls abfangen und die Werte , die zurückgegeben werden , sind Wahrscheinlichkeiten.

In Ihrer Situation, in der die Gruppen natürlich vorkommen, sollten Sie in den Daten auch die unterschiedlichen Gruppenvarianzen des zufälligen Abschnitts untersuchen. Ein Modell wie:

glmer(cbind(total_correct, total_ans-total_correct) ~
  (0 + student_type || student) + student_type,
  family = binomial, data = sub_df,
  control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE))

Möglicherweise lohnt es sich, dies zu untersuchen, da Sie derzeit davon ausgehen, dass sich die zufälligen Intercept-Varianzen nicht nach Gruppen unterscheiden. Ich habe das verwendet, ||damit lme4 nicht versucht, die beiden zufälligen Abschnitte zu korrelieren.

Ich wusste nicht, dass man für jeden Schüler einen zufälligen Abschnitt hinzufügen kann, wenn er im Wesentlichen nur eine Zeile in den Daten hat. Aber ich rationalisiere es weg, indem ich annehme, dass der Versuch gegen Fehler pro Zeile mehrere Zeilen in langer Form beträgt.

Heteroskedastischer Jim
quelle

Modell 1, Logistische Regression ohne zufälligen Effekt:

L o g i t (Pr (Y = 1)) = X β

$\mathrm{Logit}(\Pr(Y=1)) = X\beta$

Wir kennen MLE $\hat \beta$ ist asymptotisch unvoreingenommen. Aber

\hat{P} r (Y = 1) = {L o g i t}^{- 1} (X \hat{β})

$\hat Pr(Y=1) = \mathrm{Logit}^{-1}(X\hat\beta)$ ist voreingenommene Schätzung von

Pr (Y = 1)

$\Pr(Y=1)$ , weil die Nichtlinearität der Logit-Funktion. Aber asymptotisch ist es unvoreingenommen. Also für Modell 1,

\hat{P} r (Y = 1) = {L o g i t}^{- 1} (X \hat{β})

$\hat Pr(Y=1) = \mathrm{Logit}^{-1}(X\hat\beta)$ ist akzeptabel.

Modell 2, Logistische Regression mit zufälligem Schnittpunkt:

L. Ö G ich t (Pr ({Y.}_{ich j} = 1)) = {X.}_{ich j} β + γ_{ich}

$\mathrm{Logit}(\Pr(Y_{ij}=1)) = X_{ij}\beta + \gamma_i$

γ \sim N. (0, σ^{2})

$\gamma \sim N(0,\sigma^2)$

In dieser Situation ist der häufigste Fehler

Pr (Y. = 1 | X.) = E. (Y. | X.) = E. ({L. Ö G ich t}^{- - 1} (X. β + γ_{ich})) = {L. Ö G ich t}^{- - 1} (E. (X. β + γ_{ich})) = {L. Ö G ich t}^{- - 1} (X. β)

$\Pr(Y=1|X) = \mathrm{E}(Y|X) = \mathrm{E}(\mathrm{Logit}^{-1}(X\beta + \gamma_i)) =\mathrm{Logit}^{-1}( \mathrm{E}(X\beta + \gamma_i)) = \mathrm{Logit}^{-1}(X\beta)$

In diesem Prozess wird die Nichtlinearität der Logit-Funktion vollständig berücksichtigt.

Die folgende Ausgabe ist also sehr irreführend. Es gibt einen Eindruck, dass für Schüler der Gruppe 1 die Wahrscheinlichkeit einer Korrektur 4,5% beträgt.

    student_type            prob          SE df  asymp.LCL asymp.UCL
    student_group_1   0.04491007 0.004626728 NA 0.03666574 0.0549025
    student_group_2   0.06503249 0.015117905 NA 0.04097546 0.1017156

Der obige Fehler kann nicht zu einer logistischen Regression mit gemischten Effekten beitragen, sondern sollte auf eine Fehlinterpretation der Ergebnisse zurückzuführen sein.

Schauen wir uns an, wie man den Grenzmittelwert von korrekt ableitet $Y$ .

Pr (Y. = 1 | X.) = E. (Y. | X.) = E. ({L. Ö G ich t}^{- - 1} (X. β + γ_{ich})) = \int_{- - \infty}^{\infty} {L. Ö G ich t}^{- - 1} (z) ϕ (z) d z

$\Pr(Y=1|X) = \mathrm{E}(Y|X) = \mathrm{E}(\mathrm{Logit}^{-1}(X\beta + \gamma_i)) =\int_{-\infty}^{\infty} \mathrm{Logit}^{-1}(z)\phi(z)dz$ wo

z = X β + γ_{i}

$z=X\beta + \gamma_i$ und

ϕ

$\phi$ ist pdf von

N (X β, σ^{2})

$N(X\beta, \sigma^2)$

Für Schüler der Gruppe 1 zeigen die Ergebnisse dies $X\beta= -3.0571$ , $\sigma = 1.881$ . Mit der Gauß-Hermite-Quadraturmethode mit 20 Punkten bekam ich $\Pr(Y=1|X=0) = 0.1172$ . ähnlich, $\Pr(Y=1|X=0) = 0.1492$ .

user158565
quelle