Pequeno estudo de caso fazendo uso de modelos de regressão linear simples bivariada.

Para fazer essa análise, vou usar o dataset Fish quem contém diferentes medidas de peixes registrados em um mercado.

Meu objetivo aqui vai ser escolher uma variável de resposta e a melhor variável que explique seu comportamento.

Primeiro vou dar uma olhada no dataset, já carregando alguns pacotes que vou utilizar:

library(rmarkdown)
library(tidyverse) 
## ── Attaching packages ─────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.2     ✓ purrr   0.3.4
## ✓ tibble  3.0.3     ✓ dplyr   1.0.2
## ✓ tidyr   1.1.2     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.5.0
## ── Conflicts ────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(broom)
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
## The following object is masked from 'package:purrr':
## 
##     some
library(ggExtra)


paged_table( 
fish <- read_csv("https://www.dropbox.com/s/n45vml0ayoq0omx/fish.csv?dl=1")
)
## Parsed with column specification:
## cols(
##   Species = col_character(),
##   Weight = col_double(),
##   Length1 = col_double(),
##   Length2 = col_double(),
##   Length3 = col_double(),
##   Height = col_double(),
##   Width = col_double()
## )

Os nomes de coluna não são muito intuitivos, então vou descrever cada um:

  1. Species: nome da espécie mensurada.
  2. Weight: Peso do peixe em gramas.
  3. Length1: comprimento vertical em cm.
  4. Length2: comprimento diagonal em cm.
  5. Length3: comprimento cruzado em cm.
  6. Height: altura em cm.
  7. Width: largura diagonal em cm.

Vou definir como minha variável de resposta o peso do peixe, esse é o valor que quero prever em função das demais variáveis explicativas. Para dar início à análise, entao, vou plotar a correlacao de minha varável de resposta com todas as demais variáveis explanatórias e calcular o \(R^2\) de cada caso.

Para plotar tudo de uma vez preciso alterar um pouco a tabela:

fish %>% # alterando a disposição de colunas para a plotagem
  pivot_longer(cols=names(fish)[3:length(names(fish))],names_to = "explanatory") %>% 
  # plotagem
  ggplot(aes(y=Weight,x=value)) + geom_point() + geom_smooth(method="lm",linetype=2) +
  theme_bw() +
  ggtitle("Peso vs. altura, medidas de comprimento e largura.") +
  ylab("Peso - variável de resposta") +
  xlab("Variáveis explanatórias, em escala livre.") +
  facet_wrap(~explanatory, scales="free")
## `geom_smooth()` using formula 'y ~ x'

Dá para perceber que a correlação do peso com as demais medidas é positiva em todos casos, em alguns mais e outros menos, pela inclinação das retas de regressão (tracejadas azuis). Entretanto, fica difícil entender a robustez dos modelos (na forma de \(R^2\)) visualmente. Portanto, vou computar esses valores em uma tabela e ordenar em ordem decrescente:

Essa tabela me dá os valores de \(R^2\) para cada modelo simples de peso e cada uma das outras variáveis independentemente (não sendo covariáveis).

As três medidas de comprimento são as que melhor descrevem o comportamento do peso. A altura é a medida menos robusta e a largura fica no meio do caminho. Os valores altos dos modelos que levam em consideração os diferentes tipos de comprimento apontam para uma possível colinearidade entre essas variáveis. Isso significa que, em um modelo multivariado, essas medidas provavelmente seriam redundantes.

Vou desenvolver isso um pouco mais a frente, mas por hora, como vou fazer uma regressão simples, vou usar a variável Length3 como explanatória, dado que ela resultou no maior valor de \(R^2\).

Abaixo, os coeficientes do modelo:

lm(Weight~Length3,data=fish)
## 
## Call:
## lm(formula = Weight ~ Length3, data = fish)
## 
## Coefficients:
## (Intercept)      Length3  
##     -490.40        28.46

