Problems optimizing a likelihood function in R

Question

I'm trying to maximize the log-likelihood $y_1 \log(p_1) + (1-y_1)\log(p_1) + y_2\log(p_2) + (1-y_2)\log(p_2)$. I have data on success of an experiment. Y identifies if person died or not, and X identifies control or treatment group.

What's wrong with my code?

# MLE for the likelihood
y <- c(rep(1,39), rep(0,674-39), rep(1,22), rep(0,680-22))
x <- c(rep(0, 674), rep(1, 680))

binom.ll <- function(theta, y, x) {
  y[x==0]*log(theta[1]) + (1-y[x==0])*log(1-theta[1])) + y[x==1]*log(theta[2]) + (1-y[x==1])*log(1-theta[2]))
}

theta.start <- c(0, 0)
ml.res <- nlm(binom.ll, theta.start, print.level=1, y=y, x=x, hessian=T)
ml.res

iteration = 0
Step:
[1] 0 0
Parameter:
[1] 0 0
Function Value
[1] 1.797693e+308
Gradient:
[1] -Inf    0

Error in nlm(binom.ll, theta.start, print.level = 1, y = y, x = x, hessian = T): non-finite value supplied by 'nlm'
In addition: Warning messages:
1: In nlm(binom.ll, theta.start, print.level = 1, y = y, x = x, hessian = T): NA/Inf replaced by maximum positive value
2: In nlm(binom.ll, theta.start, print.level = 1, y = y, x = x, hessian = T): NA/Inf replaced by maximum positive value

Answer 1

I assume $p_1 = \mathbb P(Y=1 | X=0)$ and $p_2 = \mathbb P (Y=1|X=1)$.

The log-likelihood you want to maximize is $$ n_{10} \log p_1 + n_{00} \log(1-p_1) + n_{11} \log p_2 + n_{01} \log(1-p_2) $$ where $n_{ij}$ is the number of observations with $Y=i, X=j$.

You can do that directly, no need for numerical optimization. The optimum is $p_1 = {n_{10} \over n_{10} + n_{00}}$ and $p_2 = {n_{11} \over n_{11} + n_{01}}$.

Edit: if you really want R code, here it is.

y <- c(rep(1,39), rep(0,674-39), rep(1,22), rep(0,680-22))
x <- c(rep(0, 674), rep(1, 680))

binom.ll <- function(theta, y, x) 
{
  n <- table(y,x);
  th <- matrix( c(1-theta,theta), nrow=2, byrow=TRUE);
  return( -sum(n*log(th)) );
}


> theta.start <- c(0.5, 0.5)
> nlm(binom.ll, theta.start, y=y, x=x, hessian=T)
$minimum
[1] 246.1096

$estimate
[1] 0.05786353 0.03235302

$gradient
[1] -0.0002950512 -0.0004371658

$hessian
        [,1]     [,2]
[1,] 12323.5     0.00
[2,]     0.0 21591.73

$code
[1] 2

$iterations
[1] 21

Il y a eu 28 avis (utilisez warnings() pour les visionner)

> 39/674
[1] 0.0578635
> 22/680
[1] 0.03235294