For a paper describing the origins of the principle of Maximum Lilkelihood see Aldrich, 1997

** 1. Closed Form Solutions** The liklelihood function is the probability of the data given the parameters viewed as a function of the parameters with data fixed. Maximum likelihood Estimates (MLEs) are the values of the parameters that maximize the likelihood function.

In simple models (e.g., Poisson, Bernoully, multiple-linear regression with Gaussian IID error terms) we can find a closed-form for the ML estimates. In these cases MLEs are obtained by:

• Taking the 1st derivative of the log-likelihood with respect to each of the parameters
• (First Order Conditions, FOC) Set all the derivatives equal to zero, this gives a stationary point and renders a system with as many equations as parameters.
• Solve for the parameter values that satisfy the FOC.
• Check, based on the sign of the 2nd derivatives that the function is concave.

However, in the vast majority of the cases we cannot find a closed-form for the ML estimates. In these cases we use numerical methods. There are many numerical methdos for optimizing a function, we will briefly consider three approaches: grid search, Newton-Rapson and general purpouse optimization algorithm (functions optim and optimize in R).

1. Grid search

In it's simplest form a grid search requires:

• Implementing a function for evaluating the log-likelihood
• Defining a grid of values for the parameters
• Evaluating the log-likelihood for every value in the grid
• Finding the value within the grid that yields the largest value of the likelihood

The following R-code illustrates this for the case of Poisson data

NOTE: In this case the MLE has a closed form (it is simply the sample mean)

## Simulating data
 trueLambda=10
 y=rpois(lambda=trueLambda,n=100)

## A function to evaluate the Poisson log-likelihood
 logLikPois=function(y,lambda){     
      logLik=log(lambda)*sum(y)-length(y)*lambda
      return(logLik)
 }

## A grid
myGrid=seq(from=.01,to=100,length=10000)

## Search
logLik=rep(NA,length(myGrid))
for(i in 1:length(myGrid)){
 logLik[i]=logLikPois(y=y,lambda=myGrid[i])
 print(i)
}
plot(logLik~myGrid,type='l',col=4)
abline(v=mean(y),lty=2,col=2)

If we have a parameter vector instead of a single parameter our grid will be multi-dimensional (as many dimensions as parameters) but the princple of the grid search is the same.

Grid Search is computationally intensive, most of the values of the grid could be ruled out y looking at the derivatives of the function.

The above approach very crude. If we know that our function is concave, the we can set a very sparse grid, then find three points where the interior point gives higher likelihood than the flanking points. This interval must contain the maximum. We can then refine our grid by making a denser grid within that interval, and interate in this fashion untlil our estimate does not change more than a pre-established precision.

2 Newton Rapson (NR) Method

This method can be used to find a stationary point (either a minima or maxima) of a twice differentiable function. Starting from an intial guess theta_0, the NR method moves from theta_0 to tehta_1=theta_0+delta_0 where delta_0 is the negative of the ratio of the first and seconde derivatives of the objective function both evaluated at theta_0. The method can be motivated using a 2nd order Taylor expansion to the objective function. We will discuss this further in class.

The following example illustrates the aplictaion of the NR methods for estimating the parameter of a Poisson distribution.

## Data
 ## Simulating data
  trueLambda=10
  y=rpois(lambda=trueLambda,n=100)


 ## Derivatives
  firstD=function(lambda,y){  
     sum(y)/lambda -length(y)
  }
  
  secondD=function(lambda,y){
     -sum(y)/(lambda^2)
  }
  
  
  theta=3
  tol=1e-5
  ready=FALSE
  counter=0
  while(!ready){
  	print(theta)
  	delta=- firstD(theta,y)/secondD(theta,y)
	new_theta=theta+delta
	ready=(abs(theta-new_theta)<tol)
	counter=counter+1
	print(paste0(counter,'  ', new_theta))
	theta=new_theta
  }

3 General Purpose Optimization Functions

In this section we illustrate how to obtain MLEs using optim(), this is a general purpose optimizaton function that can be used to minimize an objective function with respect to a set of parameters. The inputs for optim are: (i) the function to be minimized (e.g., negative-log-likelihood), (ii) a parameter vector (all the parameters to be estimated must be provided in a single vector), and (iii) additional arguments to the function that are to be kept fixed during the optimization procedure (e.g., data, covariates, etc.). We illustrate the use optim in the context of logistic regression and censored regression.

3.1. Logistic Regression

