Jak obliczyć, czy moja regresja liniowa ma statystycznie istotną różnicę od znanej linii teoretycznej?


14

Mam pewne dane, które pasują do z grubsza liniowej linii:

wprowadź opis zdjęcia tutaj

Kiedy wykonuję regresję liniową tych wartości, otrzymuję równanie liniowe:

y=0,997x-0,0136

W idealnym świecie równanie powinno wynosić .y=x

Oczywiście moje wartości liniowe są zbliżone do tego ideału, ale nie do końca. Moje pytanie brzmi: jak mogę ustalić, czy ten wynik jest statystycznie istotny?

Czy wartość 0,997 znacznie różni się od 1? Czy -0,01 znacznie różni się od 0? Czy też są statystycznie takie same i mogę stwierdzić, że przy pewnym rozsądnym poziomie ufności?y=x

Jakiego testu statystycznego mogę użyć?

Dzięki


1
Możesz obliczyć, czy istnieje różnica istotna statystycznie, ale należy pamiętać, że nie oznacza to, czy nie ma różnicy. Możesz mieć pewność co do znaczenia, kiedy fałszujesz hipotezę zerową, ale jeśli nie fałszujesz hipotezy zerowej, może to być albo (1), że rzeczywiście hipoteza zerowa jest poprawna (2) test nie był skuteczny z powodu niskiej liczby próbek (3) test nie był skuteczny z powodu złej hipotezy alternatywnej (3b) fałszywej miary istotności statystycznej z powodu błędnego przedstawienia niedeterministycznej części modelu.
Sextus Empiricus

Dla mnie twoje dane nie wyglądają jak y = x + biały szum. Czy możesz powiedzieć o tym więcej? (test na założenie, że występuje taki hałas, może nie „dostrzec” znaczącej różnicy, bez względu na to, jak duża jest próbka, nawet jeśli istnieje ogromna różnica między danymi a linią y = x, tylko dlatego, że jesteś tylko porównanie z innymi liniami y = a + bx, co może nie być właściwym i najmocniejszym porównaniem)
Sextus Empiricus

Ponadto, jaki jest cel określenia znaczenia. Widzę, że wiele odpowiedzi sugeruje użycie poziomu alfa 5% (95% przedziały ufności). Jest to jednak bardzo arbitralne. Bardzo trudno jest postrzegać istotność statystyczną jako zmienną binarną (obecną lub nieobecną). Odbywa się to za pomocą takich zasad, jak standardowe poziomy alfa, ale jest to arbitralne i prawie bez znaczenia. Jeśli podasz kontekst, to użycie określonego poziomu odcięcia w celu podjęcia decyzji (zmienna binarna) na podstawie poziomu istotności ( nie zmiennej binarnej), wtedy koncepcja taka jak znaczenie binarne ma większy sens.
Sextus Empiricus

1
Jakiego rodzaju „regresję liniową” wykonujesz? Zwykle ktoś uważa, że ​​omawiasz zwykłą regresję najmniejszych kwadratów (z terminem przechwytującym), ale w takim przypadku, ponieważ oba zestawy reszt będą miały środki zerowe (dokładnie), przecięcie w regresji między resztami również powinno wynosić zero (dokładnie ). Ponieważ tak nie jest, dzieje się tutaj coś innego. Czy możesz podać informacje na temat tego, co robisz i dlaczego?
whuber

Wygląda to podobnie do problemu w pomiarze sprawdzania, czy dwa systemy dają taki sam wynik. Spróbuj spojrzeć na nijaką fabułę dla jakiegoś materiału.
mdewey

Odpowiedzi:


17

Ten typ sytuacji można rozwiązać za pomocą standardowego testu F dla modeli zagnieżdżonych . Ponieważ chcesz przetestować oba parametry na modelu zerowym ze stałymi parametrami, twoje hipotezy są następujące:

H.0:β=[01]H.ZA:β[01].

Test F obejmuje dopasowanie obu modeli i porównanie ich rezydualnej sumy kwadratów, które są:

S.S.mi0=ja=1n(yja-xja)2)S.S.miZA=ja=1n(yja-β^0-β^1xja)2)

Statystyka testu to:

fafa(y,x)=n-2)2)S.S.mi0-S.S.miZAS.S.miZA.

