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

Question

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))

It's easier to help you if you include a simple [reproducible example](https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) with sample input that can be used to test and verify possible solutions. — MrFlick, Nov 29 '22 at 18:38
You are passing values for `x`, `y` and `group` to your function. What are the free variables that you are trying to estimate? — MrFlick, Nov 29 '22 at 18:41
I am trying to calculate beta and tau and replicate lmer.mod = lmer(Reaction ~ 1 + (1|id), data = sleepstudy, REML = F) — user13064749, Nov 29 '22 at 19:16

score 1 · Answer 1 · answered Dec 02 '22 at 12:08

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)

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

1 Answers1