# Direct Standardization

## Introduction

Consider a set of observations $$(x_i,y_i)$$ drawn non-uniformly from an unknown distribution. We know the expected value of the columns of $$X$$, denoted by $$b \in {\mathbf R}^n$$, and want to estimate the true distribution of $$y$$. 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. The empirical distribution that places equal probability $$1/m$$ on each $$y_i$$ is not a good estimate.

So, we must determine the weights $$w \in {\mathbf R}^m$$ of a weighted empirical distribution, $$y = y_i$$ with probability $$w_i$$, which rectifies the skewness of the sample (Fleiss, Levin, and Paik 2003, 19.5). We can pose this problem as

$\begin{array}{ll} \underset{w}{\mbox{maximize}} & \sum_{i=1}^m -w_i\log w_i \\ \mbox{subject to} & w \geq 0, \quad \sum_{i=1}^m w_i = 1,\quad X^Tw = b. \end{array}$

Our objective is the total entropy, which is concave on $${\mathbf R}_+^m$$, and our constraints ensure $$w$$ is a probability distribution that implies our known expectations on $$X$$.

To illustrate this method, we generate $$m = 1000$$ data points $$x_{i,1} \sim \mbox{Bernoulli}(0.5)$$, $$x_{i,2} \sim \mbox{Uniform}(10,60)$$, and $$y_i \sim N(5x_{i,1} + 0.1x_{i,2},1)$$. Then we construct a skewed sample of $$m = 100$$ points that overrepresent small values of $$y_i$$, thus biasing its distribution downwards. This can be seen in Figure , where the sample probability distribution peaks around $$y = 2.0$$, and its cumulative distribution is shifted left from the populationâ€™s curve. Using direct standardization, we estimate $$w_i$$ and reweight our sample; the new empirical distribution cleaves much closer to the true distribution shown in red.

In the CVXR code below, we import data from the package and solve for $$w$$.

suppressWarnings(library(CVXR, warn.conflicts=FALSE))

## Import problem data
data(dspop)   # Population
data(dssamp)  # Skewed sample

ypop <- dspop[,1]
Xpop <- dspop[,-1]
y <- dssamp[,1]
X <- dssamp[,-1]
m <- nrow(X)

## Given population mean of features
b <- as.matrix(apply(Xpop, 2, mean))

## Construct the direct standardization problem
w <- Variable(m)
objective <- sum(entr(w))
constraints <- list(w >= 0, sum(w) == 1, t(X) %*% w == b)
prob <- Problem(Maximize(objective), constraints)

## Solve for the distribution weights
result <- solve(prob)
weights <- result$getValue(w) result$value
## [1] 4.223305

We can plot the density functions using linear approximations for the range of $$y$$.

library(ggplot2)
## Plot probability density functions
dens1 <- density(ypop)
dens2 <- density(y)
dens3 <- density(y, weights = weights)
yrange <- seq(-3, 15, 0.01)
d <- 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)
library(tidyr)
plot.data <- gather(data = d, key = "Type", value = "Estimate", True, Sample, Weighted,
factor_key = TRUE)
ggplot(plot.data) +
geom_line(mapping = aes(x = x, y = Estimate, color = Type)) +
theme(legend.position = "top")
## Warning: Removed 300 rows containing missing values (geom_path).

Followed by the cumulative distribution function.

## Return the cumulative distribution function
get_cdf <- function(data, probs, color = 'k') {
if(missing(probs))
probs <- rep(1.0/length(data), length(data))
distro <- cbind(data, probs)
dsort <- distro[order(distro[,1]),]
ecdf <- base::cumsum(dsort[,2])
cbind(dsort[,1], ecdf)
}

## Plot cumulative distribution functions
d1 <- data.frame("True", get_cdf(ypop))
d2 <- data.frame("Sample", get_cdf(y))
d3 <- data.frame("Weighted", get_cdf(y, weights))

names(d1) <- names(d2) <- names(d3) <- c("Type", "x", "Estimate")
plot.data <- rbind(d1, d2, d3)

ggplot(plot.data) +
geom_line(mapping = aes(x = x, y = Estimate, color = Type)) +
theme(legend.position = "top")

## Session Info

sessionInfo()
## R version 3.5.0 (2018-04-23)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.4
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
##
## locale:
## [1] C
##
## attached base packages:
## [1] stats     graphics  grDevices datasets  utils     methods   base
##
## other attached packages:
## [1] tidyr_0.8.1   ggplot2_2.2.1 CVXR_0.99
##
## loaded via a namespace (and not attached):
##  [1] gmp_0.5-13.1      Rcpp_0.12.17      highr_0.6
##  [4] pillar_1.2.2      compiler_3.5.0    plyr_1.8.4
##  [7] R.methodsS3_1.7.1 R.utils_2.6.0     tools_3.5.0
## [10] digest_0.6.15     bit_1.1-13        evaluate_0.10.1
## [13] tibble_1.4.2      gtable_0.2.0      lattice_0.20-35
## [16] rlang_0.2.0       Matrix_1.2-14     yaml_2.1.19
## [19] blogdown_0.6.3    xfun_0.1          Rmpfr_0.7-0
## [22] ECOSolveR_0.4     stringr_1.3.1     knitr_1.20
## [25] tidyselect_0.2.4  rprojroot_1.3-2   bit64_0.9-7
## [28] grid_3.5.0        glue_1.2.0        R6_2.2.2
## [31] rmarkdown_1.9.14  bookdown_0.7      purrr_0.2.4
## [34] magrittr_1.5      backports_1.1.2   scales_0.5.0
## [37] htmltools_0.3.6   scs_1.1-1         colorspace_1.3-2
## [40] labeling_0.3      stringi_1.2.2     lazyeval_0.2.1
## [43] munsell_0.4.3     R.oo_1.22.0

R Markdown

## References

Fleiss, J. L., B. Levin, and M. C. Paik. 2003. Statistical Methods for Rates and Proportions. Wiley-Interscience.