0

I'm trying to do MLE in R to simulate mark recapture methods. I'm trying to maximise this function;

likelihood <- function(N, p){
  likelihood <- 1
  X <- 0
  for (s in 1:S){
    X <- X + M[s]
  }
  likelihood <- likelihood * p^(X) * (1-p)^(N*S-X)
  likelihood <- likelihood * (factorial(N)/factorial((N-length(U))))
  return(likelihood)
}

I'm trying to maximise the N, all the other parameters are known or estimated and so assumed known

Have looked into nlm and optim but can't get either to do what I want...

In response to answers; Thanks for all the answers, I recoded the likelihood and managed to get something working, see the code below. To answer questions/answers specifically 1. Thanks for X <- sum(M[1:S]) nice trick, The X is counting 'all the animals ever caught (not unique captures)' 2. M is generated in the code below, based on n (captures) it is a count of number of animals marked on each survey, S is number of surveys.

 captures <- function(N, S, P){
   P <- replicate(S, P) #Capture Probability (same across animals and surveys)
   captures <- t(replicate(N, rbinom(S, 1, P))) #Generate capture data from N animals with S surveys with capture Probability P 
   return(captures)
 }

 marked <- function(N, S, captures){ #Count numbers that were marked on each survey
   M <- replicate(S, 0) #Initialise the 'marked' variable to zero
   for (s in 1:S){
     for (i in 1:N){
       if (captures[i, s] == 1){
         M[s] <- M[s] + 1
       }
     }
   }
   return(M)
 }

 unseen <- function(N, S, captures){#Count total number of animals that were unseen on any survey
   U <- array()
   seen <- replicate(N, FALSE) #All animals begin as unseen
   for (i in 1:N){
     for (s in 1:S){
       if (captures[i, s] == 1){
         seen[i] = TRUE
       }
     }
     if (seen[i] == FALSE){
       U <- c(U, i)
     }
   }  
   U <- U[-1] #fix the N/A in the first index
   return(U)
 }

 firstcaught <- function(N, S, captures){#Produce a vector containing how many animals were   caught for the first time on each survey
   fc <- 0
   seen <- replicate(N, FALSE)
     for (s in 1:S){
       fc[s] <- 0
       for (i in 1:N){
         if ((seen[i] == FALSE) && (captures[i, s] == 1)){
             fc[s] <- fc[s] + 1
             seen[i] = TRUE
         }
       }
     }
   return(fc)
 }

 ##Generate data using functions
 N <- 200
 S <- 10
 P <- 0.2
 n <- captures(N,S,P)
 M <- marked(N,S,n)
 U <- unseen(N,S,n)
 fc <- firstcaught(N,S,n) 
 Mc <- cumsum(fc) #A running total of animals that have been captured at least once


 ##Define a likelihood for model M_0
 likelihood <- function(N, P){
  likelihood <- 1
   X <- 0
  for (s in 1:S){
   X <- X + M[s]
    likelihood <- likelihood * (1-P)^(N-M[s])
    likelihood <- likelihood * choose(N-M[s],fc[s])
  }
  likelihood <- likelihood * P^(X) * (1-P)^(-X)
  return(-log(likelihood))
}

##Find the MLE
out <- nlm(f = likelihood, p = 200, P = P, check.analyticals = TRUE)
Moohan
  • 933
  • 1
  • 9
  • 27
  • 1
    You should add what you tried. Normally what you want when using optimization methods in R though is a single named input that contains all of the parameters. So you would have something like `likelihood <- function(pars, data){N <- pars[1]; p <- pars[2]; insert_other_code}`. I'm not sure what you're doing with X there but `X <- sum(M[1:S])` would accomplish the same goal. Note that you're using values defined outside the function with that though and that isn't necessarily good practice. – Dason Apr 20 '14 at 16:24
  • 3
    To make this a [reproducible example](http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example), you would need to provide at least samples of `M` and `S`. Also, are to trying to find N which maximizes likelihood for a given p, or are you trying to "maximize the N" as your question states? – jlhoward Apr 20 '14 at 16:37
  • Since you are maximizing over N only, you could use `optimize` function. Provide all the values you use (S, U, M ...) so we can better help you. – JACKY88 Apr 20 '14 at 19:48

0 Answers0