% When using TeXShop on the Mac, let it know the root document. The following must be one of the first 20 lines.
%!TEX root = ../design.tex

\chapter[Regression]{Regression}
\begin{moduleinfo}
\item[Authors] {Rahul Iyer and Hai Qian}
\item[History]
	\begin{modulehistory}
    \item[v0.3] Added section on Clustered Sandwich Estimators
    \item[v0.2] Added section on Marginal Effects
		\item[v0.1] Initial version, including background of regularization
	\end{modulehistory}
\end{moduleinfo}

% Abstract. What is the problem we want to solve?

Regression analysis is a statistical tool for the investigation of
relationships between variables. Usually, the investigator seeks to ascertain
the causal effect of one variable upon another-the effect of a price increase
upon demand, for example, or the effect of changes in the money supply upon the
inflation rate. More specifically, regression analysis helps one understand how
the typical value of the dependent variable changes when any one of the
independent variables is varied, while the other independent variables are held
fixed.

Regression models involve the following variables:
\begin{enumerate}
    \item The unknown parameters, denoted as $\beta$, which may represent a scalar or a vector.
    \item The independent variables, $x$
    \item The dependent variables, $y$
\end{enumerate}

\section{Linear Methods for Regression} % (fold)
\label{sub:linear_methods_for_regression}

% subsection linear_methods_for_regression (end)

\section{Regularization} % (fold)
\label{sub:regularization}

Usually, $y$ is the result of measurements contaminated by small errors
(noise). Frequently, ill-conditioned or singular systems also arise in the iterative solution of nonlinear systems or optimization problems. In all such situations, the vector $x = {A}^{-1}y$ (or in the full rank overdetermined
case $A^+ y$, with the pseudo inverse $A^+ = (A^T A)^{-1}A^T X)$, if it exists at all, is usually a meaningless bad approximation to x.

Regularization techniques are needed to obtain meaningful solution estimates 
for such ill-posed problems, and in particular when the number of parameters 
is larger than the number of available measurements, so that standard least 
squares techniques break down.

\subsection{Linear Ridge Regression}
Ridge regression is the most commonly used method of regularization of
ill-posed problems. Mathematically, it seeks to minimize

\begin{equation}
Q\left(\vec{w},w_0;\lambda\right)\equiv \min_{\vec{w},w_0}\left[ \frac{1}{2N} \sum_{i=1}^{N} \left( y_i - w_0 -
    \vec{w} \cdot \vec{x}_i \right)^2
  +\frac{\lambda}{2}\|\vec{w}\|_2^2 \right]\ ,
\end{equation}
for a given value of $\lambda$, where $\vec{w}$ and $w_0$ are the fitting coefficients, and $\lambda$
is a non-negative regularization parameter. $\vec{w}$ is a vector in
$d$ dimensional space, and
\begin{equation}
\|\vec{w}\|_2^2 = \sum_{j=1}^{d}w_j^2 = \vec{w}^T\vec{w}\ .
\end{equation}
When $\lambda = 0$, $Q$ is
the mean squared error of the fitting.

The intercept term $w_0$ is not regularized, because this term is
fully decided by the mean values of $y_i$ and $\vec{x}_i$ and the
values of $\vec{w}$, and does not affect the model's complexity.

$Q\left(\vec{w},w_0;\lambda\right)$ is a quadratic function of $\vec{w}$ and
  $w_0$, and thus can be solved analytically
\begin{equation}
\vec{w}_{ridge}=\left(\lambda\vec{I}_d +
  \vec{X}^T\vec{X}\right)^{-1}\vec{X}^T\vec{y}\ .
\end{equation}
By using the available Newton method (Sec. 6.2.4), the above quantity can be easily
calculated from one single step of the Newton method.

Many packages for Ridge regularization actually regularize the fitting
coefficients not for the fitting model for the original data but for
the data that has be scaled. MADlib also provides this option. When
the normalization parameter is set to be True, which is by default
False, the data will be first converted to the following before
applying the Ridge regularization.
\begin{equation}
  x_i' \leftarrow \frac{x_i - \langle x_i \rangle}{\langle (x_i -
    \langle x_i \rangle)^2\rangle} \ ,
\end{equation}
\begin{equation}
y_i \leftarrow y_i - \langle y_i \rangle \ ,
\end{equation}
where $\langle \cdot \rangle = \sum_{i=1}^{N} \cdot / N$.

