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?

r regression interaction linear-model

— Dail
źródło

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:

\hat{y} = b_{0} + b_{1} x + b_{2)} sol

$\hat{y} = b_0 + b_1x + b_2g$

\hat{y} = b_{0} + b_{1} x + b_{2)} sol + b_{3)} x sol

$\hat{y} = b_0 + b_1x + b_2g + b_3xg$

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

— Gavin Simpson
źródło

bardzo dziękuję za to bardzo dobre wytłumaczenie. Moim celem jest zrozumienie, czy sloper jest mniej więcej taki sam, więc myślę, że użyję ANOVA do przetestowania go.

— Dail

jeśli mam dwa wektory rozszczepienia i chciałbym porównać ich nachylenie, ale nie mam y (lm (x ~ y), jak mogę użyć ANOVA? Próbowałem anova (lm (x ~ 1), lm (y ~ 1)), ale dostaję ostrzeżenie

— Dail

Co rozumiesz przez wektory tutaj? W sensie R czy matematycznym? Jest to bardzo różni się od pytania, które postawione, więc należy rozpocząć nowe pytanie - czy nie zmieniać ten - to jest niemożliwe do przeprowadzenia follow-up o tak szerokim charakterze w komentarzach.

— Gavin Simpson

nie czekaj, muszę porównać dwa modele z ANOVA ... ok, ale jeśli utworzę model z tą formułą: x ~ 1 i inny model z y ~ 1, dostanę ostrzeżenie. Mówię w sensie R. Jak mogę to zrobić?

— Dail

@Dail, jeśli dopasowałeś dwie regresje, aby uzyskać dwa nachylenia / linie, masz dane xiy dla obu zestawów danych. Jak powiedziałem w mojej odpowiedzi, jeśli xs i ys są porównywalne w dwóch zestawach danych, możesz po prostu połączyć wszystkie dane i dodać zmienną grupującą. Mój przykład pokazuje, jak to zrobić, używając danych fikcyjnych, ale masz już dane xiy, są to dane użyte do dopasowania oddzielnych regresji.

— Gavin Simpson

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

y = a_{1} x + b_{1}

$y=a_1x+b_1$

y = {za}_{2)} x + b_{2)}

$y=a_2x+b_2$

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

α = \arctan {za}_{1}

$\alpha=\arctan a_1$

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

α = \arctan {za}_{1} \cdot \frac{360}{2) π} .

$\alpha=\arctan a_1\cdot \frac{360}{2\pi}.$

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:

α = 180 - \arctan {za}_{1} \cdot \frac{360}{2) π} .

$\alpha=180-\arctan a_1\cdot \frac{360}{2\pi}.$

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.

— mpiktas
źródło

Dziękuję Ci! jak uzyskać nachylenie z regresji? czy muszę uzyskać współczynnik i przechwycić?

— Dail

może regresja liniowa zwraca stopnie bezpośrednio z jakąś funkcją?

— Dail

mówiąc, że degress = +45 i degress = -315 nie są tym samym wierszem? kogo nie mówią o tej samej linii?

— Dail

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

— Roman Luštrik
źródło

dziękuję bardzo za przykład, w tym przypadku stopnie wynoszą -200?

— Dail

Tak. Jeśli uważasz, że problem został rozwiązany, zaznacz odpowiedź @mpiktas jako poprawną, dzięki.

— Roman Luštrik,

@ RomanLuštrik, przegapiłeś podział według

2 π

$2\pi$

— mpiktas

@ RomanLuštrik, poprawiłem kod i dodałem poprawne dane wyjściowe. Usuń korektę.

— mpiktas