Toy simulation

 x=rnorm(1000)
 risk=2+x*2
 prob=exp(risk)/(1+exp(risk))
 y=rbinom(prob=prob,n=length(x),size=1)

Logistic Regression using GLM

 fm0=glm(y~x,family=binomial)
 fm0$coef

Logistic regression using optim

 fm1=optim(fn=negLogLikLogistic,y=y,X=cbind(1,x),par=c(log(mean(y)/(1-mean(y))),0))
 fm1$par

Inference

Under regularity conditons, the large-sample distribution of ML estimates is Multivariate Normal with mean equal to the true parmeter values (i.e., MLs are asymptotically un-biased) and (co)variance matrix equal to the negative of the inverse of Fisher's Information matrix which is the expected value of the matrix of second derivatieves of the log-likelihood function with respect to the parameters. This matrix often depends on the true parameter values which are uknown. Thus, in practice, Expected Information is replaced with the Observed Information which is the inverse of the matrix of the second derivatives of the log-likelihood evaluated at the ML estimates. We can obtain this matrix from optim() by setting hessian=TRUE. The followin example illsutrates this.

 fm1=optim(fn=negLogLikLogistic,y=y,X=cbind(1,x),par=c(log(mean(y)/(1-mean(y))),0),hessian=T)
 OI=solve(fm1$hessian)
 cbind(fm1$par,sqrt(diag(OI)))
 summary(fm0)$coef# GLM

3.2. Censored Regression

Biomedical and engeneering data is usually censored. Censored means that for some subjects we do not know the actual value of the response but we have information (bounds) about it. There are three types of censoring:

• Right censored, we know that y is greater than some constant, y>c.
• Left censored, we know that y is smaler than some constant, y<c.
• Interval censored, means that we know that y lies in an interval, c1<y<c2.

As always, the likelihood function is the probability of the data given the parameters viewed as a function of the parameters with data fixed. For un-censored data points this is just the density function. For right-censored data point the likelihood is the Survival function (i.e., S(c)=p(y>c)=1-CDF(c)). For left-censored data points the likelihood is the CDF, p(y<c)=CDF(c), and for interval censored the likelihood is the difference in the CDF function evaluated at the right and left-bounds of the interval, p(c1<y<c2)=CDF(c2)-CDF(c1). The following example illustrates how to compute ML estimates with censored data assuming that data follows a Gamma distribution.

Suppose we have gamma data

  x=rgamma(shape=8,rate=.5,n=100)

ML estimation with the complete data

 negLogLikGamma=function(x,theta){
	logLik=dgamma(x=x,shape=exp(theta[1]),rate=exp(theta[2]),log=TRUE)
	return(-sum(logLik))
 }

 fm=optim(fn=negLogLikGamma,x=x,par=log(c(2,2)))
 exp(fm$par)
 

Suppose now that a proportion of the data are right-censored

 event=runif(length(x))<.7
 y=x
 y[!event]=x[!event]-runif(n=sum(!event),min=1,max=3)

Log-likelihood function for right-censored data

 negLogLikGammaRightCen=function(time,theta){
 	logLik=pgamma(q=time, shape=exp(theta[1]), rate = exp(theta[2]), lower.tail = FALSE, log.p = TRUE)
 	 return(-sum(logLik))
 }

 logLikGammaCen=function(time,event,theta){
 	
 		time_observed=time[event]
 		time_right_cen=time[!event]
 		
 		negLogLik_1=negLogLikGamma(x=time_observed,theta)
 		negLogLik_2=negLogLikGammaRightCen(time=time_right_cen,theta=theta)
 		return((negLogLik_1+negLogLik_2))
 }
 
 fm2=optim(fn=logLikGammaCen,time=x,event=event==T,par=log(c(2,2)))
 
  exp(fm2$par)

The survaival R-package

negLogLikRightCenNormal=function(y,d,theta){
   mu=theta[1]
   V=exp(theta[2])
   yObserved=y[d==1]
   yCensored=y[d==0]

   logLik_observed=sum(dnorm(x=yObserved,mean=mu,sd=sqrt(V),log=TRUE))
   logLik_censored=sum(pnorm(q=yCensored,mean=mu,sd=sqrt(V),lower.tail=FALSE,log=TRUE))
   return(-(logLik_observed+logLik_censored))
}


y=rnorm(n=1000,sd=2,mean=10)

d=ifelse(y>11,0,1)
time=y

time[d==0]=11

fm=optim(fn=negLogLikRightCenNormal,y=time,d=d,par=c(mean(time),log(var(time))))

Back to Outline