Modelos multinomiais e quase-verossimilhança

Abstract

Tópicos abordados: (1) Realizar uma análise exploratória dos dados e ajuste um modelo de regressão multinomial para estimar a regi’ao de origem de um vinho a partir de suas características sensorias e aromáticas. (2) Realizar uma análise exploratória da correlação entre a variável “FEV(y)” e as demais variáveis preditoras e procure ajuste o modelo mais apropriado para prever o valor médio da variável resposta FEV(y)”.

Leitura do banco de dados

\[\\\]

dados1 = read_csv("wine.csv")

dados2  = read.csv("respiratory.dat", sep="")


# Transformações dados 1 e 2

colnames(dados1)[6]  =  "Region"
dados1$Region = as.factor(dados1$Region)
levels(dados1$Region) <- c("R1","R2","R3")

colnames(dados2)[2]  =  "y"
dados2$Gender = as.factor(dados2$Gender)
dados2$Smoke = as.factor(dados2$Smoke)

Analise Exploratoria - “Wine”

\[\\\]

par(mfrow = c(3,2))
  boxplot(Clarity ~ Region, data=dados1, col  = "yellow")
  boxplot(Aroma ~ Region, data=dados1, col  = "purple")
  boxplot(Body ~ Region, data=dados1, col  = "orange")
  boxplot(Flavor ~ Region, data=dados1, col  = "dark blue")
  boxplot(Quality~ Region, data=dados1, col  = "dark red")
  mx1 = cor(dados1[,c("Clarity","Aroma","Body","Flavor","Oakiness","Quality")],method = "spearman")
    corrplot(mx1,method = 'circle', type="upper", 
          order="AOE", col= RColorBrewer::brewer.pal(n=8, name="RdYlBu"),diag=FALSE)

Análise de dados de Wine

\[\\\] - Observando os Box - Plots, é possível observar que a região denominada como “R3” fabrica vinhos de maior personalidade, ou seja, quanto maior as notas ou score associadas as características sensorias, maior são as observações dos vinhos dessa região.

• Com relação a variável qualidade, nota-se maior concentração de vinhos de R3 em notas ou scores mais altos e vinhoas de R2 concentrados em notas inferiores a R1 e R3. Para as análises futuras, um vinho com baixa nota em qualidade terá maior probabilidade de ser classificado como R2.

• Uma característica do vinho que chama atenção é ‘Clarity’. As notas dos vinhos de R1 estão muito bem disttribuídas entre o primeiro quartil e mediana. Os vinhos de R2 tem variação maior nessa variável, o que provavelmente indica um menor controle de qualidade durante o processo produtivo.

• É possível observar que há alta correlação entre as variáveis Aroma e Flavor, Flavor e Quality, Aroma e Quality, Aroma e Body e Quality e Body. Isso poderá ter implicações na construção do modelo uma vez que há indicação de que as variáveis podem explicar umas as outras (multicolinearidade).

• Pensando no contexto da produção de vinhos, as variáveis Body (corpo), Aroma e Flavor (sabor) de fato tem alta relação entre si. Vinhos encorpados tem como costume apresentar sabores intensos e aromas muito característicos e a depender da região produtora, casta da uva e denominação de origem controlada, essa relação é intensificada.

Ajuste do modelo para Wine

 # Criação das variáveisw dummy
 dados1$YR1 <- as.numeric(dados1$Region == "R1")
 dados1$YR2 <- as.numeric(dados1$Region == "R2")
 dados1$YR3 <- as.numeric(dados1$Region == "R3")

 #(dados1)
 
 ## Coloca os dados no formato da funcao "mlogit"
 mldata       <- mlogit.data(dados1, choice="Region", shape="wide") # banco de dados
 
 # Flavor + Aroma + Body + Quality
 mlogit.model <- mlogit(Region ~ 1|Aroma + Flavor + Quality, 
                        data = mldata, reflevel="R3")

  # mlogit.model <- mlogit(Region ~ 1|Clarity + Aroma + Body + Flavor + Oakiness, 
  #                       data = mldata, reflevel="R1")
  
   
 summary(mlogit.model)

