Search code examples
rmathequation

Converting a likelihood function into R code

I am needing to use the following likelihood equation to solve for "r" in a dose response experiment. In this equation, for i = 1,2,.....j, Ni is the number of individuals exposed to dose Di, and Yi the number that became infected.

Maximum likelihood equation

I have managed to construct the dataframe for Ni, Di and Yi, but as I am new to MLE, I am struggling to convert the equation into R code.

Would anyone be able to assist in the code conversion part?

Thanks

PS. Following the suggestion of Rui Barradas below, here is the dataset the paper which developed the equation used:

Di <- c(4.2,8,10.6,11.9,12.2,12.6,13,13.2,28.3,28.7,68.2,69,71.1,83,91.7,95.5,106.5,136.5,171.2,186.9,309.3,1557.1)
Ni <- c(1,1,1,1,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1)
Yi <- c(0,0,1,1,2,1,0,1,0,1,1,0,1,1,1,1,1,1,1,1,1,1)
df <- data.frame(Di, Ni, Yi)
Solution

  • You can optimize for the r parameter with the following code :

    log_Lik <- function(r, D, N, Y)
{
   return(sum(Y * log(1 - exp(-r * D))) - r * sum(D * (N - Y)))
}

optimise(f = log_Lik, lower = 0, upper = 5, maximum = TRUE, D = D, N = N, Y = Y)