Comment calculer la différence de deux pentes?

11

Existe-t-il une méthode pour comprendre si deux lignes sont (plus ou moins) parallèles? J'ai deux lignes générées à partir de régressions linéaires et je voudrais comprendre si elles sont parallèles. En d'autres termes, je voudrais obtenir les différentes pentes de ces deux lignes.

Existe-t-il une fonction R pour calculer cela?

EDIT: ... et comment puis-je obtenir la pente (en degrés) d'une droite de régression linéaire?

Dail
la source

Réponses:

23

Je me demande si je manque quelque chose d'évident, mais ne pourriez-vous pas le faire statistiquement en utilisant ANCOVA? Un problème important est que les pentes des deux régressions sont estimées avec erreur. Ce sont des estimations des pentes de l'ensemble des populations. Si la préoccupation est de savoir si les deux droites de régression sont parallèles ou non dans la population, il n'est pas logique de comparer directement a1 avec a2 pour une équivalence exacte; ils sont tous les deux sujets à des erreurs / incertitudes qui doivent être prises en compte.

Si nous y réfléchissons d'un point de vue statistique, et que nous pouvons combiner les données sur et pour les deux ensembles de données de manière significative (c'est-à-dire que et dans les deux ensembles sont tirés des deux populations avec des plages similaires pour les deux variables c'est juste la relation entre elles qui sont différentes dans les deux populations), alors on peut adapter les deux modèles suivants:y x yxyxy

y^=b0+b1x+b2g

et

y^=b0+b1x+b2g+b3xg

Où sont les coefficients du modèle et est une variable / facteur de regroupement, indiquant à quel ensemble de données appartient chaque observation. gbig

Nous pouvons utiliser une table ANOVA ou un rapport F pour tester si le deuxième modèle, plus complexe, correspond mieux aux données que le modèle plus simple. Le modèle plus simple indique que les pentes des deux lignes sont les mêmes ( ) mais les lignes sont décalées l'une de l'autre d'un montant .b 2b1b2

Le modèle plus complexe comprend une interaction entre la pente de la ligne et la variable de regroupement. Si le coefficient de ce terme d'interaction est significativement différent de zéro ou si le rapport ANOVA / F indique que le modèle plus complexe correspond mieux aux données, nous devons rejeter l'hypothèse Null selon laquelle deux lignes sont parallèles.

Voici un exemple dans R utilisant des données fictives. Tout d'abord, les données avec des pentes égales:

set.seed(2)
samp <- factor(sample(rep(c("A","B"), each = 50)))
d1 <- data.frame(y = c(2,5)[as.numeric(samp)] + (0.5 * (1:100)) + rnorm(100),
                 x = 1:100,
                 g = samp)
m1 <- lm(y ~ x * g, data = d1)
m1.null <- lm(y ~ x + g, data = d1)
anova(m1.null, m1)

Qui donne

> anova(m1.null, m1)
Analysis of Variance Table

Model 1: y ~ x + g
Model 2: y ~ x * g
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     97 122.29                           
2     96 122.13  1   0.15918 0.1251 0.7243

Indiquant que nous ne parvenons pas à rejeter l'hypothèse nulle de pentes égales dans cet échantillon de données. Bien sûr, nous voudrions nous assurer que nous avions suffisamment de puissance pour détecter une différence s'il y en avait vraiment une afin que nous ne soyons pas amenés à échouer par erreur à rejeter le zéro parce que notre taille d'échantillon était trop petite pour l'effet attendu.

Maintenant avec différentes pentes.

set.seed(42)
x <- seq(1, 100, by = 2)
d2 <- data.frame(y = c(2 + (0.5 * x) + rnorm(50),
                       5 + (1.5 * x) + rnorm(50)),
                 x = x,
                 g = rep(c("A","B"), each = 50))
m2 <- lm(y ~ x * g, data = d2)
m2.null <- lm(y ~ x + g, data = d2)
anova(m2.null, m2)

Qui donne:

> anova(m2.null, m2)
Analysis of Variance Table

