---

---

Carregamento de pacotes:

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

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

# vcd -- para trabalhar com tabelas de contingencia
if (!require('vcd')) install.packages('vcd'); library('vcd')

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

# knitr -- para trabalhar com tabelas no markdown
if (!require('knitr')) install.packages('knitr'); library('knitr')

Carregamento dos dados:

A variável resposta corresponde a ideologia do entrevistado, segundo graus de liberalismo e conservadorismo. Os níveis são, em ordem: “Very Liberal” (Muito Liberal), “Slightly Liberal” (Levemente Liberal), “Moderate” (Moderado), “Slightly Conservative” (Levemente Conservador) e “Very Conservative” (Muito Conservador). As variáveis explicativas são o partido ao qual o entrevistado pertence, cujos níveis são “Republicanos” e “Democratas”. Padronizou-se o valor 0 e 1 respectivamente para o partido. A outra variável explicativa é o gênero, separado em “Female” (Feminino) e “Male” (Masculino). Embora tal variável tenha sido inserida, não foi considerada para a modelagem dos dados.

ideologia <- factor(c(rep("Very Liberal", 44+18+36+12), rep("Slightly Liberal", 47+28+34+18),
               rep("Moderate", 118+86++53+62), rep("Slightly Conservative", 23+39+18+45),
               rep("Very Conservative", 32+48+23+51)), ord=T)

ideologia <- ordered(ideologia, levels=c("Very Liberal", #Ordenamento dos fatores
                            "Slightly Liberal", 
                            "Moderate", "Slightly Conservative", 
                            "Very Conservative"))

genero <- as.factor(c(rep("Female", 44+18), rep("Male", 36+12),
                      rep("Female", 47+28), rep("Male", 34+18), 
                      rep("Female", 118+86), rep("Male", 53+62), 
                      rep("Female", 23+39), rep("Male", 18+45), 
                      rep("Female", 32+48), rep("Male", 23+51)))

partido_politico <- as.factor(c(rep("Democratic", 44), rep("Republican", 18),
                                rep("Democratic", 36), rep("Republican", 12),
                                rep("Democratic", 47), rep("Republican", 28),
                                rep("Democratic", 34), rep("Republican", 18),
                                rep("Democratic", 118), rep("Republican", 86),
                                rep("Democratic", 53), rep("Republican", 62),
                                rep("Democratic", 23), rep("Republican", 39),
                                rep("Democratic", 18), rep("Republican", 45),
                                rep("Democratic", 32), rep("Republican", 48),
                                rep("Democratic", 23), rep("Republican", 51)))

partido_politico <- factor(partido_politico, levels=c("Republican", "Democratic")) #Ordenamento dos fatores
#para ajuste conforme o exemplo do livro.

dadosideo <- data.frame(ideologia, genero, partido_politico) #Consolidação dos dados

Tabela de Contingência:

