MLGExemplo11
Fernando L Pérez
20/07/2021
O que são modelos de regressão de Poisson?
Os modelos de regressão Poisson são usados para modelar eventos onde os resultados são contagens: dados discretos com valores inteiros não negativos que contam algo, como o número de vezes que um evento ocorre durante um determinado período de tempo ou o número de pessoas na fila do supermercado.
Os dados de contagem também podem ser expressos como dados de taxa, uma vez que o número de vezes que um evento ocorre dentro de um período de tempo pode ser expresso como uma contagem bruta, ou seja, em um dia, comemos três refeições.
A regressão Poisson nos ajuda a analisar os dados de contagem e os dados de taxa, permitindo-nos determinar quais variáveis explicativas têm efeito em uma determinada variável de resposta \(Y\), a contagem ou uma taxa. Por exemplo, a regressão Poisson pode ser aplicada por um supermercado para entender e prever melhor o número de pessoas em uma fila.
A distribuição Poisson modela a probabilidade do evento ou eventos \(Y\) ocorrerem dentro de um período de tempo específico, assumindo que as ocorrências \(Y\) não são afetadas pelo tempo das ocorrências anteriores de \(Y\). Isso pode ser expresso matematicamente usando a seguinte expressão: \[\begin{equation*} P(Y=y)=\dfrac{e^{\mu t} (\mu t)^y}{y!}, \end{equation*}\] onde \(y=0,1,2,\cdots\).
Aqui, \(\mu\) é o número médio de vezes que um evento pode ocorrer por unidade de exposição. Também é chamado de parâmetro de distribuição Poisson. A exposição pode ser tempo, espaço, tamanho da população, distância ou área, mas geralmente é tempo, denotado como \((t)\). Se o valor da exposição não for fornecido, será considerado igual a 1.
Vamos visualizar isso criamos um gráfico de distribuições Poisson para diferentes valores de \(\mu\).
Primeiro, vamos criar um vetor de 6 cores:
cores <- c("Red", "Blue", "Gold", "Black", "Pink", "Green")A seguir, criaremos uma lista para a distribuição que terá diferentes valores para \(\mu\):
poisson.dist <- list()
Em seguida, vamos criar um vetor de valores para \(\mu\) e fazer um ciclo sobre os valores de \(\mu\) cada um com intervalo de quantis 0-20, armazenando os resultados em uma lista:
# Vetor de valores de mu
a <- c(1, 2, 3, 4, 5, 6)
for (i in 1:6) {
poisson.dist[[i]] <- c(dpois(0:20, a[i])) # Armazenando o vetor de distribuição para cada valor correspondente de mu
}Finalmente, vamos representar graficamente os pontos usando plot() a qual é uma função gráfica base no R.
plot(unlist(poisson.dist[1]), type = "o", xlab="y", ylab = "P(y)")
for (i in 1:6) {
lines(unlist(poisson.dist[i]), type = "o", col = cores[i], pch = 19)
}
grid()
legend("topright", legend = a, inset = 0.08, cex = 1.0, fill = cores, title = expression(paste("Valores de ",mu)))
Modelos de regressão Poisson
Modelos Lineares Generalizados (MLG) são modelos nos quais as variáveis de resposta seguem uma distribuição diferente da distribuição normal. Isso está em contraste com os modelos de regressão linear, nos quais as variáveis de resposta seguem uma distribuição normal. Isso ocorre porque os Modelos Lineares Generalizados têm variáveis de resposta que são categóricas, como Sim, Não; ou Grupo A, Grupo B e, portanto, não variam de \(-\infty\) a \(+\infty\).
Portanto, a relação entre a resposta e as variáveis preditoras pode não ser linear. Nos modelos lineares generalizados: \[\begin{equation*} g(\mu_i)=\alpha+\beta_1 x_{1i}+\beta_2 x_{2i}+\cdots+\beta_p x_{pi}, \qquad i=1,2,\cdots,p, \end{equation*}\] sendo \(g(\cdot)\) a função de ligação escolhida.
Um modelo de regressão Poisson é um modelo linear generalizado que é usado para modelar dados de contagem e tabelas de contingência. A saída \(Y\) (contagem) é um valor que segue a distribuição de Poisson. Ele assume o logaritmo dos valores esperados (média) que podem ser modelados em uma forma linear por alguns parâmetros desconhecidos.
Para transformar o relacionamento não linear para a forma linear, uma função de ligação é usada, que é o logaritmo da regressão Poisson. Por esse motivo, um modelo de regressão Poisson também é chamado de modelo log-linear. A forma matemática geral do modelo de regressão de Poisson é: \[\begin{equation*} \log(\mu_i)=\alpha+\beta_1 x_{1i}+\beta_2 x_{2i}+\cdots+\beta_p x_{pi}, \qquad i=1,2,\cdots,p, \end{equation*}\] onde \(\mu_i\) são as esperanças das variáveis resposta \(Y_i\).
Os coeficentes \(\alpha\) e \(\beta\) são numéricos, sendo \(\alpha\) o intercepto, às vezes \(\alpha\) também é representado por \(\beta_0\) e \(x\) são as variáveis preditoras ou explicativas.
Considere uma equação com uma variável preditora e uma variável de resposta: \[\begin{equation*} \log(\mu) = \alpha + \beta x\cdot \end{equation*}\]
Isso é equivalente a: \[\begin{equation*} \mu = e^{\alpha+\beta x} = e^\alpha \times e^{\beta x}\cdot \end{equation*}\] Em modelos de regressão Poisson, as variáveis preditoras ou explicativas podem ter uma mistura de valores numéricos ou categóricos.
Uma das características mais importantes da distribuição Poisson e da regressão Poisson é a equidispersão, o que significa que a média e a variância da distribuição são iguais.
Digamos que a média \(\mbox{E}(Y)=\mu\), para a regressão Poisson, a média e a variância estão relacionadas como: \[\begin{equation*} \mbox{Var}(Y) = \sigma^2 \mbox{E}(Y), \end{equation*}\] onde \(\sigma^2\) é o parâmetro de dispersão. Uma vez que \(\mbox{var}(Y) = \mbox{E}(Y)\) (variância = média) deve valer para que o modelo de Poisson seja completamente ajustado, \(\sigma^2\) deve ser igual a 1.
Quando a variância é maior que a média, isso é chamado de sobredispersão e \(\sigma^2\) é maior que 1. Se \(\sigma^2\) for menor que 1, então é conhecido como subdispersão.
Exemplo
library(datasets)
dados <- warpbreaks
columns <- names(dados) # Extraindo os nomes das colunas do arquivo de dados
columns # mostrando as colunas## [1] "breaks" "wool" "tension"ls.str(dados)## breaks : num [1:54] 26 30 54 25 70 52 51 26 67 18 ...
## tension : Factor w/ 3 levels "L","M","H": 1 1 1 1 1 1 1 1 1 2 ...
## wool : Factor w/ 2 levels "A","B": 1 1 1 1 1 1 1 1 1 1 ...De acima, podemos ver os tipos e níveis presentes nos dados. Agora vamos trabalhar com o arquivo de dados. Lembre-se, com um modelo de distribuição Poisson, estamos tentando descobrir como algumas variáveis preditoras afetam uma variável de resposta. Aqui, breaks (quebras) é a variável de resposta e wool (lã) e tension (tensão) são variáveis de previsão.
Podemos ver a possível relação entras as variáveis explicativas e a variável dependente breaks construindo gráficos e dispersão.
Estudo descritivo
library(ggplot2)
ggplot(dados, aes(x = tension, y = breaks)) + geom_point()
# Mudando o tamanho dos pontos
ggplot(dados, aes(x = wool, y = breaks)) + geom_point(size=2, shape=19)
Tudo indica que há diferenças entre os níveis de tension, mas não parece-me tão claro se há diferenças nos níveis de wool.
Vamos verificar a média e variância da variável dependente:
mean(dados$breaks)## [1] 28.14815var(dados$breaks)## [1] 174.2041A variância é muito maior do que a média, o que sugere que teremos superdispersão no modelo.
Modelos
modelo1 = glm(breaks~tension+wool, data = dados, family = poisson(link = "log"))
modelo2 = glm(breaks~tension+wool, data = dados, family = poisson(link = "identity"))
modelo3 = glm(breaks~tension+wool, data = dados, family = poisson(link = "sqrt"))Vamos escolher qual a função de ligação mais adequada.
AIC(modelo1,modelo2,modelo3)## df AIC
## modelo1 4 493.0560
## modelo2 4 497.3612
## modelo3 4 495.3462Escolhemos a ligação canônica (menor AIC).
summary(modelo1)##
## Call:
## glm(formula = breaks ~ tension + wool, family = poisson(link = "log"),
## data = dados)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.6871 -1.6503 -0.4269 1.1902 4.2616
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.69196 0.04541 81.302 < 2e-16 ***
## tensionM -0.32132 0.06027 -5.332 9.73e-08 ***
## tensionH -0.51849 0.06396 -8.107 5.21e-16 ***
## woolB -0.20599 0.05157 -3.994 6.49e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 297.37 on 53 degrees of freedom
## Residual deviance: 210.39 on 50 degrees of freedom
## AIC: 493.06
##
## Number of Fisher Scoring iterations: 4anova(modelo1, test="Chisq")## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: breaks
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 53 297.37
## tension 2 70.942 51 226.43 3.938e-16 ***
## wool 1 16.039 50 210.39 6.206e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1library(car)
residualPlots(modelo1)
influencePlot(modelo1)
## StudRes Hat CookD
## 1 -2.476785 0.08274043 0.12222528
## 2 -1.741130 0.08274043 0.06279697
## 5 4.490911 0.08274043 0.54693078
## 9 4.071458 0.08274043 0.44260695Observamos nos resíduos suspeita de observações com excesso de influência na resposta e a relação entre o preditor linear e os residuos não é linear.
Interpretando o modelo de Poisson
Acabamos de receber muitas informações, agora precisamos interpretá-las. A primeira coluna denominada Estimate são os valores dos coeficientes de \(\alpha\), o intercepto, \(\beta_1\) e assim por diante. A seguir está a interpretação para as estimativas dos parâmetros:
- \(\exp(\alpha)\) = efeito na média \(\mu\), quando \(X = 0\).
- \(\exp(\beta)\) = com cada aumento de unidade em \(X\), a variável preditora tem efeito multiplicativo de \(\exp(\beta)\) na média de \(Y\), ou seja, \(\mu\).
- Se \(\beta = 0\), então \(\exp(\beta) = 1\), e a contagem esperada é \(\exp(\alpha)\) e, \(Y\) e \(X\) não estão relacionados.
- Se \(\beta> 0\), então \(\exp(\beta)> 1\), e a contagem esperada é \(\exp(\beta)\) vezes maior do que quando \(X = 0\).
- Se \(\beta <0\), então \(exp(\beta) <1\), e a contagem esperada é \(\exp(\beta)\) vezes menor do que quando \(X = 0\).
R trata as variáveis categóricas como variáveis dummy. Variáveis categóricas, também chamadas de variáveis indicadoras, são convertidas em variáveis dummy atribuindo aos níveis na variável alguma representação numérica. A regra geral é que se houver \(k\) categorias em uma variável de fator, a saída de glm() terá \(k-1\) categorias com 1 restante como a categoria base.
Podemos ver no resumo acima que para lã, ‘A’ foi feito a base e não é mostrado no resumo. Da mesma forma, para a tensão ‘L’ foi feita a categoria de base.
Para ver quais variáveis explicativas têm efeito na variável de resposta, examinaremos os p-valores. Se o p-valor for menor que 0.05 ou 0.10, a variável tem efeito na variável de resposta. No resumo acima, podemos ver que todos os p-valores são menores que 0.06, portanto, ambas as variáveis explicativas (lã e tensão) têm efeito significativo nas rupturas. Observe como a saída do R usou *** no final de cada variável. O número de estrelas significa significância.
Antes de começar a interpretar os resultados, vamos verificar se o modelo tem superdispersão ou subdispersão. Se o desvio residual (residual deviance) for maior do que os graus de liberdade, existe uma superdispersão. Isso significa que as estimativas estão corretas, mas os erros padrão (desvio padrão) estão errados e não são contabilizados pelo modelo.
O desvio nulo (null deviance) mostra quão bem a variável de resposta é prevista por um modelo que inclui apenas o intercepto (grande média) enquanto residual com a inclusão de variáveis independentes. Acima, podemos ver que a adição de 3 (53-50 = 3) variáveis independentes diminuiu o desvio para 210.39 de 297.37. Maior diferença de valores significa um ajuste ruim.
Modelos superdispersos
Então, para ter um erro padrão mais correto, podemos usar um modelo quase-poisson:
modelo4 <- glm(breaks ~ tension+wool, data = dados, family = quasipoisson(link = "log"))
summary(modelo4)##
## Call:
## glm(formula = breaks ~ tension + wool, family = quasipoisson(link = "log"),
## data = dados)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.6871 -1.6503 -0.4269 1.1902 4.2616
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.69196 0.09374 39.384 < 2e-16 ***
## tensionM -0.32132 0.12441 -2.583 0.012775 *
## tensionH -0.51849 0.13203 -3.927 0.000264 ***
## woolB -0.20599 0.10646 -1.935 0.058673 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 4.261537)
##
## Null deviance: 297.37 on 53 degrees of freedom
## Residual deviance: 210.39 on 50 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4Comparando os modelos
Agora que temos dois modelos diferentes, vamos compará-los para ver qual é o melhor. Primeiro, instalaremos a biblioteca arm porque ela contém uma função de que precisamos:
# carregando a biblioteca arm que contém a função se.coef()
library(arm)
# extraindo os coeficientes do primeiro modelo usando coef()
coef1 = coef(modelo1)
# extraindo os coeficientes do segundo modelo usando coef()
coef2 = coef(modelo4)
# extraindo os erros padrão do primeiro modelo usando se.coef()
se.coef1 = se.coef(modelo1)
# extraindo os erros padrão do segundo modelo usando se.coef()
se.coef2 = se.coef(modelo4)
# usando cbind() para combinar valores em um arquivo de dados
ambos.modelos = cbind(coef1, se.coef1, coef2, se.coef2, exponent = exp(coef1))
# mostrando o arquivode dados
ambos.modelos## coef1 se.coef1 coef2 se.coef2 exponent
## (Intercept) 3.6919631 0.04541069 3.6919631 0.09374352 40.1235380
## tensionM -0.3213204 0.06026580 -0.3213204 0.12440965 0.7251908
## tensionH -0.5184885 0.06395944 -0.5184885 0.13203462 0.5954198
## woolB -0.2059884 0.05157117 -0.2059884 0.10646089 0.8138425Na saída acima, podemos ver que os coeficientes são os mesmos, mas os erros padrão são diferentes.
Tendo esses pontos em mente, vamos ver a estimativa para lã (woolB). Seu valor é -0.2059884 e o expoente de -0.2059884 é 0.8138425.
1-ambos.modelos[4,5]## [1] 0.1861575Isso mostra que a mudança da lã do tipo A para a lã do tipo B resulta em uma diminuição nas rupturas 0.8138425 vezes o intercepto, porque a estimativa -0.2059884 é negativa. Outra maneira de dizer isso é se mudarmos o tipo de lã de A para B, o número de quebras cairá 18.6%, assumindo que todas as outras variáveis permaneçam iguais.
residualPlots(modelo4)
influencePlot(modelo4)
## StudRes Hat CookD
## 1 -1.2132935 0.08274043 0.02868103
## 2 -0.8464275 0.08274043 0.01473576
## 5 2.2770892 0.08274043 0.12834120
## 9 2.0457195 0.08274043 0.10386087Predizendo a partir do modelo
Uma vez que o modelo é feito, podemos usar a predict(model, data, type) para prever resultados usando novos arquivos de dados contendo dados diferentes dos dados de treinamento. Vejamos um exemplo.
# preparando um nvo arquivo com novos dados
novos.dados = data.frame(wool = "B", tension = "M")
# usando 'predict()' para executar o modelo em novos dados
predict(modelo2, newdata = novos.dados, type = "response")## 1
## 24.38907Nosso modelo está prevendo que haverá cerca de 24 rupturas com lã tipo B e nível de tensão M.
Visualizando descobertas
Quando você está compartilhando sua análise com outras pessoas, as tabelas geralmente não são a melhor maneira de chamar a atenção das pessoas. Plotagens e gráficos ajudam as pessoas a compreender suas descobertas mais rapidamente. A maneira mais popular de visualizar dados em R é provavelmente ggplot2 (que é ensinado no curso de visualização de dados do Dataquest), também vamos usar um pacote R incrível chamado jtools que inclui ferramentas para resumir e visualizar modelos de regressão especificamente. Vamos usar jtools para visualizar poisson.model2.
library("jtools"); library("ggstance"); library(broom); library(broom.mixed); library(interactions)O pacote jtools fornece as funções plot_summs() e plot_coefs() para visualizar o resumo do modelo e também nos permite comparar diferentes modelos com ggplot2.
# mostrando os coeficientes de regressão para modelo2
plot_summs(modelo1, escala = TRUE, exp = TRUE)
# mostrando os coeficientes de regressão para modelo4 e modelo1
plot_summs(modelo1, modelo4, escala = TRUE, exp = TRUE)
No código acima, plot_summs(modelo2, scale = TRUE, exp = TRUE) mostra o segundo modelo usando a família quase-poisson em glm.
- O primeiro argumento em plot_summs() é o modelo de regressão a ser usado, pode ser um ou mais de um.
- scale ajuda com o problema de escalas diferentes das variáveis.
- exp é definido como TRUE porque para a regressão Poisson é mais provável que estejamos interessados em valores exponenciais das estimativas em vez de lineares.
Também podemos visualizar a interação entre variáveis preditoras. interactions fornece funções diferentes para diferentes tipos de variáveis. Por exemplo, se todas as variáveis forem categóricas, poderíamos usar cat_plot() para entender melhor as interações entre elas. Para variáveis contínuas, interact_plot() é usado.
Neste arquivo de dados temos variáveis somente preditoras categóricas, então usaremos cat_plot() para visualizar a interação entre elas, fornecendo argumentos que especificam qual modelo gostaríamos de usar, a variável preditora que estamos olhando e a outra variável preditora com a qual ele se combina para produzir o resultado.
# modelo com interações
modelo5 <- glm(breaks ~ tension*wool, data = dados, family = quasipoisson(link = "log"))
summary(modelo5)##
## Call:
## glm(formula = breaks ~ tension * wool, family = quasipoisson(link = "log"),
## data = dados)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.3383 -1.4844 -0.1291 1.1725 3.5153
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.79674 0.09688 39.189 < 2e-16 ***
## tensionM -0.61868 0.16374 -3.778 0.000436 ***
## tensionH -0.59580 0.16253 -3.666 0.000616 ***
## woolB -0.45663 0.15558 -2.935 0.005105 **
## tensionM:woolB 0.63818 0.23699 2.693 0.009727 **
## tensionH:woolB 0.18836 0.25201 0.747 0.458436
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 3.76389)
##
## Null deviance: 297.37 on 53 degrees of freedom
## Residual deviance: 182.31 on 48 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4cat_plot(modelo5, pred = wool, modx = tension)
# argumento 1: modelo de regressão
# pred: a variável categórica que aparecerá no eixo x
# modx: variável do moderador que tem um efeito em combinação para prever o resultadoObservemos a vantagem de termos considerado um modelo com interações.
residualPlots(modelo5)
influencePlot(modelo5)
## StudRes Hat CookD
## 1 -1.658321 0.1111111 0.04811942
## 4 -1.762615 0.1111111 0.05344570
## 5 1.988905 0.1111111 0.09048121
## 24 1.899906 0.1111111 0.08626919Podemos fazer a mesma coisa para observar a tensão:
cat_plot(modelo5, pred = tension, modx = wool)
Acima, vemos como as três categorias diferentes de tensão (L, M e H) para cada uma afetam as rupturas em cada tipo de lã. Por exemplo, as quebras tendem a ser maiores com baixa tensão e lã tipo A.
Também podemos definir o tipo de gráfico criado por cat_plot() usando o parâmetro geom. Este parâmetro melhora a interpretação do gráfico. Podemos usá-lo dessa forma, passando geom como um argumento adicional para cat_plot:
cat_plot(modelo5, pred = tension, modx = wool, geom = "line")
Também podemos incluir observações no gráfico adicionando plot.points = TRUE:
cat_plot(modelo5, pred = tension, modx = wool, geom = "line", plot.points = TRUE)
Existem muitas outras opções de design, incluindo estilo de linha, cor, etc., que nos permitem personalizar a aparência dessas visualizações.