<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/cvxr-a-direct-standardization-example/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://rviews.rstudio.com/2018/07/20/cvxr-a-direct-standardization-example/"> R Views</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>In our <a rel="nofollow" href="https://rviews.rstudio.com/2017/11/27/introduction-to-cvxr/" target="_blank">first blog post</a>, we introduced <code>CVXR</code>, an R package for disciplined convex optimization, and showed how to model and solve a non-negative least squares problem using its interface. This time, we will tackle a non-parametric estimation example, which features new atoms as well as more complex constraints.</p>
<div>
<h2>Direct Standardization</h2>
<p>Consider a set of observations <span>\((x_i,y_i)\)</span> drawn non-uniformly from an unknown distribution. We know the expected value of the columns of <span>\(X\)</span>, denoted by <span>\(b \in {\mathbf R}^n\)</span>, and want to estimate the true distribution of <span>\(y\)</span>. This situation may arise, for instance, if we wish to analyze the health of a population based on a sample skewed toward young males, knowing the average population-level sex, age, etc.</p>
<p>A naive approach would be to simply take the empirical distribution that places equal probability <span>\(1/m\)</span> on each <span>\(y_i\)</span>. However, this is not a good estimation strategy when our sample is unbalanced. Instead, we will use the method of <strong>direct standardization</strong> (Fleiss, Levin, and Paik 2003, 19.5): we solve for weights <span>\(w \in {\mathbf R}^m\)</span> of a weighted empirical distribution, <span>\(y = y_i\)</span> with probability <span>\(w_i\)</span>, which rectifies the skewness of the sample. This can be posed as the convex optimization problem</p>
<p><span>\[<br>
\begin{array}{ll} \underset{w}{\mbox{maximize}} &amp; \sum_{i=1}^m -w_i\log w_i \\<br>
\mbox{subject to} &amp; w \geq 0, \quad \sum_{i=1}^m w_i = 1,\quad X^Tw = b.<br>
\end{array}<br>
\]</span></p>
<p>Our objective is the total entropy, which is concave on <span>\({\mathbf R}_+^m\)</span>. The constraints ensure <span>\(w\)</span> is a probability vector that induces our known expectations over the columns of <span>\(X\)</span>, i.e., <span>\(\sum_{i=1}^m w_iX_{ij} = b_j\)</span> for <span>\(j = 1,\ldots,n\)</span>.</p>
</div>
<div>
<h2>An Example with Simulated Data</h2>
<p>As an example, we generate <span>\(m = 1000\)</span> data points <span>\(x_{i,1} \sim \mbox{Bernoulli}(0.5)\)</span>, <span>\(x_{i,2} \sim \mbox{Uniform}(10,60)\)</span>, and <span>\(y_i \sim N(5x_{i,1} + 0.1x_{i,2},1)\)</span>. We calculate <span>\(b_j\)</span> to be the mean over <span>\(x_{.,j}\)</span> for <span>\(j = 1,2\)</span>. Then we construct a skewed sample of <span>\(m = 100\)</span> points that over-represent small values of <span>\(y_i\)</span>, thus biasing its distribution downwards.</p>
<p>Using <code>CVXR</code>, we construct the direct standardization problem. We first define the variable <span>\(w\)</span>.</p>
<pre><code>w &lt;- Variable(m)</code></pre>
<p>Then, we form the objective function by combining <code>CVXR</code>’s library of operators and atoms.</p>
<pre><code>objective &lt;- Maximize(sum(entr(w)))</code></pre>
<p>Here, <code>entr</code> is the element-wise entropy atom; the S4 object <code>entr(w)</code> represents an <span>\(m\)</span>-dimensional vector with entries <span>\(-w_i\log(w_i)\)</span> for <span>\(i=1,\ldots,m\)</span>. The <code>sum</code> operator acts exactly as expected, forming an expression that is the sum of the entries in this vector. (For a full list of atoms, see the <a rel="nofollow" href="http://cvxr.rbind.io/post/cvxr_functions/" target="_blank">function reference</a> page).</p>
<p>Our next step is to generate the list of constraints. Note that, by default, the relational operators apply over all entries in a vector or matrix.</p>
<pre><code>constraints &lt;- list(w >= 0, sum(w) == 1, t(X) %*% w == b)</code></pre>
<p>Finally, we are ready to formulate and solve the problem.</p>
<pre><code>prob &lt;- Problem(objective, constraints)
result &lt;- solve(prob)
weights &lt;- result$getValue(w)</code></pre>
<p>Using our optimal <code>weights</code>, we can then re-weight our skewed sample and compare it to the population distribution. Below, we plot the density functions using linear approximations for the range of <span>\(y\)</span>.</p>
<pre><code>## Approximate density functions
dens1 &lt;- density(ypop)
dens2 &lt;- density(y)
dens3 &lt;- density(y, weights = weights)
yrange &lt;- seq(-3, 15, 0.01)
d &lt;- data.frame(x = yrange,
                True = approx(x = dens1$x, y = dens1$y, xout = yrange)$y,
                Sample = approx(x = dens2$x, y = dens2$y, xout = yrange)$y,
                Weighted = approx(x = dens3$x, y = dens3$y, xout = yrange)$y)

