R Different Results Calculating LLike From LM and Optim

Question

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.

EDIT

Thanks the suggestions provided by Tim and the solution from here it tunes out there were two mistakes:

1) in the llike function

 sigma_sq <- par[1+n]
 #replace by
 sigma_sq <- par[1+n]^2

2) the optimization algorithm:

Change from "BFGS" to "L-BFGS-B" and obtain the desired results.

Answer 1

Few comments:

1. You don't need to implement the normal density by hand, you can just use build-in function. Using it is a good way to make sure that the results are not due to typo in the code.

llike <- function(par, y, X) {
    n <- ncol(X)
    beta <- par[1:n]
    sigma <- sqrt(par[1+n])
    -sum(dnorm(y, X %*% beta, sigma, log=TRUE))
}

2. As you could see from the above code, likelihood is not distribution for the errors, but distribution for the data. In linear regression model, we predict mean of the normal distribution given some features.
3. If you used the above code, you'd also see that the estimates are not exactly the same. Moreover, if you used different optimization algorithm, they could change again. lm function uses ordinary least squares for estimating the parameters, you used BFSG algorithm, there is no guarantee that they will give you exactly the same results.
4. For maximizing the likelihood you need only the unnormalized log-densities, i.e. $-(x-\mu)/2\sigma^2$.

Check also the Maximum Likelihood Estimation (MLE) in layman terms thread.