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.
