<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/convex-regression-model/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://freakonometrics.hypotheses.org/53384"> R-english – Freakonometrics</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>This morning during the lecture on nonlinear regression, I mentioned (very) briefly the case of convex regression. Since I forgot to mention the codes in R, I will publish them here. Assume that <img alt="y_i=m(\mathbf{x}_i)+\varepsilon_i" title="y_i=m(\mathbf{x}_i)+\varepsilon_i" src="https://s0.wp.com/latex.php?latex=y_i%3Dm%28%5Cmathbf%7Bx%7D_i%29%2B%5Cvarepsilon_i&bg=ffffff&fg=000&s=0" class="latex"> where <img alt="m:\mathbb{R}^d\rightarrow \mathbb{R}" title="m:\mathbb{R}^d\rightarrow \mathbb{R}" src="https://s0.wp.com/latex.php?latex=m%3A%5Cmathbb%7BR%7D%5Ed%5Crightarrow+%5Cmathbb%7BR%7D&bg=ffffff&fg=000&s=0" class="latex"> is some convex function. </p>
<p>Then <img alt="m" title="m" src="https://s0.wp.com/latex.php?latex=m&bg=ffffff&fg=000&s=0" class="latex"> is convex if and only if <img alt="\forall\mathbf{x}_1,\mathbf{x}_2\in\mathbb{R}^d" title="\forall\mathbf{x}_1,\mathbf{x}_2\in\mathbb{R}^d" src="https://s0.wp.com/latex.php?latex=%5Cforall%5Cmathbf%7Bx%7D_1%2C%5Cmathbf%7Bx%7D_2%5Cin%5Cmathbb%7BR%7D%5Ed&bg=ffffff&fg=000&s=0" class="latex">, <img alt="\forall t\in[0,1]" title="\forall t\in[0,1]" src="https://s0.wp.com/latex.php?latex=%5Cforall+t%5Cin%5B0%2C1%5D&bg=ffffff&fg=000&s=0" class="latex">, <img alt="m(t\mathbf{x}_1+[1-t]\mathbf{x}_2) \leq tm(\mathbf{x}_1)+[1-t]m(\mathbf{x}_2)" title="m(t\mathbf{x}_1+[1-t]\mathbf{x}_2) \leq tm(\mathbf{x}_1)+[1-t]m(\mathbf{x}_2)" src="https://s0.wp.com/latex.php?latex=m%28t%5Cmathbf%7Bx%7D_1%2B%5B1-t%5D%5Cmathbf%7Bx%7D_2%29+%5Cleq+tm%28%5Cmathbf%7Bx%7D_1%29%2B%5B1-t%5Dm%28%5Cmathbf%7Bx%7D_2%29&bg=ffffff&fg=000&s=0" class="latex"><a rel="nofollow" href="https://www.jstor.org/stable/2281132?seq=1#page_scan_tab_contents" target="_blank">Hidreth (1954)</a> proved that if<img alt=" m^\star=\underset{m \text{ convex}}{\text{argmin}}\left\lbrace\sum_{i=1}^n \big(y_i-m(\mathbf{x_i})\big)^2\right\rbrace" title=" m^\star=\underset{m \text{ convex}}{\text{argmin}}\left\lbrace\sum_{i=1}^n \big(y_i-m(\mathbf{x_i})\big)^2\right\rbrace" src="https://s0.wp.com/latex.php?latex=+m%5E%5Cstar%3D%5Cunderset%7Bm+%5Ctext%7B+convex%7D%7D%7B%5Ctext%7Bargmin%7D%7D%5Cleft%5Clbrace%5Csum_%7Bi%3D1%7D%5En+%5Cbig%28y_i-m%28%5Cmathbf%7Bx_i%7D%29%5Cbig%29%5E2%5Cright%5Crbrace&bg=ffffff&fg=000&s=0" class="latex">then <img alt="\mathbf{\theta}^\star=(m^\star(\mathbf{x_1}),\cdots,m^\star(\mathbf{x_n}))" title="\mathbf{\theta}^\star=(m^\star(\mathbf{x_1}),\cdots,m^\star(\mathbf{x_n}))" src="https://s0.wp.com/latex.php?latex=%5Cmathbf%7B%5Ctheta%7D%5E%5Cstar%3D%28m%5E%5Cstar%28%5Cmathbf%7Bx_1%7D%29%2C%5Ccdots%2Cm%5E%5Cstar%28%5Cmathbf%7Bx_n%7D%29%29&bg=ffffff&fg=000&s=0" class="latex"> is unique.</p>
<p>Let <img alt="\mathbf{y}=\mathbf{\theta}+\mathbf{\varepsilon}" title="\mathbf{y}=\mathbf{\theta}+\mathbf{\varepsilon}" src="https://s0.wp.com/latex.php?latex=%5Cmathbf%7By%7D%3D%5Cmathbf%7B%5Ctheta%7D%2B%5Cmathbf%7B%5Cvarepsilon%7D&bg=ffffff&fg=000&s=0" class="latex">, then <img alt="\mathbf{\theta}^\star=\underset{\mathbf{\theta}\in \mathcal{K}}{\text{argmin}}\left\lbrace\sum_{i=1}^n \big(y_i-\theta_i)\big)^2\right\rbrace" title="\mathbf{\theta}^\star=\underset{\mathbf{\theta}\in \mathcal{K}}{\text{argmin}}\left\lbrace\sum_{i=1}^n \big(y_i-\theta_i)\big)^2\right\rbrace" src="https://s0.wp.com/latex.php?latex=%5Cmathbf%7B%5Ctheta%7D%5E%5Cstar%3D%5Cunderset%7B%5Cmathbf%7B%5Ctheta%7D%5Cin+%5Cmathcal%7BK%7D%7D%7B%5Ctext%7Bargmin%7D%7D%5Cleft%5Clbrace%5Csum_%7Bi%3D1%7D%5En+%5Cbig%28y_i-%5Ctheta_i%29%5Cbig%29%5E2%5Cright%5Crbrace&bg=ffffff&fg=000&s=0" class="latex">where<img alt="\mathcal{K}=\{\mathbf{\theta}\in\mathbb{R}^n:\exists m\text{ convex },m(\mathbf{x}_i)=\theta_i\}" title="\mathcal{K}=\{\mathbf{\theta}\in\mathbb{R}^n:\exists m\text{ convex },m(\mathbf{x}_i)=\theta_i\}" src="https://s0.wp.com/latex.php?latex=%5Cmathcal%7BK%7D%3D%5C%7B%5Cmathbf%7B%5Ctheta%7D%5Cin%5Cmathbb%7BR%7D%5En%3A%5Cexists+m%5Ctext%7B+convex+%7D%2Cm%28%5Cmathbf%7Bx%7D_i%29%3D%5Ctheta_i%5C%7D&bg=ffffff&fg=000&s=0" class="latex">. I.e. <img alt="\mathbf{\theta}^\star" title="\mathbf{\theta}^\star" src="https://s0.wp.com/latex.php?latex=%5Cmathbf%7B%5Ctheta%7D%5E%5Cstar&bg=ffffff&fg=000&s=0" class="latex"> is the projection of <img alt="\mathbf{y}" title="\mathbf{y}" src="https://s0.wp.com/latex.php?latex=%5Cmathbf%7By%7D&bg=ffffff&fg=000&s=0" class="latex"> onto the (closed) convex cone <img alt="\mathcal{K}" title="\mathcal{K}" src="https://s0.wp.com/latex.php?latex=%5Cmathcal%7BK%7D&bg=ffffff&fg=000&s=0" class="latex">. The projection theorem gives existence and unicity.</p>
<p>For convenience, in the application, we will consider the real-valued case, <img alt="m:\mathbb{R}\rightarrow \mathbb{R}" title="m:\mathbb{R}\rightarrow \mathbb{R}" src="https://s0.wp.com/latex.php?latex=m%3A%5Cmathbb%7BR%7D%5Crightarrow+%5Cmathbb%7BR%7D&bg=ffffff&fg=000&s=0" class="latex">, i.e. <img alt="y_i=m(x_i)+\varepsilon_i" title="y_i=m(x_i)+\varepsilon_i" src="https://s0.wp.com/latex.php?latex=y_i%3Dm%28x_i%29%2B%5Cvarepsilon_i&bg=ffffff&fg=000&s=0" class="latex">. Assume that observations are ordered <img alt="x_1\leq x_2\leq\cdots \leq x_n" title="x_1\leq x_2\leq\cdots \leq x_n" src="https://s0.wp.com/latex.php?latex=x_1%5Cleq+x_2%5Cleq%5Ccdots+%5Cleq+x_n&bg=ffffff&fg=000&s=0" class="latex">. Here <img alt="\mathcal{K}=\left\lbrace\mathbf{\theta}\in\mathbb{R}^n:\frac{\theta_2-\theta_1}{x_2-x_1}\leq \frac{\theta_3-\theta_2}{x_3-x_2}\leq \cdots \leq \frac{\theta_n-\theta_{n-1}}{x_n-x_{n-1}}\right\rbrace" title="\mathcal{K}=\left\lbrace\mathbf{\theta}\in\mathbb{R}^n:\frac{\theta_2-\theta_1}{x_2-x_1}\leq \frac{\theta_3-\theta_2}{x_3-x_2}\leq \cdots \leq \frac{\theta_n-\theta_{n-1}}{x_n-x_{n-1}}\right\rbrace" src="https://s0.wp.com/latex.php?latex=%5Cmathcal%7BK%7D%3D%5Cleft%5Clbrace%5Cmathbf%7B%5Ctheta%7D%5Cin%5Cmathbb%7BR%7D%5En%3A%5Cfrac%7B%5Ctheta_2-%5Ctheta_1%7D%7Bx_2-x_1%7D%5Cleq+%5Cfrac%7B%5Ctheta_3-%5Ctheta_2%7D%7Bx_3-x_2%7D%5Cleq+%5Ccdots+%5Cleq+%5Cfrac%7B%5Ctheta_n-%5Ctheta_%7Bn-1%7D%7D%7Bx_n-x_%7Bn-1%7D%7D%5Cright%5Crbrace&bg=ffffff&fg=000&s=0" class="latex"></p>
<p>Hence, quadratic program with <img alt="n-2" title="n-2" src="https://s0.wp.com/latex.php?latex=n-2&bg=ffffff&fg=000&s=0" class="latex"> linear constraints.</p>
<p><img alt="m^\star" title="m^\star" src="https://s0.wp.com/latex.php?latex=m%5E%5Cstar&bg=ffffff&fg=000&s=0" class="latex"> is a piecewise linear function (interpolation of consecutive pairs <img alt="(x_i,\theta_i^\star)" title="(x_i,\theta_i^\star)" src="https://s0.wp.com/latex.php?latex=%28x_i%2C%5Ctheta_i%5E%5Cstar%29&bg=ffffff&fg=000&s=0" class="latex">).</p>
<p>If <img alt="m" title="m" src="https://s0.wp.com/latex.php?latex=m&bg=ffffff&fg=000&s=0" class="latex"> is differentiable, <img alt="m" title="m" src="https://s0.wp.com/latex.php?latex=m&bg=ffffff&fg=000&s=0" class="latex"> is convex if <img alt="m(\mathbf{x})+ \nabla m(\mathbf{x})^{\text{T}}\cdot[\mathbf{y}-\mathbf{x}] \leq m(\mathbf{y})" title="m(\mathbf{x})+ \nabla m(\mathbf{x})^{\text{T}}\cdot[\mathbf{y}-\mathbf{x}] \leq m(\mathbf{y})" src="https://s0.wp.com/latex.php?latex=m%28%5Cmathbf%7Bx%7D%29%2B+%5Cnabla+m%28%5Cmathbf%7Bx%7D%29%5E%7B%5Ctext%7BT%7D%7D%5Ccdot%5B%5Cmathbf%7By%7D-%5Cmathbf%7Bx%7D%5D+%5Cleq+m%28%5Cmathbf%7By%7D%29&bg=ffffff&fg=000&s=0" class="latex"></p>
<p>More generally, if <img alt="m" title="m" src="https://s0.wp.com/latex.php?latex=m&bg=ffffff&fg=000&s=0" class="latex"> is convex, then there exists <img alt="\xi_{\mathbf{x}}\in\mathbb{R}^n" title="\xi_{\mathbf{x}}\in\mathbb{R}^n" src="https://s0.wp.com/latex.php?latex=%5Cxi_%7B%5Cmathbf%7Bx%7D%7D%5Cin%5Cmathbb%7BR%7D%5En&bg=ffffff&fg=000&s=0" class="latex"> such that <img alt="m(\mathbf{x})+ \xi_{\mathbf{x}}^{\text{ T}}\cdot[\mathbf{y}-\mathbf{x}] \leq m(\mathbf{y})" title="m(\mathbf{x})+ \xi_{\mathbf{x}}^{\text{ T}}\cdot[\mathbf{y}-\mathbf{x}] \leq m(\mathbf{y})" src="https://s0.wp.com/latex.php?latex=m%28%5Cmathbf%7Bx%7D%29%2B+%5Cxi_%7B%5Cmathbf%7Bx%7D%7D%5E%7B%5Ctext%7B+T%7D%7D%5Ccdot%5B%5Cmathbf%7By%7D-%5Cmathbf%7Bx%7D%5D+%5Cleq+m%28%5Cmathbf%7By%7D%29&bg=ffffff&fg=000&s=0" class="latex"><br>
<img alt="\xi_{\mathbf{x}}" title="\xi_{\mathbf{x}}" src="https://s0.wp.com/latex.php?latex=%5Cxi_%7B%5Cmathbf%7Bx%7D%7D&bg=ffffff&fg=000&s=0" class="latex"> is a subgradient of <img alt="m" title="m" src="https://s0.wp.com/latex.php?latex=m&bg=ffffff&fg=000&s=0" class="latex"> at <img alt="{\mathbf{x}}" title="{\mathbf{x}}" src="https://s0.wp.com/latex.php?latex=%7B%5Cmathbf%7Bx%7D%7D&bg=ffffff&fg=000&s=0" class="latex">. And then <img alt="\partial m(\mathbf{x})=\big\lbrace m(\mathbf{x})+ \xi^{\text{ T}}\cdot[\mathbf{y}-\mathbf{x}] \leq m(\mathbf{y}),\forall \mathbf{y}\in\mathbb{R}^n\big\rbrace" title="\partial m(\mathbf{x})=\big\lbrace m(\mathbf{x})+ \xi^{\text{ T}}\cdot[\mathbf{y}-\mathbf{x}] \leq m(\mathbf{y}),\forall \mathbf{y}\in\mathbb{R}^n\big\rbrace" src="https://s0.wp.com/latex.php?latex=%5Cpartial+m%28%5Cmathbf%7Bx%7D%29%3D%5Cbig%5Clbrace+m%28%5Cmathbf%7Bx%7D%29%2B+%5Cxi%5E%7B%5Ctext%7B+T%7D%7D%5Ccdot%5B%5Cmathbf%7By%7D-%5Cmathbf%7Bx%7D%5D+%5Cleq+m%28%5Cmathbf%7By%7D%29%2C%5Cforall+%5Cmathbf%7By%7D%5Cin%5Cmathbb%7BR%7D%5En%5Cbig%5Crbrace&bg=ffffff&fg=000&s=0" class="latex"></p>
<p>Hence, <img alt="\mathbf{\theta}^\star" title="\mathbf{\theta}^\star" src="https://s0.wp.com/latex.php?latex=%5Cmathbf%7B%5Ctheta%7D%5E%5Cstar&bg=ffffff&fg=000&s=0" class="latex"> is solution of <img alt="\text{argmin}\big\lbrace\|\mathbf{y}-\mathbf{\theta}\|^2\big\rbrace" title="\text{argmin}\big\lbrace\|\mathbf{y}-\mathbf{\theta}\|^2\big\rbrace" src="https://s0.wp.com/latex.php?latex=%5Ctext%7Bargmin%7D%5Cbig%5Clbrace%5C%7C%5Cmathbf%7By%7D-%5Cmathbf%7B%5Ctheta%7D%5C%7C%5E2%5Cbig%5Crbrace&bg=ffffff&fg=000&s=0" class="latex"><img alt="\text{subject to }\theta_i+\xi_i^{\text{ T}}[\mathbf{x}_j-\mathbf{x}_i]\leq\mathbf{\theta}_j,~\forall i,j" title="\text{subject to }\theta_i+\xi_i^{\text{ T}}[\mathbf{x}_j-\mathbf{x}_i]\leq\mathbf{\theta}_j,~\forall i,j" src="https://s0.wp.com/latex.php?latex=%5Ctext%7Bsubject+to+%7D%5Ctheta_i%2B%5Cxi_i%5E%7B%5Ctext%7B+T%7D%7D%5B%5Cmathbf%7Bx%7D_j-%5Cmathbf%7Bx%7D_i%5D%5Cleq%5Cmathbf%7B%5Ctheta%7D_j%2C%7E%5Cforall+i%2Cj&bg=ffffff&fg=000&s=0" class="latex"> and <img alt="\xi_1,\cdots,\xi_n\in\mathbb{R}^n" title="\xi_1,\cdots,\xi_n\in\mathbb{R}^n" src="https://s0.wp.com/latex.php?latex=%5Cxi_1%2C%5Ccdots%2C%5Cxi_n%5Cin%5Cmathbb%7BR%7D%5En&bg=ffffff&fg=000&s=0" class="latex">. Now, to do it for real, use <a rel="nofollow" href="https://cran.r-project.org/web/packages/cobs/index.html" target="_blank">cobs</a> package for constrained (b)splines regression,</p>
<div><table><tr><td><pre>1
</pre></td><td><pre><span>library</span><span>(</span>cobs<span>)</span></pre></td></tr></table></div>
<p>To get a convex regression, use</p>
<div><table><tr><td><pre>1
2
3
4
5
</pre></td><td><pre><span>plot</span><span>(</span><span>cars</span><span>)</span>
x <span>=</span> <span>cars</span>$speed
y <span>=</span> <span>cars</span>$dist
rc <span>=</span> conreg<span>(</span>x,y,convex<span>=</span>TRUE<span>)</span>
<span>lines</span><span>(</span>rc, <span>col</span> <span>=</span> <span>2</span><span>)</span></pre></td></tr></table></div>
<p><img alt="" src="https://i2.wp.com/f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2018/07/convex1.png?w=456&ssl=1" data-recalc-dims="1"><br>
Here we can get the values of the knots</p>
<div><table><tr><td><pre>1
2
3
4
5
6
7
</pre></td><td><pre>rc
 
