我试图在R中从头开始计算线性混合模型的最大似然函数

Question

不运行，但存在此错误

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

使用optim来最小化最大似然函数

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

Answer 1

在R帮助文件中可以找到很多智慧。这里提供的最简单的例子是

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)

，我们看到fr接受单个输入(x) -它包含两个参数

你的函数不是一元的，但像这样;因此在盒子外是行不通的;看到:

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

可以将函数重写为一元类型;还是保持我们的好函数;然后为它写一个包装器函数，它只做x打包部分:

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