1

I want to calculate the loglikelihood for multivariate linear regression. I'm not sure whether this code is true or not.

I’ve been calculated the log likelihood using dmvnorm function in mvtnorm r package.

sdmvn_mle <- function(obj){
  sdmvn_mle_1 <- apply(obj$residuals^2,2,mean)
  sdmvn_mle_2 <- mean(residuals(obj)[,1] * residuals(obj)[,2])
  return(matrix(c(sdmvn_mle_1[1], sdmvn_mle_2, sdmvn_mle_2, sdmvn_mle_1[2]), nrow = 2))
} 


llmvn <- function(obj, sd){
  lr <- c()
  for( i in 1: nrow(obj$fitted.values)){
    lr <- c(lr, mvtnorm::dmvnorm(model.response(model.frame(obj))[i,], mean=fitted(obj)[i,], sigma=sd, log=TRUE))
  }
  return(sum(lr))
}

Y <- as.matrix(mtcars[,c("mpg","disp")])
(mvmod <- lm(Y ~ hp + drat + wt, data=mtcars))
# Call:
# lm(formula = Y ~ hp + drat + wt, data = mtcars)

# Coefficients:
#             mpg        disp     
# (Intercept)   29.39493   64.52984
# hp            -0.03223    0.66919
# drat           1.61505  -40.10238
# wt            -3.22795   65.97577

llmvn(mvmod, sdmvn_mle(mvmod))
# [1] -238.7386

I’m not sure the result is correct or not.

Additionally, Please let me know if there is another strategies for calculating log likelihood for multivariate linear regression.

JungHyun jung
  • 11
  • 3
  • 1
    Maybe the question is better vor CrossValidated. –  Jul 11 '19 at 13:57
  • The usual approach is to use the `optim` function or similar. Start with a simpler example, standard normal, checking results produced using optim correspond with lm, then think how to extend to multivariate normal case, taking into account issues with local optima on the likelihood surface and why this may require a non-default method to be specified. – JonMinton Jul 11 '19 at 14:36
  • This may offer the solution you're looking for: https://gist.github.com/cheuerde/f6ad53b4aadcce23df35 However, please start by learning how optim works for simpler likelihood methods, in particular the need to put all parameters into a 'par' vector for optim to work. – JonMinton Jul 11 '19 at 14:39

0 Answers0