Note that Ridge regressions for scaled data and un-scaled data are not equivalent.

\subsection{Elastic Net Regularization} % (fold)
\label{ssub:elastic_net_regularization}
As a continuous shrinkage method, ridge regression achieves its better prediction performance through a bias-variance trade-off. However, ridge regression cannot produce a parsimonious model, for it always keeps all the predictors in the model~\cite{zou2005}. Best subset selection in contrast produces a sparse model, but it is extremely variable because of its inherent discreteness. 

A promising technique called the lasso was proposed by Tibshirani (1996). The 
lasso is a penalized least squares method imposing an L1-penalty on the 
regression coefficients. Owing to the nature of the L1-penalty, the lasso does 
both continuous shrinkage and automatic variable selection simultaneously.

Although the lasso has shown success in many situations, it has some
limitations. Consider the following three scenarios:
\begin{enumerate}

    \item In the Number of features ($p$) >> Number of observations ($n$) case, the
    lasso selects at most $n$ variables before it saturates, because of the
    nature of the convex optimization problem. This seems to be a limiting
    feature for a variable selection method. Moreover, the lasso is not well
    defined unless the bound on the $L_1$-norm of the coefficients is smaller
    than a certain value.

    \item If there is a group of variables among which the pairwise correlations are
    very high, then the lasso tends to select only one variable from the group
    and does not care which one is selected.

    \item For usual n>p situations, if there are high correlations between predictors,
    it has been empirically observed that the prediction performance of the
    lasso is dominated by ridge regression.

\end{enumerate}

These scenarios make lasso an inappropriate variable selection method in some situations.

Hui Zou and Trevor Hastie [42] introduce a new regularization
technique called the `elastic net'. Similar to the lasso, the elastic net
simultaneously does automatic variable selection and continuous
shrinkage, and it can select groups of correlated variables. It is like a
stretchable fishing net that retains `all the big fish'.

The elastic net regularization minimizes the following target function
\begin{equation} \label{eq:target}
\min_{\vec{w} \in R^N}L(\vec{w}) + \lambda \left[\frac{1-\alpha}{2}\|\vec{w}\|_2^2 +
  \lambda\alpha \|\vec{w}\|_1\right]\ ,
\end{equation}
where $\|\vec{w}\|_1 = \sum_{i=1}^N|w_i|$ and $N$ is the number of features.

For the elastic net regularization on linear models,
\begin{equation}
L(\vec{w}) = \frac{1}{2M}\sum_{m=1}^M\left(y_m - w_0 - \vec{w} \cdot
  \vec{x}_m\right)^2\ ,
\end{equation}
where the sum is over all observations and $M$ is the total number of
observation.

For the elastic net regularization on logistic models,
\begin{equation}
L(\vec{w}) = \sum_{m=1}^M\left[y_m \log\left(1 + e^{-(w_0 +
      \vec{w}\cdot\vec{x}_m)}\right) + (1-y_m) \log\left(1 + e^{w_0 +
      \vec{w}\cdot\vec{x}_m}\right)\right]\ ,
\end{equation}
where $y_m \in \{0,1\}$.

\subsubsection{Optimizer Algorithms}
Right now, we support two algorithms for optimizer. The default one is
FISTA, and the other is IGD.

\paragraph{FISTA}

Fast Iterative Shrinkage Thresholding Algorithm (FISTA) with {\bf
  backtracking step size} [4]:
\vspace{0.2in}
\hrule
\vspace{0.2in}
{\bf Step $0$}: Choose $\delta>0$ and $\eta > 1$, and
$\vec{w}^{(0)} \in R^N$. Set $\vec{v}^{(1)}=\vec{w}^{(0)}$ and
$t_1=1$.

