Récupérer les estimations de pente par modalité depuis une ANCOVA

Postez ici vos questions, réponses, commentaires ou suggestions - Les sujets seront ultérieurement répartis dans les archives par les modérateurs

Modérateur : Groupe des modérateurs

Chrystelle Delord
Messages : 33
Enregistré le : 09 Avr 2013, 12:50

Récupérer les estimations de pente par modalité depuis une ANCOVA

Messagepar Chrystelle Delord » 31 Juil 2018, 17:52

Bonjour à tous,

Je dispose de 18 petits jeux de données correspondant à 18 espèces différentes de poissons, échantillonnés sur un maximum de 12 sites.

Chacun de ces datasets contient 9 à 12 observations, qui correspondent à 9 à 12 sites identiques d'une espèce à l'autre (aux valeurs manquantes près, d'où le fait d'avoir moins de 12 observations pour certaines espèces !)

Par ailleurs, je dispose de deux prédicteurs continus pour décrire mes sites: appelons-les X1 et X2.
Mon souhait est d'estimer la taille d'effet de X1, tout en tenant compte de X2; sur une variable réponse Y, et ce pour chacune des mes espèces. Je souhaite ainsi rapporter les effets relatifs de X1 et X2 sur Y grâce à l'estimateur de pente (slope) de chaque prédicteur et son intervalle de confiance (étant donné le faible nombre d'observations, je souhaite garder un maximum d'information). Initialement pour cela, j'avais effectué un modèle linéaire multiple pour chacune des espèces séparément. La fonction 'summ()' du package R jtools m'a permis de récupérer estimations de pentes + ICs (commande ci-dessous), et j'ai ensuite pu tracer un forest plot pour chacune des espèces (18 graphes x 2 estimations + CI).

Code : Tout sélectionner

for (i in unique(Tab$Espece)) {
    print(cat(paste(" ", "Modèle linéaire multiple pour l'espèce :", i, sep='\n')))
    print(summ(lm(Y~X1+X2, data=Tab[Tab$Species==i,]), scale=T, confint=T, digits=4))
}


(s'ensuivent ensuite mes 18 outputs)

L'avantage de cette méthode pour moi résidait dans la possibilité d'avoir les tailles d'effet relatives de chaque prédicteur pour chaque espèce, sans m'embarrasser d'une variable qualitative quelconque. Ainsi, chaque pente était simplement caractérisée par sa valeur et l'incertitude associée. Cependant, lorsque vient le besoin de comparer réellement les pentes d'une espèce à l'autre, il devient bien sûr souhaitable d'ajouter le facteur Espèce comme covariable et donc d'effectuer une ANCOVA. Comme tout cela n'est pas encore très R-technique, voici ma question réelle:

Puis-je effectuer une ANCOVA pour tester les interactions entre X1 et Espèce (tout en contrôlant toujours par X2 pour chaque Espèce), avec des commandes qui permettraient à la fois de (1) comparer bien entendu les pentes de X1, compte tenu de X2, d'une espèce à l'autre et surtout (2) de pouvoir continuer à comparer les tailles d'effet de X1 et X2 par espèce - donc de récupérer pentes + ICs pour chacune ?

Mon souci principal étant que, quelle que soit la façon de régler les contrastes, je ne sais pas comment - ni même si je peux - récupérer les estimations brutes de pentes de X1 et X2 par espèce depuis l'output de la commande

Code : Tout sélectionner

lm(Y~(X1+X2+Espece)^2, data=Tab)

comme je le faisais dans mes modèles multiples séparés avec 'summ'. Or, outre le fait de tester les différences de pentes entre Espèces (ce que je peux bien sûr faire avec la commande ci-dessus suivie de tests post-hoc), j'ai quand même besoin de récupérer les estimations brutes elle-même (les pentes), car il me faudra ces valeurs pour des analyses complémentaires. Comme 'lm()' se base sur la première modalité (contr.treatment) ou la moyenne globale (contr.sum) de la variable catégorielle pour fournir les coefficients, je ne sais comment procéder.

En conclusion: il me semble actuellement que l'ANCOVA, si elle me permettra de tester les différences de pentes de X1 entre Espèces compte tenu de X2, ne me permettra pas, en revanche, de simplement récupérer les estimations per se de pentes de X1 compte tenu de X2 que je pouvais obtenir en faisant des modèles multiples séparés par Espèce, lorsque je ne cherchais pas forcément à tester si les pentes étaient significativement différentes ou non. Et je me dis qu'en réalité j'ai probablement simplement loupé une fonction très simple qui me permettrait de récupérer ces estimations de pente directement depuis l'ANCOVA. Et donc, de tout faire d'un seul coup: tester l'effet de X1 et X2 par Espèce, comparer les pentes de X1 compte tenu de X2 entre Espèces, et récupérer l'ensemble des coefficients de pentes et CIs pour des interprétations complémentaires que je ne peux intégrer au modèle (e.g. voir si la pente évolue en fonction d'un autre critère relatif à l'espèce elle-même, comme sa taille)

J'espère avoir été claire et cependant pas trop longue !

Un grand merci par avance de votre réponse.

Chrystelle

Florent Aubry
Messages : 324
Enregistré le : 25 Juin 2010, 10:21

Re: Récupérer les estimations de pente par modalité depuis une ANCOVA

Messagepar Florent Aubry » 01 Aoû 2018, 08:29

Va voir du côté du package emmeans. Il possède les fonctions qui répondent à ta question : récupérations des estimations et tests post-hoc.

Chrystelle Delord
Messages : 33
Enregistré le : 09 Avr 2013, 12:50

Re: Récupérer les estimations de pente par modalité depuis une ANCOVA

Messagepar Chrystelle Delord » 01 Aoû 2018, 12:31

Merci Florent, c'était effectivement aussi simple que cela !

Les fonctions emmeans, lsmeans ou encore lstrends me fournissent toutes les informations nécessaires.

Bonne journée à tous,
Chrystelle

Chrystelle Delord
Messages : 33
Enregistré le : 09 Avr 2013, 12:50

Re: Récupérer les estimations de pente par modalité depuis une ANCOVA

Messagepar Chrystelle Delord » 02 Aoû 2018, 07:58

Re-bonjour,

Me revoici déjà :)
J'ai pu récupérer les coefficients de régression (B) entre Y et X1 ou X2 pour chaque niveau de mon facteur espèce grâce à la fonction lstrends du package R lsmeans. Cependant ils ont des ordres de grandeurs très différents (de l'ordre de 10²) entre X1 et X2, et j'aimerais donc travailler plutôt sur les coefficients standardisés (bêtas) pour reporter leurs effets relatifs sur ma variable réponse par espèce.

