C'est de la data et mon expeRtise afin d'en tirer le maximum

Equatiomatic : un package bien sympathique pour vos équations !

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 :

extraire l'équation en latex

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

equatiomatic avec r markdown

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

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 du modèle

Voici le rendu pdf :

equatiomatic et rendu pdf

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 equatiomatic

Ainsi que ce message d’erreur :

equatiomatic message word

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) 
Equationmatic, rendu 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) 
équation 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 

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) 
equation anova 2 avec interaction
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 

Régression linéaire multiple

mod5 <- lm(Sepal.Length~Sepal.Width+Petal.Length+Species, data=iris)
extract_eq(mod5, use_coefs=TRUE, wrap = TRUE) 
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) 
equatiomatic avec GLM
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 

GLM avec distribution quasi poisson

mod7 <- glm(breaks~wool*tension , data=warpbreaks, family = quasipoisson)
extract_eq(mod7, use_coefs=TRUE, wrap = TRUE) 
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)
 
equation 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 
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) 
equation d'un modèle mixte

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 an charge par equatiomatic

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

liste des méthodes gérées par le package equatiomatic

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 Tipeee.

Poursuivez votre lecture

Image par Ralf Designs de Pixabay

3 réponses

Laisser un commentaire

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

Bonjour !

vous venez souvent ?

Identifiez-vous pour avoir accès à toutes les fontionnalités !

Aide mémoire off'R ;)

Enregistrez vous pour recevoir gratuitement mes fiches “aide mémoire” (ou cheat sheets) qui vous permettront de réaliser facilement les principales analyses biostatistiques avec le logiciel R et pour être informés des mises à jour du site.