{\bf Step $k$}: ($k \ge 1$) Find the smallest nonnegative integers
$i_k$ such that with $\bar{\delta} = \delta_{k-1}/\eta^{i_k-1}$
\begin{equation}
F(p_{\bar{\delta}}(\vec{v}^{(k)})) \le
Q_{\bar{\delta}}(p_{\bar{\delta}}(\vec{v}^{(k)}), \vec{v}^{k})\ .
\end{equation}
Set $\delta_k = \delta_{k-1}/\eta^{i_k-1}$ and compute
\begin{equation}
\vec{w}^{(k)}  =  p_{\delta_k}\left(\vec{v}^{(k)}\right)\ ,
\end{equation}
\begin{equation}
t_{k+1} = \frac{1 + \sqrt{1 + 4t_k^2}}{2}\ ,
\end{equation}
\begin{equation}
\vec{v}^{(k+1)} = \vec{w}^{(k)} +
\frac{t_k-1}{t_{k+1}}\left(\vec{w}^{(k)} - \vec{w}^{(k-1)}\right)\ .
\end{equation}
\vspace{0.2in}
\hrule
\vspace{0.2in}
Here,
\begin{equation}
F(\vec{w}) = f(\vec{w}) + g(\vec{w})\ ,
\end{equation}
where $f(\vec{w})$ is the differentiable part of
Eq. (\ref{eq:target}) and $g(\vec{w})$ is the non-differentiable part. For linear models,
\begin{equation}
f(\vec{w}) = \frac{1}{2M}\sum_{m=1}^M\left(y_m - w_0 - \vec{w} \cdot
  \vec{x}_m\right)^2 + \frac{\lambda(1-\alpha)}{2}\|\vec{w}\|_2^2\ ,
\end{equation}
and for logistic models,
\begin{equation}
f(\vec{w}) = \sum_{m=1}^M\left[y_m \log\left(1 + e^{-(w_0 +
      \vec{w}\cdot\vec{x}_m)}\right) + (1-y_m) \log\left(1 + e^{w_0 +
      \vec{w}\cdot\vec{x}_m}\right)\right] + \frac{\lambda(1-\alpha)}{2}\|\vec{w}\|_2^2\ .
\end{equation}
And for both types of models,
\begin{equation}
g(\vec{w}) = \lambda\alpha\sum_{i=1}^N|w_i|\ .
\end{equation}

And
\begin{equation}
Q_{\delta}(\vec{a}, \vec{b}) := f(\vec{b}) + \langle \vec{a} -
\vec{b}, \nabla f(\vec{b})\rangle +
\frac{1}{2\delta}\|\vec{a} - \vec{b}\|^2 + g(\vec{a})\ ,
\end{equation}
where $\langle\cdot\rangle$ is just the usual vector dot product.

And the proxy function is defined as
\begin{equation}
p_\delta (\vec{v}) := \underset{\vec{w}}{\operatorname{arg\,min}} \left[g(\vec{w}) +
  \frac{1}{2\delta}\left\|\vec{w} - \left(\vec{v} - \delta\,\nabla f(\vec{v})\right)\right\|^2  \right]
\end{equation}
For our case, where $g(\vec{w}) = \lambda\alpha\sum_{i=1}^N|w_i|$, the
proxy function is simply equal to the soft-thresholding function

\begin{equation}
p_\delta (v_i) = \left\{ \begin{array}{ll}
v_i - \lambda\alpha\delta u_i\ , \quad  & \mbox{if } v_i > \lambda\alpha\delta
u_i \\
0\ , \quad & \mbox{otherwise} \\
v_i + \lambda\alpha\delta u_i\ , \quad & \mbox{if } v_i < - \lambda\alpha\delta u_i
\end{array}
\right\}
\end{equation}
where,
\begin{equation}
\vec{u} = \vec{v} - \delta\,\nabla f(\vec{v})\ .
\end{equation}

{\bf Active set method} is used in our implementation for FISTA to
speed up the computation. Considerable speedup is obtained by
organizing the iterations around the active set of features; those
with nonzero coefficients. After a complete cycle through all the
variables, we iterate on only the active set till convergence. If
another complete cycle does not change the active set, we are done,
otherwise the process is repeated.

{\bf Warm-up method} is also used to speed up the computation. When
the option parameter $warmup$ is set to be $True$, a series of lambda
values, which is strictly descent and ends at the lambda value that
the user wants to calculate, will be used. The larger lambda gives
very sparse solution, and the sparse solution again is used as the
initial guess for the next lambda's solution, which will speed up the
computation for the next lambda. For larger data sets, this can
sometimes accelerate the whole computation and might be faster than
computation on only one lambda value.

{\bf Note:} Our implementation is a little bit different from the
original FISTA. In the original FISTA, during the backtracking
procedure, the algorithm is searching for a non-negative integer $i_k$
and the new step size is $\delta_k = \delta_{k-1}/\eta^{i_k}$. Thus
the step size is non-increasing. Here, we allow the step size to
increase by using $\delta_k = \delta_{k-1}/\eta^{i_k-1}$ so that
larger step sizes can be tried by the algorithm. Tests show that this
is actually faster.

