---

---

Carregamento de pacotes:

Neste exemplo são utilizados os seguintes pacotes: tibble, ggplot2, plotly, 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')

# ggplot2 -- para os plots
if (!require('ggplot2')) install.packages('ggplot2'); library('ggplot2')

# plotly -- para os plots
if (!require('plotly')) install.packages('plotly'); library('plotly')

# 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 comp corresponde ao comprimento do aligator, a variável comida o tipo de comida que ele ingere. A primeira é variável explicativa, a segunda resposta. A respeito da comida, “I” significa invertebrado, “F” peixe, (do inglês “Fish”) e “O” significa outros, (do inglês “Other”).

comp <- c(1.24,1.45,1.63,1.78,1.98,2.36,2.79,3.68,
          1.30,1.45,1.65,1.78,2.03,2.39,2.84,3.71,
          1.30,1.47,1.65,1.78,2.03,2.41,3.25,3.89,
          1.32,1.47,1.65,1.80,2.16,2.44,3.28,
          1.32,1.50,1.65,1.80,2.26,2.46,3.33,
          1.40,1.52,1.68,1.85,2.31,2.56,3.56,
          1.42,1.55,1.70,1.88,2.31,2.67,3.58,
          1.42,1.60,1.73,1.93,2.36,2.72,3.66)
comida <- as.factor(c("I","I","I","I","I","F","F","O",
            "I","O","O","I","F","F","F","F",
            "I","I","I","O","F","F","O","F",
            "F","F","F","I","F","F","O",
            "F","I","F","F","F","F","F",
            "F","I","F","F","F","O","F",
            "I","I","I","I","F","F","F",
            "F","I","O","I","F","I","F"))
alligator <- data.frame(comida, comp)
tibble::as.tibble(alligator)

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. Como a saída do SAS utiliza como base a categoria “O”, utilizamos a função relevel(comida, ref="O") para que o modelo no R também tenha como base a mesma categoria para a variável comida.

alligator <- data.frame(relevel(comida, ref="O"), comp)

modall <- VGAM::vglm(comida~comp,VGAM::multinomial,alligator)
VGAM::summary(modall) 

## 
## Call:
## VGAM::vglm(formula = comida ~ comp, family = VGAM::multinomial, 
##     data = alligator)
## 
## 
## Pearson residuals:
##                       Min      1Q  Median     3Q   Max
## log(mu[,1]/mu[,3]) -2.330 -0.5075  0.5538 0.6836 1.452
## log(mu[,2]/mu[,3]) -2.687 -0.4821 -0.1653 0.7093 3.439
## 
## Coefficients: 
##               Estimate Std. Error z value Pr(>|z|)   
## (Intercept):1   1.6177     1.3073   1.237  0.21591   
## (Intercept):2   5.6974     1.7937   3.176  0.00149 **
## comp:1         -0.1101     0.5171  -0.213  0.83137   
## comp:2         -2.4654     0.8996      NA       NA   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Number of linear predictors:  2 
## 
## Names of linear predictors: log(mu[,1]/mu[,3]), log(mu[,2]/mu[,3])
## 
## Residual deviance: 98.3412 on 114 degrees of freedom
## 
## Log-likelihood: -49.1706 on 114 degrees of freedom
## 
## Number of iterations: 5 
## 
## Warning: Hauck-Donner effect detected in the following estimate(s):
## 'comp:2'
## 
## Reference group is level  3  of the response

VGAM::AIC(modall) #AIC

## [1] 106.3412

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

## [1] 98.34124

Preparação dos dados para plotagem:

