3

Eu gostaria tirar uma dúvida a respeito de modelos lineares generalizados mistos. Sei que a estrutura hierárquica das observações, a independência delas e a existência de amostras desbalanceadas entre outros pode me impedir de utilizar o GLM no lugar do GLMM, mas 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.

Partindo deste, fico com receio de ser inquirido da não possibilidade desta utilização da Logística ou na verdade o erro da logística seria somente nas estimativas dos Betas o que não me impossibilitaria de utilizar as variáveis pré-selecionadas por ela?

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

1 Resposta 1

5

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

inserir a descrição da imagem aqui

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

inserir a descrição da imagem aqui

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))

inserir a descrição da imagem aqui

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

  • H_0: os modelos modelo_treatment.visit.interacao e modelo_treatment.visit são iguais

  • H_1: os modelos modelo_treatment.visit.interacao e modelo_treatment.visit não são iguais

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.

Esta não é a resposta que você está procurando? Pesquise outras perguntas com a tag .