Regresja liniowa z ograniczeniem nachylenia

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.$y = ax + b$$a$

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.$(y - ax)$

(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, nlsmoż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 pclsfunkcja z pakietu „mgcv”, jak i nnlsfunkcja z pakietu „nnls” działają.

-

Edytuj, aby odpowiedzieć na kolejne pytanie -

Chciałem pokazać wam, jak go używać, nlminbponieważ jest on dostarczany z R, ale zdałem sobie sprawę, że nlsjuż 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.$a$$b$

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.

$y-bx$

 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:

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)

$a$

$a$

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

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