GRUPO: Rafaela Lorena; Eduardo Moreira; Ana Thiciane e Adriana Jardim
30 de outubro de 2019
Pré-requisitos
Nós vamos partir do pressuposto que você (leitor) tem pelo menos o conhecimento básico sobre os seguintes tópicos:
- Noções de regressão linear
- O que são variável dependente e independente
- O que é inferência frequentista e inferência multimodelo
- Noções do programa estatístico R e ambiente RStudio
Observação: Caso não tenha os conhecimentos básicos supracitado ou sinta dificuldade de entender algum tópico desse tutorial, nós sugerimos que você dê uma olhadinha nos links abaixos antes de prosseguir:
1. [Noções sobre regressão linear]: https://pt.khanacademy.org/search?page_search_query=regress%C3%A3o%20linear
2. [Variável dependente e independente]: http://ecologyandevolution.org/statsdocs/online-stats-manual-chapter3.html
3. [Inferência frequentista e inferência multimodelo. Livro: Princípios de Estatística em Ecologia. Gotelli & Ellison, 2011. Veja algumas páginas disponíveis deste livro em]: https://books.google.com.br/books?id=BFJCDQAAQBAJ&printsec=frontcover&dq=princ%C3%ADpios+de+estat%C3%ADstica+em+ecologia+gotelli&hl=pt-BR&sa=X&ved=0ahUKEwjI_M_OrrflAhWqK7kGHTOTDTMQ6AEIKTAA#v=onepage&q=princ%C3%ADpios%20de%20estat%C3%ADstica%20em%20ecologia%20gotelli&f=false
4. [Caso necessite fazer o download do programa R]: https://nbcgib.uesc.br/mirrors/cran/
5. [Caso necessite fazer o download do RStudio]: https://rstudio.com/
6. [Introdução ao uso do programa R]: https://cran.r-project.org/doc/contrib/Landeiro-Introducao.pdf
7. [Tutorial R e do RStudio]: https://edisciplinas.usp.br/pluginfile.php/2996937/mod_resource/content/1/Tutorial.pdf
Conteúdo
Neste tutorial sobre regressão logística, nós vamos inicialmente definir uma regressão logística, mostrar sua equação e as principais diferenças entre a regressão logística e a regressão linear. Em seguida, por meio de um estudo de caso ecológico será construído uma regressão logística binária utilizando o RStudio, onde nós vamos explorar com mais detalhes os seguintes tópicos:
A Função de ligação (logit)
Modelo estatístico
A Estimação dos parâmetros do modelo logístico por máxima verossimilhança
O teste de significância e intervalo de confiança
A qualidade de ajuste
A Inferência frequentista usando Análise de desvio (ANDEVA)
A Inferência multimodelo usando AIC
Regressão logística
Pode-se considerar a regressão logística como sendo um caso especial de regressão linear onde a variável dependente (eixo Y do gráfico) é categórica, ou seja, possue categorias/níveis
Na regressão logística, a varíavel categórica pode ser de natureza:
dicotômica e binária: quando a categoria estabelecida tem 2 níveis. Por exemplo: podemos classificar a visitação das flores de uma planta por abelhas em 2 níveis: planta visitadas por abelhas e planta não visitadas por abelhas. Quando as categorias são representadas por números (1: evento de interesse e 0: ausência de efeito de interesse), chamamos de variável binária
ordinal: quando existe uma ordenação entre 3 ou mais categorias/níveis (politômicas). Por exemplo: mês de observação dos dados coletados que podem ser janeiro, fevereiro, março…dezembro
nominal: quando não existe uma ordenação entre 3 ou mais categorias/níveis (politômicas). Por exemplo: a cor dos olhos que tem-se azul, preto, verde, etc.
Chamamos de regressão logística ordinal quando a variável dependente é categórica de natureza ordinal politômica
Chamamos de regressão logística nominal quando a variável dependente é categórica de natureza nominal politômica
Chamamos de regressão logística binária quando a variável dependente é categórica de natureza dicotômica ou binária
Já, a variável independente (eixo X do gráfico), na regressão logística, pode ser tanto categórica e/ou quantitativa. Se só tiver 1 variável independente é chamado de regressão logística simples, agora se tiver mais de 1 variável independente chamamos de regressão logística múltipla
Observação: Uma variável quantitativa pode ser categorizada e assim ser passível de realização de uma regressão logística. Por exemplo: temperatura é uma variável quantitativa que pode ser categorizada em níveis como: temperatura alta, média e baixa, sendo os níveis justificados de acordo com seu grupo de interesse
Veja a figura 1 (abaixo) que resume de forma esquemática os diferentes tipos de vaiáveis empregadas na regressão logística e exemplos gráficos para demonstrar uma regressão logística simples e múltipla utilizando uma variável dependente binária (regressão logística binária)
Equação de regressão logística
A regressão logística consiste em uma generalização de modelos lineares para casos onde os modelos lineares tradicionais não seriam adequados (variáveis dependentes categóricas). Isso é possível graças a uma série de estratégias matemáticas que permitem justamente adaptar de uma maneira simples o arcabouço dos modelos lineares a predição de variáveis categóricas.
Enquanto que a regressão linear descreve como o incremento em uma variável preditora (ou variável independente) afeta o incremento em uma variável resposta (ou variável dependente), a regressão logística descreve como o incremento em uma variável preditora afeta a probabilidade de um evento acontecer. Desta maneira, a variável resposta em uma regressão logística é a probabilidade de um determinado evento (p(x)), podendo ser representada pela expressão:
\[ p(x) = Pr(Y = 1|X = x) \]
que nada mais é do que uma probabilidade condicional do evento Y=1 (evento de interesse) dado o valor de X.
Como probabilidades (p) são razões entre o número de sucessos/evento de interesse (r) e o número total de eventos (n), ou seja p = r/n, elas se comportam de uma maneira muito distinta de uma relação linear:
Então, como podemos usar uma regressão linear para prever a probabilidade condicional de um evento acontecer?
A resposta mais óbvia seria deixar a probabilidade do evento ser uma função linear de x. O problema nessa abordagem é que a probabilidade varia entre 0 e 1 (0≤p(x)≥1), enquanto que modelos lineares são irrestritos nas duas direções (-∞, ∞). Isso significa que simplesmente aplicar um modelo linear para prever probabilidades acarretaria na predição de valores absurdos como probabilidades negativas ou maiores do que 100%;
Uma segunda alternativa seria permitir que o logaritmo da probabilidade seja uma função linear de x, para que a alteração de uma variável de entrada multiplique a probabilidade por um valor fixo. O problema é que os logaritmos são ilimitados em apenas uma direção e as funções lineares não. O que implica em valores sem sentido prático.
Para resolver os dois problemas acima os estatísticos utilizaram o logaritmo das chances do evento acontecer (log(odds)) também conhecido como função logit. Esta função possue as vantagens de ser irrestrita em ambas as direções e possuir uma relação direta com as probabilidades, permitindo seu cálculo através da aplicação de uma fórmula de conversão simples.
A chance (odds) de um evento acontecer pode ser calculada através da razão entre a probabilidade desse evento acontecer (no nosso caso p(x)) e a probabilidade do evento não acontecer (1-p(x)). Sendo assim a função logit para nossa probabilidade condicionals fica:
\[ logit~=~log\frac{p(x)}{1~-~p(x)}~=~β_0~+~β_1*X~~~~~~~~(equação~1) \]
Resolvendo essa equação chegamos a função que converte os valores obtidos na escala logit (logaritmo de odds) nos valores da probabilidade do evento Y=1 acontecerem dado um determinado valor de x, ou seja:
\[ p(x)~=~\frac{e^{β_0~+~β_1*X}}{1~+~e^{β_0~+~β_1*X}}~~~~~~~~~~~~(equação~2) \] sendo:
β0 e β1 os parâmetros do modelo linear
- enquanto que a função logit (equação 1) permite a aplicação de um modelo linear para estimar o log(odds) de um evento dada uma determinada variável preditora, a equação 2 permite a a conversão dos valores em escala logit, que varia de -∞ a +∞, para valores de probabilidade que variam entre 0 e 1.
Observação: Para obtermos as equações de regressão logística estimada (para amostra), ou seja, Ŷ = estimação de p só é substituir na equação 2 o β0 por b0 e β1 por b1…
Então, perceba que a curva de regressão logística é descrita em função de 2 parâmetros:
β0 ou b0 que significa a probabilidade de um resultado de Y quando X = 0 e
β1 ou b1: que determina se a curva em forma de S é mais acentuada ou menos acentuada. Quando β1 ou b1 é positivo tem uma curva em forma de S e quando β1 ou b1 é negativo tem-se uma curva em forma de S invertida
Agora, vamos relembrar a equação de regressão linear simples (equação 2) e regressão linear múltipla (equação 3), sendo Y uma variável dependente contínua, f(x) é a média ou valor esperado de Y para determinado valor de X. X é a variável independente e m é o número de variáveis independentes
\[
f(x)~=~β_0~+~β_1*X~~~~~~~~~~~~(equação~3)
\]
\[
f(x)~=~β_0~+~β_1*X~+~β_2*X_2~+~β_3*X_3~+...+~β_m*X_m~~~~~~~~~~~~(equação~4)
\]
Voltando a função logit agora fica fácil ver como a mesma lógica se aplica… para migrar de uma regressão logística simples para uma regressão logística múltipla basta acrescentar os parâmetros ncessários nos expoentes da função logit:
- A equação da regressão logística múltipla:
\[
p(x)~=~\frac{e^{β_0~+~β_1*X_1~+~β_2*X_2~+~β_3*X_3~+...+~β_m*X_m}}{1~+~e^{β_0~+~β_1*X_1~+~β_2*X_2~+~β_3*X_3~+...+~β_m*X_m}}~~~~~~~~~~~~(equação~5)
\]
Ampliando esse princípio, fica claro que uma vez que função logit é aplicada, praticamente qualquer caso de uma regressão linear clássica pode ser aplicado a uma regressão logística, incluindo combinações de efeitos de variáveis preditoras quantitativas, categóricas e interações.
Diferenças entre a regressão logística e regressão linear
Abaixo um quadro que resume as principais diferenças entre a regressão logística e regressão linear
Como realizar e interpretar uma regressão logística binária múltipla usando o R
Pelo fato da regressão logística binária ser o tipo mais simples e um dos mais comumente empregado em Ecologia, então nós vamos implementar este tipo de regressão no R e interpretar as saídas padrão desse tipo de análise usando um exemplo simulado de dados ecológicos.
Caso de estudo ecológico hipotético sobre a polinização da macieira
Primeiro vamos apresentar rapidamente o contexto do nosso estudo de caso: A polinização é a transferência dos grãos de pólen das anteras de uma flor para o receptor feminino (estigma). Esta transferência é realizada por exemplo por meio do vento ou por animais. A maioria das plantas com flores depende desse processo em algum nível para a formação de frutos e sementes. Por esse motivo, cultivos agrícolas como o da macieira que é totalmente dependentes dos polinizadores, principalmente das abelhas para a produção de seus frutos. Nesse sentido, a polinização é um fator limitante da produção de maçãs, assim como outros fatores mais amplamente conhecidos como por exemplo a disponibilidade de água. Diante deste contexto, pretende-se saber se as variáveis independentes grãos de pólen no estígma (variável contínua) e irrigação (variável binária) afetam a probablidade das flores de maçã estudadas formarem frutos ou não (variável dependente binária).
As nossas hipóteses de estudo serão:
Quanto maior for a quantidade de pólen no estígma da flor de maçã maior será a probabilidade da mesma formar um fruto;
Flores em plantas irrigadas terão maior probabilidade de formar frutos do que flores em plantas não irrigadas;
A irrigação aumenta o efeito do número de grãos de pólen sobre a probabilidade das flores de maçã gerarem frutos.
Geração dos dados referente ao estudo ecológico hipotético apresentado
Para exemplificar a aplicação da regressão logística nós vamos usar dados simulados que refletem nossa hipótese de estudo. Para simular esses dados usamos o código abaixo no R:
# Caso seja necessário instale os pacotes abaixo que serão utilizados para a realização das análises
## install.packages('corrplot')
## install.packages('RColorBrewer')
## install.packages('MuMIn')
# Carregando os pacotes que serão utilizados para a realização das análises
library(corrplot)## corrplot 0.84 loadedlibrary(RColorBrewer)
library(MuMIn)
# Gerando os vetores que vão receber os dados simulados
# O vetor de pólen contém os dados de número de grãos de pólen no estígma das flores de macieira e consiste em números inteiros gerados a partir de uma distribuição normal com média 6 e desvio padrão 2.5
pollen <- as.integer(sqrt(rnorm(n=200, mean=6, sd=2.5)^2))
# O vetor irrigation informa se as flores estão em plantas irrigadas ou não. Esses dois níveis estão representados pelos valores 1 e 0 respectivamente
irrigation <- sample(c(0,1), 200, replace = TRUE)
# O vetor fruit guarda os valores da nossa variável dependente, que consiste na informação de se as flores se desenvolveram em frutos ou não, representando as duas situações pelos valores 1 e 0 respectivamente. Aqui estamos criando o vetor vazio do mesmo tamanho do vetor pollen criado previamente. No próximo passo vamos preencher esse vetor de acordo com critérios que simularão o efeito das duas variáveis explicativas criadas anteriormente (i.e. pollen e irrigation)
fruit <- vector(length=length(pollen))
# O código a seguir atribuirá valores 0 e 1 ao vetor fruit de acordo com o seguinte critério:
# Em cada caso estudado (flores), se a flor recebeu quantidade suficiente de grãos de pólen (pollen > 5) e foi irrigado (irrigation == 1), a flor formou um fruto, se não foi irrigado (irrigation == 0) a flor tem 50% de chances de gerar um fruto ou não. Dessa maneira estabelecemos que na ausência de outro fator limitante as flores que receberam uma quantidade suficiente de pólen gerarão frutos e que em casos onde as plantas estão sujeitas a contingência de outro fator limitante (água), existe uma probabilidade de que os frutos não se formem. Nos casos onde as flores não receberam quantidade suficiente de pólen para garantir que os frutos se formarão na ausência de outro fator limitante (irrigation == 1) a flor tem uma probabilidade de 50% de formar um fruto ou não. Nos casos onde não há pólen suficiente e não há irrigação as flores não formam frutos
for (i in 1:length(pollen)) {
if (pollen[i] > 5) {
if (irrigation[i] == 1) {
fruit[i] <- 1
} else {
fruit[i] <- sample(c(0, 1), 1)
}
}
else {
if (irrigation[i] == 1) {
fruit[i] <- sample(c(0, 1), 1)
} else {
fruit[i] <- 0
}
}
}
# Cada vez que o código acima for executado os valores dos dados simulados serão diferentes...
# Agora agrupamos os três vetores no dataframe data e escrevemos nosso dataframe em um arquivo .CSV no nosso diretório de trabalho
# Vendo o diretório de trabalho atual
## getwd()
# Caso queira trocar o diretório de trabalho para gravar o arquivo em .CSV
## setwd()
data <- data.frame("FRUIT"=fruit, "POLLEN"=pollen, "IRRIGATION"=irrigation)
# write.csv(data, file = "./pollination_test_data.csv")
# Como queremos explorar com você as propriedades da regressão logística disponibilizamos a nossa tabela de dados usadada para os códigos abaixo no arquivo .CSV através do link: XXXXX
# Para começar, vamos explorar as propriedades dos nossos dados. O primeiro passo é carregar os dados simulados para rodar o restante do script
# Utilize o arquivo .CSV que nós compartilhamos com você. Veja o link: XXXX
ourdata <- read.csv("./pollination_test_data.csv", sep=",") #Troque o diretório
# Nosso dados simulados
head(ourdata)## X FRUIT POLLEN IRRIGATION
## 1 1 0 2 0
## 2 2 1 9 0
## 3 3 0 4 0
## 4 4 0 3 0
## 5 5 1 6 1
## 6 6 1 8 1# Resumo dos nossos dados simulados
table(ourdata$POLLEN)##
## 0 1 2 3 4 5 6 7 8 9 10 11
## 5 4 19 14 28 35 28 22 18 16 9 2table(ourdata$IRRIGATION)##
## 0 1
## 92 108table(ourdata$FRUIT)##
## 0 1
## 100 100# Histograma dos dados de pólen: Variável contínua
histograma<-hist(ourdata$POLLEN, col="lightblue", border="black", xlab="POLLEN", ylab="FREQUÊNCIA", main="Histograma dos Dados de POLLEN", breaks=10, xlim=c(0, 12)) #breaks: muda o número de barras do gráfico
abline(v=mean(ourdata$POLLEN), col="red", lwd=2) #Adicionando uma linha vertical vermelha representando a média dos dados
# Acrescentando a curva normal ao gráfico
xfit <- seq(min(ourdata$POLLEN), max(ourdata$POLLEN))
yfit <- dnorm(xfit, mean=mean(ourdata$POLLEN), sd=sd(ourdata$POLLEN))
yfit <- yfit*diff(histograma$mids[1:2]) * length(ourdata$POLLEN)
lines(xfit, yfit, col="black", lwd=2)
# Adicionando a legenda
legend(x="topright", #Posição da legenda
c("Média", "Curva normal"), #Nome da legenda
col=c("red", "black"), #Cor
lty=c(1, 1), #Estilo da linha
lwd=c(2, 2)) #Grossura das linhas
# Frequência dos dados de irrigação: Variável binária
barplot(data.frame(table(ourdata$IRRIGATION))[,2], names.arg = c(0, 1), col="lightblue", border="black", xlab="IRRIGATION",
ylab="FREQUÊNCIA", main ="Frequência de IRRIGATION")
# Frequência dos dados de fruit: Variável binária
barplot(data.frame(table(ourdata$FRUIT))[,2], names.arg = c(0, 1), col="lightpink", border = "black", xlab="FRUIT",
ylab="FREQUÊNCIA", main="Frequência de FRUIT")
Podemos começar nossas análises explorando a associação entre as variáveis através da representação gráfica da matriz de correlações de Pearson
M <-cor(ourdata[,2:4])
corrplot(M, title="Correlação entre as Variáveis: FRUIT, POLLEN e IRRIGATION", mar=c(0, 0, 2, 0), type="full", order="hclust", col=brewer.pal(n=8, name="RdYlBu"))
Podemos notar no gráfico (acima) que as duas variáveis preditoras não estão correlacionadas, o que já nos dá uma segurança razoável de que multicolinearidade não será um problema na parametrização do nosso modelo. Já a formação de frutos parece estar correlacionada com ambas as variáveis explicativas como esperávamos.
Realizando a regressão logística binária múltipla referente ao estudo ecológico hipotético apresentado
- Usando a escala logit (log(odds) ou função logit)
# Caso seja necessário instale os pacotes abaixo que serão utilizados para a realização das análises
## install.packages('visreg')
## install.packages('MASS')
## install.packages('ggplot2')
## install.packages('GGally')
# Carregando os pacotes que serão utilizados para a realização das análises
library(visreg)
library(MASS)
library(ggplot2)
library(GGally)## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2# Agora ajustamos nosso modelo logístico global
logit_fruit <- glm(factor(FRUIT) ~ POLLEN * factor(IRRIGATION), data=ourdata, family=binomial(link='logit'))
summary(logit_fruit)$coef## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.59982645 1.1905811 -4.7034396 2.558147e-06
## POLLEN 0.66440756 0.1635975 4.0612333 4.881416e-05
## factor(IRRIGATION)1 4.00700339 1.3405466 2.9890818 2.798173e-03
## POLLEN:factor(IRRIGATION)1 -0.09870833 0.2103757 -0.4692003 6.389265e-01Como podemos observar na tabela summary dos resultados da nossa regressão logística binária múltipla, no item Coefficients estão os valores estimados para os parâmetros na coluna Estimate em escala logit, seguido do Std. Error que junto com os valores de Estimate são usados para calcular os z value de cada parâmetro, que nada mais é do que uma estatística padronizada que pode ser usada no teste da hipótese de que os Estimates são significativamente diferentes de zero e na última coluna temos os valores de p para cada valor de z.
Estimação dos parâmetros do modelo
- Os parâmetros (b0, b1 e b2) são estimados a partir de um processo iterativo onde a cada passo novos valores são atribuídos a esses parâmetros, gerando um modelo alternativo a partir do qual são calculados os valores de probabilidade de acordo com a equação 5 para os dois eventos P(Y=1|X1|X2) e P(Y=0|X1|X2). Em posse das probabilidades dos dados de acordo com O modelo se calcula a soma dos logaritmos dessas probabilidades, obtendo-se assim o logaritmo da verossimilhança do modelo. Esse processo é repetido um certo número de vezes até que se possa estimar qual conjunto de valores dos parâmetros possui a máxima verossimilhança. Para otimizar esse processo a função
glm()usa como método de otimização padrão o algoritmo Iterative Weighted Least Squares (IWLS) é uma implementação do Fisher Scoring Algorithm. No summary do modelo também é informado o número de iterações necessárias para esse algorítmo encontrar o modelo com máxima verossimilhança:
summary(logit_fruit)$iter## [1] 5No nosso caso, foram apenas 5.
Observação: [Caso deseje se aprofundar sobre o assunto verossimilhança e máxima verossimilança, veja esse link]: http://cmq.esalq.usp.br/BIE5781/lib/exe/fetch.php?media=leituras:verossim.pdf
Agora nos concentramos no nosso modelo estatístico estimado: vendo a coluna do Estimate no summary e substituindo na equação 5
\[ Ŷ_i \to Binomial~(p_i,~N´)~independente \]
\[ Ŷ_i~=~estimação~de~p_i~=~\frac{exp(-5.59983~+~0.66441*POLLEN_i~_+~4.00700*IRRIGATION_i~)}{1~+~exp(-5.59983~+~0.66441*POLLEN_i~+~4.00700*IRRIGATION_i)} \]
\[
i~=~ 1,2,…,200
\]
Também podemos explorar os detalhes do teste de significância: vendo as colunas z value (valor de z) e Pr(>|z|) (valor de significância) do summary
Hipótese nula (H0) = b1 = b2
Hipótese alternativa (Ha) = um ou todos os parâmetros são ≠ 0
z value = Estimate / Std. Error:
Intercepto: -5.59983 / 1.19058 = -4.703
POLLEN: 0.66441 / 0.16360 = 4.061
factor(IRRIGATION)1: 4.00700 / 1.34055 = 2.989
POLLEN:factor(IRRIGATION)1: -0.09871 / 0.21038 = -0.469
Pr(>|z|): valor de significância para a variável POLLEN
H0 = b1 dado b2 = 0
Ha = b1 dado b2 ≠ 0
Resultado: 4.88e-05 rejeitou a H0 mostrando que POLLEN tem um efeito positivo sobre o fruit
Pr(>|z|): valor de significância para a variável factor(IRRIGATION)1:
H0 = b2 dado b1 = 0
Ha = b2 dado b1 ≠ 0
Resultado: 0.0028 *** rejeitou a H0 mostrando que IRRIGATION tem um efeito positivo sobre o fruit
Então, podemos concluir com base nos nossos resultados que número dos grãos de pólen tem um efeito positivo sobre a probabilidade da formação de frutos, juntamente com a irrigação e que a interação entre essas variáveis explicativas (ou variáveis independentes) não possui um efeito significativo sobre essas relações.
Agora podemos visualizar esse modelo de maneira gráfica usando a função visreg():
# Visualizamos os valores preditos na escala logit (log(odds))
visreg(logit_fruit, "POLLEN", "IRRIGATION", gg=T, overlay=T, xlab="POLLEN", ylab="F(POLLEN)") +
ggtitle("Valores Preditos na Escala log(Odds) ou logit") + theme(plot.title = element_text(hjust = 0.5)) +
theme(plot.title = element_text(face = "bold"))
A função visreg() plota os valores preditos de acordo com o modelo. Como podemos ver o eixo Y está na escala logit, e por isso os valores preditos possuem uma relação linear com a variável preditora POLLEN. Além disso podemos ver os efeitos positivos do número de grãos de pólen (inclinação das retas) e da irrigação (diferença entre as duas retas). Também fica evidente que a interação entre as variáveis foi muito pequeno, uma vez que o coeficiente de inclinação das duas retas (com irrigação em azul e sem irrigação em vermelho) é praticamente idêntico. As linhas coloridas representam os valores preditos e as regiões em torno dessas linhas os intervalos de confiança.
Também podemos visualizar nossos coeficientes em escala de probabilidade aplicando a fórmula de conversão descrita na equação número 5. No R basta usar a função exp sobre a tabela de coeficientes do summary do modelo:
exp(summary(logit_fruit)$coefficients)## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.003698506 3.288992 0.009064046 1.000003
## POLLEN 1.943338872 1.177740 58.045855683 1.000049
## factor(IRRIGATION)1 54.981864508 3.821132 19.867430890 1.002802
## POLLEN:factor(IRRIGATION)1 0.906006924 1.234142 0.625502296 1.894446Como podemos ver, agora os valores dos parâmetros estimados foram convertidos para a escala de probabilidades. Contudo, as demais medidas também foram afetadas e perderam o sentido.
Uma alternativa é usar a função summ()do pacote jtools que automatiza essa conversão de uma maneira muito prática, bastando agregar o argumento exp = T:
# Caso seja necessário instale os pacotes abaixo que serão utilizados para a realização das análises
## install.packages('jtools')
# Carregando o pacote
library(jtools)
summ(logit_fruit, exp = T)## MODEL INFO:
## Observations: 200
## Dependent Variable: factor(FRUIT)
## Type: Generalized linear model
## Family: binomial
## Link function: logit
##
## MODEL FIT:
## <U+03C7>²(3) = 114.33, p = 0.00
## Pseudo-R² (Cragg-Uhler) = 0.58
## Pseudo-R² (McFadden) = 0.41
## AIC = 170.93, BIC = 184.12
##
## Standard errors: MLE
## ----------------------------------------------------------------------------
## exp(Est.) 2.5% 97.5% z val. p
## -------------------------------- ----------- ------ -------- -------- ------
## (Intercept) 0.00 0.00 0.04 -4.70 0.00
## POLLEN 1.94 1.41 2.68 4.06 0.00
## factor(IRRIGATION)1 54.98 3.97 760.84 2.99 0.00
## POLLEN:factor(IRRIGATION)1 0.91 0.60 1.37 -0.47 0.64
## ----------------------------------------------------------------------------Além de fornecer algumas medidas adicionais de qualidade de ajuste como Pseudo-R² e executar um teste chi quadrado, a função summ() retorna os valores estimados de cada parâmetro na escala de probabilidades, associado aos respectivos intervalos de confiança, valores de z e p.
Também podemos visualizar o modelo em escala de probabilidade com a funão visreg(), bastando agregar o argumento scale = “response”, que a própria função aplicará a equação 5 a todos os valores preditos.
# Visualizamos os valores preditos na escala de probabilidades
visreg(logit_fruit, "POLLEN", "IRRIGATION", gg=T, overlay=T, xlab="POLLEN", ylab="F(POLLEN)", scale="response") +
ggtitle("Valores Preditos na Escala de Probabilidades") + theme(plot.title = element_text(hjust = 0.5)) +
theme(plot.title = element_text(face = "bold"))
Agora podemos ver a forma sigmóide (curva em forma de S) da nossa função de probabilidades condicionais, com os valores variando entre 0 e 1 como esperado. Também podemos ver a diferença entre as flores com irrigação em relação as flores sem irrigação e o efeito positivo do pólen.
Qualidade do ajuste do modelo
- Agora vamos olhar para o summary completo do modelo:
summary(logit_fruit)##
## Call:
## glm(formula = factor(FRUIT) ~ POLLEN * factor(IRRIGATION), family = binomial(link = "logit"),
## data = ourdata)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.93655 -0.60277 0.02671 0.71434 1.89475
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.59983 1.19058 -4.703 2.56e-06 ***
## POLLEN 0.66441 0.16360 4.061 4.88e-05 ***
## factor(IRRIGATION)1 4.00700 1.34055 2.989 0.0028 **
## POLLEN:factor(IRRIGATION)1 -0.09871 0.21038 -0.469 0.6389
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 277.26 on 199 degrees of freedom
## Residual deviance: 162.93 on 196 degrees of freedom
## AIC: 170.93
##
## Number of Fisher Scoring iterations: 5Como medidas de bondade de ajuste podemos usar a relação entre as medidas de deviance (chamado em português de desvio) e AIC (critério de informação de Akaike), comparadas com a medida de deviance do modelo nulo (sem as variáveis explicativas), ou seja, Null deviance (desvio correspondente ao modelo nulo).
\[ Null~Deviance~=~-2(logLik(modelo~nulo)~-~logLik(modelo~saturado~ou~completo)) \] \[ ou \]
\[
Null~Deviance~=~-2(logLik(modelo~nulo)~-~0)
\]
\[ Residual~deviance~=~-2(logLik(modelo~proposto)~-~logLik(modelo~saturado~ou~completo)) \]
\[ ou \]
\[ Residual~deviance~=~-2~(logLik(modelo~proposto)~-~0) \]
Note que logLik é o logaritmo de máxima verossimilhança;
Graus de liberdade (degrees of freedom = tamanho da amostra menos o número de parâmetros estatísticos);
O modelo saturado é um modelo com um parâmetro para cada observação, ou seja, é o modelo com o menor valor de deviance possível;
O modelo proposto é um modelo que assume que os dados podem ser explicados com m parâmetros mais um intercepto.
Conclusão olhando o summary do modelo gerado: Quando o modelo inclui as variáveis POLLEN e IRRIGATION, o desvio residual (residual deviance) foi menor 162.93 do que o desvio correspondente ao modelo nulo (Null deviance) que foi de 277.26. Assim, um valor mais baixo do desvio residual indica que o modelo ficou melhor quando inclui as duas variáveis independentes. Já o AIC também é baseado no desvio (deviance) e é útil quando queremos comparar modelos e seu valor por si só não é interpretável.
Contribuição de cada variável independente da regressão logística binária múltipla gerada
Para interpretarmos a contribuição de cada variável independente vamos utilizar o modelo com a transformação logit e vamos utilizar:
Inferência frequentista usando Análise de desvio (ANDEVA) e
Inferência multimodelo com AIC
Podemos avaliar a contribuição de cada variável através da análise comparativa dos desvios do modelo global em relação aos os modelos aninhados com menos parâmetros. Essa técnica é denominada de Analysis of Deviance ou ANDEVA
Podemos fazer isso aplicando a função anova() ao modelo global especificando que queremos o teste chi quadrado:
# Modelo na escala logit
anova(logit_fruit, test="Chisq")## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: factor(FRUIT)
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 199 277.26
## POLLEN 1 35.304 198 241.96 2.82e-09 ***
## factor(IRRIGATION) 1 78.805 197 163.15 < 2.2e-16 ***
## POLLEN:factor(IRRIGATION) 1 0.223 196 162.93 0.6369
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Como podemos ver esse teste compara os valores de desvio dos modelos aninhados acrescentando um parâmetro por vez de maneira sequencial partindo do modelo nulo. Em deviance temos o valor de redução no residual deviance em relação ao modelo anterior e na coluna Pr(>Chi) temos o valor de p (significância) de acordo com o teste chi quadrado dos valores da coluna deviance. Esse teste nos permite concluir que os grãos de pólen e a irrigação contribuem significativamente com a verossimilhança do modelo.
Outra opção usar inferência multimodelo, ou seleção de modelos, baseada em algum dos critérios de informação para comparar os modelos aninhados descritos acima. Nesse exemplo vamos usar a função dredge() e o critério de informação de Akaike (AIC)
options(na.action = "na.fail")
dredge(logit_fruit)## Fixed term is "(Intercept)"## Global model call: glm(formula = factor(FRUIT) ~ POLLEN * factor(IRRIGATION), family = binomial(link = "logit"),
## data = ourdata)
## ---
## Model selection table
## (Int) fct(IRR) POL fct(IRR):POL df logLik AICc delta weight
## 4 -5.199 + 0.6071 3 -81.575 169.3 0.00 0.717
## 8 -5.600 + 0.6644 + 4 -81.463 171.1 1.86 0.283
## 2 -1.346 + 2 -107.589 219.2 49.97 0.000
## 3 -2.076 0.3800 2 -120.977 246.0 76.74 0.000
## 1 0.000 1 -138.629 279.3 110.01 0.000
## Models ranked by AICc(x)Essa função gera a tabela acima com os modelos ordenados pelo menor valor de AIC. O AIC é o logLik penalizado pelo número de parâmetros do modelo. O AICc nada mais é do que o AIC corrigido para amostras pequenas. Quanto menor o valor de AIC, melhor o modelo. No nosso caso, o melhor modelo é o que contem as variáveis irrigação e pólen.
Se compararmos os resultados que obtivemos a partir da ANDEVA e da inferência multimodelos podemos concluir que chegaremos a conclusões muito similares. A principal diferença é sobre o papel da interação entre as variáveis explicativas no modelo. Na análise do AIC vemos que o modelo com a interação poderia ser considerado como igualmente plausível ao modelo sem interação (delta < 2), já na ANDEVA o valor de p está acima do valor de significância. Essa diferença pode estar associada ao fato da ANDEVA testar o efeito dos parâmetros de maneira sequêncial, como em uma ANOVA do tipo I.