Carregamento de pacotes:

Neste exemplo são utilizados os seguintes pacotes: tibble, VGAM e lmtest. Para a confecção do documento, é necessário o pacote rmarkdown.

# tibble -- para facilitar manipulacao de dados
if (!require('tibble')) install.packages('tibble'); library('tibble')

# VGAM -- para trabalhar com Regresão Logística Politômica
if (!require('VGAM')) install.packages('VGAM'); library('VGAM')

# lmtest -- para os testes de razao de verossimilhancas e para o teste de wald
if (!require('lmtest')) install.packages('lmtest'); library('lmtest')

Carregamento dos dados:

A variável resposta corresponde a deterioração mental do paciente, com graus variando da seguinte forma: “Well” (Bom), “Mild” (Ameno), “Moderate” (Moderado) e “Impaired” (Deteriorado). As variáveis explicativas são: “SES”, que corresponde ao nível social do paciente e “life events”, que corresponde à contagem dos eventos importantes na vida do paciente.

#Entrada de dados
mental <- c(rep("Well", 12), rep("Mild", 12), rep("Moderate", 7),
            rep("Impaired", 9))

mental <- ordered(mental, levels=c("Well", "Mild", "Moderate", "Impaired"))

SES <- as.integer(c(1,1,1,1,0,1,0,1,1,1,0,0,1,0,1,0,1,1,0,1,1,0,1,1,0,
         1,0,0,1,0,0,1,1,1,0,0,0,1,0,0))

life_events <- as.integer(c(1,9,4,3,2,0,1,3,3,7,1,2,5,6,3,1,8,2,5,5,9,
                 3,3,1,0,4,3,9,6,4,3,8,2,7,5,4,4,8,8,9))

dadosmental <- data.frame(mental, SES, life_events)

tibble::as.tibble(dadosmental)

Ajuste do Modelo:

