# EM Algorithm and Mixtures of Generalized Linear Models

## Generalized Linear Models (GLM)

### Component of GLM&#x20;

GLM assumes the observations follows an exponential distribution&#x20;

$$
\begin{equation\*}
f(y\_t) = \exp \left(w\_t \frac{y\_t\theta\_t - b(\theta\_t)}{\theta} + c(y\_t, \phi / w\_t)\right)
\end{equation\*}
$$

* $$\theta\_t, \phi$$ : location and scale parameters&#x20;
  * often $$\phi$$ is a known constant, such as for Binomial and Poisson distribution it equals 1. If it is unkown, it is estimated by method of moments. Note that $$\theta\_t$$ and $$\phi$$ are orthogonal parameters, so their estimation does not affect each other
* $$w\_t$$: known weights
* $$b, c$$: known functions&#x20;

The mean and variance are&#x20;

$$
\begin{align\*}
& E\[y\_t] = \mu\_t = b'(\theta\_t) & Var(y\_t) = \frac{\phi}{w\_t}b''(\theta\_t)
\end{align\*}
$$

The distribution can also depend on predictors $$x\_j$$through a linear function such that the mean is assumed to be a smooth function of the linear combination of predictors&#x20;

$$
\begin{align\*}
& \mu\_t = g(\eta\_t) = g(\beta\_0 + \beta\_1x\_{1t} + \cdots + \beta\_k x\_{kt}) \\
& \eta\_t = g^{-1}(\mu\_t) = h(\mu\_t)
\end{align\*}
$$

* $$h$$ is also called the link function
* If $$h=(b')^{-1}$$, we call this the canonical link function and $$\theta\_t = \eta\_t$$

### Model Likelihood

The log-likelihood for $$T$$ independent observations from an exponential family distribution is&#x20;

$$
\begin{equation\*}
l(\theta\_{1:T}, \phi | y\_{1:T}) = \sum\_{t=1}^T w\_t \frac{y\_t\theta\_t - b(\theta\_t)}{\phi} + \sum\_{t=1}^Tc(y\_t, \phi/w\_t)
\end{equation\*}
$$

* to note again, likelihood function is not pdf. It is a function w\.r.t parameters.&#x20;

### Estimate GLM Through Iteratively Weighted Least Squares&#x20;

GLM model beta coeffcients can be estimated through Iteratively Weighted Least Squares, and it can be shown that this is equivalent to maximum likelihood estimates.&#x20;

Start with an initial guess of $$\hat{\eta\_t}$$ and $$\hat{\mu\_t} = h(\hat{\eta\_t})$$

1. Obtain working response $$z\_t = \hat{\eta}\_t + \frac{y\_t - \hat{\mu}\_t}{\frac{\partial\mu\_t}{\partial\eta\_t}}$$
2. Obtain working weight $$W\_t = \frac{w\_t}{b''(\hat{\theta}\_t)}(\frac{\partial\mu\_t}{\partial\eta\_t})^2$$
3. Estimate $$\hat{\beta} = (X^TWX)^{-1}X^TWz$$
   * $$X$$: model matrix&#x20;
   * $$W$$: diagonal matrix with the working weights&#x20;
   * $$z$$: vector of working response
4. Obtain new estiamte $$\hat{\eta}=X\hat{\beta}$$

Repeat the steps until convergence&#x20;

### Connection to EM Algorithm and Mixture Model&#x20;

Summary: We can use IWLS to maximize individual component density numerically if each component density of mixture model is in GLM family.&#x20;

Recall in mixture model, we can define joint log-likelihood as a function of the data, hidden states, and parameters. The expectation of complete data log-likelihood over all possible hidden states is&#x20;

$$
\begin{align\*}
Q(\theta, \theta') & =\sum\_{t=1}^T \sum\_{i=1}^N \gamma\_t(i) \log P(s\_t = i|\theta\_{pr}) + \sum\_{t=1}^T \sum\_{i=1}^N \gamma\_t(i)\log f(y\_t | s\_t = i, \theta\_{obs}) \\
\gamma\_t(i) & = p(s\_t=i|y\_t, \theta')
\end{align\*}
$$

Part of the EM process is then to obtain the parameters estiamte that maximizes this expected complete data log likelihood.&#x20;

$$
\begin{equation\*} \hat{\theta}*{obs} = \overset{\text{argmax}}{\theta*{obs}} \sum\_{t=1}^T \sum\_{i=1}^N \gamma\_t(i)\log f(y\_t|s\_t = i,\theta\_{obs}) \end{equation\*}
$$

If each component density has a separate set of parameters, then they can be estimated separately. So the parameter for component $$i$$ can be estimated via maximizing the weighted likelihood:&#x20;

$$
\begin{equation\*} \sum\_{t=1}^T \gamma\_t(i)\log f\_i(y\_t|s\_t = i,\theta\_{obs}^{(i)}) \end{equation\*}
$$

Notice the tricky part of the EM for Mixture model is that it is not always possible to analytically obtain the actual maximum. But if the density function for component $$i$$ is from the generalized linear model family, we can rewrite the formula as&#x20;

$$
\begin{align\*}
& \sum\_{t=1}^T w'*t \frac{y\_t\theta\_t - b(\theta\_t)}{\phi} + \sum*{t=1}^T\gamma\_t(i)c(y\_t, \phi/w\_t)\\
& w'\_t =\gamma\_t(i)w\_t
\end{align\*}
$$

Note that only the first part impacts the estimation of $$\beta$$, and subsequently $$\theta$$, we can thus rely on IWLS routine using weights $$w'\_t$$.&#x20;

Further note that we cannot use IWLS to estimate $$\phi$$ since the wieghts are not the same ($$w'\_t$$ vs. $$w\_t$$). But generally $$\phi$$ is estimated through a method of meoments anyway.&#x20;