Odpowiednia wartość p wynosi:

pp(y,x)=F(y,x)F-Dist(r|2,n2) dr.


Implementacja w R: Załóżmy, że twoje dane są w ramce danych wywoływanej DATAzmiennymi o nazwie yi x. Test F można wykonać ręcznie za pomocą następującego kodu. W symulowanych próbnych danych, które wykorzystałem, możesz zobaczyć, że szacowane współczynniki są zbliżone do tych w hipotezie zerowej, a wartość p testu nie wykazuje istotnych dowodów na fałszowanie hipotezy zerowej, że prawdziwą funkcją regresji jest funkcja tożsamości.

#Generate mock data (you can substitute your data if you prefer)
set.seed(12345);
n    <- 1000;
x    <- rnorm(n, mean = 0, sd = 5);
e    <- rnorm(n, mean = 0, sd = 2/sqrt(1+abs(x)));
y    <- x + e;
DATA <- data.frame(y = y, x = x);

#Fit initial regression model
MODEL <- lm(y ~ x, data = DATA);

#Calculate test statistic
SSE0   <- sum((DATA$y-DATA$x)^2);
SSEA   <- sum(MODEL$residuals^2);
F_STAT <- ((n-2)/2)*((SSE0 - SSEA)/SSEA);
P_VAL  <- pf(q = F_STAT, df1 = 2, df2 = n-2, lower.tail = FALSE);

#Plot the data and show test outcome
plot(DATA$x, DATA$y,
     main = 'All Residuals',
     sub  = paste0('(Test against identity function - F-Stat = ',
            sprintf("%.4f", F_STAT), ', p-value = ', sprintf("%.4f", P_VAL), ')'),
     xlab = 'Dataset #1 Normalized residuals',
     ylab = 'Dataset #2 Normalized residuals');
abline(lm(y ~ x, DATA), col = 'red', lty = 2, lwd = 2);

Dane summarywyjściowe i plotdla tych danych wyglądają następująco:

summary(MODEL);

Call:
lm(formula = y ~ x, data = DATA)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.8276 -0.6742  0.0043  0.6703  5.1462 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.02784    0.03552  -0.784    0.433    
x            1.00507    0.00711 141.370   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.122 on 998 degrees of freedom
Multiple R-squared:  0.9524,    Adjusted R-squared:  0.9524 
F-statistic: 1.999e+04 on 1 and 998 DF,  p-value: < 2.2e-16

F_STAT;
[1] 0.5370824

P_VAL;
[1] 0.5846198

wprowadź opis zdjęcia tutaj


x

1
Tak, dobrze zauważony. Symulowane dane nie wykorzystują standardowej regresji liniowej homoskedastycznej. Użyłem heteroscedastyczności w symulacji, aby z grubsza naśladować wzór danych na wykresie pokazanym przez OP. (I myślę, że wykonałem całkiem cholernie dobrą robotę!) W takim przypadku dopasowuję standardowy homoskedastyczny model liniowy do danych symulowanych, które nie zostały wygenerowane z tego modelu. Nadal jest to uzasadnione - można symulować dane z jednego modelu, a następnie dopasować go do innego, aby zobaczyć, co się pojawi.
Przywróć Monikę

1
sd = 2/sqrt(1+abs(x))yxy=xxy=xy=x+mi
Sextus Empiricus

1
To prawda, ale przenosi cię na terytorium modeli z błędami w zmiennych, co czyni to bardziej skomplikowanym. Myślę, że OP chce po prostu użyć standardowej regresji liniowej w tym przypadku.
Przywróć Monikę

Zgadzam się, że jest to sidena, ale mimo to ważna. Prostota pytania mnie intryguje (w różnych punktach), a także martwi mnie, ponieważ może to być zbyt prosta reprezentacja. Oczywiście zależy to od tego, co faktycznie się chce osiągnąć („wszystkie modele są błędne ...”), ale ta prosta reprezentacja może stać się standardem, a złożone dodatkowe pytania, o których należy pamiętać, zostaną zapomniane, a nawet jedno nigdy nie zaczyna o tym myśleć (odwoływanie się do 95% CI w innych odpowiedziach jest przykładem takiego standardu, którego ludzie ślepo przestrzegają).
Sextus Empiricus

5