Utilizou-se o pacote VGAM (https://cran.r-project.org/web/packages/VGAM/index.html) para modelagem dos dados. Tal utilização se fez necessária tendo em vista que a função “glm”, nativa do R, não possui suporte à Regressão Logística Politômica.

modmental <- VGAM::vglm(formula=mental~life_events+SES,
                      family=VGAM::cumulative(link="logit", parallel=T, reverse=F),data=dadosmental)

VGAM::summary(modmental)
## 
## Call:
## VGAM::vglm(formula = mental ~ life_events + SES, family = VGAM::cumulative(link = "logit", 
##     parallel = T, reverse = F), data = dadosmental)
## 
## 
## Pearson residuals:
##                   Min      1Q  Median     3Q   Max
## logit(P[Y<=1]) -1.568 -0.7048 -0.2102 0.8070 2.713
## logit(P[Y<=2]) -2.328 -0.4666  0.2657 0.6904 1.615
## logit(P[Y<=3]) -3.688  0.1198  0.2039 0.4194 1.892
## 
## Coefficients: 
##               Estimate Std. Error z value Pr(>|z|)   
## (Intercept):1  -0.2819     0.6231  -0.452  0.65096   
## (Intercept):2   1.2128     0.6511   1.863  0.06251 . 
## (Intercept):3   2.2094     0.7171   3.081  0.00206 **
## life_events    -0.3189     0.1194  -2.670  0.00759 **
## SES             1.1112     0.6143   1.809  0.07045 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Number of linear predictors:  3 
## 
## Names of linear predictors: 
## logit(P[Y<=1]), logit(P[Y<=2]), logit(P[Y<=3])
## 
## Residual deviance: 99.0979 on 115 degrees of freedom
## 
## Log-likelihood: -49.5489 on 115 degrees of freedom
## 
## Number of iterations: 5 
## 
## No Hauck-Donner effect found in any of the estimates
## 
## Exponentiated coefficients:
## life_events         SES 
##   0.7269742   3.0380707
VGAM::AIC(modmental) #AIC
## [1] 109.0979
-2*logLik(modmental) #-2*loglik
## [1] 99.0979

Intervalos de confiança e testes de hipóteses:

Seguem, respectivamente, testes de hipóteses para testar se \(\beta = 0\) e os intervalos de confiança dos coeficientes ajustados. Procedeu-se aos testes de Razão de Verossimilhança e Wald.

#Testes para beta
VGAM::lrtest(modmental) #Razão de verossimilhanças
## Likelihood ratio test
## 
## Model 1: mental ~ life_events + SES
## Model 2: mental ~ 1
##   #Df  LogLik Df  Chisq Pr(>Chisq)   
## 1 115 -49.549                        
## 2 117 -54.521  2 9.9442   0.006929 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
lmtest::waldtest(modmental) #Wald.
#Intervalos de confiança
confint.default(modmental) #Para os parâmetros
##                     2.5 %      97.5 %
## (Intercept):1 -1.50305171  0.93927761
## (Intercept):2 -0.06336929  2.48897799
## (Intercept):3  0.80388808  3.61487537
## life_events   -0.55293773 -0.08479079
## SES           -0.09272235  2.31516768
exp(confint.default(modmental)) #Para a chance
##                   2.5 %     97.5 %
## (Intercept):1 0.2224503  2.5581328
## (Intercept):2 0.9385968 12.0489557
## (Intercept):3 2.2342109 37.1467162
## life_events   0.5752574  0.9187045
## SES           0.9114465 10.1266209

Qualidade do ajuste:

O pacote VGAM não apresenta em sua saída padrão o teste score para “odds” proporcionais como ocorre no SAS. Sendo assim, procedeu-se ao cálculo “manual” deste.

#Ajustamento - Pelo teste score para chances proporcionais
modmental2 <- VGAM::vglm(formula=mental~life_events+SES, 
                         family=VGAM::cumulative(link="logit", parallel=F, reverse=F),data=dadosmental)
# alteramos o parametro parallel para FALSE

pchisq(deviance(modmental)-deviance(modmental2),
       df=df.residual(modmental)-df.residual(modmental2),lower.tail=FALSE)
## [1] 0.6718069
#O modelo mais complexo não se ajusta melhor, devendo o mais simples ser preferido.

Modelo com Interação

Ajuste do Modelo:

modmental_int <- VGAM::vglm(formula=mental~life_events*SES,
                      family=VGAM::cumulative(link="logit", parallel=T, reverse=F),data=dadosmental)

VGAM::summary(modmental_int)
## 
## Call:
## VGAM::vglm(formula = mental ~ life_events * SES, family = VGAM::cumulative(link = "logit", 
##     parallel = T, reverse = F), data = dadosmental)
## 
## 
## Pearson residuals:
##                   Min      1Q  Median     3Q   Max
## logit(P[Y<=1]) -1.393 -0.7139 -0.2172 0.9084 2.262
## logit(P[Y<=2]) -2.758 -0.4862  0.2781 0.7218 1.797
## logit(P[Y<=3]) -3.364  0.1347  0.2062 0.3795 2.344
## 
## Coefficients: 
##                 Estimate Std. Error z value Pr(>|z|)   
## (Intercept):1    0.09807    0.81102   0.121  0.90375   
## (Intercept):2    1.59248    0.83717   1.902  0.05714 . 
## (Intercept):3    2.60660    0.90966   2.865  0.00416 **
## life_events     -0.42045    0.19031  -2.209  0.02715 * 
## SES              0.37090    1.13022   0.328  0.74279   
## life_events:SES  0.18131    0.23611   0.768  0.44255   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Number of linear predictors:  3 
## 
## Names of linear predictors: 
## logit(P[Y<=1]), logit(P[Y<=2]), logit(P[Y<=3])
## 
## Residual deviance: 98.5044 on 114 degrees of freedom
## 
## Log-likelihood: -49.2522 on 114 degrees of freedom
## 
## Number of iterations: 5 
## 
## No Hauck-Donner effect found in any of the estimates
## 
## Exponentiated coefficients:
##     life_events             SES life_events:SES 
##       0.6567529       1.4490350       1.1987822
VGAM::AIC(modmental_int) #AIC
## [1] 110.5044
-2*logLik(modmental_int) #-2*loglik
## [1] 98.50444

Testes de hipóteses:

Seguem, respectivamente, testes de hipóteses para testar se \(\beta = 0\) e os intervalos de confiança dos coeficientes ajustados. Procedeu-se aos testes de Razão de Verossimilhança e Wald.

#Testes para beta
VGAM::lrtest(modmental_int) #Razão de verossimilhanças
## Likelihood ratio test
## 
## Model 1: mental ~ life_events * SES
## Model 2: mental ~ 1
##   #Df  LogLik Df  Chisq Pr(>Chisq)  
## 1 114 -49.252                       
## 2 117 -54.521  3 10.538    0.01451 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
lmtest::waldtest(modmental_int) #Wald.

Referências Bibliográficas

AGRESTI, ALAN. An Introduction to Categorical Data Analysis. John Wiley & Sons, second edition, 2007.

NOTAS DE AULA. Análise de Dados Categorizados. Curso de Graduação em Estatística, UnB, 2018.

R CORE TEAM. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.