Je ne sais pas comment procéder, car appliquer un 'scale' sur X1 et X2 en amont et/ou un 'scale' sur la variable réponse également, sur l'ensemble du jeu de données toutes espèces comprises, ne change rien au problème. Lorsque je travaillais sur chaque espèce séparément, j'utilisais l'argument 'scale=TRUE' de la fonction summ, voici un exemple ci-dessous:

Coefficients non-standardisés:

Code : Tout sélectionner

> summ(mod.Sp1, [b]scale=F[/b], confint=TRUE, vifs=TRUE, digits=4)
MODEL INFO:
Observations: 10 (2 missing obs. deleted)
Dependent Variable: Fst.p
Type: OLS linear regression

MODEL FIT:
F(2,7) = 13.1896, p = 0.0042
R² = 0.7903
Adj. R² = 0.7304

Standard errors: OLS
                 Est.      2.5%     97.5%  t val.   p           VIF   
(Intercept) 2.9570  1.1173 4.7967 3.8007 0.0067   <NA> **
X1             0.0000 -0.0001 0.0002 0.4073 0.6960 1.6084   
X2             0.0542  0.0203 0.0880 3.7866 0.0068 1.6084 **


Coefficients standardisés:

Code : Tout sélectionner

> summ(mod.Sp1, [b]scale=T[/b], confint=TRUE, vifs=TRUE, digits=4)
MODEL INFO:
Observations: 10 (2 missing obs. deleted)
Dependent Variable: Fst.p
Type: OLS linear regression

MODEL FIT:
F(2,7) = 13.1896, p = 0.0042
R² = 0.7903
Adj. R² = 0.7304

Standard errors: OLS
              Est.    2.5%  97.5%  t val.      p    VIF   
(Intercept) 0.0405  0.0342 0.0468 15.2361 0.0000   <NA> ***
X1      0.0014 -0.0070 0.0098  0.4073 0.6960 1.6084   
X2      0.0135  0.0051 0.0219  3.7866 0.0068 1.6084  **

Continuous predictors are mean-centered and scaled by 1 s.d.


On peut constater les échelles plus proches des coef de X1 et X2 et de leurs IC dans la seconde version. Par défaut, summ ne centre et réduit que les prédicteurs. lstrends ne propose de transformer que la variable réponse. Y aurait-il une façon plus judicieuse de procéder ? Est-ce que le fait d'intégrer l'Espèce comme covariable dans un modèle globale explique le problème ? Un test de Breusch-Pagan pour détecter l'hétéroscédasticité des résidus entre espèces n'a pas détecté d'effet significatif.

Encore merci !

Florent Aubry
Messages : 324
Enregistré le : 25 Juin 2010, 10:21

Re: Récupérer les estimations de pente par modalité depuis une ANCOVA

Messagepar Florent Aubry » 02 Aoû 2018, 09:04

