nos modelos que tenho testado em meu trabalho, o que apresenta o melhor resultado (AIC BIC etc) não é o GLMM que partiu de todas as variáveis que possuo e sim o GLMM que parte das variáveis mais estatisticamente significativas oriundas da regressão logística.
Atenção: menor AIC ou BIC não significa, necessariamente, que o modelo ajustado é o melhor. A fórmula do AIC é dada por
em que k é o número de parâmetros estimados no modelo e L-chapéu é o valor da função de máxima verossimilhança estimada. Como desejamos o modelo com menor AIC, os GLMM, por terem mais parâmetros sendo estimados, sempre serão prejudicados em relação aos GLM.
A lógica é similar para o BIC, cuja fórmula é dada por
Ou seja, em ambos os casos, quanto maior for o número de parâmetros, maior será o valor de 2*k (ou log(n)*k) para os GLMM por, justamente, eles possuírem mais parâmetros sendo estimados. Em estatística, sempre desejamos o modelo mais simples possível para descrever nossos dados.
E apenas o valor do AIC ou BIC não diz muita coisa a respeito do modelo ajustado. Afinal, só dizer que a média de uma amostra é maior do que a outra, sem levar em conta o desvio padrão, por exemplo, não nos garante nada sobre a média populacional. Sempre temos que fazer um teste de hipóteses para confirmar isto. Portanto, o ideal é rodar um teste de hipóteses sobre o AIB ou BIC (mais sobre isso mais tarde, nesta mesma resposta).
Agradeço qualquer orientação que puder me dar para que possa justificar a utilização da regressão logística na seleção das variáveis para o GLMM
Como diria Padre Quevedo, isso non ecziste. Não é possível utilizar resultados de regressão logística (um dos modelos da classe GLM) para selecionar variáveis de um GLMM. São classes de modelos diferentes, com hipóteses diferentes. O que tu pode fazer é encontrar o melhor modelo para os teus dados utilizando apenas regressão logística ou apenas GLMM.
Particularmente, o que sugiro fazer quando ajustamos um GLMM a um conjunto de dados é começar com o modelo mais complexo possível, com todas as variáveis preditoras definidas no experimento, e ir simplificando o modelo a partir disso. Minha filosofia de ajuste de modelos, calcada basicamente na minha experência e no que aprendi durante a minha formação, é sempre ter ao final o modelo mais simples possível. A maneira que vou descrever abaixo faz justamente isto.
Vou fazer uma análise de dados a respeito de resultados em uma tratamento para fungos nas unhas do pé. Estes dados podem ser obtidos a partir do próprio R
, no data frame HSAUR2::toenail
:
library(HSAUR2)
str(toenail)
'data.frame': 1908 obs. of 5 variables:
$ patientID: Factor w/ 294 levels "1","2","3","4",..: 1 1 1 1 1 1 1 2 2 2 ...
$ outcome : Factor w/ 2 levels "none or mild",..: 2 2 2 1 1 1 1 1 1 2 ...
$ treatment: Factor w/ 2 levels "itraconazole",..: 2 2 2 2 2 2 2 1 1 1 ...
$ time : num 0 0.857 3.536 4.536 7.536 ...
$ visit : int 1 2 3 4 5 6 7 1 2 3 ...
As variáveis consideradas neste exemplo foram:
outcome: o quão separada estava a unha do dedo, com dois níveis (variável resposta)
treatment: qual remédio foi utilizado no tratamento (efeito fixo)
visit: número da visita do paciente para o tratamento (efeito fixo)
patientID: número de identificação do paciente (efeito aleatório)
Em primeiro lugar, vou fazer uma análise exploratória dos dados. Limitei a apenas 12 pacientes, de modo que a visualização ficasse melhor:
library(ggplot2)
library(dplyr)
toenail %>%
filter(patientID %in% head(unique(toenail$patientID), 12)) %>%
ggplot(., aes(x=visit, y=outcome, colour=treatment)) +
geom_point() +
facet_wrap(~ patientID) +
labs(x="Visita", y="Resultado", colour="Tratamento") +
scale_x_continuous(breaks=c(1, 3, 5, 7))
Claramente, a variável resposta deve ser considerada como binomial. Entretanto, este experimento possui medidas repetidas, pois cada paciente fez até sete visitas para tratamento. Assim o melhor é ajustar um GLMM a estes dados.
Em primeiro lugar, vamos ajustar um modelo com interação entre treatment
e visit
, considerando patientID
como efeito aleatório:
library(lme4)
modelo_treatment.visit.interacao <- glmer(outcome ~ treatment*visit + (1|patientID), data=toenail, family=binomial)
Em seguida, vamos ajustar um modelo sem interação entre treatment
e visit
:
modelo_treatment.visit <- glmer(outcome ~ treatment + visit + (1|patientID), data=toenail, family=binomial)
Como estes modelos são aninhados, podemos comparar eles através de um teste de razão de verossimilhança:
anova(modelo_treatment.visit.interacao, modelo_treatment.visit, test="Chisq")
Data: toenail
Models:
modelo_treatment.visit: outcome ~ treatment + visit + (1 | patientID)
modelo_treatment.visit.interacao: outcome ~ treatment * visit + (1 | patientID)
Df AIC BIC logLik deviance Chisq
modelo_treatment.visit 4 1260.3 1282.6 -626.17 1252.3
modelo_treatment.visit.interacao 5 1258.8 1286.5 -624.39 1248.8 3.5627
Chi Df Pr(>Chisq)
modelo_treatment.visit
modelo_treatment.visit.interacao 1 0.05909 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
O comando acima está testando as hipóteses
Como o p-valor para a comparação entre os modelos foi maior do que 0,05 (foi 0.05909), podemos assumir que estes modelos são iguais. Assim, a interação não é necessária neste caso, pois é preferível ter um modelo mais simples. A seguir, precisamos descobrir se treatment
ou visit
são necessários para o modelo.
Até o momento, o modelo escolhido possui treatment
e visit
, sem interação. Vamos ajustar um modelo sem treatment
e compará-lo com nosso modelo atual:
modelo_visit <- glmer(outcome ~ visit + (1|patientID), data=toenail, family=binomial)
anova(modelo_treatment.visit, modelo_visit, test="Chisq")
Data: toenail
Models:
modelo_visit: outcome ~ visit + (1 | patientID)
modelo_treatment.visit: outcome ~ treatment + visit + (1 | patientID)
Df AIC BIC logLik deviance Chisq Chi Df
modelo_visit 3 1259.4 1276.0 -626.69 1253.4
modelo_treatment.visit 4 1260.3 1282.6 -626.17 1252.3 1.0427 1
Pr(>Chisq)
modelo_visit
modelo_treatment.visit 0.3072
Novamente, o p-valor foi maior do que 0,05. Portanto, os modelos não são diferentes. Não necessitamos de treatment
para explicar a variável resposta. Vamos testar agora o modelo apenas com treatment
e ver o que acontece:
modelo_treatment <- glmer(outcome ~ treatment + (1|patientID), data=toenail, family=binomial)
anova(modelo_treatment.visit, modelo_treatment, test="Chisq")
Data: toenail
Models:
modelo_treatment: outcome ~ treatment + (1 | patientID)
modelo_treatment.visit: outcome ~ treatment + visit + (1 | patientID)
Df AIC BIC logLik deviance Chisq Chi Df
modelo_treatment 3 1572.6 1589.2 -783.29 1566.6
modelo_treatment.visit 4 1260.3 1282.6 -626.17 1252.3 314.25 1
Pr(>Chisq)
modelo_treatment
modelo_treatment.visit < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Agora vemos que há diferença entre o modelo com treatment
e visit
e o modelo com apenas treatment
.
Juntando os dois últimos resultados, optamos por ficar com o modelo com apenas visit
.
Por fim, basta descobrir se visit
é significante. Para isto, vamos ajusta o modelo apenas com o intercepto e fazer o mesmo que fizemos em todos os casos anteriores:
modelo_intercepto <- glmer(outcome ~ 1 + (1|patientID), data=toenail, family=binomial)
anova(modelo_visit, modelo_intercepto, test="Chisq")
Data: toenail
Models:
modelo_intercepto: outcome ~ 1 + (1 | patientID)
modelo_visit: outcome ~ visit + (1 | patientID)
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
modelo_intercepto 2 1571.5 1582.6 -783.74 1567.5
modelo_visit 3 1259.4 1276.0 -626.69 1253.4 314.09 1 < 2.2e-16
modelo_intercepto
modelo_visit ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Como o p-valor é menor do que 0,05, concluímos que há diferença entre os modelos ajustados. Portanto, visit
é significante. Ou seja, o que importa, neste caso, é ir visitar regularmente o médico, não importando qual remédio utilizar.
O método que ilustrei acima, com os testes de razão de verossimilhança (likelihood ratio test - LRT -, em inglês) servem para a comparação de quaisquer tipos de modelos, de regressão linear múltipla a GLMM. Recomendo procurar o livro Mixed Effects Models and Extensions in Ecology with R, Zuur et. al. (2009) para ver uma referência mais aprofundada sobre isso. O livro inclusive discute o fato dos p-valores obtidos com este método não serem mito confiáveis (na verdade, quase nenhum p-valor no contexto de GLMM é confiável).
A própria discussão a respeito de como ajustar modelos lineares generalizados mistos é bem sofisticada. Tem gente, como Douglas Bates, que defende não calcular p-valores. Em alguns casos, os resultados do pacote lme4
não apresentam p-valores para os testes realizados. Isto se dá porque um p-valor é, em essência, uma probabilidade associada à distribuição de alguma variável aleatória. No caso do ANOVA tradicional, esta distribuição é a F. No caso de um modelo misto, é provável que não saibamos a distribuição da estatística do teste. Portanto, não faz sentido, em alguns ajustes de modelos com dados faltantes ou desbalanceados, calcular p-valores para os testes realizados.
No teu caso específico, caso o editor da revista exija que p-valores sejam reportados, eu faria uma pesquisa a respeito dos últimos artigos publicados na área e veria como os resultados foram reportados neles. Eu tentaria reproduzir a análise realizada nestes artigos, nem que fosse pedindo o conjunto de dados original para os autores, e veria como aplicar os resultados destes artigos no contexto da tua pesquisa.
Uma boa fonte gratuita para aprender um pouco mais sobre o assunto é este FAQ do Ben Bolker. Há alguns livros que tratam do assunto, mas todos os que li são muito teóricos, demandando um conhecimento em estatística bastante avançado.