Skip to main content
1 of 2

R Different Results Calculating LLike From LM and Optim

I am trying to implement a multivariate linear regression model by ML estimation, however I ran into some discrepancies.

An assumption of the model is that the residual is a standard normal variable.

I followed this(page 6) approach, but the procedure should be clear anyway:

$$ \text{pdf}_X (x,\mu,\sigma) = \frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{(x-\mu)^2}{2\sigma^2}} $$

$$ \text{pdf}_X (x_1,...,x_n,\mu,\sigma) = \prod_{i=1}^n\frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{(x_i-\mu)^2}{2\sigma^2}} = (\frac{1}{\sqrt{2\pi\sigma^2}})^ne^{-\frac{\sum_{i=1}^n(x_i-\mu)^2}{2\sigma^2}} $$

Taking the log gives us the log likelihood:

$$ log((\frac{1}{\sqrt{2\pi\sigma^2}})^ne^{-\frac{\sum_{i=1}^n(x_i-\mu)^2}{2\sigma^2}}) = -\frac{n}{2}log(2\pi) -\frac{n}{2}log(\sigma^2) - \frac{\sum_{i=1}^n(y - X\beta)^2}{2\sigma^2} $$

Translating this into R Code:

set.seed(1)

# generate random data 
X <- cbind(rep(1, 100/4), matrix(rnorm(100), ncol=4, byrow=T))
#true beta 
beta <- c(25, 1,3,5,7)
# linear model
y <- X%*%beta + rnorm(25)

llike <- function(par, y, X){
  m <- nrow(X)
  n <- ncol(X)
  beta <- par[1:n]
  sigma_sq <- par[1+n]
  e <- y - X%*%beta
  loglike <- -(m/2)*log(2*pi) - (m/2)*log(sigma_sq) - ((t(e)%*%e) / (2*sigma_sq))
  return(-loglike)
}

res <- optim(par=c(rep(1, 6)), llike, method="BFGS", hessian=T, y=y, X=X)

res$par
24.189 1.209 3.729 5.489 7.102 8788.176

The last term is the variance it seems excessive, but keep in mind that the "rnorm" function generates normal distributed data, with mean = 0 and standard deviation = 1. Furthermore we have 25 of those so the estimated standard deviation would actually be:

 sqrt(res$par[6])/25
 3.749

Now when I calculate the logliklihood:

-llike(c(res$par[1:5], sqrt(res$par[6])/25), y, X)
-45.56091

compared to:

logLik(lm(y ~ X-1))
-42.238

So there are two issues that I have:

  1. The estimates between the optim and lm are different
  2. the loglikelihood is differes.

Am I missing something with the implementation? I really need to know if there is a mistake, since this model is the bases for Lasso and Ridge Regression which are supposed to follow.