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 :
Et si vous incluez cette commande au sein d’un chunk d’un document rmarkdown, comme ceci :
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 :
Voici le 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 :
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) 
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) 
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) 
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) 
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)
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) 
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 :
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
4 Responses
Super article claire merci pour ce travail excellent !!!!!
Euh… mais c’est incroyable ce truc!
Bonjour,
ah oui, là j’avoue, votre article tombe à point nommé !! Merci beaucoup !
Toutefois, y a-t-il moyen d’obtenir les équations autrement que sur un document markdown ?
Et autre question, si je souhaite par exemple seulement déterminer l’équation de la courbe utilisée avec ggplot geom_smooth(method = “lm”), pensez-vous que cela fonctionne ?
Merci beaucoup pour votre aide,
merci pour l article