2

The code is this:

lm.loglik3 <-function(y, x, theta){N <-length(y)# theta contains our parameters to be estimated
beta0 <- theta[1]
beta1 <- theta[2]
sigma2 <- theta[3]# x contains our covariates
x1 <- x
logl <--N/2* log(2*pi*sigma2)-1/(2*sigma2)* sum((beta0+beta1*x1-y)^2)
return(-logl)}
# Starting Values
stval <-c(0, 0, 0)

# Optim
res3 <-optim(stval, lm.loglik3, y = df2$v64, x = df2$auth_val, hessian = TRUE)

I'm trying to do a maximum likelihood estimation using as parameters the degree of authoritarianism (numerical value are given) and the income. The database contains 6000 sample. Now the problem is that when i run the code it tells me "Error in optim: function cannot be evaluated at initial parameters". Some help?

Luke
  • 21
  • 5

2 Answers2

1

Your problem is the starting values of the parameters and the log likelihood

logl <--N/2* log(2*pi*sigma2)-1/(2*sigma2)* sum((beta0+beta1*x1-y)^2)

This breaks down if sigma2 becomes zero for multiple reasons: you take the log zero (not nice) and divide by zero (not at all nice). In other words you need to ensure that your sigma2 is positive.

One simple (and classical) way of doing exactly that is to exp() the parameter input, so replace

sigma2 <- theta[3]# x contains our covariates

with 

sigma2 <- exp(theta[3])# x contains our covariates

This allows the optim() function to use the full range of values but transforms the real line to the positive line so the likelihood makes sense. Just remember that the parameter estimate for sigma2 returned by the optim() function will be the logged value.

Update with full code and output

The full code (with the change shown above) is

lm.loglik3 <-function(y, x, theta){
  N <-length(y)# theta contains our parameters to be estimated
  beta0 <- theta[1]
  beta1 <- theta[2]
  sigma2 <- exp(theta[3]) # x contains our covariates
  x1 <- x
  logl <- -N/2* log(2*pi*sigma2)-1/(2*sigma2)* sum((beta0+beta1*x1-y)^2)
  return(-logl)
}

When initialized with

stval <-c(0, 0, 0)

X <- c(8, 2, 4, 1, 1, 13, 2, 9, 2, 1, 7, 2, 1, 2, 4, 8, 3, 2, 2, 2, 
       3, 1, 14, 2, 1, 5, 2, 2, 4, 5)

Y <- c(18, 17, 13, 13, 20, 11, 11, 9, 8, 14, 8, 7, 15, 13, 12, 17, 
       10, 16, 14, 10, 13, 10, 12, 16, 18, 9, 9, 15, 13, 18)

we get 

> res3
$par
[1] 13.4598331 -0.1292182  2.4549731

$value
[1] 79.3809

$counts
function gradient 
     178       NA 

$convergence
[1] 0

$message
NULL

$hessian
           [,1]       [,2]        [,3]
[1,] 2.57596526  9.8745335  0.00558974
[2,] 9.87453350 67.7478863  0.00621170
[3,] 0.00558974  0.0062117 14.98815391
ekstroem
  • 5,957
  • 3
  • 22
  • 48
  • Show us the data you are using so we can reproduce the error. What is x and y? – ekstroem Dec 25 '19 at 14:10
  • i've put the data – Luke Dec 26 '19 at 10:01
  • The code work exactly as desired with the data you provided. I've updated my answer and shown the full code and output – ekstroem Dec 26 '19 at 11:06
  • it works with the 30 samples. But i should do it on 6000 and it doesn't work giving me the same exact error. I could use only 30 samples maybe but than when I try to plot it it says to me: Error: Aesthetics must be either length 1 or the same as the data (6000): x, y – Luke Dec 26 '19 at 14:31
  • Now i've tried to put: Y<-df2$auth_val[1:6000] X<-df2$v64[1:6000]; finally doesn't give me the error with the MLE o_0, but now when I tried to plot it it says; Error in -ggplot(df2, aes(x = X, y = Y)) : the argument is not valid for the operator..... GGplotcode: ggplot(df2, aes(x = X, y = Y)) + geom_point() + geom_smooth(method = "glm", se = FALSE, size = 1, color = "purple") + xlab("Authoritarianism") + ylab("income") + ggtitle("Regression Plot") + theme(plot.title = element_text(hjust = 0.5)) – Luke Dec 26 '19 at 14:51
  • That is another question related to ggplot and problems surrounding that and has nothing to do with the MLE question you asked. You don't even have any ggplot code in your question that we could try to fix. If you have another question then ask it as a new question – ekstroem Dec 26 '19 at 20:16
  • Yes by the way i need to do it with 6000 samples, and it doesn't work giving me the same error, so the question remained – Luke Dec 27 '19 at 08:51
  • We cannot help you without you providing the data that gives the error. – ekstroem Dec 27 '19 at 17:30
0
lm.loglik2 <- function(y, x, theta){
   N <- length(y)

   # theta contains our parameters to be estimated
   beta0 <- theta[1]
   beta1 <- theta[2]
   sigma2 <- exp(theta[3])


   logl <- - N / 2 * log(2 * pi * sigma2) - 
     1 / (2 * sigma2) * sum((beta0 + beta1 * x - y)^2)
   return(-logl)
 }

 # Starting Values
 stval <- c(0, 0, 1)
> 
> # Optim
> res3 <- optim(stval, 
               lm.loglik2,
               x=X
               y = Y,                
              hessian = TRUE)
Error in optim(stval, lm.loglik2, y = Y, x = X, hessian = TRUE) : 
  the function cannot be evaluated for its initial parameters

sigma2 coefficient: 0.46 standard error:0.07

X<- df2$v64

dput(df2$v64[1:30]) 
c(8, 2, 4, 1, 1, 13, 2, 9, 2, 1, 7, 2, 1, 2, 4, 8, 3, 2, 2, 2, 
3, 1, 14, 2, 1, 5, 2, 2, 4, 5)

Y<- df2$auth_val

dput(df2$auth_val[1:30]) 
c(18, 17, 13, 13, 20, 11, 11, 9, 8, 14, 8, 7, 15, 13, 12, 17, 
10, 16, 14, 10, 13, 10, 12, 16, 18, 9, 9, 15, 13, 18)
Luke
  • 21
  • 5