Jak obliczyć różnicę dwóch stoków?

Czy istnieje metoda zrozumienia, czy dwie linie są (mniej więcej) równoległe? Mam dwie linie wygenerowane z regresji liniowych i chciałbym zrozumieć, czy są one równoległe. Innymi słowy, chciałbym uzyskać inne nachylenie tych dwóch linii.

Czy istnieje funkcja R do obliczenia tego?

EDYCJA: ... i jak mogę uzyskać nachylenie (w stopniach) linii regresji liniowej?

Zastanawiam się, czy brakuje mi czegoś oczywistego, ale czy nie można tego zrobić statystycznie za pomocą ANCOVA? Ważną kwestią jest to, że nachylenia w dwóch regresjach są szacowane z błędem. Są to szacunkowe wartości stoków w całej populacji. Jeśli chodzi o to, czy dwie linie regresji są równoległe w populacji , czy nie, porównanie nie ma sensu$a_1$ z $a_2$bezpośrednio w celu uzyskania dokładnej równoważności; oba są obarczone błędem / niepewnością, którą należy wziąć pod uwagę.

Jeśli pomyślimy o tym ze statystycznego punktu widzenia i możemy połączyć dane $x$ i $y$ dla obu zbiorów danych w jakiś znaczący sposób (tj $x$ i $y$ w obu zestawach są narysowane z dwóch populacji o podobnych zakresach dla dwóch zmiennych, tylko relacja między nimi różni się w dwóch populacjach), następnie możemy dopasować następujące dwa modele:

i

Gdzie $b_i$ są współczynnikami modelu, oraz $g$ to zmienna / czynnik grupujący wskazujący, do którego zbioru danych należy każda obserwacja.

Możemy użyć tabeli ANOVA lub współczynnika F do przetestowania, czy drugi, bardziej złożony model lepiej pasuje do danych niż model prostszy. Prostszy model stwierdza, że ​​nachylenia dwóch linii są takie same ($b_1$), ale linie są przesunięte względem siebie o kwotę $b_2$.

Bardziej złożony model obejmuje interakcję między nachyleniem linii a zmienną grupującą. Jeśli współczynnik tego składnika interakcji jest znacząco różny od zera lub stosunek ANOVA / F wskazuje, że bardziej złożony model lepiej pasuje do danych, to musimy odrzucić hipotezę zerową, że dwie linie są równoległe.

Oto przykład w R przy użyciu danych fikcyjnych. Po pierwsze, dane o równych nachyleniach:

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)

Co daje

> 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

Wskazując, że nie odrzucamy hipotezy zerowej o równych nachyleniach w tej próbce danych. Oczywiście chcielibyśmy się upewnić, że dysponujemy wystarczającą mocą, aby wykryć różnicę, jeśli tak naprawdę istnieje, abyśmy nie doprowadzili do błędnego odrzucenia wartości zerowej, ponieważ nasza próbka była zbyt mała dla oczekiwanego efektu.

Teraz z różnymi stokami.

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)

Co daje:

> 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

Mamy tutaj istotne dowody przeciwko hipotezie zerowej, a zatem możemy ją odrzucić na korzyść alternatywy (innymi słowy, odrzucamy hipotezę, że nachylenia dwóch linii są równe).

Warunki interakcji w dwóch modelach, które dopasowałem ($b_3xg$) podać szacunkową różnicę nachyleń dla dwóch grup. W przypadku pierwszego modelu oszacowanie różnicy nachyleń jest niewielkie (~ 0,003)

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

i a $t$-test na tym nie odrzuciłby hipotezy zerowej, że ta różnica nachyleń wynosi 0:

> 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 

Jeśli przejdziemy do modelu dopasowanego do drugiego zestawu danych, w którym zmieniliśmy nachylenia dla dwóch grup, zauważymy, że szacowana różnica nachyleń dwóch linii wynosi ~ 1 jednostka.

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

Nachylenie dla grupy „A” wynosi ~ 0,49 ( xna powyższym wyjściu), podczas gdy aby uzyskać nachylenie dla grupy „B” musimy dodać nachylenia różnicowe (podać przez słowo interakcji pamiętaj) do nachylenia grupy „A” ; ~ 0,49 + ~ 1 = ~ 1,49. Jest to dość zbliżone do podanego nachylenia dla grupy „B” wynoszącej 1,5. ZA$t$-test na tej różnicy nachyleń wskazuje również, że oszacowanie różnicy jest ograniczone od 0:

> 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

Pierwsze pytanie dotyczy w rzeczywistości geometrii. Jeśli masz dwie linie formularza:

to są równoległe, jeśli $a_1=a_2$. Jeśli więc nachylenia są równe, wówczas linie są równoległe.

W przypadku drugiego pytania skorzystaj z faktu, że $\tan \alpha=a_1$, gdzie $\alpha$ to kąt, pod jakim tworzy się linia $x$-osi i $a_1$jest nachyleniem linii. Więc

i żeby przeliczyć na stopnie, przypomnijcie sobie $2\pi=360$. Tak więc odpowiedź w stopniach będzie

Funkcja R dla $\arctan$nazywa się atan.

Przykładowy kod 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 

Ostatnia linia to stopnie.

Aktualizacja. W przypadku ujemnych wartości nachylenia konwersja na stopnie powinna przebiegać według innej reguły. Zauważ, że kąt z osią x może otrzymać wartości od 0 do 180, ponieważ zakładamy, że kąt jest powyżej osi x. Tak dla ujemnych wartości$a_1$formuła jest następująca:

Uwaga. Chociaż zabawnie było mi przywoływać trygonometrię ze szkoły średniej, naprawdę użyteczną odpowiedzią jest ta udzielona przez Gavina Simpsona. Ponieważ nachylenia linii regresji są zmiennymi losowymi, do ich porównania należy zastosować ramy hipotez statystycznych.

... w odpowiedzi na odpowiedź @mpiktas, oto jak wyodrębnić nachylenie z lmobiektu i zastosować powyższą formułę.

# 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