## 
## Call:
## mlogit(formula = Region ~ 1 | Aroma + Flavor + Quality, data = mldata, 
##     reflevel = "R3", method = "nr")
## 
## Frequencies of alternatives:choice
##      R3      R1      R2 
## 0.31579 0.44737 0.23684 
## 
## nr method
## 10 iterations, 0h:0m:0s 
## g'(-H)^-1g = 1.11E-05 
## successive function values within tolerance limits 
## 
## Coefficients :
##                Estimate Std. Error z-value Pr(>|z|)  
## (Intercept):R1 101.5348    59.7317  1.6998  0.08916 .
## (Intercept):R2 114.5962    60.1153  1.9063  0.05661 .
## Aroma:R1        -8.2962     4.9840 -1.6646  0.09600 .
## Aroma:R2        -8.9477     5.0562 -1.7697  0.07678 .
## Flavor:R1        4.9747     3.4412  1.4456  0.14828  
## Flavor:R2        8.2033     3.7566  2.1837  0.02898 *
## Quality:R1      -6.3337     3.8357 -1.6513  0.09868 .
## Quality:R2      -8.5649     3.9399 -2.1739  0.02971 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log-Likelihood: -13.669
## McFadden R^2:  0.66225 
## Likelihood ratio test : chisq = 53.602 (p.value = 8.8736e-10)

\[\\\]

Análise do modelo ajustado

\[\\\] Para o ajuste do modelo foi necessário realizar previamente uma análise qualitativa das variáveis. Ao rodar o modelo com todas as variáveis, foi possível observar que nenhuma apresentava p-valor menor do que o nível de significância, ou seja, não rejeitando a hipótese nula de que o coeficiente é zero (sem efeito)

Dessa forma, para selecionar as variáveis que melhor se ajustariam e explicariam as variações para classificação, foram feitos ajustes com cada variável isolada. Vale ressaltar que, a região de referência escolhida foi ‘R3’, já que proporcionou melhores interpretações do modelo.

Foram selecionadas as variáveis Flavor e Aroma e Quality.

Para fins preditivos, verificou-se que a seleção de apenas duas variáveis explicativas (Aroma e Flavor) não seriam suficientes para prever vinhos da região ‘R2’, portanto selecionamos para modelagem a variável Quality. Na análise descritiva, ela foi a que conseguiu realizar uma melhor distinção entre ‘R1’ e ‘R2’

\[\\\]

Interpretação do modelo

\[\\\] Considerando que temos como referência a região R3, temos a seguinte interpretação

• Para a variável Aroma, quanto maior for o valor da variável, menor será a probabilidade do vinho ser reconhecido como das regiões R1 ou R2

• A variável Flavor e Quality não parecem ter influências significativas em relação a classificação do vinho quando combinada com Aroma. Entretanto, nas previsões podemos observar o comportamento do modelo

\[\\\]

Previsão de probabilidaes para modelo de ‘Wine’

\[\\\] Para demonstração da previsão de probabilidades do modelo, vamos utilizar os valores médios das variáveis preditoras.

saida_p <- predict(mlogit.model, 
                 newdata=data.frame(Flavor=rep(mean(dados1$Flavor),3),
                                    Quality=rep(mean(dados1$Quality),3),
                                    Aroma=rep(mean(dados1$Aroma, 3))))

prev = c("R3","R1","R2")[which.max(saida_p)]

cat("\nPara os casos médios das notas dos vinhos, ele seria classificado como:",prev)

## 
## Para os casos médios das notas dos vinhos, ele seria classificado como: R1

Graficos das probabilidades

\[\\\] Abaixo, temos o gráfico de como as probabilidades de classificação do vinho para as regiões R1, R2 e R3 se alteram quando fixamos as variáveis ‘Quality’e ’Flavor’ em suas respectivas médias e realizamos o incremento dos valores da variável Aroma por 0.1 inciando de 3.3:

 dt <- data.frame(Aroma = seq(3.3, 7.7, by=0.1))
 dt$pR3 <- rep(NA, nrow(dt))
 dt$pR1 <- rep(NA, nrow(dt))
 dt$pR2 <- rep(NA, nrow(dt))
 
  for(i in 1:nrow(dt)){
   saida <- predict(mlogit.model, 
                    newdata=data.frame(Flavor=rep(mean(dados1$Flavor), 3),
                                       Quality=rep(mean(dados1$Quality), 3),
                                       Aroma=rep(dt$Aroma[i], 3)))
   
   dt$pR3[i] <- saida[1]
   dt$pR2[i] <- saida[2]
   dt$pR1[i] <- saida[3]
  }
 
 plot(pR1 ~ Aroma, data=dt, col="blue", type="l", ylim=c(0,1), ylab="probs")
 lines(pR2 ~ Aroma, data=dt, col="red", type="l", ylim=c(0,1))
 lines(pR3 ~ Aroma, data=dt, col="dark green", type="l", ylim=c(0,1))
 legend("topleft", c("prob(R1)","prob(R2)","prob(R3)"),
        fill=c("blue","red","dark green"), cex=0.7, bty="n")

\[\\\]

