Introdução Link para o cabeçalho
Recentemente tomei conhecimento do pacote equatiomatic do R. Ele é capaz de reportar automaticamente as equações de modelos lineares ajustados a dados.
O leitor pode se perguntar se funções como summary já fazem isso. A diferença é que o pacote equatiomatic é capaz de tomar o resultado de funções como lm, glm e outras e transformá-los em equações já formatadas em LaTeX. O objetivo deste post é ilustrar como isso pode ser feito.
Instalação Link para o cabeçalho
Antes de mais nada é preciso instalar o pacote. Isso pode ser feito de duas maneiras. Caso o usuário queira usar a versão oficial mais recente do pacote, basta rodar
install.packages("equatiomatic")
para resolver os seus problemas. Por outro lado, aqueles que como eu desejam ter acesso à versão mais recente, com as mais novas funcionalidades, devem rodar
devtools::install_github("datalorax/equatiomatic")
para isso.
Regressão Linear Link para o cabeçalho
Vou começar a ilustrar o uso do pacote com uma regressão linear simples. O conjunto de dados mpg possui diversas 11 informações a respeito de 234 carros de passeio. Em particular, ele posssui informações sobre consumo desses carros na estrada (hwy) e na cidade (cty).
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.0 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
theme_set(theme_bw())
ggplot(mpg, aes(x = cty, y = hwy)) +
geom_point()
Se quisermos prever através de uma regressão linear simples qual será o consumo desses carros na estrada a partir do consumo na cidade, podemos ajustar o seguinte modelo:
ajuste_lm <- lm(hwy ~ cty, data = mpg)
summary(ajuste_lm)
##
## Call:
## lm(formula = hwy ~ cty, data = mpg)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.3408 -1.2790 0.0214 1.0338 4.0461
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.89204 0.46895 1.902 0.0584 .
## cty 1.33746 0.02697 49.585 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.752 on 232 degrees of freedom
## Multiple R-squared: 0.9138, Adjusted R-squared: 0.9134
## F-statistic: 2459 on 1 and 232 DF, p-value: < 2.2e-16
Embora já seja fácil ver quais são os coeficientes do modelo acima, o comando equatiomatic::extract_eq deixa isso mais explícito:
library(equatiomatic)
extract_eq(ajuste_lm)
$$ \operatorname{hwy} = \alpha + \beta_{1}(\operatorname{cty}) + \epsilon $$
Ainda é possível utilizar o argumento use_coefs = TRUE para substituir os coeficientes genéricos mostrados acima por seus valores estimados:
extract_eq(ajuste_lm, use_coefs = TRUE)
$$ \operatorname{\widehat{hwy}} = 0.89 + 1.34(\operatorname{cty}) $$
Regressão Logística Link para o cabeçalho
Em um post de 2018, intitulado Teria sido possível evitar o desastre da Challenger?, eu analiso dados referentes ao acidente do ônibus espacial Challenger. Vou retornar àquele exemplo para mostrar como o pacote equatiomatic é capaz de
temp = c(66, 70, 69, 68, 67, 72, 73, 70, 57, 63, 70, 78, 67, 53, 67, 75, 70, 81, 76, 79, 75, 76, 58)
ring = c(0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1)
# converter Fahrenheit para Celsius
temp <- (temp-32)*5/9
dados <- data.frame(temperatura = temp, falha = ring)
ajuste_glm <- glm(falha ~ temperatura,
data = dados,
family = binomial)
summary(ajuste_glm)
##
## Call:
## glm(formula = falha ~ temperatura, family = binomial, data = dados)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 7.6137 3.9334 1.936 0.0529 .
## temperatura -0.4179 0.1948 -2.145 0.0320 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 28.267 on 22 degrees of freedom
## Residual deviance: 20.315 on 21 degrees of freedom
## AIC: 24.315
##
## Number of Fisher Scoring iterations: 5
A equação resultante do ajuste é essa:
extract_eq(ajuste_glm)
$$ \log\left[ \frac { P( \operatorname{falha} = \operatorname{1} ) }{ 1 - P( \operatorname{falha} = \operatorname{1} ) } \right] = \alpha + \beta_{1}(\operatorname{temperatura}) $$
Agora o resultado com os coeficientes estimados:
extract_eq(ajuste_glm, use_coefs = TRUE)
$$ \log\left[ \frac { \widehat{P( \operatorname{falha} = \operatorname{1} )} }{ 1 - \widehat{P( \operatorname{falha} = \operatorname{1} )} } \right] = 7.61 - 0.42(\operatorname{temperatura}) $$
Modelos Mistos Link para o cabeçalho
E não é necessário parar por nos modelos lineares generalizados. É possível, com a versão de desenvolvimento do pacote, explicitar equações de modelos mistos.
Iremos analisar aqui os resultados de um experimento de tempo de resposta. Tome, por exemplo, o conjunto de dados sleepstudy do pacote lme4. 18 sujeitos foram submetidos a privação de sono e, em seguida, foi medido o tempo de reação médio (em milissegundo) destes sujeitos a uma série de estímulos. O dia 0 significa que o sujeito teve uma noite de sono normal. Cada dia de 1 a 9 significa que o sujeito dormiu no máximo três horas durante a noite anterior aos testes.
library(lme4)
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
ggplot(sleepstudy, aes(x = Days, y = Reaction)) +
geom_point() +
facet_wrap(~ Subject) +
scale_x_continuous(breaks = seq(0, 9, 3))
Vamos ajustar um modelo misto a esses dados, utilizando Days como efeito fixo e Subject como efeito aleatório. Os resultados obtidos estão a seguir.
ajuste_lmer <- lmer(Reaction ~ Days + (1|Subject),
data = sleepstudy)
summary(ajuste_lmer)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Reaction ~ Days + (1 | Subject)
## Data: sleepstudy
##
## REML criterion at convergence: 1786.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2257 -0.5529 0.0109 0.5188 4.2506
##
## Random effects:
## Groups Name Variance Std.Dev.
## Subject (Intercept) 1378.2 37.12
## Residual 960.5 30.99
## Number of obs: 180, groups: Subject, 18
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 251.4051 9.7467 25.79
## Days 10.4673 0.8042 13.02
##
## Correlation of Fixed Effects:
## (Intr)
## Days -0.371
extract_eq(ajuste_lmer)
$$
\begin{aligned} \operatorname{Reaction}_{i} &\sim N \left(\alpha_{j[i]} + \beta_{1}(\operatorname{Days}), \sigma^2 \right) \\ \alpha_{j} &\sim N \left(\mu_{\alpha_{j}}, \sigma^2_{\alpha_{j}} \right) \text{, for Subject j = 1,} \dots \text{,J} \end{aligned}
$$
extract_eq(ajuste_lmer, use_coefs = TRUE)
$$
\begin{aligned} \operatorname{\widehat{Reaction}}_{i} &\sim N \left(251.41_{\alpha_{j[i]}} + 10.47_{\beta_{1}}(\operatorname{Days}), \sigma^2 \right) \\ \alpha_{j} &\sim N \left(0, 37.12 \right) \text{, for Subject j = 1,} \dots \text{,J} \end{aligned}
$$
Conclusão Link para o cabeçalho
O pacote equatiomatic veio cobrir uma lacuna muito importante na pesquisa científica. Ele é uma excelente ajuda para estudantes que estão recém aprendendo como modelagem paramétrica funciona. É possível ver o output tradicional da função summary e compará-lo com o resultado apresentado pelo pacote equatiomatic.
Além disso, pesquisadores experientes podem agilizar seu fluxo de trabalho, já criando automaticamente as equações dos modelos ajustados em R Markdown.
Note que estes modelos que listei aqui não são os únivos que o pacote equatiomatic é capaz de suportar. A vignette do pacote exibe muito mais aplicações possíveis para ele.