[R-sig-ME] evaluating fixed effects variance-covariance matrix at specified values of theta

Jim Hughes jphughe@ @end|ng |rom uw@edu
Fri Jun 13 23:52:01 CEST 2025


Hello

I would like to compute the first derivative of (a simple function of) the fixed effects variance-covariance matrix from a g/lmer model with respect to the random effects parameters (the theta's), evaluated at \hat{theta}. [Note: the function of the variance-covariance matrix I am interested in will produce a scalar, ie t(lvec)%*%vcov(obj)%*%lvec where lvec is a p x 1 matrix of constants]. I was wondering if there is a way to modify the approach outlined in a previous post by Ben Bolker:

>> fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
>> fm1Fun <- update(fm1,devFunOnly=TRUE)
>> library(numDeriv)
>> fm1_thpar <- getME(fm1,"theta")
>> h <- hessian(fm1Fun,fm1_thpar)
>>
>>  and similarly for the gradient.

In this case fm1Fun is a function that computes the deviance (as a function of the theta parameters) so calls to hessian and grad produce the second and first derivatives of the deviance, respectively, with respect to the theta parameters. I would like a function
that essentially computes vcov(fm1) as a function of the theta parameters.

I did try the following:

  fm1Fun = function(newtheta,obj,lvec) {
    beta=fixef(obj)
    dat=model.frame(obj)
    fam=family(obj)
    form=formula(obj)
    newobj = suppressWarnings(glmer(formula=form,family=fam,data=dat,start=list(theta=newtheta,fixef=beta),
                                                                              control=glmerControl(optCtrl=list(maxfun=1))))
    as.numeric(t(lvec)%*%vcov(newobj)%*%lvec)
  }
#
  library(numDeriv)
  rslt = glmer(response.var ~ ftime + ftimeOnTx + (1 | fcluster), family=binomial, data=sdat)
lvec = matrix(c(0,0,0,0,1/3,1/3,1/3),7,1)
  theta = getME(rslt,"theta")
  g = grad(fm1Fun,theta,obj=rslt,lvec=lvec)

My thought was that the glmer call in fm1Fun would evaluate everything at the starting values I provided and then stop, but I don't know for sure if that is what happens when you specify maxfun=1?

Any thoughts or suggestions?

Jim
--------------
Jim Hughes



	[[alternative HTML version deleted]]



More information about the R-sig-mixed-models mailing list