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

[Jim Hughes]

Thank you ... I like that you were about to use update to do this!

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


-----Original Message-----
From: R-sig-mixed-models <r-sig-mixed-models-bounces using r-project.org> On Behalf Of Ben Bolker
Sent: Friday, June 13, 2025 4:47 PM
To: r-sig-mixed-models using r-project.org
Subject: Re: [R-sig-ME] evaluating fixed effects variance-covariance matrix at specified values of theta

   This seems reasonable at first glance.

   A  few details:

*  you might want to shut off gradient/Hessian computations by specifying checkConv = FALSE
* whether maxfun or some other control parameter is appropriate depends
(sigh) on which optimizer is being used. Both Nelder-Mead and bobyqa (the default setting is c("bobyqa", "Nelder_Mead"), meaning that "bobyqa" is used for a first-stage optimization and Nelder-Mead for the second.  I might consider setting nAGQ0initStep = FALSE and optimizer="bobyqa" (skip the first optimization phase, use bobyqa for the second)
* a slightly more compact form of the function (the only caveat is that
update() can occasionally get messed up if you create data, objects, etc. in enough different disconnected environments, as it relies on re-evaluating the original function call ...)


ctrl <- glmerControl(optCtrl = list(maxfun = 1),
                      nAGQ0initStep = FALSE, checkConv = FALSE,
                      optimizer = "bobyqa") fm1Fun = function(newtheta,obj,lvec) {
    beta <- fixef(obj)
    newobj <- update(obj, control = ctrl,
        start=list(theta=newtheta,fixef=beta))
    as.numeric(t(lvec)%*%vcov(newobj)%*%lvec)
}

   It probably doesn't matter, but emulator::quad.tform offers a slightly more efficient way to compute the quadratic form (it's also a one-liner, so you can copy it if you don't want to depend on the package)

   cheers
    Ben Bolker





On 2025-06-13 5:52 p.m., Jim Hughes wrote:
> 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]]
> 
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list 
> https://urldefense.com/v3/__https://stat.ethz.ch/mailman/listinfo/r-si
> g-mixed-models__;!!K-Hz7m0Vt54!iuIQhaS-yBaHD-1S8rzX92WgHLHwpbaY8xtfjaa
> uCCWdlOF5tM09ZrM4T6CuMpnPFdt6Z7H1gs8A0k8$

--
Dr. Benjamin Bolker
Professor, Mathematics & Statistics and Biology, McMaster University Director, School of Computational Science and Engineering  > E-mail is sent at my convenience; I don't expect replies outside of working hours.

_______________________________________________
R-sig-mixed-models using r-project.org mailing list https://urldefense.com/v3/__https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models__;!!K-Hz7m0Vt54!iuIQhaS-yBaHD-1S8rzX92WgHLHwpbaY8xtfjaauCCWdlOF5tM09ZrM4T6CuMpnPFdt6Z7H1gs8A0k8$