Chcę wykonać bardzo prostą regresję liniową w R
. Formuła jest tak prosta, jak . Chciałbym jednak, aby nachylenie ( ) znajdowało się w przedziale, powiedzmy, między 1,4 a 1,6.
Jak można to zrobić?
źródło
Chcę wykonać bardzo prostą regresję liniową w R
. Formuła jest tak prosta, jak . Chciałbym jednak, aby nachylenie ( ) znajdowało się w przedziale, powiedzmy, między 1,4 a 1,6.
Jak można to zrobić?
Chcę wykonać ... regresję liniową w R. ... Chciałbym, aby nachylenie było w przedziale, powiedzmy, między 1,4 a 1,6. Jak można to zrobić?
(i) Prosty sposób:
dopasuj regresję. Jeśli jest w granicach, skończyłeś.
Jeśli nie ma go w granicach, ustaw nachylenie do najbliższej granicy i
oszacuj punkt przecięcia jako średnią ze wszystkich obserwacji.
(ii) Bardziej złożony sposób: rób najmniejsze kwadraty z ograniczeniami pola na zboczu; wiele procedur optymalizacyjnych implementuje ograniczenia pola, np. nlminb
(które jest dostarczane z R).
Edycja: faktycznie (jak wspomniano w poniższym przykładzie), w wanilii R, nls
można wykonywać ograniczenia pola; jak pokazano w przykładzie, jest to naprawdę bardzo łatwe do zrobienia.
Możesz użyć regresji ograniczonej bardziej bezpośrednio; Myślę, że zarówno pcls
funkcja z pakietu „mgcv”, jak i nnls
funkcja z pakietu „nnls” działają.
-
Edytuj, aby odpowiedzieć na kolejne pytanie -
Chciałem pokazać wam, jak go używać, nlminb
ponieważ jest on dostarczany z R, ale zdałem sobie sprawę, że nls
już używa tych samych procedur (procedur PORT) do implementacji ograniczonych najmniejszych kwadratów, więc mój przykład poniżej to robi.
Uwaga: w poniższym przykładzie to punkt przecięcia, a b to nachylenie (bardziej powszechna konwencja w statystykach). Po umieszczeniu tutaj zdałem sobie sprawę, że zacząłeś na odwrót; Ale zostawię przykład „zacofany” w stosunku do twojego pytania.
Najpierw skonfiguruj niektóre dane z „prawdziwym” nachyleniem w zakresie:
set.seed(seed=439812L)
x=runif(35,10,30)
y = 5.8 + 1.53*x + rnorm(35,s=5) # population slope is in range
plot(x,y)
lm(y~x)
Call:
lm(formula = y ~ x)
Coefficients:
(Intercept) x
12.681 1.217
... ale oszacowanie LS jest daleko poza tym, spowodowane tylko przypadkową zmiennością. Wykorzystajmy więc regresję ograniczoną w nls
:
nls(y~a+b*x,algorithm="port",
start=c(a=0,b=1.5),lower=c(a=-Inf,b=1.4),upper=c(a=Inf,b=1.6))
Nonlinear regression model
model: y ~ a + b * x
data: parent.frame()
a b
9.019 1.400
residual sum-of-squares: 706.2
Algorithm "port", convergence message: both X-convergence and relative convergence (5)
Jak widzisz, masz nachylenie tuż przy granicy. Jeśli przekażesz dopasowany model summary
, wygeneruje on nawet standardowe błędy i wartości t, ale nie jestem pewien, czy są one znaczące / interpretowalne.
b=1.4
c(a=mean(y-x*b),b=b)
a b
9.019376 1.400000
To ta sama ocena ...
Na poniższym wykresie niebieska linia to najmniej kwadratów, a czerwona linia to najmniejszych kwadratów:
nls
do tego celu.Druga metoda Glen_b, polegająca na użyciu najmniejszych kwadratów z ograniczeniem pola, może być łatwiej zaimplementowana poprzez regresję grzbietu. Rozwiązanie regresji grzbietu można postrzegać jako Lagrangian dla regresji z ograniczeniem wielkości normy wektora ciężaru (a tym samym jej nachylenia). Zatem zgodnie z sugestią Whubera poniżej podejściem byłoby odjęcie trendu (1,6 + 1,4) / 2 = 1,5, a następnie zastosowanie regresji grzbietu i stopniowe zwiększanie parametru grzbietu, aż wielkość nachylenia będzie mniejsza lub równa 0,1.
Zaletą tego podejścia jest to, że nie są wymagane żadne fantazyjne narzędzia do optymalizacji, wystarczy regresson ridge, który jest już dostępny w R (i wielu innych pakietach).
Jednak proste rozwiązanie Glen_b (i) wydaje mi się rozsądne (+1)
źródło
Ten wynik nadal da wiarygodne przedziały parametrów będących przedmiotem zainteresowania (oczywiście znaczenie tych przedziałów będzie oparte na racjonalności twoich wcześniejszych informacji o nachyleniu).
źródło
Innym podejściem może być przeformułowanie regresji jako problemu optymalizacji i użycie optymalizatora. Nie jestem pewien, czy można to przeformułować w ten sposób, ale pomyślałem o tym pytaniu, czytając ten blog na optymalizatorach R:
http://zoonek.free.fr/blosxom/R/2012-06-01_Optimization.html
źródło