summ ne fait que présenter différemment les résultats d'une analyse et le fait de normaliser ou non les valeurs ne change rien aux résultats du test, sauf les valeurs elles-mêmes, mais pas les F, t, p et R2 ajusté ou non et le rapport entre la pente normalisée (coefficient beta) et celle non normalisée est égal à la déviation standard du régresseur dans tous les cas.

Les différences proviennent des modèles. En effet, d'un côté tu testes un modèle bilinéaire dont la formule est Y ~ X1 + X2, donc considéré comme deux régresseurs indépendants et d'un autre un modèle avec 3 prédicteurs ayant entre eux des interactions d'ordre 2, dont la formule est : Y ~ (X1 + X2 + Espece)^2. Développé, le modèle devient : Y ~ X1 + X2 + Espece + X1:X2 + X1:Espece + X2:Espece. Ce modèle suppose donc deux choses : i) les régresseurs ne sont pas orthogonaux et il existe un effet lié à leur interaction, ii) que les pentes dépendent de l'espèce. Si au moins une de ces conditions est vraie, il est normal que les résultats diffèrent largement entre ces deux modèles.

Chrystelle Delord
Messages : 33
Enregistré le : 09 Avr 2013, 12:50

Re: Récupérer les estimations de pente par modalité depuis une ANCOVA

Messagepar Chrystelle Delord » 02 Aoû 2018, 12:11

Merci Florent, il est vrai que dans mon modèle global, j'ai omis d'ôter l'interaction entre X1 et X2 qui de toute façon n'est pas significative puisque les deux régresseurs sont bel et bien indépendants, comme le sous-entendaient les régressions linéaires multiples effectuées précédemment.

Je réécris donc:

Code : Tout sélectionner

mod.1 <- lm(formula = Y ~ X1+X2+Espece+X1:Espece+X2:Espece, data = Tab)
lstrends(mod.1, "Espece", var="X1")

et j'obtiens pour X1:

Code : Tout sélectionner

 > lstrends(mod.1, "Espece", var="X1")
 Species  X1.trend           SE  df      lower.CL     upper.CL
 Sp1     2.452719e-05 2.888757e-05 131 -3.261933e-05 8.167370e-05
 Sp2     1.901338e-05 2.825246e-05 131 -3.687673e-05 7.490348e-05
 Sp3     4.781648e-06 2.075202e-05 131 -3.627080e-05 4.583410e-05
 [...]

... et pour X2:

Code : Tout sélectionner

> lstrends(mod.1, "Espece", var="X2")
 Species    X2.trend          SE  df     lower.CL    upper.CL
 Sp1     0.0541616421 0.006860727 131  0.040589487 0.067733797
 Sp2     0.0222920812 0.005461422 131  0.011488086 0.033096077
 Sp3     0.0031760890 0.005511610 131 -0.007727190 0.014079368
[...]


Si je reprends ma régression multiple sur Sp1 par exemple, avec summ et scale=FALSE (la même que sur mon message précédent à ceci près que je prends jusqu'à la 5è décimale):

Code : Tout sélectionner

> summ(mod.Sp1, scale=F, confint=TRUE, vifs=TRUE, digits=5)
MODEL INFO:
Observations: 10 (2 missing obs. deleted)
Dependent Variable: Fst.p
Type: OLS linear regression

MODEL FIT:
F(2,7) = 13.18960, p = 0.00422
R² = 0.79029
Adj. R² = 0.73037

Standard errors: OLS
               Est.     2.5%   97.5%  t val.       p     VIF   
(Intercept) 2.95702  1.11731 4.79672 3.80074 0.00671    <NA> **
X1      0.00002 -0.00012 0.00017 0.40725 0.69599 1.60840   
X2      0.05416  0.02034 0.08798 3.78658 0.00683 1.60840 **


Ici on voit que les coefficients B sont quasi identiques entre ce que donne summ pour Sp1 seule, et ce que donne lstrends pour Sp1 via le modèle global: B=0.00002 pour X1 et B=0.05416 pour X2. Cela me semblait logique puisque, bien que les interactions X1:Espece et X2:Espece soient bien significatives, lstrends me donne normalement chaque pente spécifique. Du coup, je crois que je ne comprends toujours pas bien pourquoi les coefficients standardisés, eux, varient entre les deux types de modèle.

Vous souhaitiez peut-être souligner, dans votre réponse, que standardiser les coefficients ne donnerait pas la même chose selon qu'une covariable soit ajoutée ou pas; dans la mesure où elle interagit avec les prédicteurs continus; et ce même si les moyennes marginales de pentes restent similaires entre les deux approches. C'est peut-être cela que j'ai du mal à me représenter mathématiquement ! En bref, sois je saisis effectivement mal cette subtilité statistique des coefficients standardisés, soit j'ai mal procédé pour standardiser mes coefficients du modèle global.

J'aurais pu me contenter de garder mes coefficients B, mais la différence d'échelle est telle que c'est un cauchemar à rapporter graphiquement (forest plot) ! Encore une fois merci de prendre tout ce temps pour me répondre :)

Florent Aubry
Messages : 324
Enregistré le : 25 Juin 2010, 10:21

Re: Récupérer les estimations de pente par modalité depuis une ANCOVA

Messagepar Florent Aubry » 03 Aoû 2018, 08:55

1) Le nombre de décimales ne diffèrent qu'au niveau de l'affichage et pas au niveau des valeurs calculées. Cela provient que les deux procédures n'utilisent pas le même arrondi pour l'affichage. Exemple :