Oto fajna metoda graficzna, którą przytoczyłem z doskonałej książki Juliana Faraway'a „Modele liniowe z R (wydanie drugie)”. To równoczesne 95% przedziały ufności dla punktu przecięcia i nachylenia, wykreślone jako elipsa.

Dla ilustracji stworzyłem 500 obserwacji ze zmienną „x” o rozkładzie N (średnia = 10, sd = 5), a następnie zmienną „y”, której rozkład wynosi N (średnia = x, sd = 2). Daje to korelację nieco powyżej 0,9, która może nie być tak ścisła jak twoje dane.

Możesz sprawdzić elipsę, aby zobaczyć, czy punkt (punkt przecięcia = 0, nachylenie = 1) mieści się w tym lub przed tym przedziałem ufności.

library(tidyverse)
library(ellipse)
#> 
#> Attaching package: 'ellipse'
#> The following object is masked from 'package:graphics':
#> 
#>     pairs

set.seed(50)
dat <- data.frame(x=rnorm(500,10,5)) %>% mutate(y=rnorm(n(),x,2))

lmod1 <- lm(y~x,data=dat)
summary(lmod1)
#> 
#> Call:
#> lm(formula = y ~ x, data = dat)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -6.9652 -1.1796 -0.0576  1.2802  6.0212 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)  0.24171    0.20074   1.204    0.229    
#> x            0.97753    0.01802  54.246   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 2.057 on 498 degrees of freedom
#> Multiple R-squared:  0.8553, Adjusted R-squared:  0.855 
#> F-statistic:  2943 on 1 and 498 DF,  p-value: < 2.2e-16

cor(dat$y,dat$x)
#> [1] 0.9248032

plot(y~x,dat)
abline(0,1)


confint(lmod1)
#>                  2.5 %    97.5 %
#> (Intercept) -0.1526848 0.6361047
#> x            0.9421270 1.0129370

plot(ellipse(lmod1,c("(Intercept)","x")),type="l")
points(coef(lmod1)["(Intercept)"],coef(lmod1)["x"],pch=19)

abline(v=confint(lmod1)["(Intercept)",],lty=2)
abline(h=confint(lmod1)["x",],lty=2)

points(0,1,pch=1,size=3)
#> Warning in plot.xy(xy.coords(x, y), type = type, ...): "size" is not a
#> graphical parameter

abline(v=0,lty=10)
abline(h=0,lty=10)

Utworzono 21.01.2019 przez pakiet reprezentx (v0.2.1)


1

Można obliczyć współczynniki za pomocą n próbek próbnych. Spowoduje to prawdopodobnie uzyskanie normalnych wartości współczynników rozkładu (centralne twierdzenie graniczne). Dzięki temu możesz następnie skonstruować (np. 95%) przedział ufności z wartościami t (n-1 stopni swobody) wokół średniej. Jeśli twój CI nie zawiera 1 (0), jest statystycznie istotny różny, a dokładniej: Możesz odrzucić hipotezę zerową o równym nachyleniu.


Jak to tutaj sformułowałeś, testuje tylko dwie hipotezy osobno, ale potrzebujesz wspólnego testu.
kjetil b halvorsen

0

β0=0β1=1


1
Ale potrzebny jest wspólny test, jak w przypadku innych odpowiedzi.
kjetil b halvorsen

@kjetilbhalvorsen Zdałem sobie sprawę, że dzisiaj rano się myliłem, czytając inne odpowiedzi. Usunę to.
RScrlli

0

Powinieneś dopasować regresję liniową i sprawdzić 95% przedziały ufności dla dwóch parametrów. Jeżeli CI nachylenia obejmuje 1, a CI przesunięcia obejmuje 0, to dwustronny test jest nieznaczny w przybliżeniu. na poziomie (95%) ^ 2 - ponieważ stosujemy dwa osobne testy, ryzyko typ-I wzrasta.

Za pomocą R:

fit = lm(Y ~ X)
confint(fit)

lub używasz

summary(fit)

i samodzielnie oblicz 2 przedziały sigma.

Korzystając z naszej strony potwierdzasz, że przeczytałeś(-aś) i rozumiesz nasze zasady używania plików cookie i zasady ochrony prywatności.
Licensed under cc by-sa 3.0 with attribution required.