Introdução
Aqui, vamos aprender o que é e como ajustar um Modelo Linear Generalizado (GLM) no R. Um GLM é uma extensão do modelo linear clássico, permitindo ajustar modelos a dados com distribuições diferentes da distribuição normal. Um exemplo disso poderia ser uma variável binária, como presença e ausência que é muito utilizada na ecologia. Para isso, vamos utilizar um conjunto de dados simulados da seleção de habitats de aranhas.
Conjunto de dados
Vamos começar carregando e explorando o conjunto de dados. Suponha que estamos interessados em prever a probabilidade de um evento ocorrer (variável binária) com base em uma variável preditora contínua.
# Pacotes nescessários####
#install.packages("ggplot2")
#install.packages("dplyr")
#install.packages("MASS")
require(ggplot2)#Pacote para gerar plots melhores
require(ggpubr)#Para paineis gráficos
require(ggthemes)#temas do ggplot
require(dplyr)#Para manipulação de dados
require(MASS)#Para GLM
# Gerando um conjunto de dados
set.seed(7) # Para padronizar o clock do processador
dados = data.frame(
densidade_bromelias = rnorm(100, mean = 10, sd = 3),
numero_folhas = rpois(100, lambda = 15),
altura_planta = rnorm(100, mean = 100, sd = 20),
altura_chao = rnorm(100, mean = 30, sd = 5),
raio_planta = rnorm(100, mean = 50, sd = 10),
dossel = rnorm(100, mean = 80, sd = 15),
presenca_aranha = rbinom(100, size = 1, prob = 0.4)
)
# Visualizando os primeiros registros
head(dados)## densidade_bromelias numero_folhas altura_planta altura_chao raio_planta
## 1 16.861741 17 109.45093 30.43845 65.66575
## 2 6.409685 17 88.82444 33.83762 41.55155
## 3 7.917122 14 124.84129 26.93421 60.97649
## 4 8.763121 10 99.69483 30.76251 53.00880
## 5 7.087980 14 84.16423 27.04146 42.23034
## 6 7.158160 13 91.96763 25.19648 49.61566
## dossel presenca_aranha
## 1 69.64818 0
## 2 52.41526 0
## 3 84.24557 0
## 4 46.21115 1
## 5 91.01525 1
## 6 82.18709 0Avaliando distribuição das variáveis
Agora vamos explorar a relação entre cada variável preditora e a ocorrência de aranhas (variável resposta). Vamos plotar um gráfico de dispersão para cada par de variáveis.
# Gráficos de dispersão
ggplot(dados, aes(x = densidade_bromelias, y = presenca_aranha)) +
geom_point(color="#006699") +
labs(title = "Densidade de Bromélias e Presença de Aranhas", x = "Densidade de Bromélias", y = "Presença de Aranhas") + theme_classic() + theme(plot.title = element_text(hjust = 0.5))
ggplot(dados, aes(x = numero_folhas, y = presenca_aranha)) +
geom_point(color="#006699") +
labs(title = "Número de Folhas e Presença de Aranhas", x = "Número de Folhas", y = "Presença de Aranhas") + theme_classic() + theme(plot.title = element_text(hjust = 0.5))
ggplot(dados, aes(x = altura_planta, y = presenca_aranha)) +
geom_point(color="#006699") +
labs(title = "Altura da Planta e Presença de Aranhas", x = "Altura da Planta", y = "Presença de Aranhas") + theme_classic() + theme(plot.title = element_text(hjust = 0.5))
ggplot(dados, aes(x = altura_chao, y = presenca_aranha)) +
geom_point(color="#006699") +
labs(title = "Altura do Chão e Presença de Aranhas", x = "Altura do Chão", y = "Presença de Aranhas") + theme_classic() + theme(plot.title = element_text(hjust = 0.5))
ggplot(dados, aes(x = raio_planta, y = presenca_aranha)) +
geom_point(color="#006699") +
labs(title = "Raio da Planta e Presença de Aranhas", x = "Raio da Planta", y = "Presença de Aranhas") + theme_classic() + theme(plot.title = element_text(hjust = 0.5))
ggplot(dados, aes(x = dossel, y = presenca_aranha)) +
geom_point(color="#006699") +
labs(title = "Dossel e Presença de Aranhas", x = "Dossel", y = "Presença de Aranhas") + theme_classic() + theme(plot.title = element_text(hjust = 0.5))
Outra forma de visualizar é usando boxplots
ggplot(dados, aes(x = factor(presenca_aranha), y = densidade_bromelias)) +
geom_boxplot(fill = c("#006666", "#006699")) +
labs(title = "Variação da Densidade de Bromélias entre Presença e Ausência de Aranhas",
x = "Presença de Aranhas", y = "Densidade de Bromélias") + theme_classic() + theme(plot.title = element_text(hjust = 0.5))
ggplot(dados, aes(x = factor(presenca_aranha), y = numero_folhas)) +
geom_boxplot(fill = c("#006666", "#006699")) +
labs(title = "Variação do Número de Folhas entre Presença e Ausência de Aranhas",
x = "Presença de Aranhas", y = "Número de Folhas") + theme_classic() + theme(plot.title = element_text(hjust = 0.5))
ggplot(dados, aes(x = factor(presenca_aranha), y = altura_planta)) +
geom_boxplot(fill = c("#006666", "#006699")) +
labs(title = "Variação da Altura da Planta entre Presença e Ausência de Aranhas",
x = "Presença de Aranhas", y = "Altura da Planta") + theme_classic() + theme(plot.title = element_text(hjust = 0.5))
ggplot(dados, aes(x = factor(presenca_aranha), y = altura_chao)) +
geom_boxplot(fill = c("#006666", "#006699")) +
labs(title = "Variação da Altura do Chão entre Presença e Ausência de Aranhas",
x = "Presença de Aranhas", y = "Altura do Chão") + theme_classic() + theme(plot.title = element_text(hjust = 0.5))
ggplot(dados, aes(x = factor(presenca_aranha), y = raio_planta)) +
geom_boxplot(fill = c("#006666", "#006699")) +
labs(title = "Variação do Raio da Planta entre Presença e Ausência de Aranhas",
x = "Presença de Aranhas", y = "Raio da Planta") + theme_classic() + theme(plot.title = element_text(hjust = 0.5))
ggplot(dados, aes(x = factor(presenca_aranha), y = dossel)) +
geom_boxplot(fill = c("#006666", "#006699")) +
labs(title = "Variação do Dossel entre Presença e Ausência de Aranhas",
x = "Presença de Aranhas", y = "Dossel") + theme_classic() + theme(plot.title = element_text(hjust = 0.5))
Ajustando o modelo GLM
Agora vamos ajustar um modelo linear generalizado para prever a presença de aranhas com base nas variáveis preditoras. Neste exemplo, usaremos a distribuição binomial e a função de ligação logit. A função logit no GLM é uma transformação matemática que mapeia a probabilidade de sucesso em um preditor linear. Em outras palavras, é uma forma de “forçar” que a variavél resposta seja ajustada em relação a um preditor linear. Essa função de ligação é essencial para o ajuste adequado do modelo e para a interpretação dos resultados no contexto de um modelo de regressão logística binária.
# Ajustando o modelo GLM
modelo_glm = glm(presenca_aranha ~ densidade_bromelias + numero_folhas + altura_planta + altura_chao + raio_planta + dossel,
data = dados, family = binomial(link = "logit"))
# Visualizando os resultados do modelo
summary(modelo_glm)##
## Call:
## glm(formula = presenca_aranha ~ densidade_bromelias + numero_folhas +
## altura_planta + altura_chao + raio_planta + dossel, family = binomial(link = "logit"),
## data = dados)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.4265 -0.8617 -0.6202 1.1362 1.7244
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.80960 2.74227 -0.660 0.509
## densidade_bromelias -0.08216 0.08198 -1.002 0.316
## numero_folhas 0.09309 0.06188 1.504 0.133
## altura_planta 0.01312 0.01061 1.236 0.216
## altura_chao 0.06968 0.04741 1.470 0.142
## raio_planta -0.03412 0.02306 -1.480 0.139
## dossel -0.01494 0.01670 -0.895 0.371
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 126.84 on 99 degrees of freedom
## Residual deviance: 116.41 on 93 degrees of freedom
## AIC: 130.41
##
## Number of Fisher Scoring iterations: 4Interpretando os resultados
O modelo GLM fornece números importantes para entender como as características das plantas (densidade de bromélias, número de folhas, altura, raio e dossel) afetam a presença de aranhas. Esses números mostram se as características aumentam ou diminuem a probabilidade de encontrar aranhas. Para facilitar a compreensão, usamos uma fórmula matemática chamada “função logística” para transformar esses números em chances reais de encontrar aranhas com base nas características das plantas. Isso nos ajuda a entender onde podemos encontrar aranhas com base na caracteristica da bromelia e densidade de plantas.
# Obtendo as probabilidades preditas
dados$prob_predita = predict(modelo_glm, type = "response")
# Visualizando as probabilidades preditas
ggplot(dados, aes(x = densidade_bromelias, y = prob_predita)) +
geom_point(color="#006699") +
labs(title = "Probabilidade Preditiva da Presença de Aranhas", x = "Densidade de Bromélias", y = "Probabilidade") + theme_classic() + theme(plot.title = element_text(hjust = 0.5))
ggplot(dados, aes(x = numero_folhas, y = prob_predita)) +
geom_point(color="#006699") +
labs(title = "Probabilidade Preditiva da Presença de Aranhas", x = "Número de Folhas", y = "Probabilidade") + theme_classic() + theme(plot.title = element_text(hjust = 0.5))
ggplot(dados, aes(x = altura_planta, y = prob_predita)) +
geom_point(color="#006699") +
labs(title = "Probabilidade Preditiva da Presença de Aranhas", x = "Altura da Planta", y = "Probabilidade") + theme_classic() + theme(plot.title = element_text(hjust = 0.5))
ggplot(dados, aes(x = altura_chao, y = prob_predita)) +
geom_point(color="#006699") +
labs(title = "Probabilidade Preditiva da Presença de Aranhas", x = "Altura do Chão", y = "Probabilidade") + theme_classic() + theme(plot.title = element_text(hjust = 0.5))
ggplot(dados, aes(x = raio_planta, y = prob_predita)) +
geom_point(color="#006699") +
labs(title = "Probabilidade Preditiva da Presença de Aranhas", x = "Raio da Planta", y = "Probabilidade") + theme_classic() + theme(plot.title = element_text(hjust = 0.5))
ggplot(dados, aes(x = dossel, y = prob_predita)) +
geom_point(color="#006699") +
labs(title = "Probabilidade Preditiva da Presença de Aranhas", x = "Dossel", y = "Probabilidade") + theme_classic() + theme(plot.title = element_text(hjust = 0.5))
Para facilitar a comparação entre as variáveis podemos usar o painéis:
#Primeiro vamos criar um objeto para cada um dos plots anteriores
plot1 = ggplot(dados, aes(x = densidade_bromelias, y = prob_predita)) +
geom_point(color="#006699") +
labs( x = "Densidade de Bromélias", y = "Probabilidade") +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5))
plot2 = ggplot(dados, aes(x = numero_folhas, y = prob_predita)) +
geom_point(color="#006699") +
labs( x = "Número de Folhas", y = "Probabilidade") +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5))
plot3 = ggplot(dados, aes(x = altura_planta, y = prob_predita)) +
geom_point(color="#006699") +
labs( x = "Altura da Planta", y = "Probabilidade") +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5))
plot4 = ggplot(dados, aes(x = altura_chao, y = prob_predita)) +
geom_point(color="#006699") +
labs( x = "Altura do Chão", y = "Probabilidade") +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5))
plot5 = ggplot(dados, aes(x = raio_planta, y = prob_predita)) +
geom_point(color="#006699") +
labs( x = "Raio da Planta", y = "Probabilidade") +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5))
plot6 = ggplot(dados, aes(x = dossel, y = prob_predita)) +
geom_point(color="#006699") + geom_smooth(color="#006666", method = "glm", method.args = list(family = binomial(link = "logit"))) +
labs( x = "Dossel", y = "Probabilidade") +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5))
#Agora vamos usar o ggarrange() para dispor os gráficos em uma matriz de 2x3 com um tema chique
ggarrange(plot1, plot2, plot3, plot4, plot5, plot6, ncol = 3, nrow = 2) + theme_minimal() 
Conclusão
Este é apenas um exemplo didático para introduzir o conceito de GLM. Em casos reais, é necessário considerar a adequação do modelo, diagnósticos e outras análises adicionais, mas já é útil para você começar a aprender sobre Modelos Lineares Generalizados!