Call<span>:</span>  conreg<span>(</span>x <span>=</span> x, y <span>=</span> y, convex <span>=</span> TRUE<span>)</span> 
Convex regression<span>:</span> From <span>19</span> separated x<span>-</span>values, using <span>5</span> inner <span>knots</span>,
     <span>7</span>,    <span>8</span>,    <span>9</span>,   <span>20</span>,   <span>23</span>.
<span>RSS</span> <span>=</span>  <span>1356</span><span>;</span> R<span>^</span><span>2</span> <span>=</span> <span>0.8766</span><span>;</span>
 needed <span>(</span><span>5</span>,<span>0</span><span>)</span> iterations</pre></td></tr></table></div>
<p>and actually, if we use them in a linear-spline regression, we get the same output here</p>
<div><table><tr><td><pre>1
2
3
4
</pre></td><td><pre>reg <span>=</span> <span>lm</span><span>(</span><span>dist</span>~bs<span>(</span>speed,degree<span>=</span><span>1</span>,<span>knots</span><span>=</span><span>c</span><span>(</span><span>4</span>,<span>7</span>,<span>8</span>,<span>9</span>,,<span>20</span>,<span>23</span>,<span>25</span><span>)</span><span>)</span>,<span>data</span><span>=</span><span>cars</span><span>)</span>
u <span>=</span> <span>seq</span><span>(</span><span>4</span>,<span>25</span>,<span>by</span><span>=</span>.1<span>)</span>
v <span>=</span> <span>predict</span><span>(</span>reg,newdata<span>=</span><span>data.<span>frame</span></span><span>(</span>speed<span>=</span>u<span>)</span><span>)</span>
<span>lines</span><span>(</span>u,v,<span>col</span><span>=</span><span>&quot;green&quot;</span><span>)</span></pre></td></tr></table></div>
<p>Let us add vertical lines for the knots</p>
<div><table><tr><td><pre>1
</pre></td><td><pre><span>abline</span><span>(</span>v<span>=</span><span>c</span><span>(</span><span>4</span>,<span>7</span>,<span>8</span>,<span>9</span>,<span>20</span>,<span>23</span>,<span>25</span><span>)</span>,<span>col</span><span>=</span><span>&quot;grey&quot;</span>,lty<span>=</span><span>2</span><span>)</span></pre></td></tr></table></div>
<p><img alt="" src="https://i2.wp.com/f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2018/07/convex2.png?w=456&ssl=1" data-recalc-dims="1"></p>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://freakonometrics.hypotheses.org/53384"> R-english – Freakonometrics</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/e1N6pB_HfHs" height="1"></a></div>