Utilizou-se o pacote vcd (https://cran.r-project.org/web/packages/vcd/index.html) para a criação da tabela de contingência.

tabideo <- vcd::structable(cbind(dadosideo[1:2], #Tabela de contingência
                                 partido_politico=ordered(partido_politico, 
                                                          levels=c("Democratic", "Republican"))), 
                                                     split_vertical = c(T,F,F))

knitr::kable(as.matrix(tabideo))

	Very Liberal	Slightly Liberal	Moderate	Slightly Conservative	Very Conservative
Female_Democratic	44	47	118	23	32
Female_Republican	18	28	86	39	48
Male_Democratic	36	34	53	18	23
Male_Republican	12	18	62	45	51

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.

modideo <- VGAM::vglm(formula=ideologia~partido_politico, #Ajuste do modelo
                      family=VGAM::cumulative(link="logit", parallel=T, reverse=F),data=dadosideo)


VGAM::summary(modideo) #Principais resultados

## 
## Call:
## VGAM::vglm(formula = ideologia ~ partido_politico, family = VGAM::cumulative(link = "logit", 
##     parallel = T, reverse = F), data = dadosideo)
## 
## 
## Pearson residuals:
##                    Min      1Q  Median      3Q   Max
## logit(P[Y<=1]) -1.0462 -0.3161 -0.2172 -0.1442 3.309
## logit(P[Y<=2]) -0.9327 -0.6669 -0.2964  0.5826 2.628
## logit(P[Y<=3]) -2.4630 -0.5430  0.3121  0.7122 1.024
## logit(P[Y<=4]) -2.6380  0.1444  0.2589  0.3761 1.236
## 
## Coefficients: 
##                            Estimate Std. Error z value Pr(>|z|)    
## (Intercept):1              -2.46899    0.13182 -18.730  < 2e-16 ***
## (Intercept):2              -1.47454    0.10909 -13.517  < 2e-16 ***
## (Intercept):3               0.23712    0.09485   2.500   0.0124 *  
## (Intercept):4               1.06954    0.10457  10.228  < 2e-16 ***
## partido_politicoDemocratic  0.97451    0.12905   7.551 4.31e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Number of linear predictors:  4 
## 
## Names of linear predictors: 
## logit(P[Y<=1]), logit(P[Y<=2]), logit(P[Y<=3]), logit(P[Y<=4])
## 
## Residual deviance: 2474.985 on 3335 degrees of freedom
## 
## Log-likelihood: -1237.493 on 3335 degrees of freedom
## 
## Number of iterations: 4 
## 
## No Hauck-Donner effect found in any of the estimates
## 
## Exponentiated coefficients:
## partido_politicoDemocratic 
##                   2.649882

VGAM::AIC(modideo) #AIC

## [1] 2484.985

-2*logLik(modideo) #-2*loglik

## [1] 2474.985

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(modideo) #Razão de verossimilhanças

## Likelihood ratio test
## 
## Model 1: ideologia ~ partido_politico
## Model 2: ideologia ~ 1
##    #Df  LogLik Df  Chisq Pr(>Chisq)    
## 1 3335 -1237.5                         
## 2 3336 -1266.8  1 58.645  1.888e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

lmtest::waldtest(modideo) #Wald.

#Intervalos de confiança
confint.default(modideo) #Para os parâmetros

##                                  2.5 %     97.5 %
## (Intercept):1              -2.72735262 -2.2106313
## (Intercept):2              -1.68834585 -1.2607401
## (Intercept):3               0.05122329  0.4230262
## (Intercept):4               0.86458202  1.2744934
## partido_politicoDemocratic  0.72157568  1.2274543

exp(confint.default(modideo)) #Para a chance

##                                 2.5 %    97.5 %
## (Intercept):1              0.06539218 0.1096314
## (Intercept):2              0.18482500 0.2834442
## (Intercept):3              1.05255790 1.5265743
## (Intercept):4              2.37401359 3.5768890
## partido_politicoDemocratic 2.05767290 3.4125311

Qualidade do ajuste:

Por todas as variável se tratarem de dados categóricos, é possível calcular-se os resíduos de Pearson e resíduos Deviance, para realizar teste de hipótese para a qualidade do ajustamento. Não se encontrou nos pacotes utilizados função para a realização direta destes testes. Sendo assim, procedeu-se ao cálculo dos coeficientes e probabilidades acumuladas e em cada nível segundo o modelo ajustado. Com isto, foi possível calcular os resíduos de Pearson e Deviance, realizando-se em seguida um teste qui-quadrado para qualidade do ajustamento.

#Coeficientes
alfa1 <- VGAM::coefvgam(modideo)[1]
alfa2 <- VGAM::coefvgam(modideo)[2]
alfa3 <- VGAM::coefvgam(modideo)[3]
alfa4 <- VGAM::coefvgam(modideo)[4]
beta <- VGAM::coefvgam(modideo)[5]

#Probabilidades acumuladas - Democratas (1)
dpmenorq1 <- exp(alfa1+beta)/(1+exp(alfa1+beta)) #P(Y <= 1)
dpmenorq2 <- exp(alfa2+beta)/(1+exp(alfa2+beta)) #P(Y <= 2)
dpmenorq3 <- exp(alfa3+beta)/(1+exp(alfa3+beta)) #P(Y <= 3)
dpmenorq4 <- exp(alfa4+beta)/(1+exp(alfa4+beta)) #P(Y <= 4)

#Probabilidades - Democratas(1)
dpigual1 <- dpmenorq1 #P(Y = 1)
dpigual2 <- dpmenorq2-dpmenorq1 #P(Y = 2)
dpigual3 <- dpmenorq3-dpmenorq2 #P(Y = 3)
dpigual4 <- dpmenorq4-dpmenorq3 #P(Y = 4)
dpigual5 <- 1-dpmenorq4 #P(Y = 5)

probd <- c(dpigual1,dpigual2,dpigual3,dpigual4,dpigual5)
sum(c(dpigual1,dpigual2,dpigual3,dpigual4,dpigual5)) #OK - Soma das prob. = 1

## [1] 1

#Probabilidades acumuladas - Republicanos (0)
rpmenorq1 <- exp(alfa1)/(1+exp(alfa1)) #P(Y <= 1)
rpmenorq2 <- exp(alfa2)/(1+exp(alfa2)) #P(Y <= 2)
rpmenorq3 <- exp(alfa3)/(1+exp(alfa3)) #P(Y <= 3)
rpmenorq4 <- exp(alfa4)/(1+exp(alfa4)) #P(Y <= 4)

#Probabilidades - Republicanos(0)
rpigual1 <- rpmenorq1 #P(Y = 1)
rpigual2 <- rpmenorq2-rpmenorq1 #P(Y = 2)
rpigual3 <- rpmenorq3-rpmenorq2 #P(Y = 3)
rpigual4 <- rpmenorq4-rpmenorq3 #P(Y = 4)
rpigual5 <- 1-rpmenorq4 #P(Y = 5)

probr <- c(rpigual1,rpigual2,rpigual3,rpigual4,rpigual5)
sum(c(rpigual1,rpigual2,rpigual3,rpigual4,rpigual5)) #OK - Soma das prob. = 1

## [1] 1

#Ajustamento

#Valor esperado e residuo de Pearson - Democratas
observadosd <- apply(as.matrix(tabideo)[c(1,3),], 2, sum)
esperadosd <- sum(observadosd)*probd
pearsond <- sum((observadosd-esperadosd)^2/esperadosd)

#Valor esperado e residuo de Pearson - Republicanos
observadosr <- apply(as.matrix(tabideo)[c(2,4),], 2, sum)
esperadosr <- sum(observadosr)*probr
pearsonr <- sum((observadosr-esperadosr)^2/esperadosr)

#Residuo de Pearson e teste para o ajuste
pearson <- pearsond+pearsonr
1-pchisq(pearson, df=3) #Bom ajuste pelos residuos de Pearson

## [1] 0.3002341

#Residuos Deviance e teste para o ajuste
devianced <- 2*sum(observadosd*log(observadosd/esperadosd))
deviancer <- 2*sum(observadosr*log(observadosr/esperadosr))
deviance <- devianced+deviancer

1-pchisq(deviance, df=3) #Bom ajuste pelos residuos Deviance

## [1] 0.2972241

---

---