Segundo esse modelo, existe um incremento de cerca de 28 gramas no peso de um peixe conforme a largura cruzada do peixe aumenta em uma unidade, ou seja, 1 centímetro.

Abaixo, vou plotar o gráfico da regressão, já atribuindo cores diferentes aos pontos para as diferentes espécies da amostra:

library(paletteer) # pacote com paletas de cores

fish %>% ggplot(aes(y=Weight,x=Length3)) + 
  geom_point(alpha=0.7,aes(color=Species)) + 
  geom_smooth(method="lm", linetype=2) +
  theme_bw() +
  ylab("Peso (g)") +
  xlab("Comprimento cruzado (cm)") +
  labs(color='Espécie') +
  ggtitle(" Relação entre peso e comprimento cruzado,\n com cores referentes às espécies das observações.") +
  scale_color_paletteer_d("pals::alphabet") +
  theme(plot.title = element_text(size = 10))
## `geom_smooth()` using formula 'y ~ x'

Um fator interessante de se notar é a provável alta influência dos pontos associados à espécie Pike (em roxo), na inclinação e erro associado da reta, devido ao alto valor tanto de peso quanto de comprimento cruzado de alguns de seus pontos. Essa espécie provavelmente tem maiores valores, em média, de distância de cook para suas observações, o que vai influenciar no comportamento da reta. Podemos checar isso sumarizando o modelo e esses valores, por espécie:

mod <- lm(Weight~Length3,data=fish)

augment(mod,data=fish) %>% 
  group_by(Species) %>% summarise(mean_cooksd=mean(.cooksd)) %>% 
  arrange(desc(mean_cooksd))
## `summarise()` ungrouping output (override with `.groups` argument)

Como previsto, a espécie Pike de fato tem a maior média para os valores de distância de cook associados seus pontos, seguida pela espécie Smelt (verde escuro), que no gráfico tem uma concentração de pontos próximos de zero. Esses pontos, em ambas espécies, se distanciam bastante da reta de regressão, explicando o comportamento das médias de distância de cooks associadas a elas.

Vou checar mais a fundo a influência que estas espécies tem na reta e no valor de \(R^2\).

Lembrando que o modelo tem \(R^2\) de 0.85, vou verificar qual valor que \(R^2\) assume retirando essas espécies dele:

fish %>% filter(!(Species %in% c("Pike","Smelt"))) %>% 
  lm(Weight~Length3,data=.) %>% summary() %>% {.$r.squared}
## [1] 0.9109848

Retirando as espécies, a robustez do modelo aumenta significativamente. A inclinação da reta também se altera, conforme podemos ver na mudança do coeficiente associado à variável explanatória de comprimento cruzado:

fish %>% filter(!(Species %in% c("Pike","Smelt"))) %>% 
  lm(Weight~Length3,data=.)
## 
## Call:
## lm(formula = Weight ~ Length3, data = .)
## 
## Coefficients:
## (Intercept)      Length3  
##     -656.48        34.14

Para enxergar essa alteração melhor ainda, podemos sobrepor as retas de cada modelo:

df_alt <- fish %>% filter(!(Species %in% c("Pike","Smelt")))
mod2 <- lm(formula = Weight ~ Length3, data = df_alt)


df_alt %>% ggplot(aes(y=Weight,x=Length3)) + geom_point(alpha=0) +
  geom_abline(intercept = -656.48, slope = 34.14, linetype=2,color="blue") +
  geom_abline(intercept = -490.40, slope = 28.46, linetype=2,color="red") +
  ylab("Peso (g)") +
  xlab("Comprimento cruzado (cm)") +
  ggtitle(" Gráfico de pontos com reta de regressão (em azul) sem as espécies Pike e Smelt.\n A reta em vermelho é derivada do modelo com todas espécies.") +
  theme_bw() +
  theme(plot.title = element_text(size = 10))

