Likelihood Function Core code

Generating some data with a know probability of germination -> myP

myP<-0.67
myData<-rbinom(n=10,size=1,prob=myP)
myData

 [1] 1 1 1 0 0 0 1 1 1 1

mean(myData)

[1] 0.7

Here is the form of the Likelihood Functions that we have been using to fit models. This is specific to the Binomial probability function.

nllBinomial<-function(parmEstimate,seedlingDat){

  nllik<- -sum(dbinom(seedlingDat,size=1,parmEstimate,log=TRUE))

  cat("nllik= ",nllik,sep=" ",fill=T)
  return(nllik)
}

nllBinomial(0.6,myData)

nllik=  6.324652
[1] 6.324652

nllBinomial(0.1,myData)

nllik=  16.43418
[1] 16.43418

We would then minimize the function using optim().

parVec<-c(0.1) # Initial parameter values 
outLogReg<-optim(par=parVec,fn=nllBinomial,method="L-BFGS-B",lower=c(0.01),upper=c(0.99),seedlingDat=myData)

nllik=  16.43418
nllik=  16.36786
nllik=  16.5012
nllik=  13.88586
nllik=  13.88586
nllik=  13.60701
nllik=  10.07198
nllik=  10.05039
nllik=  10.09367
nllik=  8.576073
nllik=  8.560991
nllik=  8.591218
nllik=  6.609054
nllik=  6.602835
nllik=  6.61531
nllik=  6.11829
nllik=  6.117376
nllik=  6.119248
nllik=  6.108907
nllik=  6.109091
nllik=  6.108772
nllik=  6.108643
nllik=  6.108661
nllik=  6.108673
nllik=  6.108643
nllik=  6.108667
nllik=  6.108667
nllik=  6.108643
nllik=  6.108667
nllik=  6.108667

outLogReg$par # 

[1] 0.6999994

The core portion of the likelihood function is this line of code:

nllik<- -sum(dbinom(seedlingDat,size=1,parmEstimate,log=TRUE))

I want to step through this line of code so that its clear what its doing.

dataVector<-c(1,0,1)

where 1 indicates a seed produced a seedling and 0 indicates no seedling emerged.

So the likelihood is the likelihood of observing a 1, 0, and 1 given some germination probability,G.

Likelihood = P(1|G) * P(0|G) * P(1|G) where the P() is the binomial distribution. For the binomial distribution, we have to specify the size of the trials. For each of our ‘trials’, we plant a single seed and observe if a seedling emerges, So the size of our trial is 1, because we plant a single seed.

We next choose a value of G, then calculate the likelihood of observing our data as follows:

G<-0.10
dbinom(1,size=1,G) * dbinom(0,size=1,G) * dbinom(1,size=1,G) 

[1] 0.009

G<-0.7
dbinom(1,size=1,G) * dbinom(0,size=1,G) * dbinom(1,size=1,G) 

[1] 0.147

G<-0.9
dbinom(1,size=1,G) * dbinom(0,size=1,G) * dbinom(1,size=1,G) 

[1] 0.081

By taking the logs of each of these probabilities, we add them together rather than multiply them. This helps avoid numerical issues associated with getting very small numbers from multiplying long sequences of probabilities together. Eventually, the computer will just say that the final product is 0.

So we modify the above code by setting log=T and then adding the terms together.

G<-0.10
dbinom(1,size=1,G, log=T) + dbinom(0,size=1,G, log=T) + dbinom(1,size=1,G, log=T) 

[1] -4.710531

G<-0.7
dbinom(1,size=1,G, log=T) + dbinom(0,size=1,G, log=T) + dbinom(1,size=1,G, log=T)

[1] -1.917323

G<-0.9
dbinom(1,size=1,G, log=T) + dbinom(0,size=1,G, log=T) + dbinom(1,size=1,G, log=T)

[1] -2.513306

Finally, we want to be able to compute the likelihood over a set of alternative values of G. The value of G that produces the highest likelihood is the most likely value or maximum likelihood estimate (mle).

In the code below, I create a sequence of values of G, then calculate the likelihood associated with each of these values.

Note that I am using the map_dbl function from the purrr library. The map_dbl function applies the myLike function to each value of G.

library(purrr)

G<-seq(from=0, to=1,by=0.10)

myLike<-function(G) {dbinom(1,size=1,G, log=T) + dbinom(0,size=1,G, log=T) + dbinom(1,size=1,G, log=T) }

map_dbl(G,myLike)

 [1]      -Inf -4.710531 -3.442019 -2.764621 -2.343407 -2.079442 -1.937942 -1.917323 -2.055725 -2.513306      -Inf

Since the dbinom function can be vectorized, we can simplify this calculation as follows

G<-0.1
dataVector<-c(1,0,1)
sum(dbinom(dataVector,size=1, G, log=T))

[1] -4.710531

so that our code above becomes

library(purrr)

dataVector<-c(1,0,1)
G<-seq(from=0, to=1,by=0.10)

myLike2<-function(G) {sum(dbinom(dataVector,size=1, G, log=T)) }

map_dbl(G,myLike2)

 [1]      -Inf -4.710531 -3.442019 -2.764621 -2.343407 -2.079442 -1.937942 -1.917323 -2.055725 -2.513306      -Inf

Finally, we add a minus sign in front of the sum to calculate the negative log likelihood, so that we can minimize the function using optim() rather than maximize it.

library(purrr)

dataVector<-c(1,0,1)
G<-seq(from=0, to=1,by=0.10)

myLike2<-function(G) {-sum(dbinom(dataVector,size=1, G, log=T)) }

nllVector<-map_dbl(G,myLike2)
nllVector

 [1]      Inf 4.710531 3.442019 2.764621 2.343407 2.079442 1.937942 1.917323 2.055725 2.513306      Inf

Now we can plot the likelihood function. Its not smooth because we have a small number of values of G that we are using. The more values of G that we use, the smoother the function will become.

plot(G,nllVector,typ='l')

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Cmd+Shift+Enter.

plot(cars)

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Cmd+Option+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Cmd+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.