Estou tentando estimar a pontuação média para dois grupos de estudantes. Eu uso um modelo de regressão binomial. Essa total_ans
é a pergunta total que eles responderam, que podem ser diferentes para diferentes alunos.
O modelo 1 estima diretamente
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
No modelo 2, uso um efeito aleatório para explicar melhor as variações individuais.
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
Estou surpreso que exista uma diferença tão grande entre os resultados nos dois grupos. Qual pode ser a razão disso?
mais informações: o grupo 1 tem 1434 alunos, o grupo 2 tem 119 alunos. estes são grupos que ocorrem naturalmente
emmeans
. A sintaxe provavelmente será quase idêntica.Respostas:
Há informações suficientes na pergunta para esclarecer isso. lsmeans é simplesmente usar os coeficientes para obter as probabilidades previstas pelo grupo.
Portanto, para o GLM, o modelo implícito do OP é:
onde é um indicador de participação no grupo 2. Portanto, as probabilidades previstas são: ed (1+e−(−1.9684))−1 (1+e−(−1.9684+0.2139))−1 para os grupos 1 e 2, respectivamente. Isso resulta em probabilidades previstas de aproximadamente12.25% e 14.75% respectivamente.
Para o modelo multinível (ou GLMM), o modelo implícito do OP é:
Ondeu^ é a interceptação aleatória assumida como normal normal. As probabilidades previstas de lsmeans assumem um valor de interceptação aleatória igual a zero(u^=0) resultando em: (1+e−(−3.0571))−1 e (1+e−(−3.0571+0.3915))−1 para os grupos 1 e 2, respectivamente. Isso resulta em probabilidades previstas de aproximadamente4.49% e 6.50% respectivamente. Estes são os resultados lsmeans do GLMM.
Um problema é que a interceptação no GLMM é a probabilidade de log esperada de sucesso para alguém que é "médio" em relação a outros. Portanto, usar isso como base para relatar todo o modelo é um problema. Sobre o outro coeficiente, uma sugestão de por que o coeficiente de diferença de grupo aumenta é que a qualidade do modelo é melhor, para que o coeficiente aumente, procure recolhível neste site ou consulte A meta-análise de odds ratio é essencialmente sem esperança? .
Tornar os resultados do GLMM lsmeans comparáveis aos resultados do GLM lsmeans. Temos que usar os valores observados da interceptação aleatória,u^ . Pode-se simular interceptações aleatórias para corresponder ao modelo específico do OP:
Neste exemplo simulado, esses valores estão muito mais próximos dos resultados lsmeans GLM. Se você executar algo parecido com a sintaxe abaixo em seus dados:
Onde
model
está o GLMM, você deve obter valores bem próximos dos valores lsmeans. Estou assumindo que quando você chamarfitted()
emglmer()
, ele também inclui a intercepção aleatória e os valores que são retornados são probabilidades.Na sua situação em que os grupos ocorrem naturalmente, outra coisa a explorar nos dados são as diferentes variações de grupos na interceptação aleatória, portanto, um modelo como:
pode valer a pena explorar, pois, no momento, você está assumindo que as variações de interceptação aleatória não diferem por grupo. Eu usei o
||
so lme4 não tenta correlacionar as duas interceptações aleatórias.Eu não sabia que era possível adicionar uma interceptação aleatória para cada aluno quando eles tinham apenas uma linha nos dados. Mas estou racionalizando, assumindo que o julgamento versus falhas por linha equivale a várias linhas na forma longa.
fonte
Modelo 1, regressão logística sem efeito aleatório:
Conhecemos o MLEβ^ é assintoticamente imparcial. Mas
é uma estimativa enviesada de , devido à não linearidade da função logit. Mas assintoticamente é imparcial. Portanto, para o modelo 1, é aceitável.P^r(Y=1)=Logit−1(Xβ^) Pr(Y=1) P^r(Y=1)=Logit−1(Xβ^)
Modelo 2, regressão logística com interceptação aleatória:
Nesta situação, o erro comum éPr(Y=1|X)=E(Y|X)=E(Logit−1(Xβ+γi))=Logit−1(E(Xβ+γi))=Logit−1(Xβ)
Nesse processo, a não linearidade da função logit é totalmente ignorada.
Portanto, a seguinte saída é muito enganadora. Isso indica que, para o aluno do grupo 1, a probabilidade de correção é de 4,5%.
O erro acima não pode ser contribuído para a regressão logística de efeito misto, deve vir da interpretação incorreta dos resultados.
Vejamos como derivar corretamente a média marginal deY .
Pr(Y=1|X)=E(Y|X)=E(Logit−1(Xβ+γi))=∫∞−∞Logit−1(z)ϕ(z)dz
Onde z=Xβ+γi e ϕ é pdf de N(Xβ,σ2)
Para os alunos do Grupo 1, os resultados indicam queXβ=−3.0571 , σ=1.881 . Usando o método da quadratura de Gauss-Hermite com 20 pontos, obtive Pr(Y=1|X=0)=0.1172 . similarmente,Pr(Y=1|X=0)=0.1492 .
fonte