Visualizando dessa forma fica mais claro o comportamento da nova reta em comparação com a reta associada ao primeiro modelo (tracejada vermelha). Essa inclinação indica uma correlação mais forte entre as variáveis na ausência das espécies Pike e Smelt, fato reforçado quando se calcula os índices de pearson para os dois casos:

tibble(df=c("df_original","df_sem_especies_ruido"),
       indice_pearson=c(
         cor(fish$Weight,fish$Length3),
         cor(df_alt$Weight,df_alt$Length3)
       ))

Finalmente, sem aprofundar muito, existe outra maneira prática de verificar o quanto os pontos das espécies podem estar influenciando na qualidade da regressão, que é plotando as curvas de regressão por espécie:

fish %>% ggplot(aes(y=Weight,x=Length3)) + geom_point() + geom_smooth(method="lm", linetype=2) +
  ylab("Peso (g)") +
  xlab("Comprimento cruzado (cm)") +
  coord_cartesian(ylim = c(0,1500)) + 
  ggtitle(" Peso vs. comprimento cruzado, por espécie da amostra.") +
  theme_bw() +
  theme(plot.title = element_text(size = 10)) +
  facet_wrap(~Species)
## `geom_smooth()` using formula 'y ~ x'

A espécie Pike é a única com valores de peso maiores que 1500. Já a espécie Smelt tem apenas valores próximos de zero, reforçando um pouco da influência destas na qualidade da regressão. Essa informação é redundante com o gráfico onde os pontos foram pintados de acordo com as espécies de referência, mas é uma maneira distinta de enxergar essa questão.

Esse tipo de fator levanta a questão: quais espécies escolher para estabelecer um modelo geral? Uma maneira de fazer isso seria calculando o \(R^2\) associado ao modelo escolhido por espécie, mas também indicando o número de observações por espécie existe nos dados:

# funcao para extrair os valores de R2 dos modelos
rsquared_mod <- function(x,data) {
  return(
    summary( 
    lm(paste0("Weight ~", x), data = data)
    )[["r.squared"]]
  )
}

# lista de vars explanatorias
explanatory <- names(fish)[3:length(names(fish))]

# tabela com valores de r2 dos modelos
tabela_r <- tibble(variables=explanatory,r_squared=unlist(map(explanatory,rsquared_mod,data=fish)))

# tabela ordenada

tabela_r %>% arrange(desc(r_squared))
miau <- function(sp){
  data <- fish %>% filter(Species==sp)
  return( 
  tabela_r <- tibble(Species=rep(sp,5),
                     variables=explanatory,
                     r_squared=unlist(map(explanatory,rsquared_mod,data=data)))
  )
}

auau <- map(unique(fish$Species),miau)

moo <- do.call("rbind",auau) 

left_join( 
moo %>% filter(variables=="Length3") %>% arrange(desc(r_squared)),
fish %>% group_by(Species) %>% summarise(n=n()),
by="Species"
)
## `summarise()` ungrouping output (override with `.groups` argument)

Quando analisadas em separado, a espécie Pike (que apresentava altos valores de distância cook, em média), não apresenta um \(R^2\) tão baixo. Mas é importante observar que o número de observações da maior parte das espécies da amostra é bem baixo (menor que 20). Apenas as espécies Perch e Bream tem mais de 20 observações. A espécie Whitefish, que contém o maior valor de \(R^2\) em seu modelo, contém apenas 6 observações. Se a ideia fossse utilizar um modelo generalista, seria razoável manter apenas espécies número alto de observações (20 talvez) ou aumentar o número de medidas das espécies subamostradas.

1 Colinearidade de variáveis

Como observado, todas medidas de comprimento apresentaram altos valores de \(R^2\) quando correlacionados com o peso. Recapitulando a tabela:

# funcao para extrair os valores de R2 dos modelos
rsquared_mod <- function(x) {
  return(
    summary( 
    lm(paste0("Weight ~", x), data = fish)
    )[["r.squared"]]
  )
}

