CODEWITHSUNDEEP
r
Mamba
Mamba

Reputation: 1203

Understanding different results of optim() and lm()

Given:

   set.seed(1001)
   outcome<-rnorm(1000,sd = 1)
   covariate<-rnorm(1000,sd = 1)

log-likelihood of normal pdf:

   loglike <- function(par, outcome, covariate){
               cov <- as.matrix(cbind(1, covariate))
               xb <- cov * par
              (- 1/2* sum((outcome - xb)^2))
   }

optimize:

                  opt.normal <- optim(par = 0.1,fn = loglike,outcome=outcome,cov=covariate, method = "BFGS",  control = list(fnscale = -1),hessian = TRUE)

However I get different results when running an simple OLS. However maximizing log-likelihhod and minimizing OLS should bring me to a similar estimate. I suppose there is something wrong with my optimization.

                  summary(lm(outcome~covariate))

Upvotes: 0

Views: 1044

Answers (2)

peterweethetbeter
peterweethetbeter

Reputation: 21

Maybe you will find it usefull to see the difference between these two methods from my code. I programmed it the following way.

data.matrix <- as.matrix(hprice1[,c("assess","bdrms","lotsize","sqrft","colonial")])
loglik <- function(p,z){
  beta <- p[1:5]
  sigma <- p[6]
  y <- log(data.matrix[,1])
  eps <- (y - beta[1] - z[,2:5] %*% beta[2:5])
  -nrow(z)*log(sigma)-0.5*sum((eps/sigma)^2)
}
p0 <- c(5,0,0,0,0,2)
m <- optim(p0,loglik,method="BFGS",control=list(fnscale=-1,trace=10),hessian=TRUE,z=data.matrix)


rbind(m$par,sqrt(diag(solve(-m$hessian))))

And for the lm() method I find this 

m.ols <- lm(log(assess)~bdrms+lotsize+sqrft+colonial,data=hprice1)
summary(m.ols)

Also if you would like to estimate the elasticity of assessed value with respect to the lotsize or calculate a 95% confidence interval for this parameter, you could use the following

elasticity.at.mean <- mean(hprice1$lotsize) * m$par[3]
var.coefficient <- solve(-m$hessian)[3,3]
var.elasticity <- mean(hprice1$lotsize)^2 * var.coefficient
# upper bound
elasticity.at.mean + qnorm(0.975)* sqrt(var.elasticity)
# lower bound
elasticity.at.mean + qnorm(0.025)* sqrt(var.elasticity)

A more simple example of the optim method is given below for a binomial distribution.

loglik1 <- function(p,n,n.f){
  n.f*log(p) + (n-n.f)*log(1-p)
}

m <- optim(c(pi=0.5),loglik1,control=list(fnscale=-1),
           n=73,n.f=18)
m

m <- optim(c(pi=0.5),loglik1,method="BFGS",hessian=TRUE,
           control=list(fnscale=-1),n=73,n.f=18)
m

pi.hat <- m$par

numerical calculation of s.d

rbind(pi.hat=pi.hat,sd.pi.hat=sqrt(diag(solve(-m$hessian))))

analytical

rbind(pi.hat=18/73,sd.pi.hat=sqrt((pi.hat*(1-pi.hat))/73))

Or this code for the normal distribution.

loglik1 <- function(p,z){
  mu <- p[1]
  sigma <- p[2]

  -(length(z)/2)*log(sigma^2) - sum(z^2)/(2*sigma^2) + 
    (mu*sum(z)/sigma^2) - (length(z)*mu^2)/(2*sigma^2)
}

m <- optim(c(mu=0,sigma2=0.1),loglik1,
           control=list(fnscale=-1),z=aex)

Upvotes: 0

Spacedman
Spacedman

Reputation: 94317

Umm several things... Here's a proper working likelihood function (with names x and y):

loglike =
function(par,x,y){cov = cbind(1,x); xb = cov %*% par;(-1/2)*sum((y-xb)^2)}

Note use of matrix multiplication operator.

You were also only running it with one par parameter, so it was not only broken because your loglike was doing element-element multiplication, it was only returning one value too.

Now compare optimiser parameters with lm coefficients:

opt.normal <- optim(par = c(0.1,0.1),fn = loglike,y=outcome,x=covariate, method = "BFGS",  control = list(fnscale = -1),hessian = TRUE)
opt.normal$par

[1]  0.02148234 -0.09124299

 summary(lm(outcome~covariate))$coeff

               Estimate Std. Error    t value    Pr(>|t|)
(Intercept)  0.02148235 0.03049535  0.7044466 0.481319029
covariate   -0.09124299 0.03049819 -2.9917515 0.002842011

shazam.

Helpful hints: create data that you know the right answer for - eg x=1:10; y=rnorm(10)+(1:10) so you know the slope is 1 and the intercept 0. Then you can easily see which of your things are in the right ballpark. Also, run your loglike function on its own to see if it behaves as you expect.

Upvotes: 5

Related Questions