[R-sig-ME] Questions about porting pedigreemm to the new lme4

Tue Apr 3 19:53:41 CEST 2012

Hi all:

I have created a patch for the latest lmer.R to accommodate calls by pedigreemm. I see that there is potentially a way to create a single "finalize" method instead of one for each lmer, glmer, and nlmer. I assume that nAGQ is 0 for lmer, correct? Or is it possible for glmer or nlmer to specify nAGQ=0?

Perhaps I should submit the patch into the patch database to see whether it is accepted or not. If it is accepted, then I can make the changes to pedigreemm accordingly.

Thank you,
Roby

On Apr 2, 2012, at 6:19 PM, Joehanes, Roby (NIH/NHLBI) [F] wrote:

> Hi all:
>
> Just a follow up: It seems that some changes need to be made in the new lme4 so that the pedigreemm can be modified to use it. They are:
> 1. To allow the return of all necessary internal variables required for invoking the optimizer (especially rho plus devfun in one call), without invoking the optimizer itself. Right now I did two calls, the first to get the devfun, the second to get the rest. However, I cannot obtain rho in all cases.
> #lmerc$doFit <- FALSE # call lmer without pedigree and with doFit = FALSE
> l... <- list(...)
> lmerc$devFunOnly <- TRUE
> devFun <- eval(lmerc, parent.frame())
> lmerc$devFunOnly <- FALSE
> lmf <- eval(lmerc, parent.frame())
>
> 2. Some partitioning of Zt that allows duplicate level names for multiple random factors. That is, if random factors 1 and 2 (let's call them RF1 and RF2) share the same level-string, my current hack will not work. Perhaps someone can propose a better way. My current hack is as follows:
> Zt <- getME(lmf, "Zt")
> Zt <- Zt[rownames(Zt) %in% levels(getME(lmf, "flist")[[tn]]), ]
> (The original was: Zt <- lmf$FL$trms[[tn]]$Zt)
>
> Since I cannot get rho, I cannot realize the lmer_finalize or glmer_finalize. Any help is appreciated. I am stuck right now since I have installed the lme4 SVN version, but the pedigreemm package still assumes the CRAN version.
>
> Thank you,
> Roby
>
>
> On Apr 2, 2012, at 4:24 PM, Joehanes, Roby (NIH/NHLBI) [F] wrote:
>
> Hi Ben and Andrzej:
>
> Thank you for the info. I am trying to get pedigreemm to work with the new lme4 and I am now stuck on how to translate the following line:
> lmf$FL$trms[[tn]]$Zt <- lmf$FL$trms[[tn]]$A <- relfac[[i]] %*% Zt
>
> Because the matrix A is defined in getME method as:
> PR$Lambdat %*% PR$Zt
>
> I am still trying to figure the code out. I think we can do away with the A matrix since it is apparently not essential to the computation. I think the crux of the pedigreemm call to the underlying lme4 is in the following short lines:
>
>   pnms <- names(pedigree)
>   stopifnot(all(pnms %in% names(getME(lmf, "flist"))))
>   asgn <- attr(getME(lmf, "flist"), "assign")
>   for (i in seq_along(pedigree)) {
>       tn <- which(match(pnms[i], names(getME(lmf, "flist"))) == asgn)
>       if (length(tn) > 1)
>           stop("a pedigree factor must be associated with only one r.e. term")
>       Zt <- getME(lmf, "Zt")
>       relfac[[i]] <- relfactor(pedigree[[i]], rownames(Zt))
>       lmf$FL$trms[[tn]]$Zt <- lmf$FL$trms[[tn]]$A <- relfac[[i]] %*% Zt
>   }
>   ans <- do.call(if (!is.null(lmf$glmFit)) lme4:::glmer_finalize else lme4:::lmer_finalize, lmf)
>   ans <- new("pedigreemm", relfac = relfac, ans)
>
> Any help is appreciated.
>
> Thanks,
> Roby
>
>
> On Apr 2, 2012, at 4:17 PM, Ben Bolker wrote:
>
> Andrzej T Galecki 1 <agalecki at ...> writes:
>
> Ben prepared getME(} function to extract various components of the model
> fit
> obtained using lme4.0 and new lme4.
>
> Try
> getME(lmf,"Zt")
> getME(lmf,"A")
> getME(lmf,"flist")
>
> Andrzej
>
> On Mon, 2 Apr 2012 14:27:20 -0400, "Joehanes, Roby (NIH/NHLBI) [F]"
> <roby.joehanes at ...> wrote:
> Hi:
>
> I am trying to port pedigreemm to use the new lme4. However, the new lme4
> is somewhat opaque to me and hence I need help. Specifically, how can I
> extract the following terms:
> lmf$FL$trms[[tn]]$Zt
> lmf$FL$trms[[tn]]$A
>
> The information lies within lmf <at> pp object, but it is too opaque to me to
> fish anything out.
>
> I understand that lmf$FL$fl is lmf <at> flist and that lmer_finalize calls
> should be replaced with optwrap and mkMerMod. Nevertheless, I think that
> the lmer or glmer could use some refactoring to ease the calls.
>
> I appreciate any help on these.
>
>
> A couple of questions:
>
> 1. how broadly are you using A?  Off-list, Doug Bates has commented to
> me that defining A is relatively straightforward for LMMs but not necessarily
> as simple (nor uniquely defined!) for [GN]LMMs -- for the time being, I
> am quite likely to add code that *disallows* getME(.,"A") except for LMMs --
> I can imagine that the bulk of pedigreemm use is for LMMs, but I suppose
> that people are using it for [GN]LMMs as welll ...
>
> 2. I can't quite tell from your e-mail -- would it be useful to have
> the stuff between the first optwrap() and the last mkMerMod abstracted
> into a single (exposed) function?  For [GN]LMMs it should probably
> also allow control of whether the preliminary "nAGQ=0" optimization
> is done or not ...
>
> (All comments are mine alone, not necessarily speaking for all lme4
> authors ...)
>
> Ben Bolker
>
> _______________________________________________
> R-sig-mixed-models at r-project.org<mailto:R-sig-mixed-models at r-project.org> mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: lmer-patch-for-pedigreemm.txt
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20120403/ea362fa5/attachment-0002.txt>