[R-sig-ME] Patch for pedigreemm 0.3-3 plus questions

[Ben Bolker]

  This is great -- thank you! -- but it would be better to e-mail the
maintainer (try maintainer("pedigreemm")) directly -- they may or may
not actually read this mailing list ...

On 17-02-08 02:20 PM, Roby Joehanes wrote:
> Dear all:
> 
> It looks like lme4::lFormula or lme4::glFormula works. The time
> savings is about 40% now. Here is the code. If this code gets the
> blessings of the maintainers, I wonder if this code could be
> incorporated into the official pedigreemm.
> 
> Sincerely, Roby
> 
> # Modified pedigreemm source # Based on pedigreemm 0.3-3 # By Roby
> Joehanes
> 
> # What is modified: # 1. Allowing optimizer options to pass through #
> 2. Allow direct specification of pedigree matrix # 3. Shortcut to
> reconstruct necessary matrices instead of calling lmer/glmer
> 
> pedigreemm <- function(formula, data, family = NULL, REML = TRUE,
> pedigree = list(), control = list(), start = NULL, verbose = FALSE, 
> subset, weights, na.action, offset, contrasts = NULL, model = TRUE, x
> = TRUE, ...) { gaus <- FALSE if (is.null(family)) { gaus <- TRUE }
> else { ## copied from glm() if (is.character(family)) family <-
> get(family, mode = "function", envir = parent.frame()) if
> (is.function(family)) family <- family() if (!inherits(family,
> "family")) stop("unknown family type") gaus <- family$family ==
> "gaussian" && family$link == "identity" } mc <- match.call() lmerc <-
> mc                         # create a call to lmer lmerc[[1]] <- if
> (gaus) as.name("lmer") else as.name("glmer") lmerc$pedigree <- NULL 
> if (!gaus) lmerc$REML <- NULL  if (!length(pedigree))              #
> call [g]lmer instead return(eval.parent(lmerc))  
> stopifnot(is.list(pedigree),        # check the pedigree argument 
> length(names(pedigree)) == length(pedigree), all(sapply(pedigree,
> function(x) is(x, "pedigree") | is(x, "matrix")))) # RJ modification 
> 

> lmerc[[1]] <- if (gaus) quote(lme4::lFormula) else
> quote(lme4::glFormula) # RJ modification lmf <- eval(lmerc,
> parent.frame())   relfac <- pedigree          # copy the pedigree
> list for relfactor pnms <- names(pedigree) #pp <- lmf at pp # RJ
> modification #resp <- lmf at resp # RJ modification fl <-
> lmf$reTrms$flist # RJ modification stopifnot(all(pnms %in%
> names(fl))) asgn <- attr(fl, "assign") #Zt <- pp$Zt # RJ
> modification for (i in seq_along(pedigree)) { tn <-
> which(match(pnms[i], names(fl)) == asgn) if (length(tn) > 1) stop("a
> pedigree factor must be associated with only one r.e. term") ind <-
> (lmf$reTrms$Gp)[tn:(tn+1L)] # RJ modification rowsi <-
> (ind[1]+1L):ind[2] if (is(relfac[[i]], "pedigree")) relfac[[i]] <-
> relfactor(pedigree[[i]], rownames(lmf$reTrms$Zt)[rowsi]) # RJ
> modification lmf$reTrms$Zt[rowsi,] <- relfac[[i]] %*%
> lmf$reTrms$Zt[rowsi,] # RJ modification } #reTrms <-
> list(Zt=Zt,theta=lmf at theta,Lambdat=pp$Lambdat,Lind=pp$Lind,  # RJ
> modification starts #		lower=lmf at lower,flist=lmf at flist,cnms=lmf at cnms,
> Gp=lmf at Gp) dfl <- list(fr=lmf$fr, X=lmf$X, reTrms=lmf$reTrms,
> start=lmf$reTrms$theta) optimizer_opt <- ifelse(gaus,
> control$optimizer, control$optimizer[[2]]); if
> (is.null(optimizer_opt)) optimizer_opt <- "Nelder_Mead"; # RJ
> modification ends if (gaus) { dfl$REML = lmf$REML > 0L # RJ
> modification devfun <- do.call(mkLmerDevfun,dfl) opt <-
> optimizeLmer(devfun, optimizer=optimizer_opt,  # RJ modification
> starts restart_edge = control$restart_edge, boundary.tol =
> control$boundary.tol, control = control$optCtrl, 
> calc.derivs=control$calc.derivs, 
> use.last.params=control$use.last.params,...) # RJ modification ends }
> else { dfl$family <- family devfun <- do.call(mkGlmerDevfun,dfl) opt
> <- optimizeGlmer(devfun, optimizer=optimizer_opt, # RJ modification
> starts restart_edge=control$restart_edge, 
> boundary.tol=control$boundary.tol, control = control$optCtrl, 
> stage=2, calc.derivs=control$calc.derivs, 
> use.last.params=control$use.last.params,...) # RJ modification ends 
> } mm <- mkMerMod(environment(devfun), opt, lmf$reTrms, lmf$fr, mc)  #
> RJ modification cls <- if (gaus) "lmerpedigreemm" else
> "glmerpedigreemm" ans <- do.call(new, list(Class=cls, relfac=relfac, 
> frame=mm at frame, flist=mm at flist, cnms=mm at cnms, Gp=mm at Gp, 
> theta=mm at theta, beta=mm at beta,u=mm at u,lower=mm at lower, 
> devcomp=mm at devcomp, pp=mm at pp,resp=mm at resp,optinfo=mm at optinfo)) 
> ans at call <- evalq(mc) ans }
> 
> 
> 
> -----Original Message----- From: R-sig-mixed-models
> [mailto:r-sig-mixed-models-bounces at r-project.org] On Behalf Of Roby
> Joehanes Sent: Wednesday, February 8, 2017 1:52 PM To:
> r-sig-mixed-models at r-project.org Subject: [R-sig-ME] Patch for
> pedigreemm 0.3-3 plus questions
> 
> Dear all:
> 
> Attached please find the patch for pedigreemm 0.3-3, with the
> following changes: 1. Control options are now passed through to the
> optimizeLmer / optimizeGlmer, including optimizer, calc.derivs, etc.
> This saves about 15% of time. 2. Allowing direct matrix specification
> on pedigree. This will allow broader application in genetics, such as
> allowing cryptic relatedness to stand in as pedigree matrix.
> 
> A few questions from me. I think pedigreemm tries to call lmer /
> glmer once to get @pp (for Zt, primarily, but also X, Lambdat, and
> Lind), @resp (only for REML), @flist and @Gp, then followed up with
> optimizeLmer or optimizeGlmer. I think it is a big waste of time
> since most of these can be constructed offline (i.e., without
> actually computing the results), like in pedigreemm version 1.1 or
> so. Is there a way in lme4 package to do that instead of calling lmer
> / glmer? Is lme4::lFormula / lme4::glFormula the right way? Did I
> overlook something?
> 
> Thank you!
> 
> Sincerely, Roby
> 
> 
> 
> ----------------------------------------------------------------------
>
> 
CONFIDENTIAL NOTICE:
> 
> This electronic mail transmission contains confidential information
> including Protected Health Information (PHI) that is legally
> privileged. If you are not the intended recipient, or designee, you
> are hereby notified that any disclosure, copying, distribution or use
> of any and all attachments to this transmission is STRICTLY
> PROHIBITED. If you have received this transmission in error, please
> notify the sender immediately to arrange for return or destruction of
> these documents. _______________________________________________ 
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> 
> ----------------------------------------------------------------------
>
> 
CONFIDENTIAL NOTICE:
> 
> This electronic mail transmission contains confidential information 
> including Protected Health Information (PHI) that is legally
> privileged. If you are not the intended recipient, or designee, you
> are hereby notified that any disclosure, copying, distribution or use
> of any and all attachments to this transmission is STRICTLY
> PROHIBITED. If you have received this transmission in error, please
> notify the sender immediately to arrange for return or destruction of
> these documents.
> 
> _______________________________________________ 
> R-sig-mixed-models at r-project.org mailing list 
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>