\paragraph{IGD}

The Incremental Gradient Descent (IGD) algorithm is a stochastic
algorithm by its nature. So it is difficult to get sparse
solutions. What we implemented is Stochastic Mirror Descent Algorithm
made Sparse (SMIDAS). The inverse p-form link function is used
\begin{equation} \label{eq:f}
h_j^{-1}(\vec{\theta}) = \frac{\mbox{sign}(\theta_j)\vert \theta_j
  \vert^{p-1}}{\|\theta\|_p^{p-2}}\ ,
\end{equation}
where
\begin{equation}
\|\theta\|_p = \left(\sum_j \vert \theta \vert ^p\right)^{1/p}\ ,
\end{equation}
and $\mbox{sign}(0) = 0$.
\vspace{0.2in}
\hrule
\vspace{0.2in}
Choose step size $\delta > 0$.

Let $p = 2\log d$ and let $h^{-1}$ be as in Eq. (\ref{eq:f})

Let $\vec{\theta} = \vec{0}$

{\bf for} $k = 1,2,\dots$

\quad $\vec{w} = h^{-1}(\vec{\theta})$

\quad $\vec{v} = \nabla f(\vec{w})$

\quad $\tilde{\vec{\theta}} = \theta - \delta \, \vec{v}$

\quad $\forall j, \theta_j = \mbox{sign}(\tilde{\theta}_j)\max (0,
\vert \tilde{\theta}_j \vert- \lambda\alpha\delta)$
\vspace{0.2in}
\hrule
\vspace{0.2in}

The resulting fitting coefficients of this algorithm is not really
sparse, but the values are very small (usually $< 10^{15}$), which can
be safely set to be zero after filtering with a threshold.

This is done as the following: (1) multiply each coefficient with the
standard deviation of the corresponding feature (2) compute the
average of absolute values of re-scaled coefficients (3) divide each
rescaled coefficients with the average, and if the resulting absolute
value is smaller than \emph{threshold}, set the original coefficient
to be zero.

IGD is in nature a sequential algorithm, and when running in a
distributed way, each segment of the data runs its own SGD model, and
the models are averaged to get a model for each iteration. This
average might slow down the convergence speed, although we acquire the
ability to process large data set on multiple machines. So this
algorithm provides an option  \emph{parallel} to let the user choose
whether to do parallel computation.

IGD also implements the {\bf warm-up method}.

{\bf Stopping Criteria} Both {\bf FISTA} and {\bf IGD} compute the average difference
between the coefficients of two consecutive iterations, and if the
difference is smaller than \emph{tolerance} or the iteration number
is larger than \emph{max\_iter}, the computation stops.

\section{Robust Variance via Huber-White Sandwich Estimators}
Given  $N$ data points, where the $i$th point is defined by a feature  vector $x_i \in \mathbb{R}^M$ and a category scalar $y_i$, where $y_i \in \mathbb{R}, y_i \in \{0,1 \}$, and $y_i \in \{1,\dots, L \}$ for linear, logistic, and multi-logistic regression respectively.  We assume that $y_i$ is drawn from an independent and identically distributed (i.i.d.) distribution determined by a $K$-dimensional parameter vector $\beta$ (which is $K\times L$ dimensional for multi-logistic regression).  

Generally, we are interested in finding the values of $\beta$ that best predict $y_i$ from $x_i$, with \textit{best} being defined as the values that maximize some log-likelihood function $l(y,x,\beta)$.  The maximization is typically solved using the derivative of the likelihood $\psi$  and the Hessian $H$.  More formally, $\psi$ is defined as 
\begin{align}
\psi(y,x, \beta) = \frac{\partial l(x,y,\beta)}{\partial \beta}
\end{align} 
and $H$ is defined as
\begin{align}
H(y,x, \beta) = \frac{\partial^2 l(x,y,\beta)}{\partial \beta^2}.
\end{align} 
Using these derivatives, one can solve the logistic or linear regression for the optimal solution, and compute the variance and other statistics of the regression.  

