Modelo de regressão para avaliar os pesos de pinguins no arquipélago Palmer
Eduardo Freire Mangabeira, Julio Cesar dos Santos Magalhães
23/01/2022
Os dados de palmerpenguins contêm medidas de tamanho para três espécies de pinguins observadas em três ilhas no Arquipélago de Palmer, na Antártida.
Esses dados foram coletados de 2007 a 2009 pela Dra. Kristen Gorman com o Palmer Station Long Term Ecological Research Program, parte da US Long Term Ecological Research Network.
Fonte: https://allisonhorst.github.io/palmerpenguins/
Modelo de regressão
No presente trabalho, propomos um modelo de regressão para avaliar a influência das variáveis sexo, ano, ilha, espécie, comprimento do bico, profundidade do bico e comprimento da nadadeira na massa corporal em gramas de um pinguim do Arquipélago Palmer, que pode ser da espécie Adelie(pinguim-de-Adélia), Gentoo ou Chinstrap(Pinguim-de-barbicha). Descobrimos que o modelo que inclui as variáveis [blablabla] teve uma performance superior comparado aos outros e que existe[ou não existe] uma relação linear entre [blablabla] e a massa corporal dos pinguins.
Váriaveis explicativas
- sex - Sexo do Pinguim. (variável binária)
- year - Ano de coleta do dado. (variável qualitativa)
- island - Ilha de origem do pinguim. (variável qualitativa)
- species - Espécie do pinguim.(variável qualitativa)
- bill_depth_mm - Profundidade do bico do pinguim em milímetros.(variável quantitativa contínua)
- flipper_length_mm - Comprimento da nadadeira do Pinguim em milímetros. (variável quantitativa contínua)
- body_mass_g - Massa corporal em gramas do pinguim.
Variável resposta
- bill_length_mm - Comprimento do bico do pinguim em milímetros. (variável quantitativa - contínua)
Dataset
# install.packages("palmerpenguins")
library(palmerpenguins)
dataset <- palmerpenguins::penguinsLimpeza do dataset
Valores faltantes
Existem 11 linhas com valores faltantes no dataset. Para não prejudicar o experimento, retiramos todas as linhas em que pelo menos uma coluna não tenha um dado.
library(tidyr)
library(dplyr)
library(purrr)
dados <- dataset %>%
tidyr::drop_na(names(dataset)) %>%
dplyr::filter(species=="Gentoo") %>%
dplyr::select(bill_length_mm, everything(), -year)Conferindo outliers
Observando os box-plots é possível ver que não há valores discrepantes para quase todas as variáveis numéricas, com exceção da variável resposta.
boxplot(dados$"body_mass_g", ylab="body_mass_g")
boxplot(dados$"bill_depth_mm", ylab="bill_depth_mm")
boxplot(dados$"flipper_length_mm", ylab="flipper_length_mm")
boxplot(dados$"bill_length_mm", ylab="bill_length_mm")
Vamos testar o modelo retirando o valor discrepante usando o teste de Grubbs para identificá-lo.
library(outliers)## Warning: package 'outliers' was built under R version 4.0.3outlier <- outliers::grubbs.test(dados$bill_length_mm)
if(outlier$p.value < 0.05){
dados <- dados %>%
dplyr::filter(bill_length_mm != max(bill_length_mm))
}
boxplot(dados$bill_length_mm, ylab="bill_length_mm")
Transformando a variável sexo em numérica
Para poder realizar operações com a variável sexo, escolhemos “feminino” como categoria de referência e transformamos a coluna em numérica.
dados$sex <- ifelse(dados$sex=="female", 1, 0)One-hot encoding das variáveis categóricas
Criaremos algumas variáveis “Dummy” para simular variáveis indicadoras para os campos espécie e ilha.
library(mltools)
library(data.table)
dados <- mltools::one_hot(data.table::as.data.table(dados))
dados## bill_length_mm species_Adelie species_Chinstrap species_Gentoo
## 1: 46.1 0 0 1
## 2: 50.0 0 0 1
## 3: 48.7 0 0 1
## 4: 50.0 0 0 1
## 5: 47.6 0 0 1
## ---
## 114: 47.2 0 0 1
## 115: 46.8 0 0 1
## 116: 50.4 0 0 1
## 117: 45.2 0 0 1
## 118: 49.9 0 0 1
## island_Biscoe island_Dream island_Torgersen bill_depth_mm
## 1: 1 0 0 13.2
## 2: 1 0 0 16.3
## 3: 1 0 0 14.1
## 4: 1 0 0 15.2
## 5: 1 0 0 14.5
## ---
## 114: 1 0 0 13.7
## 115: 1 0 0 14.3
## 116: 1 0 0 15.7
## 117: 1 0 0 14.8
## 118: 1 0 0 16.1
## flipper_length_mm body_mass_g sex
## 1: 211 4500 1
## 2: 230 5700 0
## 3: 210 4450 1
## 4: 218 5700 0
## 5: 215 5400 0
## ---
## 114: 214 4925 1
## 115: 215 4850 1
## 116: 222 5750 0
## 117: 212 5200 1
## 118: 213 5400 0Seleção do modelo
Usaremos o método de “todos os modelos possíveis” para escolher o melhor modelo.
ajusta_todos_modelos <- function(dados){
explicativas <- colnames(dados)[2:ncol(dados)]
ajustes <- list()
for (i in seq_along(explicativas)) {
combinacoes <- combn(explicativas, i)
soma_explicativas <- apply(combinacoes, 2, paste, collapse="+")
ajustes[[i]] <- paste0("bill_length_mm~", soma_explicativas)
}
ajustes <- sapply(unlist(ajustes), as.formula)
modelos <- lapply(ajustes, lm, data=dados)
return(modelos)
}Métricas de avaliação
Avaliaremos as métricas \(R^2\), \(R^2\) ajustado, estatística PRESS e informação de Akaike(AIC).
Escolhemos os seguintes critérios para tornar automática a comparação entre modelos:
Diferença entre o \(R^2\) do modelo mais simples e do mais completo menor que 0.05,
\(R^2\) do modelo mais recente for maior que o modelo antigo em 0.09.
AIC do modelo antigo maior que do modelo novo.
Estatística PRESS do modelo antigo maior que do modelo novo.
Se TODOS os critérios acima são respeitados, o modelo novo é escolhido.
complexidade <- function(modelo){
return(length(coefficients(modelo)))
}
calcula_press <- function(modelo){
influencia <- residuals(modelo)/ (1-lm.influence(modelo)$hat)
return(sum(influencia^2))
}
confere_diferenca_metricas <- function(modelo_novo, modelo_antigo){
# R2 e R2 Ajustado
s_modelo_novo = summary(modelo_novo)
s_modelo_antigo = summary(modelo_antigo)
diff_r2 <- s_modelo_novo$r.squared - s_modelo_antigo$r.squared
diff_r2_ajustado<- s_modelo_novo$adj.r.squared - s_modelo_antigo$adj.r.squared
r2_significativo <- abs(diff_r2) <= 0.03
r2_ajustado_significativo <- abs(diff_r2_ajustado) <= 0.03
r2_grande <- diff_r2 > 0 & abs(diff_r2) >= 0.07
r2_ajustado_grande <- diff_r2_ajustado > 0 & abs(diff_r2) >= 0.07
# PRESS e AIC
aic <- AIC(modelo_novo) < AIC(modelo_antigo)
press <- calcula_press(modelo_novo) < calcula_press(modelo_antigo)
diff_significativa <- r2_significativo & r2_ajustado_significativo & aic & press
if(diff_significativa){
return("significativa")
}
# & aic & press
diff_grande <- r2_grande & r2_ajustado_grande & aic & press
if(diff_grande){
return("grande")
}
return("não significante")
}
escolhe_melhor_modelo <- function(modelos){
melhor_modelo = modelos[[1]]
for(modelo in modelos){
diferenca_metricas <- confere_diferenca_metricas(modelo, melhor_modelo)
if((complexidade(melhor_modelo) > complexidade(modelo)) & diferenca_metricas == "significativa"){
melhor_modelo <- modelo
}
if(diferenca_metricas == "grande"){
melhor_modelo <- modelo
}
}
return(melhor_modelo)
}Todos os modelos possíveis
seleciona_modelo <- function(dados){
melhor_modelo <- dados %>%
ajusta_todos_modelos() %>%
escolhe_melhor_modelo()
return(melhor_modelo)
}
model <- seleciona_modelo(dados)
summary(model)##
## Call:
## FUN(formula = X[[i]], data = ..1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.8897 -1.1924 -0.0294 1.2791 4.7175
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.7099 6.8161 -1.131 0.260350
## bill_depth_mm 1.0715 0.2790 3.841 0.000201 ***
## flipper_length_mm 0.1802 0.0417 4.322 3.31e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.094 on 115 degrees of freedom
## Multiple R-squared: 0.4919, Adjusted R-squared: 0.4831
## F-statistic: 55.66 on 2 and 115 DF, p-value: < 2.2e-16Modelo escolhido
Como os modelos dos artigos de referência são modelos lineares generalizados, nossa técnica de encontrar o melhor modelo linear foi falha. Usamos o método de todos os modelos possíveis, mas chegamos à conclusão de que o modelo não reforçava o que diz a literatura e não respeitava as suposições de um modelo linar. Portanto, testamos outros modelos que usam variáveis que aparecem normalmente nas principais publicações sobre os dados apresentados.
#model <- lm(log(dados$bill_length_mm)~dados$body_mass_g + dados$sex)
model <- lm(log(dados$bill_length_mm)~dados$body_mass_g + dados$sex + dados$flipper_length_mm)
summary(model)##
## Call:
## lm(formula = log(dados$bill_length_mm) ~ dados$body_mass_g +
## dados$sex + dados$flipper_length_mm)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.109955 -0.027128 0.000716 0.025198 0.103189
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.069e+00 1.825e-01 16.819 < 2e-16 ***
## dados$body_mass_g 3.020e-05 1.455e-05 2.076 0.04016 *
## dados$sex -2.871e-02 1.380e-02 -2.081 0.03966 *
## dados$flipper_length_mm 2.992e-03 8.849e-04 3.381 0.00099 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.04296 on 114 degrees of freedom
## Multiple R-squared: 0.5195, Adjusted R-squared: 0.5069
## F-statistic: 41.09 on 3 and 114 DF, p-value: < 2.2e-16Considerando as hipóteses das variáveis explicativas serem ou não significativas, de acordo com os p-valores, ao nível de significância de 5%, temos evidências para rejeitar que as variáveis explicativas e o intercepto são não significativos.
Pela a hipótese da significância do modelo, o p-valor calculado a partir da estatística F mostra que ao nível de significância de 5%, temos evidências para rejeitar que o modelo é não significativo.
O \(R^2\) mostra que conseguimos explicar por volta de 47% dos erros do modelo, o que para este presente trabalho é satisfatório para contribuir na sustentação de nossa hipótese inicial de que há uma relação linear entre as variáveis explicativas e a variável resposta.
Interpretação dos coeficientes do modelo
Intercepto: Caso não houvessem outras variáveis influenciando, o logaritmo do tamanho base do bico seria de 3.069 mm.
Massa corporal em gramas: a cada grama, o logaritmo do comprimento do bico do Gentoo aumenta a uma taxa de 0,00003 mm.
Sexo: Como fêmeas compoem a referência da variável, o coeficiente indica que o logaritmo do bico das fêmeas é menor em comprimento em -0,02871 mm.
Comprimento da nadadeira: A cada milímetro da nadadeira o logaritmo do bico do pinguim Gentoo aumenta em 0,002992 mm.
Análise dos resíduos
Homecedasticidade
Observando o gráfico dos resíduos, é possível notar que os resíduos não aparentam seguir um padrão.
plot(fitted.values(model),rstandard(model),xlab="Valores Ajustados",ylab="Resíduos normalizados")
Normalidade
Observando o histograma, podemos ver que os resíduos aparentam seguir uma distribuição normal.
hist(model$residuals)
qqnorm(model$residuals, pch = 1, frame = FALSE)
qqline(model$residuals, col = "steelblue", lwd = 2)
## Linearidade
Usamos o teste de falta de ajuste para verificar a linearidade entre as variáveis resposta e explicativas. A conclusão foi que não há evidências ao nível de 95% de significância para rejeitar a hipótese nula de que o modelo mais simples é adequado aos dados.
full_model <- lm(log(dados$bill_length_mm) ~., data=dados)
anova(model, full_model)## Analysis of Variance Table
##
## Model 1: log(dados$bill_length_mm) ~ dados$body_mass_g + dados$sex + dados$flipper_length_mm
## Model 2: log(dados$bill_length_mm) ~ species_Adelie + species_Chinstrap +
## species_Gentoo + island_Biscoe + island_Dream + island_Torgersen +
## bill_depth_mm + flipper_length_mm + body_mass_g + sex
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 114 0.21040
## 2 113 0.20539 1 0.0050066 2.7544 0.09976 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Independência
Temos informações sobre a sequência de coleta dos dados, portanto usaremos o teste de Durbin-Watson para verificar independência entre os erros, onde a hipótese nula é de que os erros não são autocorrelacionados entre si. Como o p-valor da estatística de teste é maior que 0.05, ao nível de significância de 95% não há evidências para rejeitar a hipótese nula, e assim, o modelo respeita mais uma suposição.
library(car)
car::durbinWatsonTest(model)## lag Autocorrelation D-W Statistic p-value
## 1 0.0722902 1.845033 0.424
## Alternative hypothesis: rho != 0Detectando valores influentes
Observando a tabela com os DFBETAS é possível perceber que não há nenhum valor que ultrapasse o limiar \(\frac{2}{\sqrt{118}}\), logo não há de se preocupar com valores influentes alterando o ajuste do modelo.
limiar <- 2/sqrt(118)
dfbetas <- dfbeta(model)
dfbetas <- as.data.frame(dfbetas)
names(dfbetas) <- c("intercept", "body_mass_g", "sex", "flipper_length_mm")
head(dfbetas)## intercept body_mass_g sex flipper_length_mm
## 1 0.003494417 -3.853030e-07 4.807606e-05 -6.254392e-06
## 2 0.012420921 -6.754160e-08 -2.362247e-04 -5.579800e-05
## 3 0.019470547 -1.641739e-06 -2.598024e-04 -4.734120e-05
## 4 0.005259875 6.599578e-07 -1.298488e-04 -3.863748e-05
## 5 -0.006939316 -1.009729e-07 4.224338e-04 3.285934e-05
## 6 0.007633885 -2.705210e-07 1.302935e-04 -2.778984e-05Como os fatores de aumento de variância não são muito altos é seguro afirmar que não temos problemas graves de multicolinearidade seguindo essse modelo.
library(DAAG)## Warning: package 'DAAG' was built under R version 4.0.5## Loading required package: lattice##
## Attaching package: 'DAAG'## The following object is masked from 'package:car':
##
## vifDAAG::vif(model)## dados$body_mass_g dados$sex dados$flipper_length_mm
## 3.2975 3.0416 2.1016