Nota-se que a medida que aumentamos as notas de Aroma, as probabilidades de classificação de R1 e R2 diminuiem e a de R3 aumenta. Quando um vinho é classificado em 7.7 (nota máxima obtida na base de dados) a probabilidade de ser classificado como R3 é de aproximadamente 1.

\[\\\]

Taxas de acertos - Qualidade de ajuste do modelo (matriz de confusão)

\[\\\]

Para análise da qualidade do ajuste do modelo, fazemos a previsão para cada observação do banco de dados e a comparamos com o que de fato foi classificado:

\[\\\]

 # Estima os eventos mais provaveis (maior probabilidade)
 dados1$fit <- rep(NA, nrow(dados1))


 for(j in 1:nrow(dados1)){
   
   saida <- predict(mlogit.model, 
                    newdata=data.frame(Flavor=rep(dados1$Flavor[j],3),
                                       Quality = rep(dados1$Quality[j],3),
                                       Aroma=rep(dados1$Aroma[j],3)))
   
   dados1$fit[j] = c("R3","R1","R2")[which.max(saida)]

 }

 # Matriz de confusao
 (saida <- with(dados1, table(fit, Region)))

##     Region
## fit  R1 R2 R3
##   R1 15  3  1
##   R2  1  6  0
##   R3  1  0 11

\[\\\]

Taxas de acerto:

• R1: 88%

• R2: 66%

• R3: 91%

• Geral: 84%

É possível verificar que o modelo atingiu uma taxa de acerto considerável, principalmente para as regiões R1 e R3 e no geral. A região R2 é classficada erroneamente como R1 devido a similaridade de suas características, é possível observar na análise descritiva que é muito sutil a diferença das características dos vinhos entre as duas regiões. De posse dessas informações, considero que o modelo tem um ajuste aceitável para classificação de vinhos entre as regiões R1, R2 e R3 com as variáveis aqui expostas.

Análise de dados de ‘Respiratory’

\[\\\]

Analise Exploratória - “Respiratory”

 par(mfrow = c(2,2))
  boxplot(y ~ Age, data=dados2, col  = "grey")
  boxplot(y ~ Smoke, data=dados2, col  = "red")
  boxplot(y ~ Gender, data=dados2, col  = "blue")
  mosaicplot(Gender ~ Smoke, data=dados2, main = "Gender x Smoke")

\[\\\]

g1 = ggplot(data = dados2,aes(x = Ht, y = y))+
  geom_point(aes(col = Gender))
g2 = ggplot(data = dados2,aes(x = Ht, y = y))+
  geom_point(aes(col = Smoke))

  
  grid.arrange(g1,g2,ncol = 2, nrow = 1)

\[\\\]

• Podemos observar com essa rápida análise exploratória de dados que o índice da capacidade respiratória aumenta conforme a idade. Destaque para a idade de 17 anos com variações relevantes entre o primeiro e terceiro quartil (dados mais dispersos).

• É possível observar também que o índice da capacidade respiratória aumenta conforme conforme o inicador ‘Ht’. Ht tem alta correlação com a idade o que pode indicar que esteja associado com determinada capacidade atlética dos indivíduos. É importante observar que indivíduos com ‘Gender’ igual 1 possuem maior índice da capacidade respiratória (provavelmente do gênero masculino) e indivíduos com marcação de que fumam também. Como temos uma limitação de idade gritante no banco de dados, é possível afirmar que não há fumantes de longa data portanto o efeito do cigarro na capacidade respiratória tem pouca ou nenhuma influência. É importante salientar que temos uma maioria de não fumantes no banco de dados.

Ajuste de um modelo linear

\[\\\]

  # Ajuste do modelo linear

  modelo21 <- lm(y ~ Age + Ht + Gender + Smoke, data=dados2)

 
 summary(modelo21)

## 
## Call:
## lm(formula = y ~ Age + Ht + Gender + Smoke, data = dados2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.37656 -0.25033  0.00894  0.25588  1.92047 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.456974   0.222839 -20.001  < 2e-16 ***
## Age          0.065509   0.009489   6.904 1.21e-11 ***
## Ht           0.104199   0.004758  21.901  < 2e-16 ***
## Gender1      0.157103   0.033207   4.731 2.74e-06 ***
## Smoke1      -0.087246   0.059254  -1.472    0.141    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4122 on 649 degrees of freedom
## Multiple R-squared:  0.7754, Adjusted R-squared:  0.774 
## F-statistic:   560 on 4 and 649 DF,  p-value: < 2.2e-16

 # Gráfico dos resíduos 

 plot(modelo21$fitted.values,modelo21$residuals,
      main = "Residuals vs Fitted",
      xlab = "Fitted",
      ylab = "Residuals")
 abline(h = 0, col  ="red",lty= 2)

