# Using an EM algorithm for a random parameters model with mixed data

Suppose $$x_i \sim f(\theta,\nu_i)$$ and $$y_i \sim g(\tau,\theta,\nu_i)$$ and $$\nu_i \sim N(0,\sigma)$$. Importantly, $$f$$ is smooth in $$\theta$$, but $$g$$ is not.

We have data on $$x_i$$, $$y_i$$ and the goal is to estimate $$A=(\tau,\theta,\sigma)$$.

The marginal likelihoods are: $$L_{x_i}(\theta,\sigma) = \int f(\theta,\nu_i) \ d\Phi(\nu_{i};0,\sigma)$$ $$L_{y_i}(\tau,\theta,\sigma) = \int g(\tau,\theta,\nu_i) \ d\Phi(\nu_{i};0,\sigma)$$

However, maximizing the individual likelihoods would lead to biased estimates of $$\sigma$$ due to the correlation via $$\nu_i$$ (Is this correct?). Alternatively, the joint likelihood is then: $$L(\tau,\theta,\sigma) = \prod_i \int f(\theta,\nu_i) g(\tau,\theta,\nu_i) \ d\Phi(\nu_{i};0,\sigma)$$

To maximize this (simulated) likelihood, given $$R$$ draws of $$\nu_i$$, the EM algorithm develops a surrogate function: $$Q(A | A^m) = \sum_{i,r} w_{i,r}^m \ln f(\theta,\nu_{ir}) g(\tau,\theta,\nu_{ir}) + \sum_{i,r} w_{i,r}^m \ln \phi(\nu_{ir};0,\sigma)$$ where $$w_{ir}^m = \frac{f(\theta^m,\nu_{ir}) g(\tau^m,\theta^m,\nu_{ir})}{\sum_{r'}f(\theta^m,\nu_{ir'}) g(\tau^m,\theta^m,\nu_{ir'}) }$$ Iteratively maximizing this function with respect to parameters, converges to the maximum of the simulated joint likelihood.

Call the first sum $$Q_1(\theta,\tau; A^m)$$ and the second $$Q_2(\sigma; A^m)$$, Maximizing $$Q_2$$ is trivial, and reduces to the (weighted) sample standard deviation of $$\nu_{ir}$$.

Maximizing $$Q_1$$ over $$\tau,\theta$$ is challenging because $$g$$ is not smooth in $$\theta$$. Instead, we decompose this into $$\tilde{Q}_1$$: $$\tilde{Q}_1(\theta,\tau; A^m) = \sum_{i,r} w_{i,r}^m \ln f(\theta,\nu_{ir}) + \sum_{i,r} w_{i,r}^m \ln g(\tau,\theta^m,\nu_{ir})$$ where the $$\theta$$ in the second sum is replaced by $$\theta^m$$. Now, the first sum can be maximized because $$f$$ is smooth in $$\theta$$ and the second sum can be maximized because $$g$$ is smooth in $$\tau$$ (and no longer depends on $$\theta$$).

The intuition is that, because $$\tau,\theta$$ are fixed parameters, the data $$x_i$$ provides information to separately identify $$\theta$$ and $$y_i$$ provides information to separately identify $$\tau$$ conditional on $$\nu_i$$. Meanwhile, the estimation of $$\sigma$$, which comes from maximizing $$Q_2$$, still accounts for the joint likelihood via $$w_{ir}^m$$.

I understand that this would be a loss in efficiency, as we are essentially dropping data by not using $$y_i$$ to estimate $$\theta$$. However, it seems reasonable that the estimate would still be unbiased and consistent. I'm confused about how to argue this or if I'm missing something completely.

• Clarification question: as $\nu$ has an $i$ subscript, it looks like that you have an incidental parameter problem here. Why do you also use notation $\nu_{ir}$? What is $\Phi(0,\sigma)$ ? It depends on two parameters only? why is it evaluated at $0$? Sep 19 at 7:19
• $\nu_i$ is an unobserved random effect, with known mean zero and unknown variance $\sigma$. Its CDF is $\Phi(\nu_i; 0,\sigma)$. The $r$ subscript denotes the $r$th draw of $\nu_i$ in the EM algorithm (which takes $R$ random draws of $\nu_i$ in each iteration). So it should only appear when I define $Q$, thanks for the edit! Sep 19 at 16:47
• What is $A^m$ ? Where is $\theta^m$ in the expression of $\tilde{Q}_1$? If the dimension of $\theta$ is small, I would tackle the optimization problem using Nelder-Mead or other numerical methods which do not rely on derivatives. Sep 21 at 7:46