\subsection{Sandwich Operators}
However, we may believe that the underlying distribution is not i.i.d., (in particular, the variance is not independent of $x_i$), in which case, our variance estimates will be incorrect. 
 The Huber sandwich estimator is used to get a robust estimate of the variance even if the i.i.d. assumption is wrong.  The estimator is known as a sandwich estimator because the robust covariance matrix $S(\beta)$ of $\beta$ can be expressed in a \textit{sandwich formulation}, of the form
\begin{align}
S(\beta) = B(\beta) M(\beta) B(\beta).  
\end{align}
The $B(\beta)$ matrix is commonly called the \textit{bread}, whereas the $M(\beta)$ matrix is the \textit{meat}.  

\paragraph{The Bread}

Computing $B$ is relatively straightforward, 
\begin{align}
B(\beta) = N\left(\sum_i^N -H(y_i, x_i, \beta) \right)^{-1}
\end{align}

\paragraph{The Meat}

There are several choices for the $M$ matrix, each with different robustness properties.  In the Huber-White estimator, the matrix $M$ is defined as
\begin{align}
M_{H} =\frac{1}{N} \sum_i^N \psi(y_i,x_i, \beta)^T  \psi(y_i,x_i, \beta).
\end{align}


\subsection{Implementation}

The Huber-White sandwich estimators implemented for linear, logistic, and multinomial-logistic regression mimic the same framework as  the linear/logistic regression implementations.  In these implementations, the gradient and Hessian are computed in parallel, and a final step operates on the aggregate result.  

This framework breaks the computation into three steps: \textit{transition} step, the \textit{merge states} step, and the \textit{final} step.  The transition step computes the gradient and Hessian contribution from each row in the data matrix.  To compute this step, we need to define the derivatives for each regression.  
\subsubsection{Linear Regression Derivatives}
For linear regression, the derivative is 
\begin{align}
\frac{\partial l(x,y,\beta)}{\partial \beta} = X^T(y - X \beta)
\end{align}
and the Hessian is
\begin{align}
\frac{\partial^2 l(x,y,\beta)}{\partial^2 \beta} = -X^TX. 
\end{align}

\subsubsection{Logistic Regression Derivatives}
For logistic, the derivative is 
\begin{align}
\frac{\partial l(x,y,\beta)}{\partial \beta} = \sum_{i=1}^N \frac{-1^{y_i}\cdot \beta}{1 + \exp(-1^{y_i} \cdot  \beta^Tx)}  \exp(-1^{y_i}\cdot \beta^Tx)
\end{align}
and the Hessian is
\begin{align}
\frac{\partial^2 l(x,y,\beta)}{\partial^2 \beta} = -X^TAX = -\sum_{i=1}^N A_{ii} x_i^T x_i 
\end{align}
where $A$ is a diagonal matrix with 
\begin{align}
A_{ii} = \left( [1 + \exp(c^Tx_i) ) (1 + \exp(-c^Tx_i)] \right)^{-1}.
\end{align}

\subsubsection{Multi-Logistic Regression Derivatives}

For multi-logistic regression,  we replace $y_i$ with a vector $Y_i \in \{0,1\}^L$, where all entries of $Y_i$ are zero expect the $y_i$th entry, which is set to 1.  In addition, we choose a \textit{baseline} category, for which the odds of all the categories are measured against.  Let $J \in \{1, \dots, L \}$ be the baseline category.  

We define the variables
\begin{align}
\pi_{i,j} &= \frac{\exp\left(\sum_{k=1}^K x_{i,k} \beta_{k,j} \right)}{1 + \sum_{j\ne J} \exp \left(\sum_{k=1}^K x_{i,k} \beta_{k,j} \right)}, \ \ \forall j \ne J\\
\pi_{i,J} &= \frac{1}{1 + \sum_{j\ne J} \exp \left(\sum_{k=1}^K x_{i,k} \beta_{k,j} \right)}.
\end{align}

The derivatives are then
\begin{equation}\label{eq:first_derivative}
\frac{\partial l}{\partial \beta_{k,j}} = \sum_{i=1}^{N} Y_{i,j}x_{i,k} - \pi_{i,j}x_{i,k} \ \ \ \ \forall k \  \forall j
\end{equation}
The Hessian is then 
\begin{align}\label{eq:second_derivative}
\frac{\partial^2 l({\beta})}{\partial \beta_{k_2,j_2} \partial \beta_{k_1,j_1}} 
&= \sum_{i=1}^{N} -\pi_{i,j_2}x_{i,k_2}(1-\pi_{i,j_1})x_{i,k_1} &&j_1 = j_2 \\
&= \sum_{i=1}^{N} \pi_{i,j_2}x_{i,k_2}\pi_{i,j_1}x_{i,k_1} &&j_1 \neq j_2 
\end{align}