## Plot probability distribution functions
plot.data &lt;- gather(data = d, key = &quot;Type&quot;, value = &quot;Estimate&quot;, True, Sample, Weighted,
                    factor_key = TRUE)
ggplot(plot.data) +
    geom_line(mapping = aes(x = x, y = Estimate, color = Type)) +
    theme(legend.position = &quot;top&quot;)</code></pre>
<pre><code>## Warning: Removed 300 rows containing missing values (geom_path).</code></pre>
<div><span></span><br>
<img alt="Probability distribution functions: population, skewed sample and reweighted sample" src="https://i1.wp.com/rviews.rstudio.com/post/2018-07-20-cvxr-a-direct-standardization-example_files/figure-html/unnamed-chunk-6-1.png?w=450&ssl=1" data-recalc-dims="1">
<p>
Figure 1: Probability distribution functions: population, skewed sample and reweighted sample
</p>
</div>
<pre><code>## Return the cumulative distribution function
get_cdf &lt;- function(data, probs, color = 'k') {
    if(missing(probs))
        probs &lt;- rep(1.0/length(data), length(data))
    distro &lt;- cbind(data, probs)
    dsort &lt;- distro[order(distro[,1]),]
    ecdf &lt;- base::cumsum(dsort[,2])
    cbind(dsort[,1], ecdf)
}

## Plot cumulative distribution functions
d1 &lt;- data.frame(&quot;True&quot;, get_cdf(ypop))
d2 &lt;- data.frame(&quot;Sample&quot;, get_cdf(y))
d3 &lt;- data.frame(&quot;Weighted&quot;, get_cdf(y, weights))

names(d1) &lt;- names(d2) &lt;- names(d3) &lt;- c(&quot;Type&quot;, &quot;x&quot;, &quot;Estimate&quot;)
plot.data &lt;- rbind(d1, d2, d3)

ggplot(plot.data) +
    geom_line(mapping = aes(x = x, y = Estimate, color = Type)) +
    theme(legend.position = &quot;top&quot;)</code></pre>
<div><span></span><br>
<img alt="Cumulative distribution functions: population, skewed sample and reweighted sample" src="https://i1.wp.com/rviews.rstudio.com/post/2018-07-20-cvxr-a-direct-standardization-example_files/figure-html/unnamed-chunk-7-1.png?w=450&ssl=1" data-recalc-dims="1">
<p>
Figure 2: Cumulative distribution functions: population, skewed sample and reweighted sample
</p>
</div>
<p>As is clear from the plots, the sample probability distribution peaks around <span>\(y = 2.0\)</span>, and its cumulative distribution is shifted left from the population’s curve, a result of the downward bias in our sampled <span>\(y_i\)</span>. However, with the direct standardization weights, the new empirical distribution cleaves much closer to the true distribution shown in red.</p>
<p>We hope you’ve enjoyed this demonstration of <code>CVXR</code>. For more examples, check out our <a rel="nofollow" href="http://cvxr.rbind.io/" target="_blank">official site</a> and recent presentation <a rel="nofollow" href="https://www.youtube.com/watch?v=MyglbtnmQ8A" target="_blank">“Disciplined Convex Optimization with CVXR”</a> at useR! 2018.</p>
</div>
<p> </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://rviews.rstudio.com/2018/07/20/cvxr-a-direct-standardization-example/"> R Views</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/XOMVu_WLD1M" height="1"></a></div>