Equatiomatic : un package bien sympathique !

 

Si vous employez des modèles linéaires (ANOVA, régressions linéaires multiples, glm, modèles mixtes, etc…) pour analyser vos données, vous avez sans doute déjà eu besoin, notamment pour une publication, d’écrire l’équation du modèle obtenu.

Restituer l’équation du modèle obtenu, c’est fastidieux, et parfois même un peu casse-tête : je pense au ANOVA à deux facteurs avec interaction, par exemple.

Mais ça, c’était avant le package `equatiomatic` ! Ce package, et notamment sa fonction `extract_eq`, extrait automatiquement, l’équation en Latex, d’un modèle.

Et si vous l’utilisez dans un document rmarkdown, celui-ci va automatiquement interpréter l’équation en Latex, pour vous l’écrire directement dans le document rendu !

Dans cet article, je vous montre en pas, comment l’utiliser pour obtenir, ultra-facilement l’équation de vos modèles.

Le package equatiomatic

Installation

Ce package est disponible sur CRAN, il est donc facile de le télécharger, SAUF que il n’est pas disponible pour la dernière version de R (la 4.0.2, en tout cas au moment de la rédaction decet article).

Pour ceux qui sont à jour de la dernière version de R, il faudra donc le télécharger depuis Github. Et pour cela, il faut d’abord télécharger le package `remotes` , sauf évidemment si vous l’avez déjà !

install.packages("remotes")
library(remotes)
remotes::install_github("datalorax/equatiomatic")
library(equatiomatic)

Remarque : si vous obtenez un message vous disant que d’autres packages peuvent être mis à jour, vous pouvez dire oui pour All.

Fonctionnement

Ajustement du modèle

Je réalise ici une régression linéaire simple de la longueur des sépales en fonction de leur largeur :

mod1 <- lm(Sepal.Length~Sepal.Width, data=iris)

Obtention de l’équation

Pour obtenir l’équation générale du modèle de régression linéaire simple, en Latex, il suffit simplement d’exécuter la commande `extract_eq(mod1)` dans la console :

extraction de l'équation en latex

Et si vous incluez cette commande au sein d’un chunk d’un document rmarkdown, comme ceci :

euqatiomatic dans un document rmarkdown

Vous obtiendrez directement l’équation dans le document rendu (word, dans l’exemple ci-dessous):

Rendu de l'équation

C’est déjà pas mal, mais ce qui nous intéresse vraiment, c’est d’obtenir l’équation, avec les coefficients estimés. Pour cela, il suffit simplement d’ajouter l’argument use_coefs = TRUE dans la fonction extract_eq(), comme cela :

extraction des coefficients estimés

Voici le rendu pdf :

Rendu pdf d'equatiomatic

C’est chouette non ?

Nous pouvons vérifier que les coefficients de l’équation sont bien ceux estimés :

summary(mod1)$coefficients

##               Estimate Std. Error   t value     Pr(>|t|)
## (Intercept)  6.5262226  0.4788963 13.627631 6.469702e-28
## Sepal.Width -0.2233611  0.1550809 -1.440287 1.518983e-01

Bon… alors, pour être complètement honnête, tout ne semble pas parfait : avec le rendu word, je n’obtiens pas directement l’équation, mais son code Latex :

rendu word d'equatiomatic

Ainsi que ce message d’erreur :

Ce n’est quand même pas très grave, puisque vous pouvez obtenir cette équation avec les rendus pdf et html.

 

 

Gestion de l’apparence

Parfois, l‘équation est très longue, et dans cette situation, on peut l’étendre sur plusieurs lignes en employant l’argument wrap=TRUE.
Par exemple avec un modèle de covariance :

mod2 <- lm(Sepal.Length~Sepal.Width*Species, data=iris)
extract_eq(mod2, use_coefs=TRUE, wrap = TRUE)

rendu de l'équation sur plusieurs lignes

summary(mod2)$coefficients

##                                Estimate Std. Error   t value     Pr(>|t|)
## (Intercept)                   2.6390012  0.5715053 4.6176320 8.526118e-06
## Sepal.Width                   0.6904897  0.1657268 4.1664328 5.311035e-05
## Speciesversicolor             0.9007335  0.7987512 1.1276772 2.613317e-01
## Speciesvirginica              1.2678352  0.8161509 1.5534324 1.225149e-01
## Sepal.Width:Speciesversicolor 0.1745880  0.2598919 0.6717717 5.028054e-01
## Sepal.Width:Speciesvirginica  0.2110448  0.2557557 0.8251811 4.106337e-01

Les coefficients sont bien identiques.

Remarque :  Vous trouverez d’autres options de gestion de l’apparence de l’équation ici.

Test avec différents modèles

Les codes présentés sont à intégrer dans un chunk d’un document rmarkdown.

 

ANOVA à un facteur

mod3 <- lm(Sepal.Length~Species, data=iris)
extract_eq(mod3, use_coefs=TRUE, wrap = TRUE)

equation de l'anova

summary(mod3)$coefficients

##                   Estimate Std. Error   t value      Pr(>|t|)
## (Intercept)          5.006 0.07280222 68.761639 1.134286e-113
## Speciesversicolor    0.930 0.10295789  9.032819  8.770194e-16
## Speciesvirginica     1.582 0.10295789 15.365506  2.214821e-32

👉Cliquez ici pour accéder à un tutoriel sur l’ANOVA à un facteur avec R

 

ANOVA à deux facteurs avec interaction

library(funModeling)
mod4 <- lm(max_heart_rate~gender*has_heart_disease , data=heart_disease)
extract_eq(mod4, use_coefs=TRUE, wrap = TRUE)

équation anova à 2 facteurs

summary(mod4)$coefficients

