Como automatizar a produção de equações utilizando o pacote equatiomatic

Introdução

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

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

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 packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.5     ✓ purrr   0.3.4
## ✓ tibble  3.1.3     ✓ dplyr   1.0.7
## ✓ tidyr   1.1.3     ✓ stringr 1.4.0
## ✓ readr   2.0.1     ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
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

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)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.0611  -0.7613  -0.3783   0.4524   2.2175  
## 
## 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

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

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.