Example of E-M: Genetic Linkage Problem

Rao (1973) considers 197 animals randomly divided into four categories as follows:

Y = (y_{1}, y_{2}, y_{3}, y_{4})^{T} = (125, 18, 20, 34)^{T}

with cell probabilities

(1/2 + θ/4, (1 − θ)/4, (1 − θ)/4, θ/4)^{T}

The observed-data likelihood is

g(Y |θ) ∝ ·(2 + θ)^{y_{1}}(1 − θ)^{y_{2}+y_{3}}θ^{y_{4}}

If we “augment” the data by splitting the first category into two,

X = (x_{1}, x_{2}, x_{3}, x_{4}, x_{5})^{T}
with probabilities  (1/2, θ/4, (1 − θ)/4, (1 − θ)/4, θ/4)^{T}

the complete-data likelihood is now simpler:

f(X|θ) ∝ θ^{x_{2}+x_{5}}(1 − θ)^{x_{3}+x_{4}}

The “missing data” Z is like a coin toss for which category each observation of y_{1} gets split into:

(Z|Y, θ) ∼ Bin(y_{1}, 2/(2 + θ))

Now we compute the Q-function,

Q(θ, θ^{k}) = E[(x_{2} + x_{5}) log(θ) + (x_{3} + x_{4}) log(1 − θ) | Y, θ^{k}]
= (y_{1}θ^{k}/(2 + θ^{k}) + x_{5}) log(θ) + (x_{3} + x_{4}) log(1 − θ)

This is of the form a log(θ) + b log(1 − θ), a binomial whose MLE is θ̂ = a/(a + b); translating back to the particulars of Q(·, ·) we get

θ^{k+1} = \frac{y_{1}θ^{k}/(2 + θ^{k}) + x_{5}}{y_{1}θ^{k}/(2 + θ^{k}) + x_{5} + x_{3} + x_{4}}

Now we attempt to find the MLE of θ, two ways:

• By Newton-Raphson

     > Y _c(125, 18, 20, 34)
     > g _function(th) {
     + (2 + th)^(125) * (1-th)^(18+20) * th^(34)
     + }
     > ng _function(th) {-log(g(th))}
     > summary(ms( ~ ng(th),start=c(th=0.5)))
     Final value: -67.3841
     Solution:
             Par.

th 0.6268215
Convergence: RELATIVE FUNCTION CONVERGENCE.
Computations done:
 Iterations Function Gradient
          6        8        9
> source("infomat.s")
> info.mtx.fun(ng,0.6268215)
          [,1]
[1,] -377.5169
> sqrt(1/377.5169)
[1] 0.05146735

• By E-M

> emfun <-
function(th0, maxit = 1000, releps = 1e-09)
{
        i <- 0
        th1 <- th0
        th0 <- th1 + 1
        while((i != maxit) && (abs(th1 - th0) > releps * abs(th0))) {
                i <- i + 1
                th0 <- th1
                th1 <- ((125 * th1)/(2 + th1) + 34)/((125 * th1)/(2 + th1) +
                        34 + 18 + 20)
                print(c(th0,th1,i))
        }
        return(i, th1)
}
> emfun (0.5)
[1] 0.5000000 0.6082474 1.0000000
[1] 0.6082474 0.6243211 2.0000000
[1] 0.6243211 0.6264889 3.0000000
[1] 0.6264889 0.6267773 4.0000000
[1] 0.6267773 0.6268156 5.0000000
[1] 0.6268156 0.6268207 6.0000000
[1] 0.6268207 0.6268214 7.0000000
[1] 0.6268214 0.6268215 8.0000000
[1] 0.6268215 0.6268215 9.0000000
[1]  0.6268215  0.6268215 10.0000000
[1]  0.6268215  0.6268215 11.0000000
$i:
[1] 11
$th1:
[1] 0.6268215