O objetivo deste projeto é treinar a aplicação de modelos de regressão linear e criar gráficos de dispersão que apresentem paramêtros relacionados a qualidade do modelo. Para isso, será replicada algumas análises realizadas pelo Ladislas Nalborczyk em seu texto “Using R to make sense of the generalised linear model”
Para criação das figuras e modelos foi utilizada a base de dados Howell1 do pacote Rethinking (McElreath, 2016), que contém dados sobre 544 indivíduos, incluindo altura (centímetros), peso (quilogramas), idade (anos) e sexo (0 indicando feminino e 1 indicando masculino).
# Lendo os pacotes --------------------------------------------------------
library(rethinking)
library(ggplot2)
library(ggExtra)Após filtrar a base de dados apenas com os dados pertencentes aos indíviduos com idade igual ou maior que 18 anos (>=18), foi criado um gráfico de dispersão utilizando o ggplot2 (i.e., geom_point e geom_smooth). Um objeto foi criado com o gráfico (p). Em seguida, foi utilizada a função ggMarginal do pacote ggExtra para criar um gráfico de dispersão com os histogramas em cada um dos eixos.
# Abrindo, criando um elemento (d) e visualizando a base ------------------
data("Howell1")
# Filtrando a idade (>18yr) -----------------------------------------------
d <- Howell1 |>
dplyr::filter(age >= 18)# Criando um grafico de dispersão -----------------------------------------
p <- d |>
ggplot(mapping = aes(x = weight, y = height)) +
geom_point(pch = 21, color = "white", fill = "black", size = 2, alpha = 0.8) +
geom_smooth(method = "lm", color = "black", se = TRUE)+
theme_bw(base_size = 12)
# Usando ggExtra para fazer grafico de dispersão com histograma -----------
ggMarginal(p, type = "histogram")
Uma rápida inspeção visual do conjunto de dados revela uma relação positiva entre altura e peso. A linha da regressão traçada no grafico acima corresponde ao seguinte modelo linear assumindo uma distribuição normal:
Peso ~ Normal (µi, σ)
µi = α + β ∙ 43
que pode ser excrita com a sintaxe:
lm(formula, data, subset, weights, na.action, method = “qr”, model = TRUE, x = FALSE, y = FALSE, qr = TRUE, singular.ok = TRUE, contrasts = NULL, offset, …)
Assim, o código para esta análise e o resultado serão:
mod1 <- lm(height ~ weight, data = d)
mod1
#>
#> Call:
#> lm(formula = height ~ weight, data = d)
#>
#> Coefficients:
#> (Intercept) weight
#> 113.879 0.905O intercept (i.e., 113.879) representa a altura prevista quando o peso está em 0 (o que não faz muito sentido neste caso), enquanto o slope (i.e., 0.905) representa a mudança na altura quando o peso aumenta em uma unidade (e.g., um quilograma).
A principal vantagem de um modelo de regressão é a possibilidade de fazer previsões como, por exemplo, baseado no modelo apresentado acima, qual é a altura esperada para um indivíduo com 43 kg. No R isto é realizado criando um geom_point e marcando a linha referente ao peso que você quer fazer a predição da estatura e usando a função predict, que é igual a α + β * 43.
Vamos testar!
# Fazendo gráfico com linhas referentes ao peso que queremos preve --------
wght <- 43
d |>
ggplot(mapping = aes(x = weight, y = height)) +
geom_line(mapping = aes(y = predict(mod1)), size = 1) +
geom_point(size = 2, alpha = .2) +
geom_segment(
x = wght, xend = wght,
y = 0,
yend = predict(mod1, newdata = data.frame(weight = wght)),
linetype = 2, lwd = 0.5) +
geom_segment(
x = 0, xend = wght,
y = predict(mod1, newdata = data.frame(weight = wght)),
yend = predict(mod1, newdata = data.frame(weight = wght) ),
linetype = 2, lwd = 0.5
) +
theme_bw(base_size = 12)
Conforme o autor do texto ensina, implementar a função predict
manualmente é fácil. Vamos testar para ver como funciona.
# Aplicando a função predict ----------------------------------------------
d <- d |>
dplyr::mutate(
pred_mod1 = predict(mod1),
pred_mod1_2 = coef(mod1)[1] + coef(mod1)[2] * weight
)
head(d)
#> height weight age male pred_mod1 pred_mod1_2
#> 1 151.765 47.82561 63 1 157.1630 157.1630
#> 2 139.700 36.48581 63 0 146.9001 146.9001
#> 3 136.525 31.86484 65 0 142.7180 142.7180
#> 4 156.845 53.04191 41 1 161.8839 161.8839
#> 5 145.415 41.27687 51 0 151.2362 151.2362
#> 6 163.830 62.99259 35 1 170.8895 170.8895Como dá para perceber, a função predict está simplesmente recuperando parâmetros do modelo ajustado (neste caso, pelo intercept e slope) para fazer previsões sobre a variável de resultado, dados alguns valores do (s) preditor (es). Este modelo ajustado corresponde à segunda linha do código descrito acima (i.e.,pred_mod1_2).
É vbem comum querer predizer desfechos binários (e.g., sim e não, tem e não tem). E o modelo mais adequado para fazer isso é por meio de modelos de regressão logística.
Para realizar uma regressão logística no R utiliza-se a função glm, onde o argumento family é usado para especificar a probabilidade do modelo e a função de ligação. Para isso, utiliza-se a função glm cuja sintaxe é:
glm(formula, family = gaussian, data, weights, subset, na.action, start = NULL, etastart, mustart, offset, control = list(…), model = TRUE, method = “glm.fit”, x = FALSE, y = TRUE, singular.ok = TRUE, contrasts = NULL, …)
Para testar a função glm, segue abaixo o código sendo aplicado na base de dados que estamos usando para predizer o sexo pela estatura e o resultado.
mod2 <-
glm(male ~ height,
data = d,
family = binomial(link = "logit"))
mod2 <-
glm(male ~ height,
data = d,
family = binomial(link = "logit"))
# Aplicado o modelo para predizer o sexo pela estatura usando a predict----
# e fitted ----------------------------------------------------------------
d <-
d |>
dplyr::mutate(
pred_mod2 = predict(mod2),
fitted_mod2 = fitted(mod2)
)
# Selecionando e vendo o resultado ----------------------------------------
d |>
dplyr::select(height, weight, male, pred_mod2, fitted_mod2) |>
head()
#> height weight male pred_mod2 fitted_mod2
#> 1 151.765 47.82561 1 -1.1910139 0.233077651
#> 2 139.700 36.48581 0 -5.3387584 0.004778882
#> 3 136.525 31.86484 0 -6.4302701 0.001609421
#> 4 156.845 53.04191 1 0.5554049 0.635388646
#> 5 145.415 41.27687 0 -3.3740373 0.033116790
#> 6 163.830 62.99259 1 2.9567306 0.950580634Conforme o autor destaca, diferente do modelo linear utilizado para variável continuas, neste modelo o predic e o fitted apresentou resultados bastante diferentes. O autor explica que a saída das funções predict e fitted são diferentes quando usamos um GLM porque a função predict retorna previsões do modelo na escala do preditor linear (neste exempo, na escala log-odds), enquanto a função fitted retorna previsões na escala da resposta. De maneira generosa, o autor disponibiliza uma função criada por ele (logit_dotplot.R) para realizar um gráfico para uma regressão logística. Abaixo, segue a chamada da função mais a aplicação dela nesta base de dados.
source("logit_dotplot.R")
logit_dotplot(d$height, d$male, xlab = "height", ylab = "p(male)")
O autor deste texto comenta que um erro corriqueiro sobre as suposições do modelo linear (Gaussiano) é que o desfecho deve apresentar distribuição normal. Contudo, essa suposição diz respeito à distribuição do desfecho em torno de seu valor previsto (ou seja, a distribuição dos erros). Mas os erros são as diferenças não observadas (e não observáveis) entre o valor teórico previsto e os resultados observados e, consequentemente, não temos acesso a ele. Em vez disso, podemos trabalhar com os resíduos, que podem ser vistos como uma estimativa (da amostra) dos erros De uma forma simplicada, os erros pertencem ao processo de geração de dados, enquanto os resíduos são a diferença entre a estimativa do modelo e os resultados observados. No R, podemos obtê-los facilmente usando a função de residuals ou subtraindo a cada resultado observado o previsto.
# Calculando os residuos --------------------------------------------------
d <- d |>
dplyr::mutate(
res1 = residuals(mod1),
res2 = height - pred_mod1
)
d |>
dplyr::select(height, weight, male, pred_mod1, res1, res2) |>
head()
#> height weight male pred_mod1 res1 res2
#> 1 151.765 47.82561 1 157.1630 -5.397960 -5.397960
#> 2 139.700 36.48581 0 146.9001 -7.200111 -7.200111
#> 3 136.525 31.86484 0 142.7180 -6.193000 -6.193000
#> 4 156.845 53.04191 1 161.8839 -5.038870 -5.038870
#> 5 145.415 41.27687 0 151.2362 -5.821164 -5.821164
#> 6 163.830 62.99259 1 170.8895 -7.059520 -7.059520Abaixo, é apresetado um gráfico com os residuos onde a tranparência e o tamanho dos pontos são dependentes da distância até o valor previsto (de modo que resíduos maiores apareçam como maiores e menos transparentes). Essa distância também é representada pelo comprimento das linhas verticais.
d |>
dplyr::slice_sample(prop = 0.5) |>
ggplot(mapping = aes(x = weight,
y = height)) +
geom_line(mapping = aes(y = pred_mod1), size = 1) +
geom_point(mapping = aes(alpha = abs(res1),
size = abs(res1))) +
guides(alpha = "none", size = "none") +
geom_segment(mapping = aes(xend = weight,
yend = pred_mod1,
alpha = abs(res1))) +
theme_bw(base_size = 12)
Agora, se pegarmos os resíduos e traçar sua distribuição, podemos
verificar a suposição de que a distribuição dos resíduos apresenta
distribuição normal. Como podemos observar na figura abaixo, sim, a
distribuição dos resíduos apresenta distribuição normal.
# Distribuição dos resíduos -----------------------------------------------
d |>
ggplot(mapping = aes(x = res1)) +
geom_histogram(mapping = aes(y = ..density..), bins = 20,
alpha = 0.6) +
geom_line(aes(y = dnorm(res1, mean = 0, sd = sd(res1))),
size = 1) +
guides(fill = "none") +
theme_bw(base_size = 12) +
labs(x = "Residuals", y = "Density")
Após replicar algumas das análises e gráficos ensinados pelo autor n post replicado foi possível praticar alguns conceitos relacionados ao modelo linear generalizado utilizando funções do R para prever e ajustar o modelo de regressão e observar a distribuição dos resíduos que deve ser apresentar distrivuição Gaugassiana.