# lista de vars explanatorias
explanatory <- names(fish)[3:length(names(fish))]

# tabela com valores de r2 dos modelos
tabela_r <- tibble(variables=explanatory,r_squared=unlist(map(explanatory,rsquared_mod)))

# tabela ordenada

tabela_r %>% arrange(desc(r_squared))

Esse comportamento pode ser um indicativo de colinearidade, o que faria sentido dada a natureza similar dessas medidas, em termos de proporção. Isso pode ser conferido computando os fatores de inflação de variância (vif) para cada variável, em um modelo que leva em conta todas variáveis:

library(car)

fish %>% 
  select(-Species) %>%
  lm(Weight~.,data=.) %>%
  vif()
##    Length1    Length2    Length3     Height      Width 
## 1681.49649 2084.25783  422.46825   14.57009   12.27536

Os altos valores de VIF para as medidas de comprimento indicam que essas variáveis são, de fato, redundantes para o modelo, não fazendo diferença qual é utilizada: ambas tem boa perfomance para explicar os valores de peso.

Esse comportamento pode ser visualizado se plotarmos, por exemplo, o peso em função de dois comprimentos quaisquer:

library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
miau <- fish

auau <- lm(Weight~Length3+Length1, data=miau)

p <- plot_ly(data = miau, z = ~Weight, x = ~Length3, y = ~Length1, 
             opacity = 0.6, colorbar = list(title = "Peso previsto")) %>%
  add_markers()

cf.mod <- coef(auau)


x1.seq <- seq(min(miau$Length3),max(miau$Length3),length.out=1000)
x2.seq <- seq(min(miau$Length1),max(miau$Length1),length.out=1000)

z <- t(outer(x1.seq, x2.seq, function(x,y) cf.mod[1]+cf.mod[2]*x+cf.mod[3]*y))

p %>% # da para eu criar uma escala personalizada
  add_surface(x = ~x1.seq, y = ~x2.seq, z = ~z, 
              showscale = TRUE,colorscale="Viridis") %>% 
  layout(scene = list(xaxis = list(title = "Comprimento cruzado (cm)"), 
                      yaxis = list(title = "Comprimento vertical (cm)"),
                      zaxis = list(title = "Peso (g)")))
## Warning: `arrange_()` is deprecated as of dplyr 0.7.0.
## Please use `arrange()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
## Warning: 'scatter3d' objects don't have these attributes: 'colorbar'
## Valid attributes include:
## 'type', 'visible', 'showlegend', 'legendgroup', 'opacity', 'name', 'uid', 'ids', 'customdata', 'meta', 'hoverlabel', 'stream', 'transforms', 'uirevision', 'x', 'y', 'z', 'text', 'texttemplate', 'hovertext', 'hovertemplate', 'mode', 'surfaceaxis', 'surfacecolor', 'projection', 'connectgaps', 'line', 'marker', 'textposition', 'textfont', 'hoverinfo', 'error_x', 'error_y', 'error_z', 'xcalendar', 'ycalendar', 'zcalendar', 'scene', 'idssrc', 'customdatasrc', 'metasrc', 'xsrc', 'ysrc', 'zsrc', 'textsrc', 'texttemplatesrc', 'hovertextsrc', 'hovertemplatesrc', 'textpositionsrc', 'hoverinfosrc', 'key', 'set', 'frame', 'transforms', '_isNestedKey', '_isSimpleKey', '_isGraticule', '_bbox'

Os pontos nao se espalham quase nada. Mas nao sei se vale fazer esse grafico aqui. Vale, o fato dos pontos formarem quase uma linha mostra que x e y crescem quase do mesmo jeito, linearmente, reforcando a colinearidade.

Colocar titulo nesse grafico.

2 Conclusões

O comprimento, em geral, é um bom preditor.

Alguns peixes causam ruido.

As medidas de comprimento apresentam colinearidade. Isso indica que não vale um modelo multivariado aqui.