user13064749
user13064749

Reputation: 1

I am trying to compute max likelihood function an implement from scratch in R for linear mix models

It do not run but exist with this error

Error in fn(par, ...) : unused argument (par)

Use optim to minimalist maximum likelihood function

x=sleepstudy$Days
group=sleepstudy$id
y=sleepstudy$Reaction

mylmax=function(y,x,group){
  n=length(y)
  x=cbind(x,group)
  sigma2=1000
  tau2=1000
  beta=matrix(1:1,2,1)
  v=sigma2*diag(n)+tau2+diag(n)
  inv=(1/sigma2)*(diag(n)-tau2/(sigma2+n*tau2)*diag(n))
  reml=FALSE
  l=-0.5*(log(det(v))+t(y-x%*%beta)%*%inv%*%(y-x%*%beta))
  if(reml==TRUE){
    l=l-0.5*log(det(t(v)%*%inv%*%x))
  }
  
}

optim(par=c(0,0,0),fn=mylmax,y=y,x=x,group=group, method = "L-BFGS-B",hessian = TRUE,control = list(fnscale = -1))

Upvotes: 0

Views: 70

Answers (1)

Nir Graham
Nir Graham

Reputation: 5137

Much wisdom is to be found in consulting the R help files for ?optim. The simplest example provided there is of the form


fr <- function(x) {   ## Rosenbrock Banana function
  x1 <- x[1]
  x2 <- x[2]
  100 * (x2 - x1 * x1)^2 + (1 - x1)^2
}

optim(c(-1.2,1), fr)

where we see fr takes a single input (x) - which is packed with the two params

your function is not unary but is like so; which will therefore not work out of the box ; see :

frx <- function(x1,x2) {   ## Rosenbrock Banana function
  100 * (x2 - x1 * x1)^2 + (1 - x1)^2
}
optim(c(-1.2,1), frx)
# Error in fn(par, ...) : argument "x2" is missing, with no default

One can either rewrite ones function to be of the unary type; or keep our nice function; and write a wrapper function for it that just does the x packing part:

frxwrap <- function(x){
  frx(x1=x[1],
      x2=x[2])
}
optim(c(-1.2,1), frxwrap)

Upvotes: 1

Related Questions