\subsubsection{Merge Step and Final Step}

For the logistic and multi-logistic, the derivative and Hessian are sums in which the terms can be computed in parallel, and thus in a distributed manner.  Each transition step computes a single term in the sum.  The merge step sums two or more terms computed by the transition steps, and the final step computes the Hessian inverse, and the matrix product between the bread and meat matrices.  

We note that the Hessian and derivative both depend on $\beta$, so the variance estimates will also depend on $\beta$.  Rather than allow the user to specify a $\beta$ value, the implementation computes the optimal $\beta$ by running the appropriate regression  before computing the robust variance.  In cases where the regression has parameters (regression tolerance, max iterations), the interface allows the user to specify those parameters.

% subsubsection elastic_net_regularization (end)
% subsection regularization (end)

% TODO: Remove the pagebreak before merging to master
\pagebreak

\section{Marginal Effects} % (fold)
\label{sub:marginal_effects}
\textit{Most of the notes below are based on~\cite{diekmann2008}}

A \emph{marginal effect} (ME) or partial effect measures the effect on the
conditional mean of $y$ of a change in one of the regressors, say
$X_k$~\cite{cameron2009}. In the linear regression model, the ME equals the
relevant slope coefficient, greatly simplifying analysis. For nonlinear models,
this is no longer the case, leading to remarkably many different methods for
calculating MEs.

Let $E(y_i | x_i)$ represent the expected value of a dependent variable $y_i$
given a vector of independent variables $x_i$. In the case where $y$ is a
linear function of $(x_1, \dots, x_m) = X$ and $y$ is a continuous variable, a
linear regression model can be stated as follows:
\begin{align*}
    & y = X' \beta \\
    & \text{or} \\
    & y = \beta_0 + \beta_1 x_1 +  \dots  + \beta_l x_m.
\end{align*}
From the above equation it is straightforward to see that the marginal effect of
variable $x_k$ on the dependent variable is $\partial y / \partial x = \beta_k$.