A fim de propiciar a plotagem do gráfico das probabilidades x comprimento por meio do pacote ggplot2 (https://cran.r-project.org/web/packages/ggplot2/index.html), obteve-se os valores dos coeficientes do modelo ajustado a fim de que se pudesse calcular as probabilidades em um espectro com maior quantidade de valores para a variável explicativa, o que será importante para suavizar o gráfico.

#Preparação dos dados - Amigável pro ggplot2
xx <- 1.24 + (0:249)*(3.89-1.24)/249
coeficientes <- c(VGAM::coefvgam(modall),VGAM::coefvgam(modall)[1]-VGAM::coefvgam(modall)[2],
                  VGAM::coefvgam(modall)[3]-VGAM::coefvgam(modall)[4])
preditores1 <- coeficientes[1]+coeficientes[3]*xx
preditores2 <- coeficientes[2]+coeficientes[4]*xx
preditores3 <- coeficientes[5]+coeficientes[6]*xx

# probabilidades e proporções (para o plot)
prob1=exp(preditores1)/(1+exp(preditores1)+exp(preditores2))
prob2=exp(preditores2)/(1+exp(preditores1)+exp(preditores2))
prob3=1/(1+exp(preditores1)+exp(preditores2))

#Consolidação dos dados
dados <- data.frame(prob=c(prob1,prob2,prob3),rep(xx,3),tipo=rep(c("F","I","O"),each=length(xx)))
names(dados) <- c("prob", "x", "Alimento")

Plotagem dos gráficos:

plotly_pallette <- c('#1F77B4', '#FF7F0E', '#2CA02C', '#D62728')

#Plotagem dos gráficos
library(ggplot2)
p <- ggplot(dados,ggplot2::aes(x=dados$x, y=dados$prob, colour=Alimento)) + 
  geom_line(size=0.75)+
  ggtitle("Curvas Logísticas Ajustadas - Peixe, Invertebrados e Outros")+
  xlab("Comprimento") + ylab("Probabilidade") +
  theme_bw()+
  theme(panel.border = element_blank())+ # para ficar igual o plotly
  # guides(color=guide_legend(title=NULL))+
  scale_color_manual(values=plotly_pallette)+
  theme(plot.title = element_text(hjust=0.5))
  
plotly::ggplotly(p)

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

## Likelihood ratio test
## 
## Model 1: comida ~ comp
## Model 2: comida ~ 1
##   #Df  LogLik Df  Chisq Pr(>Chisq)    
## 1 114 -49.171                         
## 2 116 -57.571  2 16.801  0.0002248 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

lmtest::waldtest(modall) #Wald.

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

##                    2.5 %     97.5 %
## (Intercept):1 -0.9444737  4.1799367
## (Intercept):2  2.1817746  9.2131138
## comp:1        -1.1235700  0.9033519
## comp:2        -4.2286348 -0.7022574

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

##                    2.5 %       97.5 %
## (Intercept):1 0.38888418 6.536172e+01
## (Intercept):2 8.86201848 1.002777e+04
## comp:1        0.32511704 2.467861e+00
## comp:2        0.01457227 4.954656e-01

Análise de resíduos:

Calculou-se os resíduos padronizados a fim de observar o comportamento destes e realizar conjecturas acerca da qualidade do ajuste do modelo.

#Calculo dos resíduos padronizados para log(pi1/pi3)
resnorm1 <- VGAM::residuals(modall)[,1]/sd(VGAM::residuals(modall)[,1])

#Calculo dos resíduos padronizados para log(pi2/pi3)
resnorm2 <- VGAM::residuals(modall)[,2]/sd(VGAM::residuals(modall)[,2]) 

#Calculo dos resíduos padronizados para log(pi1/pi2)
pi1pi2 <- VGAM::residuals(modall)[,1]-VGAM::residuals(modall)[,2]
resnorm3 <- pi1pi2/sd(pi1pi2)

#Gráficos
p1 <- ggplot(as.data.frame(resnorm1), aes(y=resnorm1, x=alligator$comp)) + 
  geom_point(colour=plotly_pallette[1], alpha=0.7) +
  geom_hline(yintercept = 0, size=0.4, alpha=0.5)+
  ggtitle("Resíduos Padronizados - log[pi1/pi3]")+
  xlab("Comprimento") + ggplot2::ylab("Valor") + 
  theme_bw()+
  theme(panel.border = element_blank())+ # para ficar igual o plotly
  # guides(color=guide_legend(title=NULL))+
  scale_color_manual(values=plotly_pallette)+
  theme(plot.title = element_text(hjust=0.5))
  
plotly::ggplotly(p1)

p2 <- ggplot(as.data.frame(resnorm2), aes(y=resnorm2, x=alligator$comp)) + 
  geom_point(colour=plotly_pallette[4], alpha=0.7) +
  geom_hline(yintercept = 0, size=0.4, alpha=0.5)+
  ggtitle("Resíduos Padronizados - log[pi2/pi3]")+
  xlab("Comprimento") + ggplot2::ylab("Valor") + 
  theme_bw()+
  theme(panel.border = element_blank())+ # para ficar igual o plotly
  # guides(color=guide_legend(title=NULL))+
  scale_color_manual(values=plotly_pallette)+
  theme(plot.title = element_text(hjust=0.5))

plotly::ggplotly(p2)

p3 <- ggplot(as.data.frame(pi1pi2), aes(y=pi1pi2, x=alligator$comp))+
  geom_point(colour=plotly_pallette[3], alpha=0.7) +
  geom_hline(yintercept = 0, size=0.4, alpha=0.5)+
  ggtitle("Resíduos Padronizados - log[pi1/pi2]")+
  xlab("Comprimento") + ggplot2::ylab("Valor") + 
  theme_bw()+
  theme(panel.border = element_blank())+ # para ficar igual o plotly
  # guides(color=guide_legend(title=NULL))+
  scale_color_manual(values=plotly_pallette)+
  theme(plot.title = element_text(hjust=0.5))
  
plotly::ggplotly(p3)

---

---

---