Análise do ajuste linear

\[\\\]

• O modelo linear ajustado apresentaa um coeficiente de determinação ajustado igual a 77.40%, indicando que o modelo explica uma parcela considerável da variabilidade da variável resposta. O modelo apresentou p-valores abaixo de 0.05 confirmando a significância do mesmo.

• Quando avaliamos o gráfico de resisíduos vs ajuste, é possível notar pela forma da distribuição de dados, que não há homocedasticidade. Portanto o modelo linear não é adequado para esse conjunto de dados.

• Observando o comportamento dos resíduos, podemos observar que a estrtutura de variância do modelo tem indícios de crescimento proporcinal a média.

Ajustes não lineares

\[\\\] - Conforme exposto acima, buscamos um ajuste que considere a esrtutura de variância quadrática, provavelmente que varie com a esperança e que atenda as características dos dados.

• O modelo de Poisson seria um possível candidato para ajuste, dado a estrutura de dados. Entretanto o modelo de Poisson apresenta algumas restrições, uma vez que é adequado para variáveis discretas e ssua estrutura de variância é proporcional a esperança. Dessa forma, podemos usar uma quasi Poisson e averiguar os resulltados do modelo.

• Outro modelo candidato é o modelo Gama, que pode atender a estrutura de variância que buscamos \[\\\]

MODELO GAMA \[\\\]

  # Modelo Gama  
  modelo22 <- glm(y ~ Age + Ht + Gender + Smoke,
                family=Gamma(link = "log"), data=dados2)

 summary(modelo22)

## 
## Call:
## glm(formula = y ~ Age + Ht + Gender + Smoke, family = Gamma(link = "log"), 
##     data = dados2)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.58134  -0.09456   0.00209   0.08461   0.42455  
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.940198   0.076402 -25.395  < 2e-16 ***
## Age          0.022816   0.003253   7.013 5.86e-12 ***
## Ht           0.042985   0.001631  26.351  < 2e-16 ***
## Gender1      0.029409   0.011385   2.583   0.0100 *  
## Smoke1      -0.040853   0.020315  -2.011   0.0447 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.01997449)
## 
##     Null deviance: 70.791  on 653  degrees of freedom
## Residual deviance: 13.414  on 649  degrees of freedom
## AIC: 525.61
## 
## Number of Fisher Scoring iterations: 4

 p22 = 1 - pchisq(modelo22$deviance, modelo22$df.residual)
 cat("\nDeviance do modelo é igual a: ",p22,"portanto apresenta subdispersão")

## 
## Deviance do modelo é igual a:  1 portanto apresenta subdispersão

\[\\\] Modelo Gama apresenta subdispersão, portanto não é adequado para representação dos dados.

\[\\\]

MODELO QUASI POISSON

  # Modelo quasi Poisson 
  modelo23 <- glm(y ~ Age + Ht + Gender + Smoke,
                family=quasipoisson, data=dados2)

  # summary(modelo23)
   modelo25 <- glm(y ~ Age + Ht + Gender + Smoke,
               family=quasi(link = "log", 
                            variance = "mu"),
               data=dados2,
               weights = rep(1/0.0530696, nrow(dados2)))
   
   summary(modelo25)

## 
## Call:
## glm(formula = y ~ Age + Ht + Gender + Smoke, family = quasi(link = "log", 
##     variance = "mu"), data = dados2, weights = rep(1/0.0530696, 
##     nrow(dados2)))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.7451  -0.6434   0.0061   0.5990   3.5313  
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.914444   0.078726 -24.318  < 2e-16 ***
## Age          0.023453   0.003040   7.715 4.57e-14 ***
## Ht           0.042427   0.001612  26.312  < 2e-16 ***
## Gender1      0.033698   0.011727   2.873  0.00419 ** 
## Smoke1      -0.040775   0.018415  -2.214  0.02716 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasi family taken to be 1)
## 
##     Null deviance: 3427.79  on 653  degrees of freedom
## Residual deviance:  656.14  on 649  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4

  p25 = (1 - pchisq(modelo25$deviance, modelo25$df.residual))
  cat("\nDeviance do modelo é igual a: ",p25)

## 
## Deviance do modelo é igual a:  0.4145323

Análise do ajuste

\[\\\]

• Foi possível verificar que o ajuste por Gama do modelo não se mostrou adequado devido à subdispersão

• Portanto, recorreu-se a estrutura do modelo de Poisson para estimar \(\phi\) e ajusar um modelo que tivesse estrutura expenoncial, variância proporcional a esperança onde a variável resposta pudesse ser contínua.

• Após esses ajustes foi possível estimar um modelo para prever a o valor médio do índice de capacidade respiratória.