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

[Roby Joehanes]

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.