如何使用OPTIM in R修复纵向数据(和随机效果)的问题($convergence=10和$gradient=15)？

Question

我试图从纵向数据( y )中估计参数，其中随机效应被用作y均值的一部分。随机效应遵循双变量正态分布。在r中使用optim时，不同的方法会提供不同的结果。

library(JM)
library(MASS)
library(survival)
library(splines)
library(nlme)
library(mvtnorm)

数据：

data("aids")
aids<-na.omit(aids)
aids<-aids[order(aids$Time),]
aids=aids[aids$patient %in% c(1:11), c("patient", "CD4",
  "obstime","Time","death","drug")]
delta<-aids$death
tv<-aids$Time
Zmat<-cbind(rep(1,length(delta)),aids$obstime)
aids$trt=as.numeric(aids$drug=="ddI")
Xmat<-cbind(rep(1,length(delta)),aids$obstime, aids$trt*aids$obstime)
Wmat=cbind(aids$trt)
y=aids$CD4
n=length(delta)
k=length(unique(aids$patient))

已创建要在optim中使用的函数

p_y_hat=function(par,y,Xmat,Zmat){
  v.beta=par[1:3]
  v.sigma=par[4]
  mu=c(0,0)
  v.s=matrix(c(par[5],par[6],par[6],par[7]),2,2)
  b=mvrnorm(k,mu,v.s)
  m=drop(Xmat %*% v.beta + rowSums(Zmat * b[aids$patient, ]))
  return(dnorm(y,m,v.sigma))
}

m.LL1<-function(par,y,Xmat,Zmat) {
  LL1.pr<-p_y_hat(par,y,Xmat,Zmat)
  return( - sum(log(LL1.pr) ))
}

按默认optim方法

fit1<-optim(par=c(1,0.2,0.2,5,2,1,3),m.LL1,
            control=list(maxit=1000),y=y,Xmat=Xmat,Zmat=Zmat)

它提供了

> fit1

$par  
[1] 1.2277439 0.2643519 0.2702051 5.2276476 1.9062541 0.9728032 2.8033234    
$value    
[1] 102.8805    
$counts    
function gradient     
     853       NA     
$convergence
[1] 10
$message
NULL

L-BFGS-B方法

    fit1<-optim(par=c(1,0.2,0.2,5,2,1,3),m.LL1,
  method="L-BFGS-B",control=list(maxit=1000),y=y,Xmat=Xmat,Zmat=Zmat)

$par
[1] 1.07742647 0.47206512 0.04462127 5.09672419 2.03116102 1.15173075 3.06604291
$value
[1] 305.7507
$counts
function gradient 
      15       15 
$convergence
[1] 0
$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

我对正确的方法感到有点困惑。第一个给了$convergence=10 (这应该是错的)，第二个给了$convergence=0，但$gradient=15给了一些额外的消息。

从上述纵向和双变量正态分布问题中估计参数，哪种方法是合适的？

Answer 1

L-BFGS-B认为一切都很好(收敛code=0)；您在那里看到的"gradient=15“只表示梯度计算的次数。消息(“融合: REL_REDUCTION_OF_F ...")

使用默认的(奈尔德-米德)优化器，optim()报告说“奈尔德-米德单纯形的退化”(参见?optim)，即optim()用来搜索最优的参数空间区域已经收缩到较低的维度。

一如既往，你应该评估你的结果，以确保它们在你的问题的上下文中是有意义的……

注：我只是在回答这里提出的狭义问题。我实际上并没有检查你的代码，看看它是否合理。