##                                   Estimate Std. Error   t value      Pr(>|t|)
## (Intercept)                     154.027778   2.434886 63.258725 3.789051e-175
## gendermale                        7.754831   3.250923  2.385425  1.768205e-02
## has_heart_diseaseyes            -10.867778   4.796169 -2.265929  2.417054e-02
## gendermale:has_heart_diseaseyes -12.511322   5.602454 -2.233186  2.627647e-02

👉 Cliquez ici pour accéder à un tutoriel sur l’ANOVA à 2 facteurs avec R 

Régression linéaire multiple

mod5 <- lm(Sepal.Length~Sepal.Width+Petal.Length+Species, data=iris)
extract_eq(mod5, use_coefs=TRUE, wrap = TRUE)

équation régression linéaire multiple

summary(mod5)$coefficients

##                     Estimate Std. Error   t value     Pr(>|t|)
## (Intercept)        2.3903891 0.26226815  9.114294 5.942826e-16
## Sepal.Width        0.4322172 0.08138982  5.310458 4.025982e-07
## Petal.Length       0.7756295 0.06424566 12.072869 1.151112e-23
## Speciesversicolor -0.9558123 0.21519853 -4.441537 1.759999e-05
## Speciesvirginica  -1.3940979 0.28566053 -4.880261 2.759618e-06

GLM avec distribution de Poisson

str(warpbreaks)

## 'data.frame':    54 obs. of  3 variables:
##  $ breaks : num  26 30 54 25 70 52 51 26 67 18 ...
##  $ wool   : Factor w/ 2 levels "A","B": 1 1 1 1 1 1 1 1 1 1 ...
##  $ tension: Factor w/ 3 levels "L","M","H": 1 1 1 1 1 1 1 1 1 2 ...
mod6 <- glm(breaks~wool*tension , data=warpbreaks, family = poisson)
extract_eq(mod6, use_coefs=TRUE, wrap = TRUE)

équation glm poisson

summary(mod6)$coefficients

##                  Estimate Std. Error   z value     Pr(>|z|)
## (Intercept)     3.7967368 0.04993753 76.029734 0.000000e+00
## woolB          -0.4566272 0.08019202 -5.694172 1.239721e-08
## tensionM       -0.6186830 0.08440012 -7.330357 2.295399e-13
## tensionH       -0.5957987 0.08377723 -7.111702 1.146202e-12
## woolB:tensionM  0.6381768 0.12215312  5.224400 1.747203e-07
## woolB:tensionH  0.1883632 0.12989529  1.450115 1.470263e-01

👉 Cliquez ici pour accéder au Tutoriel : GLM sur données de comptage (régression de Poisson) avec R  

 

 

GLM avec distribution quasi poisson

mod7 <- glm(breaks~wool*tension , data=warpbreaks, family = quasipoisson)
extract_eq(mod7, use_coefs=TRUE, wrap = TRUE)

equation quasi poisson

summary(mod7)$coefficients


##                  Estimate Std. Error    t value     Pr(>|t|)
## (Intercept)     3.7967368 0.09688254 39.1890723 4.193343e-38
## woolB          -0.4566272 0.15557852 -2.9350269 5.104563e-03
## tensionM       -0.6186830 0.16374255 -3.7783889 4.359386e-04
## tensionH       -0.5957987 0.16253410 -3.6656845 6.160288e-04
## woolB:tensionM  0.6381768 0.23698620  2.6928860 9.726580e-03
## woolB:tensionH  0.1883632 0.25200659  0.7474534 4.584362e-01

Régression logistique

library(funModeling)
mod8 <- glm(has_heart_disease~age+max_heart_rate , data=heart_disease, family = binomial)
extract_eq(mod8, use_coefs=TRUE, wrap = TRUE)

équation régression logistique

summary(mod8)$coefficients

##                   Estimate  Std. Error   z value     Pr(>|z|)
## (Intercept)     4.83959656 1.498563177  3.229491 1.240107e-03
## age             0.01970614 0.015338849  1.284721 1.988899e-01
## max_heart_rate -0.04071090 0.006808203 -5.979684 2.235714e-09

👉 Cliquez ici pour accéder à un tutoriel dur la La régression logistique avec R

Modèles mixtes

library(lme4)
str(sleepstudy)

## 'data.frame':    180 obs. of  3 variables:
##  $ Reaction: num  250 259 251 321 357 ...
##  $ Days    : num  0 1 2 3 4 5 6 7 8 9 ...
##  $ Subject : Factor w/ 18 levels "308","309","310",..: 1 1 1 1 1 1 1 1 1 1 ...
mod9 <- lmer(Reaction~Days + (1|Subject), data=sleepstudy)
extract_eq(mod9, use_coefs=TRUE, wrap = TRUE)

rendu équation modèle mixtes

Les modèles mixtes ne sont pas encore entièrement pris en charge par ce package, mais ça semble être dans les tuyaux : https://github.com/datalorax/equatiomatic/issues/105.

Les méthodes prises en charges

Les méthodes prise en charges sont résumées dans ce tableau

liste des méthodes prises en charge par equatiomatic

D’après https://cran.r-project.org/web/packages/equatiomatic/vignettes/intro-equatiomatic.html

Conclusion

J’aime beaucoup ce package, qui non seulement simplifie la vie, mais aussi contribue à la compréhension des approches statistiques.
Et vous, qu’en pensez-vous ?

Si cet article vous a plu ou vous a été utile, n’oubliez pas de le partager ! Vous pouvez également soutenir le blog par un don libre sur la page Tipee.

soutenez le blog

Image par Ralf Designs de Pixabay

Poursuivez votre lecture:

Laisser un commentaire

Votre adresse de messagerie ne sera pas publiée. Les champs obligatoires sont indiqués avec *