Model 1: y ~ x + g
Model 2: y ~ x * g
  Res.Df     RSS Df Sum of Sq     F    Pr(>F)    
1     97 21132.0                                 
2     96   103.8  1     21028 19439 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1

Ici, nous avons des preuves substantielles contre l'hypothèse nulle et nous pouvons donc la rejeter en faveur de l'alternative (en d'autres termes, nous rejetons l'hypothèse que les pentes des deux droites sont égales).

Les termes d'interaction dans les deux modèles que j'ai ajustés ( ) donnent la différence estimée de pentes pour les deux groupes. Pour le premier modèle, l'estimation de la différence de pente est faible (~ 0,003)b3xg

> coef(m1)
(Intercept)           x          gB        x:gB 
2.100068977 0.500596394 2.659509181 0.002846393

et un test à ce sujet ne permettrait pas de rejeter l'hypothèse nulle selon laquelle cette différence de pente est de 0:t

> summary(m1)

Call:
lm(formula = y ~ x * g, data = d1)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.32886 -0.81224 -0.01569  0.93010  2.29984 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 2.100069   0.334669   6.275 1.01e-08 ***
x           0.500596   0.005256  95.249  < 2e-16 ***
gB          2.659509   0.461191   5.767 9.82e-08 ***
x:gB        0.002846   0.008047   0.354    0.724    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1 

Residual standard error: 1.128 on 96 degrees of freedom
Multiple R-squared: 0.9941, Adjusted R-squared: 0.9939 
F-statistic:  5347 on 3 and 96 DF,  p-value: < 2.2e-16 

Si nous nous tournons vers le modèle ajusté au deuxième ensemble de données, où nous avons fait différer les pentes pour les deux groupes, nous voyons que la différence estimée dans les pentes des deux lignes est d'environ 1 unité.

> coef(m2)
(Intercept)           x          gB        x:gB 
  2.3627432   0.4920317   2.8931074   1.0048653 

