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.

  • $\begingroup$ 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$? $\endgroup$
    – Bertrand
    Sep 19, 2023 at 7:19
  • $\begingroup$ $\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! $\endgroup$ Sep 19, 2023 at 16:47
  • $\begingroup$ 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. $\endgroup$
    – Bertrand
    Sep 21, 2023 at 7:46


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.