Code : Tout sélectionner

set.seed( 1)
z
round( z, 3)
round( z, 5)


2) Écriture du modèle :

Code : Tout sélectionner

Y ~ (X1 + X2) * Espece


3) lsmeans est obsolète. Il faut maintenant utiliser emmeans.
IMPORTANT : ne pas oublier d'utiliser la dernière version de R.

4) Tests des pentes en utilisant emtrends :

Code : Tout sélectionner

summary( emtrends( mod.1, pairwise ~ Espece, var="X1") , infer=TRUE) # ou var="X2"


5) Standardiser des coefficients revient à calculer à partir de régresseurs dont on a modifié les valeurs puisqu'elles sont divisées par la déviation standard du régresseur. Les centrer changent l'ordonnée à l'origine. Donc si on a :
y = ax + b
et qu'on transforme x en x0 = (x - mean( x)) / sd( x)
d'où x = sd( x) * x0 + mean( x)
alors y = a * (sd( x) * x0 + mean( x)) + b = a * sd( x) * x0 + (b + a * mean( x))
Ce qui montre bien que les pentes et les ordonnées à l'origine diffèrent.

6) jtools::summ recalcule éventuellement les coefficients qu'on obtient par la fonction summary, donc cette procédure donne le coefficient (pente) pour le régresseur indépendamment du facteur puis des différences / contrastes entre les pentes de chacun des groupes, différences qui dépendent du contraste utilisé (par défaut pour R, contr.treatment). emtrends (lstrends) donne les pentes dans chacun des groupes puis si on a donné une spécification, fait les tests post-hoc demandés. On ne peut donc pas comparer les résultats de jtools::summ (ou de summary) et de emtrends, ces procédures ne donnant pas les mêmes informations puisque emtrends travaille sur les combinaisons linéaires ad-hoc (i.e., dépendantes des contrastes) de summary.

Chrystelle Delord
Messages : 33
Enregistré le : 09 Avr 2013, 12:50

Re: Récupérer les estimations de pente par modalité depuis une ANCOVA

Messagepar Chrystelle Delord » 03 Aoû 2018, 12:15

Encore une fois, merci à vous Florent.

J'ai désormais compris l'origine de mon erreur: j'ai fait une confusion entre (1) les paramètres estimés par le modèle de type ANCOVA (pente globale des prédicteurs et contrastes entre facteurs, i.e. au niveau des ordonnées à l'origine + pentes par groupe pour chaque régresseur), et (2) les estimations des moyennes marginales par niveau de facteurs. Et du coup il apparaît logique que les coefficients betas se calculent sur les estimations des paramètres uniquement, et non pas sur les moyennes marginales de l'ANCOVA. Ca marchait donc avec summ, logique, et ça n'est pas possible avec emtrends, logique également.

Merci pour vos informations et les corrections/suggestions de commandes ! :)

Bonne journée à tous,
Chrystelle

Florent Aubry
Messages : 324
Enregistré le : 25 Juin 2010, 10:21

Re: Récupérer les estimations de pente par modalité depuis une ANCOVA

Messagepar Florent Aubry » 03 Aoû 2018, 12:31

Attention à l'interprétation des résultat de summary ou de jtools::summ. Les coefficients calculés dépendent du contraste utilisé, tandis que ce n'est pas vrai pour emtrends. Pour t'en convaincre, compare les résultats suivants :

Code : Tout sélectionner

lm.1 <- lm( Y ~ (X1 + X2) * Espece, data=Tab, contrasts=list( Espece=contr.treatment)) # defaut pour R
lm.2 <- lm( Y ~ (X1 + X2) * Espece, data=Tab, contrasts=list( Espece=contr.sum))

jtools::summ( lm.1)
jtools::summ( lm.2)

emtrends( lm.1, pairwise ~ Espece, "X1")
emtrends( lm.2, pairwise ~ Espece, "X1")

emtrends( lm.1, pairwise ~ Espece, "X2")
emtrends( lm.2, pairwise ~ Espece, "X2")


Retourner vers « Questions en cours »

Qui est en ligne

Utilisateurs parcourant ce forum : Aucun utilisateur enregistré et 1 invité