The standard approach to modeling dichotomous/binary variables
(so $y \in {0, 1}$) is to estimate a generalized linear model under the
assumption that y follows some form of Bernoulli distribution. Thus the expected
value of $y$ becomes,
\begin{equation*}
    y = G(X' \beta),
  \end{equation*}
where G is the specified binomial distribution. Here we assume to use
logistic regression and use $g$ to refer to the inverse logit function.

\subsection{Logistic regression} % (fold)
\label{sub:logistic_regression}
In logistic regression:
\begin{align*}
  P &= \frac{1}{1 + e^{-(\beta_0 + \beta_1 x_1 + \dots  \beta_m x_m)}} \\
    &= \frac{1}{1 + e^{-z}}
\end{align*}
\begin{align*}
  \implies \frac{\partial P}{\partial X_k} &= \beta_k \cdot \frac{1}{1 + e^{-z}} \cdot
              \frac{e^{-z}}{1 + e^{-z}} \\
      &= \beta_k \cdot P \cdot (1-P)
\end{align*}


There are two main methods of calculating the marginal effects for dichotomous
dependent variables.
\begin{enumerate}
  \item The first uses the average of the marginal effects at every sample
  observation. This is calculated as follows:
  \begin{gather*}
    \frac{\partial y}{\partial x_k} = \beta_k \frac{\sum_{i=1}^{n} P(y_i = 1)(1-P(y_i = 1))}{n}, \\
    \text{where, } P(y_i=1) = g(X^{(i)} \beta) \\
    \text{and, } g(z) = \frac{1}{1 + e^{-z}} \\
  \end{gather*}

  \item The second approach calculates the marginal effect for $x_k$ by taking
    predicted probability calculated when all regressors are held at their mean
    value from the same formulation with the exception of adding one unit to $x_k$.
    The derivation of this marginal effect is captured by the following:
    \begin{gather*}
      \frac{\partial y}{\partial x_k} = \quad \beta_k P(y=1|\bar{X})(1-P(y=1|\bar{X})) \\
      \text{where, } \bar{X} = \frac{\sum_{i=1}^{n}X^{(i)}}{n}
    \end{gather*}
\end{enumerate}
% subsection logistic_regression (end)

\subsection{Discrete change effect} % (fold)
\label{sub:discrete_change_effect}
Along with marginal effects we can also compute the following discrete change
effects.

\begin{enumerate}
  \item Unit change effect, which should not be confused with the
marginal effect:
\begin{align*}
   \frac{\partial y}{\partial x_k} & = P(y=1|X,x_k+1) - P(y=1|X, x_k) \\
                                     & = g(\beta_0 + \beta_1 x_1 + \dots  + \beta_k (x_k+1) + \beta_l x_l) \\
                                    & \qquad  - g(\beta_0 + \beta_1 x_1 + \dots  + \beta_k x_k + \beta_l x_l)
 \end{align*}
    \item Centered unit change:
        \begin{gather*}
          \frac{\partial y}{\partial x_k} = P(y=1|X,x_k+0.5) - P(y=1|X, x_k-0.5) \\
        \end{gather*}
    \item Standard deviation change:
        \begin{gather*}
          \frac{\partial y}{\partial x_k} = P(y=1|X,x_k+0.5\delta_k) - P(y=1|X, x_k-0.5\delta_k) \\
        \end{gather*}
    \item Min-max change:
        \begin{gather*}
          \frac{\partial y}{\partial x_k} = P(y=1|X,x_k=x_k^{max}) - P(y=1|X, x_k=x_k^{min}) \\
        \end{gather*}
\end{enumerate}
% subsection discrete_change_effect (end)

\subsection{Multilogistic regression} % (fold)
\label{sub:multilogistic_regression}
The probabilities of different outcomes for multilogistic regression are expressed as,
\begin{gather*}
  P(y=j | X)  = \frac{e^{X\beta_j}}{\sum_{l=1}^{j} e^{X\beta_l}},
\end{gather*}
with $\beta$ set to zero for one of the outcomes. The output which the
$\beta$ vector is set to zero is called the ``base outcome'' or the ``reference category''.

The odds of outcome $j$ verus outcome $m$ are
\begin{align*}
  \frac{P(y=j | X)}{P(y=m | X)} & = \frac{\frac{e^{X\beta_j}}{\sum_{l=1}^{j} e^{X\beta_l}}}
                                       {\frac{e^{X\beta_m}}{\sum_{l=1}^{j} e^{X\beta_l}}} \\
                                & = \frac{e^{X\beta_j}}{e^{X\beta_m}}
\end{align*}
Thus the partial derivative would be,
\begin{gather*}
  \frac{\partial ln(P_j/P_m)}{\partial x_k} = \beta_{kj} - \beta_{km}
\end{gather*}

Thus, if $x_k$ is increased by one unit the odds of outcome $j$ against the
base outcome changes by exp($\beta_{kj}$) (since $\beta_{km}=0$ for the base
outcome). Expanding $P_m$ we finally get,
\begin{gather*}
  \frac{\partial P(y=j|X)}{\partial x_k} = P(y=j|X)\left[ \beta_{kj} - \sum_{l=1}^{j}\beta_{kl} P(y=l|X) \right]
\end{gather*}

The discrete change effect can be computed as,
\begin{gather*}
  \frac{\Delta P(y=j|X)}{\Delta x_k} = P(y=j|X, x_k+1) - P(y=j|X)
\end{gather*}
% subsection multilogistic_regression (end)

\subsection{Standard Errors} % (fold)
\label{sub:standard_errors}
The delta method is a popular way to estimate standard errors of non-linear
functions of model parameters. While it is straightforward to calculate the
variance of a linear function of a random variable, it is not for a nonlinear
function. The delta method therefore relies on finding a linear approximation
of the function by using a first-order Taylor expansion.

We can approximate a function $f(x)$ about a value $a$ as,
\[
  f(x) \approx f(a) + (x-a)f'(a)
\]

Taking the variance and setting $a = \mu_x$,
\[
  Var(f(X)) \approx \left[f'(\mu_x)\right]^2 Var(X)
\]

Using this technique, to compute the variance of the marginal effects at the
mean observation value in \emph{logistic regression}, we obtain:
\begin{gather*}
  Var(ME_k) = \frac{\partial (\beta_k \bar{P} (1- \bar{P}))}{\partial \beta_k} Var(\beta_k),\\
  \text{where, } \bar{P} = g(\bar{X}' \beta) = \frac{1}{1 + e^{-\bar{z}}} \\
  \text{and }    \bar{z} = \beta_0 + \beta_1 \bar{x}_1 + \dots + \beta_m \bar{x}_m
\end{gather*}

Thus, using the rule for differentiating compositions of functions, we get
\begin{align*}
  Var(ME_k) & = \left(-\beta_k \bar{P} \frac{\partial \bar{P}}{\partial \beta_k} +
              \beta_k (1-\bar{P})\frac{\partial \bar{P}}{\partial \beta_k} +
              \bar{P}(1-\bar{P}) \right) Var(\beta_k) \\
          & = \left( (1-2\bar{P})\beta_k \frac{\partial \bar{P}}{\partial \beta_k} + \bar{P}(1-\bar{P}) \right) Var(\beta_k)
\end{align*}
We have,
\begin{align*}
  \frac{\partial \bar{P}}{\partial \beta_k} & = \frac{\partial (\frac{1}{1 + e^{-z}})}{\partial \beta_k} \\
                                 & = \frac{1}{(1+e^{-z})^2} e^{-z} \frac{\partial z}{\partial \beta_k} \\
                                 & = \frac{x_k e^{-z}}{(1+e^{-z})^2} \\
                                 & = x_k \bar{P} (1 - \bar{P})
\end{align*}
Replacing this in the equation for $Var(ME_k)$,

\begin{align*}
  Var(ME_k) =  \bar{P}(1-\bar{P}) \left(1 + (1-2\bar{P})\beta_k x_k \right) Var(\beta_k)
\end{align*}

Since $\beta$, is a multivariate variable, we will have to use the variance-
covariance matrix of $\beta$  to compute the variance of the marginal effects.
Thus for the vector of marginal effects the equation becomes,

\begin{gather*}
  Var(ME) = \bar{P}^2(1-\bar{P})^2 \left[I + (1-2\bar{P})\beta \bar{X}' \right] V \left[I+ (1-2\bar{P}) \bar{X} \beta' \right],
\end{gather*}
where $V$ is the estimated variance-covariance matrix of $\beta$.
% subsection standard_errors (end)
% section marginal_effects (end)

\section{Clustered Standard Errors} % (fold)
\label{sec:clustered_standard_errors}
Adjusting standard errors for clustering can be important. For
example, replicating a dataset 100 times should not increase the
precision of parameter estimates. However, performing this procedure
with the IID assumption will actually do this. Another example is in
economics of education research, it is reasonable to expect that the
error terms for children in the same class are not
independent. Clustering standard errors can correct for this.

\subsection{Overview of Clustered Standard Errors}

Assume that the data can be separated into $m$ clusters. Usually this
can be down by grouping the data table according to one or multiple
columns.

The estimator has a similar form to the usual sandwich estimator
\begin{equation}
  S(\vec{\beta}) = B(\vec{\beta}) M(\vec{\beta}) B(\vec{\beta})
\end{equation}

The bread part is the same as sandwich estimator
\begin{eqnarray}
  B(\vec{\beta}) & = & \left(-\sum_{i=1}^{n} H(y_i, \vec{x}_i,
    \vec{\beta})\right)^{-1}\\
  & = & \left(-\sum_{i=1}^{n}\frac{\partial^2 l(y_i, \vec{x}_i,
      \vec{\beta})}{\partial \beta_\alpha \partial \beta_\beta}\right)^{-1}
\end{eqnarray}
where $H$ is the hessian matrix, which is the second derivative of the
target function
\begin{equation}
  L(\vec{\beta}) = \sum_{i=1}^n l(y_i, \vec{x}_i, \vec{\beta})\ .
\end{equation}

The meat part is different
\begin{equation}
  M(\vec{\beta}) = \bf{A}^T\bf{A}
\end{equation}
where the $m$-th row of $\bf{A}$ is
\begin{equation}
  A_m = \sum_{i\in G_m}\frac{\partial
      l(y_i,\vec{x}_i,\vec{\beta})}{\partial \vec{\beta}}
\end{equation}
where $G_m$ is the set of rows that belong to the same cluster.

\subsection{Implementation}

We can compute the quantities of $B$ and $A$ for each cluster during one scan through
the data table in an aggregate function. Then sum over all clusters to
the full $B$ and $A$ in the outside of the aggregate function. At last, the matrix mulplitications
are done in a separate function on master.
% section clustered_standard_errors (end)
