Régression linéaire avec contrainte de pente

Je veux effectuer une régression linéaire très simple dans R. La formule est aussi simple que . Cependant, je voudrais que la pente ( ) soit à l'intérieur d'un intervalle, disons, entre 1,4 et 1,6.$y = ax + b$$a$

Comment cela peut-il être fait?

Je veux effectuer ... une régression linéaire dans R. ... Je voudrais que la pente soit à l'intérieur d'un intervalle, disons, entre 1,4 et 1,6. Comment cela peut-il être fait?

(i) Manière simple:

• adapter la régression. Si c'est dans les limites, vous avez terminé.

• S'il n'est pas dans les limites, définissez la pente sur la limite la plus proche, et

• estimer l'ordonnée à l'origine comme la moyenne de sur toutes les observations.$(y - ax)$

(ii) Méthode plus complexe: faire les moindres carrés avec des contraintes de boîte sur la pente; de nombreuses routines d'optimisation implémentent des contraintes de boîte, par exemple nlminb(qui vient avec R).

Edit: en fait (comme mentionné dans l'exemple ci-dessous), dans vanilla R, nlspeut faire des contraintes de boîte; comme le montre l'exemple, c'est vraiment très facile à faire.

Vous pouvez utiliser la régression contrainte plus directement; Je pense que la pclsfonction du package "mgcv" et la nnlsfonction du package "nnls" font toutes les deux.

-

Modifier pour répondre à la question de suivi -

J'allais vous montrer comment l'utiliser avec nlminbpuisque cela vient avec R, mais j'ai réalisé qu'il nlsutilise déjà les mêmes routines (les routines PORT) pour implémenter les moindres carrés contraints, donc mon exemple ci-dessous le fait.

NB: dans mon exemple ci-dessous, est l'ordonnée à l'origine et est la pente (la convention la plus courante dans les statistiques). J'ai réalisé après l'avoir mis ici que vous avez commencé l'inverse; Je vais cependant laisser l'exemple «en arrière» par rapport à votre question.$a$$b$

Tout d'abord, configurez certaines données avec la «vraie» pente à l'intérieur de la plage:

 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  

... mais l'estimation LS est bien en dehors d'elle, juste causée par une variation aléatoire. Permet donc d'utiliser la régression contrainte dans 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)

Comme vous le voyez, vous obtenez une pente à droite sur la frontière. Si vous lui passez le modèle ajusté, summarycela produira même des erreurs standard et des valeurs t, mais je ne sais pas à quel point elles sont significatives / interprétables.

Alors, comment ma suggestion (1) se compare-t-elle? (c.-à-d. régler la pente à la limite la plus proche et faire la moyenne des résidus pour estimer l'ordonnée à l'origine)$y-bx$

 b=1.4
 c(a=mean(y-x*b),b=b)
       a        b 
9.019376 1.400000

C'est la même estimation ...

Dans le graphique ci-dessous, la ligne bleue correspond aux moindres carrés et la ligne rouge aux moindres carrés contraints:

La deuxième méthode de Glen_b, utilisant les moindres carrés avec une contrainte de boîte peut être plus facilement implémentée via la régression de crête. La solution à la régression des crêtes peut être considérée comme le lagrangien pour une régression avec une borne sur l'amplitude de la norme du vecteur de poids (et donc sa pente). Donc, suivant la suggestion de whuber ci-dessous, l'approche serait de soustraire une tendance de (1,6 + 1,4) / 2 = 1,5, puis d'appliquer une régression de crête et d'augmenter progressivement le paramètre de crête jusqu'à ce que l'amplitude de la pente soit inférieure ou égale à 0,1.

L'avantage de cette approche est qu'aucun outil d'optimisation sophistiqué n'est requis, juste ridge regresson, qui est déjà disponible dans R (et de nombreux autres packages).

Cependant la solution simple de Glen_b (i) me semble sensée (+1)

$a$

$a$

Ce résultat donnera toujours des intervalles crédibles des paramètres d'intérêt (bien sûr, la signification de ces intervalles sera basée sur le caractère raisonnable de vos informations antérieures sur la pente).

Une autre approche pourrait être de reformuler votre régression en tant que problème d'optimisation et d'utiliser un optimiseur. Je ne sais pas si cela peut être reformulé de cette façon, mais j'ai pensé à cette question en lisant cet article de blog sur les optimiseurs R:

http://zoonek.free.fr/blosxom/R/2012-06-01_Optimization.html