La pente pour le groupe "A" est ~ 0,49 ( xdans la sortie ci-dessus), tandis que pour obtenir la pente pour le groupe "B", nous devons ajouter les différences de pente (données par le terme d'interaction rappelez-vous) à la pente du groupe "A" ; ~ 0,49 + ~ 1 = ~ 1,49. Ceci est assez proche de la pente indiquée pour le groupe "B" de 1,5. Un test sur cette différence de pentes indique également que l'estimation de la différence est limitée à 0:t

> summary(m2)

Call:
lm(formula = y ~ x * g, data = d2)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.1962 -0.5389  0.0373  0.6952  2.1072 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 2.362743   0.294220   8.031 2.45e-12 ***
x           0.492032   0.005096  96.547  < 2e-16 ***
gB          2.893107   0.416090   6.953 4.33e-10 ***
x:gB        1.004865   0.007207 139.424  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1 

Residual standard error: 1.04 on 96 degrees of freedom
Multiple R-squared: 0.9994, Adjusted R-squared: 0.9994 
F-statistic: 5.362e+04 on 3 and 96 DF,  p-value: < 2.2e-16
Gavin Simpson
la source
merci beaucoup pour cette très bonne explication. Mon objectif est de comprendre si le sloper est plus ou moins le même donc je pense que je vais utiliser l'ANOVA pour le tester.
Dail
si j'ai deux vecteurs distint et que je voudrais comparer leurs slops mais je n'ai pas le y (lm (x ~ y), comment puis-je utiliser ANOVA? J'ai essayé anova (lm (x ~ 1), lm (y ~ 1)) mais je reçois un avertissement
Dail
Qu'entendez-vous par vecteurs ici? Au sens R ou au sens mathématique? Ceci est très différent de la question que vous avez posée, alors veuillez commencer une nouvelle question - ne modifiez pas celle-ci - il est impossible d'effectuer des suivis d'une telle nature dans les commentaires.
Gavin Simpson,
pas d'attente, je dois comparer deux modèles avec ANOVA ... ok, mais si je crée un modèle avec cette formule: x ~ 1 et un autre modèle avec y ~ 1 je reçois l'avertissement. Je parle dans le sens R. Comment puis-je faire?
Dail
1
@Dail si vous avez ajusté deux régressions pour obtenir deux pentes / lignes, vous disposez de données x et y pour les deux ensembles de données. Comme je l'ai dit dans ma réponse, si les x et les y sont comparables dans les deux ensembles de données, vous pouvez simplement combiner toutes les données et ajouter une variable de regroupement. Mon exemple montre comment faire cela en utilisant des données factices, mais vous avez déjà des données x et y, ce sont les données que vous avez utilisées pour ajuster les régressions distinctes.
Gavin Simpson
8

La première question vient en fait de la géométrie. Si vous avez deux lignes du formulaire:

y = a 2 x + b 2

y=a1x+b1
y=a2x+b2

alors ils sont parallèles si . Donc, si les pentes sont égales, les lignes sont parallèles.a1=a2

Pour la deuxième question, utilisez le fait que , où est l'angle que la ligne fait avec l' axe des et est la pente de la ligne. Donc α x a 1tanα=a1αxa1

α=arctana1

et pour convertir en degrés, rappelez-vous que . Ainsi, la réponse dans les degrés sera2π=360

α=arctana13602π.

La fonction R pour est appelée .arctanatan

Exemple de code R:

> x<-rnorm(100)
> y<-x+1+rnorm(100)/2
> mod<-lm(y~x)
> mod$coef
    (Intercept)           x 
      0.9416175   0.9850303 
    > mod$coef[2]
        x 
0.9850303 
> atan(mod$coef[2])*360/2/pi
       x 
44.56792 

La dernière ligne correspond aux degrés.

Mise à jour. Pour les valeurs de pente négatives, la conversion en degrés doit suivre une règle différente. Notez que l'angle avec l'axe des x peut obtenir des valeurs de 0 à 180, car nous supposons que l'angle est au-dessus de l'axe des x. Ainsi, pour les valeurs négatives de , la formule est:a1

α=180arctana13602π.

Remarque. Alors que c'était amusant pour moi de rappeler la trigonométrie au lycée, la réponse vraiment utile est celle donnée par Gavin Simpson. Étant donné que les pentes des droites de régression sont des variables aléatoires, pour les comparer, un cadre d'hypothèse statistique doit être utilisé.

mpiktas
la source
Merci! comment obtenir la pente de la régression? dois-je obtenir le coefficient et l'interception?
Dail
peut-être que la régression linéaire renvoie directement les degrés avec une fonction?
Dail
dire degress = +45 et degress = -315 ne sont pas la même ligne? whare ne parle pas de la même ligne?
Dail
1

... suite à la réponse de @mpiktas, voici comment extraire la pente d'un lmobjet et appliquer la formule ci-dessus.

# prepare some data, see ?lm
ctl <- c(4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14)
trt <- c(4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69)
group <- gl(2,10,20, labels=c("Ctl","Trt"))
weight <- c(ctl, trt)

lm.D9 <- lm(weight ~ group)
# extract the slope (this is also used to draw a regression line if you wrote abline(lm.D9)
coefficients(lm.D9)["groupTrt"] 
      groupTrt 
   -0.371 
# use the arctan*a1 / (360 / (2*pi)) formula provided by mpiktas
atan(coefficients(lm.D9)["groupTrt"]) * (360/(2 * pi)) 
 groupTrt 
-20.35485 
180-atan(coefficients(lm.D9)["groupTrt"]) * (360/(2 * pi))
 groupTrt 
200.3549 
Roman Luštrik
la source
merci beaucoup pour l'exemple, dans ce cas les degrés sont -200?
Dail
Oui. Si vous pensez que votre problème a été résolu, cochez la réponse de @mpiktas comme étant la bonne, merci.
Roman Luštrik
@ RomanLuštrik, vous avez raté la division par2π
mpiktas
1
@ RomanLuštrik, j'ai corrigé le code et ajouté la sortie correcte. N'hésitez pas à supprimer la correction.
mpiktas du