From joint to marginal likelihood

Question

I am trying to calculate the marginal maximum likelihood estimate of a parameter ρ given a vector of data x for person i across n items. This is in the context of Item Response Theory. The first (usual IRT) parameter, θ, is the ability (theoretically unbounded but usually between -3 and 3). The second parameter, ρ, is the probability that the respondent answers in a non-optimal way. I want to integrate θ out of my joint likelihood equation:

likelihood http://lehrfeld.me/images/lik.jpg

where the probability function in the integrand is

probability http://lehrfeld.me/images/prob.jpg

The second term in the integrand, f(θ), will just be a numerical approximation to the probability density function of θ. The three dots in the first set of curly brackets above can take three different forms depending on certain conditions: 0 (in the "normal" condition, so person i will definitely not respond in a non-optimal way), the formula in the second set of curly brackets, or one minus the formula in the second set of curly brackets.

My problem is just in wrapping my head around how I can accomplish these calculations. I am using R, but I suppose any good explanation of how I can do this computationally would be sufficient for me to understand enough to try to implement it. It is my understanding that in order to be able to carry out these calculations, I need to express the first term in the integrand (i.e., the joint probability function) as a function of 0. Then, I multiply this function of 0 by many numerical approximations to the second term (pdf of 0). Finally, I integrate each of these products over 0. Does this mean that the integrand is entirely a function of 0, and once I integrate 0 out I essentially go from a formula to a number? At which point I would sum of the results over all n items, and then implement an optimization routine to search for the most likely value of ρ within the bounds allowed?

Answer 1

Your intuitions in your final sentences sound about right - the integrand is entirely a function of $\theta$ and the likelihood you've written there is only a function of $\rho_i$ (assuming the data is fixed). For each likelihood evaluation (i.e. for a given $\rho_i$), you will need to calculate the integrals for each $j$ and sum them.

If these integrals do not have a closed form (which, evidently, they do not), you will need to approximate them to calculate the likelihood. If you can calculate $P(x_{ij}|\theta, \rho_i) \cdot f(\theta)$ as a function of $\theta$, then you can do a numerical approximation (read about Quadrature for some methods).

If you can easily simulate from $f$ then another approach that may be easier to implement is to estimate the integral by monte carlo. This is done by generating $\tilde{\theta}_1, ..., \tilde{\theta}_N$ from $f$ and calculating $$\widehat I = \frac{1}{N} \sum_{k=1}^{N} P(x_{ij}|\tilde{\theta}_k, \rho_i)$$

which will converge to the true integral (assuming it exists) as $N \rightarrow \infty$ and, for finite $N$, the mean squared error of this approximation equals ${\rm var} \big(P(x_{ij}|\Theta, \rho_i) \big)/N$. To see why this is true, notice that the integral may be thought of as the expectation of the random variable $P(x_{ij}|\Theta, \rho_i)$, with respect to $\Theta$ and $\widehat{I}$ is the sample mean of $N$ realizations of $P(x_{ij}|\Theta, \rho_i)$.