Error in function Optim

Asked

Viewed 168 times

1

I am trying to maximize the likelihood of logistic distribution with a regression structure. The code is:

cov1 <- rep(1,115)
cov2 <- rnorm(115,0,1)
e    <- rlogis(115, 0,1)

yy <- 1*cov1 + 0*cov2 + 2*e
n  <- 115

logL <- function(par, x1, x2, y, n){
    b1 <- par[1]
    b2 <- par[2]
    l  <- -(sum(y) - n*(b1*x1+b2*x2) - 2* (sum(log(1+exp(y-(b1*x1+b2*x2))))))
    return(l)
}

optim(c(1,0), logL, method="BFGS", x1=cov1, x2=cov2, y=yy, n=115)$par

Returns the following error:

Error in optim(start, logL, method = "BFGS", x1 = cov1, x2 = cov2, y = yy,  : 
função alvo em optim retorna um objeto de comprimento 115 ao invés de 1

I can’t find what’s wrong.

r
  • 1

    Your logL function is poorly defined. It returns n values when, in reality, it should return only one. Check exactly which function to minimize and implement in R.

1 answer

3

In logistic regression the function to be maximized is as follows:

inserir a descrição da imagem aqui

If you were to implement it in the R, it would look like this (the Optim minimizes, so the - in front). 

logL <- function(par, x1, x2, y, n){
  b1 <- par[1]
  b2 <- par[2]
  l  <- -sum(y*(x1*b1 + x2*b2) - n*log(1 + exp(x1*b1 + x2*b2)))
  return(l)
}

However, your example seems to have something strange... Logistic regression serves for discrete data and not for continuous numbers. (Your variable yy is continuous). In addition, the n is not the sample size, but the n of the binomial you are estimating for each individual.

I would make an example as follows:

> set.seed(1)
> cov1 <- runif(100)
> cov2 <- runif(100)
> 
> yy <- rbinom(n = 100, prob = (cov1 + cov2)/2, size = 1)
> 
> glm(yy ~ 0 + cov1 + cov2, family = "binomial")

Call:  glm(formula = yy ~ 0 + cov1 + cov2, family = "binomial")

Coefficients:
  cov1    cov2  
0.2138  0.1638  

Degrees of Freedom: 100 Total (i.e. Null);  98 Residual
Null Deviance:      138.6 
Residual Deviance: 137.6    AIC: 141.6
> optim(par = c(1,1), fn = logL, x1 = cov1, x2 = cov2, y = yy, n = 1)
$par
[1] 0.2132318 0.1643093

$value
[1] 68.77538

$counts
function gradient 
      55       NA 

$convergence
[1] 0

$message
NULL

Note that in this case, the calculation using the function glm of R, is equal to the version using the optim.

Browser other questions tagged r

You are not signed in. Login or sign up in order to post.