1## NB: doc in ../man/*.Rd ***not*** auto generated 2## FIXME: need to document S3 methods better (can we pull from r-forge version?) 3 4##' Fit a linear mixed model (LMM) 5lmer <- function(formula, data=NULL, REML = TRUE, 6 control = lmerControl(), start = NULL 7 , verbose = 0L 8 , subset, weights, na.action, offset 9 , contrasts = NULL 10 , devFunOnly=FALSE 11 ) 12 ## , ...) 13{ 14 mc <- mcout <- match.call() 15 missCtrl <- missing(control) 16 ## see functions in modular.R for the body .. 17 if (!missCtrl && !inherits(control, "lmerControl")) { 18 if(!is.list(control)) stop("'control' is not a list; use lmerControl()") 19 ## back-compatibility kluge 20 warning("passing control as list is deprecated: please use lmerControl() instead", 21 immediate.=TRUE) 22 control <- do.call(lmerControl, control) 23 } 24 ## if (!is.null(list(...)[["family"]])) { 25 ## warning("calling lmer with 'family' is deprecated; please use glmer() instead") 26 ## mc[[1]] <- quote(lme4::glmer) 27 ## if(missCtrl) mc$control <- glmerControl() 28 ## return(eval(mc, parent.frame(1L))) 29 ## } 30 mc$control <- control ## update for back-compatibility kluge 31 32 ## https://github.com/lme4/lme4/issues/50 33 ## parse data and formula 34 mc[[1]] <- quote(lme4::lFormula) 35 lmod <- eval(mc, parent.frame(1L)) 36 mcout$formula <- lmod$formula 37 lmod$formula <- NULL 38 39 ## create deviance function for covariance parameters (theta) 40 devfun <- do.call(mkLmerDevfun, 41 c(lmod, 42 list(start=start, verbose=verbose, control=control))) 43 if (devFunOnly) return(devfun) 44 ## optimize deviance function over covariance parameters 45 if (identical(control$optimizer,"none")) 46 stop("deprecated use of optimizer=='none'; use NULL instead") 47 opt <- if (length(control$optimizer)==0) { 48 s <- getStart(start, environment(devfun)$pp) 49 list(par=s,fval=devfun(s), 50 conv=1000,message="no optimization") 51 } else { 52 optimizeLmer(devfun, optimizer = control$optimizer, 53 restart_edge = control$restart_edge, 54 boundary.tol = control$boundary.tol, 55 control = control$optCtrl, 56 verbose=verbose, 57 start=start, 58 calc.derivs=control$calc.derivs, 59 use.last.params=control$use.last.params) 60 } 61 cc <- checkConv(attr(opt,"derivs"), opt$par, 62 ctrl = control$checkConv, 63 lbound = environment(devfun)$lower) 64 mkMerMod(environment(devfun), opt, lmod$reTrms, fr = lmod$fr, 65 mc = mcout, lme4conv=cc) ## prepare output 66}## { lmer } 67 68 69##' Fit a generalized linear mixed model (GLMM) 70glmer <- function(formula, data=NULL 71 , family = gaussian 72 , control = glmerControl() 73 , start = NULL 74 , verbose = 0L 75 , nAGQ = 1L 76 , subset, weights, na.action, offset, contrasts = NULL 77 , mustart, etastart 78 , devFunOnly = FALSE) 79{ 80 if (!inherits(control, "glmerControl")) { 81 if(!is.list(control)) stop("'control' is not a list; use glmerControl()") 82 ## back-compatibility kluge 83 if (class(control)[1]=="lmerControl") { 84 warning("please use glmerControl() instead of lmerControl()", 85 immediate.=TRUE) 86 control <- 87 ## unpack sub-lists 88 c(control[!names(control) %in% c("checkConv","checkControl")], 89 control$checkControl,control$checkConv) 90 control["restart_edge"] <- NULL ## not implemented for glmer 91 } else { 92 msg <- "Use control=glmerControl(..) instead of passing a list" 93 if(length(cl <- class(control))) { 94 msg <- paste(msg, "of class", dQuote(cl[1])) 95 } 96 warning(msg, immediate.=TRUE) 97 } 98 99 control <- do.call(glmerControl, control) 100 } 101 mc <- mcout <- match.call() 102 103 ## family-checking code duplicated here and in glFormula (for now) since 104 ## we really need to redirect at this point; eventually deprecate formally 105 ## and clean up 106 if (is.character(family)) 107 family <- get(family, mode = "function", envir = parent.frame(2)) 108 if( is.function(family)) family <- family() 109 if (isTRUE(all.equal(family, gaussian()))) { 110 ## redirect to lmer (with warning) 111 warning("calling glmer() with family=gaussian (identity link) as a shortcut to lmer() is deprecated;", 112 " please call lmer() directly") 113 mc[[1]] <- quote(lme4::lmer) 114 mc["family"] <- NULL # to avoid an infinite loop 115 return(eval(mc, parent.frame())) 116 } 117 118 ## see https://github.com/lme4/lme4/issues/50 119 ## parse the formula and data 120 mc[[1]] <- quote(lme4::glFormula) 121 glmod <- eval(mc, parent.frame(1L)) 122 mcout$formula <- glmod$formula 123 glmod$formula <- NULL 124 125 ## create deviance function for covariance parameters (theta) 126 127 nAGQinit <- if(control$nAGQ0initStep) 0L else 1L 128 devfun <- do.call(mkGlmerDevfun, c(glmod, list(verbose = verbose, 129 control = control, 130 nAGQ = nAGQinit))) 131 if (nAGQ==0 && devFunOnly) return(devfun) 132 ## optimize deviance function over covariance parameters 133 134 ## FIXME: perhaps should be in glFormula instead?? 135 if (is.list(start)) { 136 start.bad <- setdiff(names(start),c("theta","fixef")) 137 if (length(start.bad)>0) { 138 stop(sprintf("bad name(s) for start vector (%s); should be %s and/or %s", 139 paste(start.bad,collapse=", "), 140 shQuote("theta"), 141 shQuote("fixef")),call.=FALSE) 142 } 143 if (!is.null(start$fixef) && nAGQ==0) 144 stop("should not specify both start$fixef and nAGQ==0") 145 } 146 147 ## FIX ME: allow calc.derivs, use.last.params etc. if nAGQ=0 148 if(control$nAGQ0initStep) { 149 opt <- optimizeGlmer(devfun, 150 optimizer = control$optimizer[[1]], 151 ## DON'T try fancy edge tricks unless nAGQ=0 explicitly set 152 restart_edge=if (nAGQ==0) control$restart_edge else FALSE, 153 boundary.tol=if (nAGQ==0) control$boundary.tol else 0, 154 control = control$optCtrl, 155 start=start, 156 nAGQ = 0, 157 verbose=verbose, 158 calc.derivs=FALSE) 159 } 160 161 if(nAGQ > 0L) { 162 163 164 ## update deviance function to include fixed effects as inputs 165 devfun <- updateGlmerDevfun(devfun, glmod$reTrms, nAGQ = nAGQ) 166 167 if (control$nAGQ0initStep) { 168 start <- updateStart(start,theta=opt$par) 169 } 170 ## if nAGQ0 was skipped 171 ## we don't actually need to do anything here, it seems -- 172 ## getStart gets called again in optimizeGlmer 173 174 if (devFunOnly) return(devfun) 175 ## reoptimize deviance function over covariance parameters and fixed effects 176 opt <- optimizeGlmer(devfun, 177 optimizer = control$optimizer[[2]], 178 restart_edge=control$restart_edge, 179 boundary.tol=control$boundary.tol, 180 control = control$optCtrl, 181 start=start, 182 nAGQ=nAGQ, 183 verbose = verbose, 184 stage=2, 185 calc.derivs=control$calc.derivs, 186 use.last.params=control$use.last.params) 187 } 188 cc <- if (!control$calc.derivs) NULL else { 189 if (verbose > 10) cat("checking convergence\n") 190 checkConv(attr(opt,"derivs"),opt$par, 191 ctrl = control$checkConv, 192 lbound=environment(devfun)$lower) 193 } 194 195 ## prepare output 196 mkMerMod(environment(devfun), opt, glmod$reTrms, fr = glmod$fr, 197 mc = mcout, lme4conv=cc) 198 199}## {glmer} 200 201##' Fit a nonlinear mixed-effects model 202nlmer <- function(formula, data=NULL, control = nlmerControl(), start = NULL, verbose = 0L, 203 nAGQ = 1L, subset, weights, na.action, offset, 204 contrasts = NULL, devFunOnly = FALSE) 205{ 206 207 vals <- nlformula(mc <- match.call()) 208 p <- ncol(X <- vals$X) 209 if ((rankX <- rankMatrix(X)) < p) 210 stop(gettextf("rank of X = %d < ncol(X) = %d", rankX, p)) 211 212 rho <- list2env(list(verbose=verbose, 213 tolPwrss=0.001, # this is reset to the tolPwrss argument's value later 214 resp=vals$resp, 215 lower=vals$reTrms$lower), 216 parent=parent.frame()) 217 rho$pp <- do.call(merPredD$new, 218 c(vals$reTrms[c("Zt","theta","Lambdat","Lind")], 219 list(X=X, n=length(vals$respMod$mu), Xwts=vals$respMod$sqrtXwt, 220 beta0=qr.coef(qr(X), unlist(lapply(vals$pnames, get, 221 envir = rho$resp$nlenv)))))) 222 rho$u0 <- rho$pp$u0 223 rho$beta0 <- rho$pp$beta0 224 ## deviance as a function of theta only : 225 devfun <- mkdevfun(rho, 0L, verbose=verbose, control=control) 226 if (devFunOnly && !nAGQ) return(devfun) 227 devfun(rho$pp$theta) # initial coarse evaluation to get u0 and beta0 228 rho$u0 <- rho$pp$u0 229 rho$beta0 <- rho$pp$beta0 230 rho$tolPwrss <- control$tolPwrss # Reset control parameter (the initial optimization is coarse) 231 232 ## set lower and upper bounds: if user-specified, select 233 ## only the ones corresponding to random effects 234 if (!is.null(lwr <- control$optCtrl$lower)) { 235 rho$lower <- lwr[seq_along(rho$lower)] 236 control$optCtrl$lower <- NULL 237 } 238 upper <- rep(Inf, length(rho$lower)) 239 if (!is.null(upr <- control$optCtrl$upper)) { 240 upper <- upr[seq_along(rho$lower)] 241 control$optCtrl$upper <- NULL 242 } 243 244 opt <- optwrap(control$optimizer[[1]], devfun, rho$pp$theta, 245 lower=rho$lower, 246 upper=upper, 247 control=control$optCtrl, 248 adj=FALSE) 249 250 rho$control <- attr(opt,"control") 251 252 if (nAGQ > 0L) { 253 ## set lower/upper to values already harvested from control$optCtrl$upper 254 rho$lower <- if(!is.null(lwr)) lwr else c(rho$lower, rep.int(-Inf, length(rho$beta0))) 255 upper <- if(!is.null(upr)) upr else c( upper, rep.int( Inf, length(rho$beta0))) 256 rho$u0 <- rho$pp$u0 257 rho$dpars <- seq_along(rho$pp$theta) 258 ## fixed-effect parameters 259 rho$beta0 <- pmin(upper[-rho$dpars], 260 pmax(rho$pp$beta0,rho$lower[-rho$dpars])) 261 if (nAGQ > 1L) { 262 if (length(vals$reTrms$flist) != 1L || length(vals$reTrms$cnms[[1]]) != 1L) 263 stop("nAGQ > 1 is only available for models with a single, scalar random-effects term") 264 rho$fac <- vals$reTrms$flist[[1]] 265 } 266 devfun <- mkdevfun(rho, nAGQ, verbose=verbose, control=control) 267 if (devFunOnly) return(devfun) 268 269 opt <- optwrap(control$optimizer[[2]], devfun, 270 par = c(rho$pp$theta, rho$beta0), 271 lower = rho$lower, 272 upper = upper, 273 control = control$optCtrl, 274 adj = TRUE, verbose=verbose) 275 276 } 277 mkMerMod(environment(devfun), opt, vals$reTrms, fr = vals$frame, mc = mc) 278}## {nlmer} 279 280## R 3.1.0 devel [2013-08-05]: This does not help yet 281if(getRversion() >= "3.1.0") utils::suppressForeignCheck("nlmerAGQ") 282if(getRversion() < "3.1.0") dontCheck <- identity 283 284## *not* exported (had help page till early 2018) 285## -> issue #92: -> also look at devfun2() in ./profile.R (which returns class!) 286##' Create a deviance evaluation function from a predictor and a response module 287##' @param rho an `environment` already containing `verbose` and tolPwrss 288##' @param nAGQ for glmer/nlmer: #{AGQ steps}; 0 <==> Laplace 289##' @param maxit maximal number of PIRLS iterations 290##' @param verbose integer specifying if outputs should be produced 291##' @param control a list as from lmerControl() etc 292mkdevfun <- function(rho, nAGQ=1L, maxit = if(extends(rho.cld, "nlsResp")) 300L else 100L, 293 verbose=0, control=list()) { 294 ## FIXME: should nAGQ be automatically embedded in rho? 295 stopifnot(is.environment(rho), ## class definition, compute and save : 296 extends(rho.cld <- getClass(class(rho$resp)), "lmResp")) 297 298 ## silence R CMD check warnings *locally* in this function 299 ## (clearly preferred to using globalVariables() !] 300 fac <- pp <- resp <- lp0 <- compDev <- dpars <- baseOffset <- tolPwrss <- 301 pwrssUpdate <- ## <-- even though it's a function below 302 GQmat <- nlmerAGQ <- NULL 303 304 ## The deviance function (to be returned, with 'rho' as its environment): 305 ff <- 306 if (extends(rho.cld, "lmerResp")) { 307 rho$lmer_Deviance <- lmer_Deviance 308 function(theta) .Call(lmer_Deviance, pp$ptr(), resp$ptr(), as.double(theta)) 309 } else if (extends(rho.cld, "glmResp")) { 310 ## control values will override rho values *if present* 311 if (!is.null(tp <- control$tolPwrss)) rho$tolPwrss <- tp 312 if (!is.null(cd <- control$ compDev)) rho$compDev <- cd 313 if (nAGQ == 0L) 314 function(theta) { 315 resp$updateMu(lp0) 316 pp$setTheta(theta) 317 p <- pwrssUpdate(pp, resp, tol=tolPwrss, GQmat=GHrule(0L), 318 compDev=compDev, maxit=maxit, verbose=verbose) 319 resp$updateWts() 320 p 321 } 322 else ## nAGQ > 0 323 function(pars) { 324 ## pp$setDelu(rep(0, length(pp$delu))) 325 resp$setOffset(baseOffset) 326 resp$updateMu(lp0) 327 pp$setTheta(as.double(pars[dpars])) # theta is first part of pars 328 spars <- as.numeric(pars[-dpars]) 329 offset <- if (length(spars)==0) baseOffset else baseOffset + pp$X %*% spars 330 resp$setOffset(offset) 331 p <- pwrssUpdate(pp, resp, tol=tolPwrss, GQmat=GQmat, 332 compDev=compDev, grpFac=fac, maxit=maxit, verbose=verbose) 333 resp$updateWts() 334 p 335 } 336 } else if (extends(rho.cld, "nlsResp")) { 337 if (nAGQ <= 1L) { 338 rho$nlmerLaplace <- nlmerLaplace 339 rho$tolPwrss <- control$tolPwrss 340 rho$maxit <- maxit 341 switch(nAGQ + 1L, 342 function(theta) 343 .Call(nlmerLaplace, pp$ptr(), resp$ptr(), as.double(theta), 344 as.double(u0), beta0, verbose, FALSE, tolPwrss, maxit), 345 function(pars) 346 .Call(nlmerLaplace, pp$ptr(), resp$ptr(), pars[dpars], 347 u0, pars[-dpars], verbose, TRUE, tolPwrss, maxit)) 348 } else { 349 stop("nAGQ > 1 not yet implemented for nlmer models") 350 rho$nlmerAGQ <- nlmerAGQ 351 rho$GQmat <- GHrule(nAGQ) 352 ## function(pars) { 353 ## .Call(nlmerAGQ, ## <- dontCheck(nlmerAGQ) should work according to docs but does not 354 ## pp$ptr(), resp$ptr(), fac, GQmat, pars[dpars], 355 ## u0, pars[-dpars], tolPwrss) 356 ##} 357 } 358 } 359 else stop("code not yet written") 360 environment(ff) <- rho 361 ff 362} 363 364## Determine a step factor that will reduce the pwrss 365## 366## The penalized, weighted residual sum of squares (pwrss) is the sum 367## of the weighted residual sum of squares from the resp module and 368## the squared length of u from the predictor module. The predictor module 369## contains a base value and an increment for the coefficients. 370## @title Determine a step factor 371## @param pp predictor module 372## @param resp response module 373## @param verbose logical value determining verbose output 374## @return NULL if successful 375## @note Typically all this is done in the C++ code. 376## The R code is for debugging and comparisons of 377## results. 378## stepFac <- function(pp, resp, verbose, maxSteps = 10) { 379## stopifnot(is.numeric(maxSteps), maxSteps >= 2) 380## pwrss0 <- resp$wrss() + pp$sqrL(0) 381## for (fac in 2^(-(0:maxSteps))) { 382## wrss <- resp$updateMu(pp$linPred(fac)) 383## pwrss1 <- wrss + pp$sqrL(fac) 384## if (verbose > 3L) 385## cat(sprintf("pwrss0=%10g, diff=%10g, fac=%6.4f\n", 386## pwrss0, pwrss0 - pwrss1, fac)) 387## if (pwrss1 <= pwrss0) { 388## pp$installPars(fac) 389## return(NULL) 390## } 391## } 392## stop("step factor reduced below ",signif(2^(-maxSteps),2)," without reducing pwrss") 393## } 394 395RglmerWrkIter <- function(pp, resp, uOnly=FALSE) { 396 pp$updateXwts(resp$sqrtWrkWt()) 397 pp$updateDecomp() 398 pp$updateRes(resp$wtWrkResp()) 399 if (uOnly) pp$solveU() else pp$solve() 400 resp$updateMu(pp$linPred(1)) # full increment 401 resp$resDev() + pp$sqrL(1) 402} 403 404##' @param pp pred module 405##' @param resp resp module 406##' @param tol numeric tolerance 407##' @param GQmat matrix of Gauss-Hermite quad info 408##' @param compDev compute in C++ (as opposed to doing as much as possible in R) 409##' @param grpFac grouping factor (normally found in environment ..) 410##' @param verbose verbosity, of course 411glmerPwrssUpdate <- function(pp, resp, tol, GQmat, compDev=TRUE, grpFac=NULL, maxit = 70L, verbose=0) { 412 nAGQ <- nrow(GQmat) 413 if (compDev) { 414 if (nAGQ < 2L) 415 return(.Call(glmerLaplace, pp$ptr(), resp$ptr(), 416 nAGQ, tol, as.integer(maxit), 417 verbose)) 418 return(.Call(glmerAGQ, pp$ptr(), resp$ptr(), 419 tol, as.integer(maxit), 420 GQmat, grpFac, verbose)) 421 } 422 oldpdev <- .Machine$double.xmax 423 uOnly <- nAGQ == 0L 424 i <- 0 425 repeat { 426 ## oldu <- pp$delu 427 ## olddelb <- pp$delb 428 pdev <- RglmerWrkIter(pp, resp, uOnly=uOnly) 429 if (verbose > 2) cat(i,": ",pdev,"\n",sep="") 430 ## check convergence first so small increases don't trigger errors 431 if (is.na(pdev)) stop("encountered NA in PWRSS update") 432 if (abs((oldpdev - pdev) / pdev) < tol) 433 break 434 ## if (pdev > oldpdev) { 435 ## ## try step-halving 436 ## ## browser() 437 ## k <- 0 438 ## while (k < 10 && pdev > oldpdev) { 439 ## pp$setDelu((oldu + pp$delu)/2.) 440 ## if (!uOnly) pp$setDelb((olddelb + pp$delb)/2.) 441 ## pdev <- RglmerWrkIter(pp, resp, uOnly=uOnly) 442 ## k <- k+1 443 ## } 444 ## } 445 if (pdev > oldpdev) stop("PIRLS update failed") 446 oldpdev <- pdev 447 i <- i+1 448 } 449 resp$Laplace(pp$ldL2(), 0., pp$sqrL(1)) ## FIXME: should 0. be pp$ldRX2 ? 450} 451 452## create a deviance evaluation function that uses the sigma parameters 453## df2 <- function(dd) { 454## stopifnot(is.function(dd), 455## length(formals(dd)) == 1L, 456## is((rem <- (rho <- environment(dd))$rem), "Rcpp_reModule"), 457## is((fem <- rho$fem), "Rcpp_deFeMod"), 458## is((resp <- rho$resp), "Rcpp_lmerResp"), 459## all((lower <- rem$lower) == 0)) 460## Lind <- rem$Lind 461## n <- length(resp$y) 462## function(pars) { 463## sigma <- pars[1] 464## sigsq <- sigma * sigma 465## sigmas <- pars[-1] 466## theta <- sigmas/sigma 467## rem$theta <- theta 468## resp$updateMu(numeric(n)) 469## solveBetaU(rem, fem, resp$sqrtXwt, resp$wtres) 470## resp$updateMu(rem$linPred1(1) + fem$linPred1(1)) 471## n * log(2*pi*sigsq) + (resp$wrss + rem$sqrLenU)/sigsq + rem$ldL2 472## } 473## } 474 475## bootMer() ---> now in ./bootMer.R 476 477 478## Methods for the merMod class 479 480## Anova for merMod objects 481## 482## @title anova() for merMod objects 483## @param a merMod object 484## @param ... further such objects 485## @param refit should objects be refitted with ML (if applicable) 486## @return an "anova" data frame; the traditional (S3) result of anova() 487anovaLmer <- function(object, ..., refit = TRUE, model.names=NULL) { 488 mCall <- match.call(expand.dots = TRUE) 489 dots <- list(...) 490 .sapply <- function(L, FUN, ...) unlist(lapply(L, FUN, ...)) 491 modp <- (as.logical(vapply(dots, is, NA, "merMod")) | 492 as.logical(vapply(dots, is, NA, "lm"))) 493 if (any(modp)) { ## multiple models - form table 494 ## opts <- dots[!modp] 495 mods <- c(list(object), dots[modp]) 496 nobs.vec <- vapply(mods, nobs, 1L) 497 if (var(nobs.vec) > 0) 498 stop("models were not all fitted to the same size of dataset") 499 500 ## model names 501 if (is.null(mNms <- model.names)) 502 mNms <- vapply(as.list(mCall)[c(FALSE, TRUE, modp)], deparse1, "") 503 504 ## HACK to try to identify model names in situations such as 505 ## 'do.call(anova,list(model1,model2))' where the model names 506 ## are lost in the call stack ... this doesn't quite work but might 507 ## be useful for future attempts? 508 ## maxdepth <- -2 509 ## depth <- -1 510 ## while (depth >= maxdepth & 511 ## all(grepl("S4 object of class structure",mNms))) { 512 ## xCall <- match.call(call=sys.call(depth)) 513 ## mNms <- .sapply(as.list(xCall)[c(FALSE, TRUE, modp)], deparse) 514 ## depth <- depth-1 515 ## } 516 ## if (depth < maxdepth) { 517 if (any(substr(mNms, 1,4) == "new(") || 518 any(duplicated(mNms)) || ## <- only if S4 objects are *not* properly deparsed 519 max(nchar(mNms)) > 200) { 520 warning("failed to find model names, assigning generic names") 521 mNms <- paste0("MODEL",seq_along(mNms)) 522 } 523 if (length(mNms) != length(mods)) 524 stop("model names vector and model list have different lengths") 525 names(mods) <- sub("@env$", '', mNms) # <- hack 526 models.reml <- vapply(mods, function(x) is(x,"merMod") && isREML(x), NA) 527 models.GHQ <- vapply(mods, function(x) is(x,"glmerMod") && getME(x,"devcomp")$dims["nAGQ"]>1 , NA) 528 if (any(models.GHQ) && any(vapply(mods, function(x) is(x,"glm"), NA))) 529 stop("GLMMs with nAGQ>1 have log-likelihoods incommensurate with glm() objects") 530 if (refit) { 531 ## message only if at least one models is REML: 532 if (any(models.reml)) message("refitting model(s) with ML (instead of REML)") 533 mods[models.reml] <- lapply(mods[models.reml], refitML) 534 } else { ## check that models are consistent (all REML or all ML) 535 if(any(models.reml) && any(!models.reml)) 536 warning("some models fit with REML = TRUE, some not") 537 } 538 ## devs <- sapply(mods, deviance) 539 llks <- lapply(mods, logLik) 540 ## Order models by increasing degrees of freedom: 541 ii <- order(npar <- vapply(llks, attr, FUN.VALUE=numeric(1), "df")) 542 mods <- mods[ii] 543 llks <- llks[ii] 544 npar <- npar [ii] 545 calls <- lapply(mods, getCall) 546 data <- lapply(calls, `[[`, "data") 547 if(!all(vapply(data, identical, NA, data[[1]]))) 548 stop("all models must be fit to the same data object") 549 header <- paste("Data:", abbrDeparse(data[[1]])) 550 subset <- lapply(calls, `[[`, "subset") 551 if(!all(vapply(subset, identical, NA, subset[[1]]))) 552 stop("all models must use the same subset") 553 if (!is.null(subset[[1]])) 554 header <- c(header, paste("Subset:", abbrDeparse(subset[[1]]))) 555 llk <- unlist(llks) 556 chisq <- 2 * pmax(0, c(NA, diff(llk))) 557 dfChisq <- c(NA, diff(npar)) 558 val <- data.frame(npar = npar, 559 ## afraid to swap in vapply here; wondering 560 ## why .sapply was needed in the first place ... 561 AIC = .sapply(llks, AIC), # FIXME? vapply() 562 BIC = .sapply(llks, BIC), # " " 563 logLik = llk, 564 deviance = -2*llk, 565 Chisq = chisq, 566 Df = dfChisq, 567 "Pr(>Chisq)" = ifelse(dfChisq==0,NA,pchisq(chisq, dfChisq, lower.tail = FALSE)), 568 row.names = names(mods), check.names = FALSE) 569 class(val) <- c("anova", class(val)) 570 forms <- lapply(lapply(calls, `[[`, "formula"), deparse1) 571 structure(val, 572 heading = c(header, "Models:", 573 paste(rep.int(names(mods), lengths(forms)), 574 unlist(forms), sep = ": "))) 575 } 576 else { ## ------ single model --------------------- 577 if (length(dots)>0) { 578 warnmsg <- "additional arguments ignored" 579 nd <- names(dots) 580 nd <- nd[nzchar(nd)] 581 if (length(nd)>0) { 582 warnmsg <- paste0(warnmsg,": ", 583 paste(sQuote(nd),collapse=", ")) 584 } 585 warning(warnmsg) 586 } 587 dc <- getME(object, "devcomp") 588 X <- getME(object, "X") 589 stopifnot(length(asgn <- attr(X, "assign")) == dc$dims[["p"]]) 590 ss <- as.vector(object@pp$RX() %*% object@beta)^2 591 names(ss) <- colnames(X) 592 terms <- terms(object) 593 nmeffects <- attr(terms, "term.labels")[unique(asgn)] 594 if ("(Intercept)" %in% names(ss)) 595 nmeffects <- c("(Intercept)", nmeffects) 596 ss <- unlist(lapply(split(ss, asgn), sum)) 597 stopifnot(length(ss) == length(nmeffects)) 598 df <- lengths(split(asgn, asgn)) 599 ## dfr <- unlist(lapply(split(dfr, asgn), function(x) x[1])) 600 ms <- ss/df 601 f <- ms/(sigma(object)^2) 602 ## No longer provide p-values, but still the F statistic (may not be F distributed): 603 ## 604 ## P <- pf(f, df, dfr, lower.tail = FALSE) 605 ## table <- data.frame(df, ss, ms, dfr, f, P) 606 table <- data.frame(df, ss, ms, f) 607 dimnames(table) <- 608 list(nmeffects, 609 ## c("npar", "Sum Sq", "Mean Sq", "Denom", "F value", "Pr(>F)")) 610 c("npar", "Sum Sq", "Mean Sq", "F value")) 611 if ("(Intercept)" %in% nmeffects) 612 table <- table[-match("(Intercept)", nmeffects), ] 613 structure(table, heading = "Analysis of Variance Table", 614 class = c("anova", "data.frame")) 615 } 616}## {anovaLmer} 617 618##' @importFrom stats anova 619##' @S3method anova merMod 620anova.merMod <- anovaLmer 621 622##' @S3method as.function merMod 623as.function.merMod <- function(x, ...) { 624 rho <- list2env(list(resp = x@resp$copy(), 625 pp = x@pp$copy(), 626 beta0 = x@beta, 627 u0 = x@u), 628 parent=as.environment("package:lme4")) 629 ## FIXME: extract verbose [, maxit] and control 630 mkdevfun(rho, getME(x, "devcomp")$dims[["nAGQ"]], ...) 631} 632 633## coef() method for all kinds of "mer", "*merMod", ... objects 634## ------ should work with fixef() + ranef() alone 635coefMer <- function(object, ...) 636{ 637 if(...length()) 638 warning('arguments named ', paste(sQuote(...names()), collapse = ", "), 639 ' ignored') 640 fef <- data.frame(rbind(fixef(object)), check.names = FALSE) 641 ref <- ranef(object, condVar = FALSE) 642 ## check for variables in RE but missing from FE, fill in zeros in FE accordingly 643 refnames <- unlist(lapply(ref,colnames)) 644 nmiss <- length(missnames <- setdiff(refnames,names(fef))) 645 if (nmiss > 0) { 646 fillvars <- setNames(data.frame(rbind(rep(0,nmiss))),missnames) 647 fef <- cbind(fillvars,fef) 648 } 649 val <- lapply(ref, function(x) 650 fef[rep.int(1L, nrow(x)),,drop = FALSE]) 651 for (i in seq_along(val)) { 652 refi <- ref[[i]] 653 row.names(val[[i]]) <- row.names(refi) 654 nmsi <- colnames(refi) 655 if (!all(nmsi %in% names(fef))) 656 stop("unable to align random and fixed effects") 657 for (nm in nmsi) val[[i]][[nm]] <- val[[i]][[nm]] + refi[,nm] 658 } 659 class(val) <- "coef.mer" 660 val 661} ## {coefMer} 662 663##' @importFrom stats coef 664##' @S3method coef merMod 665coef.merMod <- coefMer 666 667## FIXME: should these values (i.e. ML criterion for REML models 668## and vice versa) be computed and stored in the object in the first place? 669##' @importFrom stats deviance 670##' @S3method deviance merMod 671deviance.merMod <- function(object, REML = NULL, ...) { 672 ## type = c("conditional", "unconditional", "penalized"), 673 ## relative = TRUE, ...) { 674 if (isGLMM(object)) { 675 return(sum(residuals(object,type="deviance")^2)) 676 ## ------------------------------------------------------------ 677 ## proposed change to deviance function for GLMMs 678 ## ------------------------------------------------------------ 679 ## @param type Type of deviance (can be unconditional, 680 ## penalized, conditional) 681 ## @param relative Should deviance be shifted relative to a 682 ## saturated model? (only available with type == penalized or 683 ## conditional) 684 ## ------------------------------------------------------------ 685 ## ans <- switch(type[1], 686 ## unconditional = { 687 ## if (relative) { 688 ## stop("unconditional and relative deviance is undefined") 689 ## } 690 ## c(-2 * logLik(object)) 691 ## }, 692 ## penalized = { 693 ## sqrL <- object@pp$sqrL(1) 694 ## if (relative) { 695 ## object@resp$resDev() + sqrL 696 ## } else { 697 ## useSc <- unname(getME(gm1, "devcomp")$dims["useSc"]) 698 ## qLog2Pi <- unname(getME(object, "q")) * log(2 * pi) 699 ## object@resp$aic() - (2 * useSc) + sqrL + qLog2Pi 700 ## } 701 ## }, 702 ## conditional = { 703 ## if (relative) { 704 ## object@resp$resDev() 705 ## } else { 706 ## useSc <- unname(getME(gm1, "devcomp")$dims["useSc"]) 707 ## object@resp$aic() - (2 * useSc) 708 ## } 709 ## }) 710 ## return(ans) 711 } 712 if (isREML(object) && is.null(REML)) { 713 warning("deviance() is deprecated for REML fits; use REMLcrit for the REML criterion or deviance(.,REML=FALSE) for deviance calculated at the REML fit") 714 return(devCrit(object, REML=TRUE)) 715 } 716 devCrit(object, REML=FALSE) 717} 718 719REMLcrit <- function(object) { 720 devCrit(object, REML=TRUE) 721} 722 723## original deviance.merMod -- now wrapped by REMLcrit 724## REML=NULL: 725## if REML fit return REML criterion 726## if ML fit, return deviance 727## REML=TRUE: 728## if not LMM, stop. 729## if ML fit, compute and return REML criterion 730## if REML fit, return REML criterion 731## REML=FALSE: 732## if ML fit, return deviance 733## if REML fit, compute and return deviance 734devCrit <- function(object, REML = NULL) { 735 ## cf. (1) lmerResp::Laplace in respModule.cpp 736 ## (2) section 5.6 of lMMwR, listing lines 34-42 737 if (isTRUE(REML) && !isLMM(object)) 738 stop("can't compute REML deviance for a non-LMM") 739 cmp <- object@devcomp$cmp 740 if (is.null(REML) || is.na(REML[1])) 741 REML <- isREML(object) 742 if (REML) { 743 if (isREML(object)) { 744 cmp[["REML"]] 745 } else { 746 ## adjust ML results to REML 747 lnum <- log(2*pi*cmp[["pwrss"]]) 748 n <- object@devcomp$dims[["n"]] 749 nmp <- n - length(object@beta) 750 ldW <- sum(log(weights(object, method = "prior"))) 751 - ldW + cmp[["ldL2"]] + cmp[["ldRX2"]] + nmp*(1 + lnum - log(nmp)) 752 } 753 } else { 754 if (!isREML(object)) { 755 cmp[["dev"]] 756 } else { 757 ## adjust REML results to ML 758 n <- object@devcomp$dims[["n"]] 759 lnum <- log(2*pi*cmp[["pwrss"]]) 760 ldW <- sum(log(weights(object, method = "prior"))) 761 - ldW + cmp[["ldL2"]] + n*(1 + lnum - log(n)) 762 } 763 } 764} 765 766## copied from stats:::safe_pchisq 767safe_pchisq <- function (q, df, ...) { 768 df[df <= 0] <- NA 769 pchisq(q = q, df = df, ...) 770} 771 772##' @importFrom stats drop1 773##' @S3method drop1 merMod 774drop1.merMod <- function(object, scope, scale = 0, test = c("none", "Chisq", "user"), 775 k = 2, trace = FALSE, 776 sumFun=NULL, ...) { 777 evalhack <- "formulaenv" 778 test <- match.arg(test) 779 if ((test=="user" && is.null(sumFun)) || 780 ((test!="user" && !is.null(sumFun)))) 781 stop(sQuote("sumFun"),' must be specified if (and only if) test=="user"') 782 tl <- attr(terms(object), "term.labels") 783 if(missing(scope)) scope <- drop.scope(object) 784 else { 785 if(!is.character(scope)) { 786 scope <- attr(terms(getFixedFormula(update.formula(object, scope))), 787 "term.labels") 788 } 789 if(!all(match(scope, tl, 0L) > 0L)) 790 stop("scope is not a subset of term labels") 791 } 792 ns <- length(scope) 793 if (is.null(sumFun)) { 794 sumFun <- function(x,scale,k,...) 795 setNames(extractAIC(x,scale,k,...),c("df","AIC")) 796 } 797 ss <- sumFun(object, scale=scale, k=k, ...) 798 ans <- matrix(nrow = ns + 1L, ncol = length(ss), 799 dimnames = list(c("<none>", scope), names(ss))) 800 ans[1, ] <- ss 801 n0 <- nobs(object, use.fallback = TRUE) 802 env <- environment(formula(object)) # perhaps here is where trouble begins?? 803 for(i in seq_along(scope)) { ## was seq(ns), failed on empty scope 804 tt <- scope[i] 805 if(trace > 1) { 806 cat("trying -", tt, "\n", sep='') 807 flush.console() 808 } 809 ## FIXME: make this more robust, somehow? 810 ## three choices explored so far: 811 ## (1) evaluate nfit in parent frame: tests in inst/tests/test-formulaEval.R 812 ## will fail on lapply(m_data_List,drop1) 813 ## (formula environment contains r,x,y,z but not d) 814 ## (2) evaluate nfit in frame of formula: tests will fail when data specified and formula is character 815 ## (3) update with data=NULL: fails when ... 816 ## 817 if (evalhack %in% c("parent","formulaenv")) { 818 nfit <- update(object, 819 as.formula(paste("~ . -", tt)), 820 evaluate = FALSE) 821 ## nfit <- eval(nfit, envir = env) # was eval.parent(nfit) 822 if (evalhack=="parent") { 823 nfit <- eval.parent(nfit) 824 } else if (evalhack=="formulaenv") { 825 nfit <- eval(nfit,envir=env) 826 } 827 } else { 828 nfit <- update(object, 829 as.formula(paste("~ . -", tt)),data=NULL, 830 evaluate = FALSE) 831 nfit <- eval(nfit,envir=env) 832 } 833 if (test=="user") { 834 ans[i+1, ] <- sumFun(object, nfit, scale=scale, k=k, ...) 835 } else { 836 ans[i+1, ] <- sumFun(nfit, scale, k = k, ...) 837 } 838 nnew <- nobs(nfit, use.fallback = TRUE) 839 if(all(is.finite(c(n0, nnew))) && nnew != n0) 840 stop("number of rows in use has changed: remove missing values?") 841 } 842 if (test=="user") { 843 aod <- as.data.frame(ans) 844 } else { 845 dfs <- ans[1L, 1L] - ans[, 1L] 846 dfs[1L] <- NA 847 aod <- data.frame(npar = dfs, AIC = ans[,2]) 848 if(test == "Chisq") { 849 ## reconstruct deviance from AIC (ugh) 850 dev <- ans[, 2L] - k*ans[, 1L] 851 dev <- dev - dev[1L] ; dev[1L] <- NA 852 nas <- !is.na(dev) 853 P <- dev 854 P[nas] <- safe_pchisq(dev[nas], dfs[nas], lower.tail = FALSE) 855 aod[, c("LRT", "Pr(Chi)")] <- list(dev, P) 856 } else if (test == "F") { 857 ## FIXME: allow this if denominator df are specified externally? 858 stop("F test STUB -- unfinished maybe forever") 859 dev <- ans[, 2L] - k*ans[, 1L] 860 dev <- dev - dev[1L] ; dev[1L] <- NA 861 nas <- !is.na(dev) 862 P <- dev 863 P[nas] <- safe_pchisq(dev[nas], dfs[nas], lower.tail = FALSE) 864 aod[, c("LRT", "Pr(F)")] <- list(dev, P) 865 } 866 } 867 head <- c("Single term deletions", "\nModel:", deparse(formula(object)), 868 if(scale > 0) paste("\nscale: ", format(scale), "\n")) 869 if (!is.null(method <- attr(ss,"method"))) { 870 head <- c(head,"Method: ",method,"\n") 871 } 872 structure(aod, heading = head, class = c("anova", "data.frame")) 873} 874 875##' @importFrom stats extractAIC 876##' @S3method extractAIC merMod 877extractAIC.merMod <- function(fit, scale = 0, k = 2, ...) { 878 L <- logLik(refitML(fit)) 879 edf <- attr(L,"df") 880 c(edf,-2*L + k*edf) 881} 882 883##' @importFrom stats family 884##' @S3method family merMod 885family.merMod <- function(object, ...) family(object@resp, ...) 886 887##' @S3method family glmResp 888family.glmResp <- function(object, ...) { 889 # regenerate initialize 890 # expression if necessary 891 892 ## FIXME: may fail with user-specified/custom family? 893 ## should be obsolete 894 if(is.null(object$family$initialize)) 895 return(do.call(object$family$family, 896 list(link=object$family$link))) 897 898 object$family 899} 900 901##' @S3method family lmResp 902family.lmResp <- function(object, ...) gaussian() 903 904##' @S3method family nlsResp 905family.nlsResp <- function(object, ...) gaussian() 906 907##' @importFrom stats fitted 908##' @S3method fitted merMod 909fitted.merMod <- function(object, ...) { 910 xx <- object@resp$mu 911 if (length(xx)==0) { 912 ## handle 'fake' objects created by simulate() 913 xx <- rep(NA,nrow(model.frame(object))) 914 } 915 if (is.null(nm <- rownames(model.frame(object)))) nm <- seq_along(xx) 916 names(xx) <- nm 917 if (!is.null(fit.na.action <- attr(model.frame(object),"na.action"))) 918 napredict(fit.na.action, xx) 919 else 920 xx 921} 922 923##' Extract the fixed-effects estimates 924##' 925##' Extract the estimates of the fixed-effects parameters from a fitted model. 926##' @name fixef 927##' @title Extract fixed-effects estimates 928##' @aliases fixef fixed.effects fixef.merMod 929##' @docType methods 930##' @param object any fitted model object from which fixed effects estimates can 931##' be extracted. 932##' @param \dots optional additional arguments. Currently none are used in any 933##' methods. 934##' @return a named, numeric vector of fixed-effects estimates. 935##' @keywords models 936##' @examples 937##' fixef(lmer(Reaction ~ Days + (1|Subject) + (0+Days|Subject), sleepstudy)) 938##' @importFrom nlme fixef 939##' @export fixef 940##' @method fixef merMod 941##' @export 942fixef.merMod <- function(object, add.dropped=FALSE, ...) { 943 X <- getME(object,"X") 944 ff <- structure(object@beta, names = dimnames(X)[[2]]) 945 if (add.dropped) { 946 if (!is.null(dd <- attr(X,"col.dropped"))) { 947 ## restore positions dropped for rank deficiency 948 vv <- numeric(length(ff)+length(dd)) 949 all.pos <- seq_along(vv) 950 kept.pos <- all.pos[-dd] 951 vv[kept.pos] <- ff 952 names(vv)[kept.pos] <- names(ff) 953 vv[dd] <- NA 954 names(vv)[dd] <- names(dd) 955 ff <- vv 956 } 957 } 958 return(ff) 959} 960 961 962getFixedFormula <- function(form) { 963 RHSForm(form) <- nobars(RHSForm(form)) 964 form 965} 966 967##' @importFrom stats formula 968##' @S3method formula merMod 969formula.merMod <- function(x, fixed.only=FALSE, random.only=FALSE, ...) { 970 if (missing(fixed.only) && random.only) fixed.only <- FALSE 971 if (fixed.only && random.only) stop("can't specify 'only fixed' and 'only random' terms") 972 if (is.null(form <- attr(x@frame,"formula"))) { 973 if (!grepl("lmer$",deparse(getCall(x)[[1]]))) 974 stop("can't find formula stored in model frame or call") 975 form <- as.formula(formula(getCall(x),...)) 976 } 977 if (fixed.only) { 978 form <- getFixedFormula(form) 979 } 980 if (random.only) { 981 ## from predict.R 982 form <- reOnly(form,response=TRUE) 983 } 984 form 985} 986 987##' @S3method isREML merMod 988isREML.merMod <- function(x, ...) as.logical(x@devcomp$dims[["REML"]]) 989 990##' @S3method isGLMM merMod 991isGLMM.merMod <- function(x,...) { 992 as.logical(x@devcomp$dims[["GLMM"]]) 993 ## or: is(x@resp,"glmResp") 994} 995 996##' @S3method isNLMM merMod 997isNLMM.merMod <- function(x,...) { 998 as.logical(x@devcomp$dims[["NLMM"]]) 999 ## or: is(x@resp,"nlsResp") 1000} 1001 1002##' @S3method isLMM merMod 1003isLMM.merMod <- function(x,...) { 1004 !isGLMM(x) && !isNLMM(x) 1005 ## or: is(x@resp,"lmerResp") ? 1006} 1007 1008npar.merMod <- function(object) { 1009 n <- length(object@beta) + length(object@theta) + 1010 object@devcomp[["dims"]][["useSc"]] 1011 ## FIXME: this is a bit of a hack: a user *might* have specified 1012 ## negative binomial family with a known theta, in which case we 1013 ## shouldn't count it as extra. Either glmer.nb needs to set a 1014 ## flag somewhere, or we need class 'nbglmerMod' to extend 'glmerMod' ... 1015 ## We do *not* want to use the 'useSc' slot (as above), because 1016 ## although theta is in some sense a scale parameter, it's not 1017 ## one in the formal sense (and isn't stored in the 'sigma' slot) 1018 if (grepl("Negative Binomial",family(object)$family)) { 1019 n <- n+1 1020 } 1021 return(n) 1022 ## TODO: how do we feel about counting the scale parameter ??? 1023} 1024 1025##' @importFrom stats logLik 1026##' @S3method logLik merMod 1027logLik.merMod <- function(object, REML = NULL, ...) { 1028 if (is.null(REML) || is.na(REML[1])) 1029 REML <- isREML(object) 1030 val <- -devCrit(object, REML = REML)/2 1031 ## dc <- object@devcomp 1032 nobs <- nobs.merMod(object) 1033 structure(val, 1034 nobs = nobs, 1035 nall = nobs, 1036 df = npar.merMod(object), 1037 ## length(object@beta) + length(object@theta) + dc$dims[["useSc"]], 1038 class = "logLik") 1039} 1040 1041##' @importFrom stats df.residual 1042##' @S3method df.residual merMod 1043## TODO: not clear whether the residual df should be based 1044## on p=length(beta) or p=length(c(theta,beta)) ... but 1045## this is just to allow things like aods3::gof to work ... 1046## 1047df.residual.merMod <- function(object, ...) { 1048 nobs(object)-npar.merMod(object) 1049} 1050 1051##' @importFrom stats logLik 1052##' @S3method model.frame merMod 1053model.frame.merMod <- function(formula, fixed.only=FALSE, ...) { 1054 fr <- formula@frame 1055 if (fixed.only) { 1056 vars <- attr(terms(fr),"varnames.fixed") 1057 if (is.null(vars)) { 1058 ## back-compatibility: saved objects pre 1.1-15 1059 ff <- formula(formula,fixed.only=TRUE) 1060 ## thanks to Thomas Leeper and Roman Lustrik, Stack Overflow 1061 ## https://stackoverflow.com/questions/18017765/extract-variables-in-formula-from-a-data-frame 1062 vars <- rownames(attr(terms.formula(ff), "factors")) 1063 } 1064 vars <- gsub("`","",vars) ## weirdness in deparsing variable names with spaces 1065 fr <- fr[vars] 1066 } 1067 fr 1068} 1069 1070##' @importFrom stats model.matrix 1071##' @S3method model.matrix merMod 1072model.matrix.merMod <- 1073 function(object, 1074 type = c("fixed", "random", "randomListRaw"), ...) { 1075 switch(type[1], 1076 "fixed" = object@pp$X, 1077 "random" = getME(object, "Z"), 1078 "randomListRaw" = mmList(object)) 1079} 1080 1081##' Dummy variables (experimental) 1082##' 1083##' Largely a wrapper for \code{model.matrix} that 1084##' accepts a factor, \code{f}, and returns a dummy 1085##' matrix with \code{nlevels(f)-1} columns. 1086dummy <- function(f, levelsToKeep){ 1087 f <- as.factor(f) 1088 mm <- model.matrix(~ 0 + f) 1089 colnames(mm) <- levels(f) 1090 1091 # sort out levels to keep 1092 missingLevels <- missing(levelsToKeep) 1093 if(missingLevels) levelsToKeep <- levels(f)[-1] 1094 if(!any(levels(f) %in% levelsToKeep)) 1095 stop("at least some of the levels in f ", 1096 "must also be present in levelsToKeep") 1097 if(!all(levelsToKeep %in% levels(f))) 1098 stop("all of the levelsToKeep must be levels of f") 1099 mm <- mm[, levelsToKeep, drop=FALSE] 1100 1101 ## # communicate that some usages are unlikely 1102 ## # to help with readibility, which is the 1103 ## # whole purpose of dummy() 1104 ## if((!missingLevels)&&(ncol(mm) > 1)) 1105 ## message("note from dummy: explicitly specifying more than one ", 1106 ## "level to keep may do little to improve readibility") 1107 return(mm) 1108} 1109 1110 1111##' @importFrom stats nobs 1112##' @S3method nobs merMod 1113nobs.merMod <- function(object, ...) nrow(object@frame) 1114 1115## used in summary.merMod(): 1116ngrps <- function(object, ...) UseMethod("ngrps") 1117 1118ngrps.default <- function(object, ...) stop("Cannot extract the number of groups from this object") 1119 1120ngrps.merMod <- function(object, ...) vapply(object@flist, nlevels, 1) 1121 1122ngrps.factor <- function(object, ...) nlevels(object) 1123 1124##' @importFrom nlme ranef 1125##' @export ranef 1126NULL 1127 1128##' Extract the modes of the random effects 1129##' 1130##' A generic function to extract the conditional modes of the random effects 1131##' from a fitted model object. For linear mixed models the conditional modes 1132##' of the random effects are also the conditional means. 1133##' 1134##' If grouping factor i has k levels and j random effects per level the ith 1135##' component of the list returned by \code{ranef} is a data frame with k rows 1136##' and j columns. If \code{condVar} is \code{TRUE} the \code{"postVar"} 1137##' attribute is an array of dimension j by j by k. The kth face of this array 1138##' is a positive definite symmetric j by j matrix. If there is only one 1139##' grouping factor in the model the variance-covariance matrix for the entire 1140##' random effects vector, conditional on the estimates of the model parameters 1141##' and on the data will be block diagonal and this j by j matrix is the kth 1142##' diagonal block. With multiple grouping factors the faces of the 1143##' \code{"postVar"} attributes are still the diagonal blocks of this 1144##' conditional variance-covariance matrix but the matrix itself is no longer 1145##' block diagonal. 1146##' @name ranef 1147##' @aliases ranef ranef.merMod 1148##' @param object an object of a class of fitted models with random effects, 1149##' typically an \code{"\linkS4class{merMod}"} object. 1150##' @param condVar an optional logical argument indicating if the conditional 1151##' variance-covariance matrices of the random effects should be added as an attribute. 1152##' @param postVar a (deprecated) synonym for \code{condVar} 1153##' @param drop an optional logical argument indicating components of the return 1154##' value that would be data frames with a single column, usually a column 1155##' called \sQuote{\code{(Intercept)}}, should be returned as named vectors. 1156##' @param whichel an optional character vector of names of grouping factors for 1157##' which the random effects should be returned. Defaults to all the grouping 1158##' factors. 1159##' @param \dots some methods for this generic function require additional 1160##' arguments. 1161##' @return A list of data frames, one for each grouping factor for the random 1162##' effects. The number of rows in the data frame is the number of levels of 1163##' the grouping factor. The number of columns is the dimension of the random 1164##' effect associated with each level of the factor. 1165##' 1166##' If \code{condVar} is \code{TRUE} each of the data frames has an attribute 1167##' called \code{"postVar"} which is a three-dimensional array with symmetric 1168##' faces. 1169##' 1170##' When \code{drop} is \code{TRUE} any components that would be data frames of 1171##' a single column are converted to named numeric vectors. 1172##' @note To produce a \dQuote{caterpillar plot} of the random effects apply 1173##' \code{\link[lattice:xyplot]{dotplot}} to the result of a call to 1174##' \code{ranef} with \code{condVar = TRUE}. 1175##' @examples 1176##' fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy) 1177##' fm2 <- lmer(Reaction ~ Days + (1|Subject) + (0+Days|Subject), sleepstudy) 1178##' fm3 <- lmer(diameter ~ (1|plate) + (1|sample), Penicillin) 1179##' ranef(fm1) 1180##' str(rr1 <- ranef(fm1, condVar = TRUE)) 1181##' dotplot(rr1) ## default 1182##' ## specify free scales in order to make Day effects more visible 1183##' dotplot(rr1,scales = list(x = list(relation = 'free')))[["Subject"]] 1184##' if(FALSE) { ##-- condVar=TRUE is not yet implemented for multiple terms -- FIXME 1185##' str(ranef(fm2, condVar = TRUE)) 1186##' } 1187##' op <- options(digits = 4) 1188##' ranef(fm3, drop = TRUE) 1189##' options(op) 1190##' @keywords models methods 1191##' @method ranef merMod 1192##' @export 1193ranef.merMod <- function(object, condVar = TRUE, drop = FALSE, 1194 whichel = names(ans), postVar = FALSE, ...) 1195{ 1196 1197 if (length(L <- list(...))>0) { 1198 warning(paste("additional arguments to ranef.merMod ignored:", 1199 paste(names(L),collapse=", "))) 1200 } 1201 if (!missing(postVar) && missing(condVar)) { 1202 warning(sQuote("postVar")," is deprecated: please use ", 1203 sQuote("condVar")," instead") 1204 condVar <- postVar 1205 } 1206 ans <- object@pp$b(1) ## not always == c(matrix(unlist(getME(object,"b")))) 1207 if (!is.null(object@flist)) { 1208 ## evaluate the list of matrices 1209 levs <- lapply(fl <- object@flist, levels) 1210 asgn <- attr(fl, "assign") 1211 cnms <- object@cnms 1212 nc <- lengths(cnms) ## number of terms 1213 ## nb <- nc * lengths(levs)[asgn] ## number of cond modes per term 1214 nb <- diff(object@Gp) 1215 nbseq <- rep.int(seq_along(nb), nb) 1216 ml <- split(ans, nbseq) 1217 for (i in seq_along(ml)) 1218 ml[[i]] <- matrix(ml[[i]], ncol = nc[i], byrow = TRUE, 1219 dimnames = list(NULL, cnms[[i]])) 1220 ## create a list of data frames corresponding to factors 1221 ans <- lapply(seq_along(fl), 1222 function(i) { 1223 m <- ml[asgn == i] 1224 b2 <- vapply(m,nrow,numeric(1)) 1225 ub2 <- unique(b2) 1226 if (length(ub2)>1) 1227 stop("differing numbers of b per group") 1228 ## if number of sets of modes != number of levels (e.g. Gaussian process/phyloglmm), 1229 ## generate numeric sequence for names 1230 1231 rnms <- if (ub2==length(levs[[i]])) levs[[i]] else seq(ub2) 1232 data.frame(do.call(cbind, m), 1233 row.names = rnms, 1234 check.names = FALSE) 1235 }) 1236 names(ans) <- names(fl) 1237 # process whichel 1238 stopifnot(is(whichel, "character")) 1239 whchL <- names(ans) %in% whichel 1240 ans <- ans[whchL] 1241 1242 if (condVar) { 1243 sigsqr <- sigma(object)^2 1244 rp <- rePos$new(object) 1245 if(any(lengths(rp$terms) > 1L)) { 1246 ## use R machinery here ... 1247 vv <- arrange.condVar(object,condVar(object, scaled=TRUE)) 1248 } else { 1249 vv <- .Call(merPredDcondVar, object@pp$ptr(), as.environment(rp)) 1250 vv <- lapply(vv, "*", sigsqr) 1251 } 1252 for (i in names(ans)) { 1253 attr(ans[[i]], "postVar") <- vv[[i]] 1254 } 1255 } 1256 if (drop) 1257 ans <- lapply(ans, function(el) 1258 { 1259 if (ncol(el) > 1) return(el) 1260 pv <- drop(attr(el, "postVar")) 1261 el <- drop(as.matrix(el)) 1262 if (!is.null(pv)) 1263 attr(el, "postVar") <- pv 1264 el 1265 }) 1266 class(ans) <- "ranef.mer" 1267 } 1268 ans 1269}## ranef.merMod 1270 1271print.ranef.mer <- function(x, ...) { 1272 print(unclass(x), ...) 1273 if(any(has.pv <- vapply(x, function(el) 1274 !is.null(attr(el, "postVar")), NA))) 1275 cat('with conditional variances for', 1276 paste(dQuote(names(x)[has.pv]), sep=", "), "\n") 1277 invisible(x) 1278} 1279 1280## try to redo refit by calling modular structure ... 1281refit2.merMod <- function(object, 1282 newresp=NULL) { 1283 ## the idea is to steal as much structure as we can from the 1284 ## previous fit, including 1285 ## * starting parameter values 1286 ## * random-effects structure 1287 ## * fixed-effects structure 1288 ## * model frame 1289 ## and jump into the modular structure at an appropriate place; 1290 ## essentially, this should merge with a smart-as-possible 1291 ## version of 'update' ... 1292 1293} 1294 1295## FIXME DRY: much of copy'n'paste from lmer() etc .. ==> become more modular (?) 1296refit.merMod <- function(object, 1297 newresp=NULL, 1298 ## formula=NULL, weights=NULL, 1299 rename.response=FALSE, 1300 maxit = 100L, ...) 1301{ 1302 1303 l... <- list(...) 1304 1305 ctrl.arg <- NULL 1306 if("control" %in% names(l...)) ctrl.arg <- l...$control 1307 1308 if(!all(names(l...) %in% c("control", "verbose"))) { 1309 warning("additional arguments to refit.merMod ignored") 1310 } 1311 ## TODO: not clear whether we should reset the names 1312 ## to the new response variable. Maybe not. 1313 1314 ## retrieve name before it gets mangled by operations on newresp 1315 newrespSub <- substitute(newresp) 1316 1317 ## for backward compatibility/functioning of refit(fit,simulate(fit)) 1318 if (is.list(newresp)) { 1319 if (length(newresp)==1) { 1320 na.action <- attr(newresp,"na.action") 1321 newresp <- newresp[[1]] 1322 attr(newresp,"na.action") <- na.action 1323 } else { 1324 stop("refit not implemented for 'newresp' lists of length > 1: ", 1325 "consider ", sQuote("lapply(object,refit)")) 1326 } 1327 } 1328 1329 ## oldresp <- object@resp$y # need to set this before deep copy, 1330 ## # otherwise it gets reset with the call 1331 ## # to setResp below 1332 1333 ## somewhat repeated from profile.merMod, but sufficiently 1334 ## different that refactoring is slightly non-trivial 1335 ## "three minutes' thought would suffice ..." 1336 control <- 1337 if (!is.null(ctrl.arg)) { 1338 if (length(ctrl.arg$optCtrl) == 0) { ## use object's version: 1339 obj.control <- object@optinfo$control 1340 ignore.pars <- c("xst", "xt") 1341 if (any(ign <- names(obj.control) %in% ignore.pars)) 1342 obj.control <- obj.control[!ign] 1343 ctrl.arg$optCtrl <- obj.control 1344 } 1345 ctrl.arg 1346 } 1347 else if (isGLMM(object)) 1348 glmerControl() 1349 else 1350 lmerControl() 1351 if (object@optinfo$optimizer == "optimx") { 1352 control$optCtrl <- object@optinfo$control 1353 } 1354 ## we need this stuff defined before we call .glmerLaplace below ... 1355 pp <- object@pp$copy() 1356 dc <- object@devcomp 1357 nAGQ <- dc$dims["nAGQ"] # possibly NA 1358 nth <- dc$dims[["nth"]] 1359 verbose <- l...$verbose; if (is.null(verbose)) verbose <- 0L 1360 if (!is.null(newresp)) { 1361 ## update call and model frame with new response 1362 rcol <- attr(attr(model.frame(object), "terms"), "response") 1363 if (rename.response) { 1364 attr(object@frame,"formula")[[2]] <- object@call$formula[[2]] <- 1365 newrespSub 1366 names(object@frame)[rcol] <- deparse(newrespSub) 1367 } 1368 if (!is.null(na.act <- attr(object@frame,"na.action")) && 1369 is.null(attr(newresp,"na.action"))) { 1370 ## will only get here if na.action is 'na.omit' or 'na.exclude' 1371 ## *and* newresp does not have an 'na.action' attribute 1372 ## indicating that NAs have already been filtered 1373 newresp <- if (is.matrix(newresp)) 1374 newresp[-na.act, ] 1375 else newresp[-na.act] 1376 } 1377 object@frame[,rcol] <- newresp 1378 1379 ## modFrame <- model.frame(object) 1380 ## modFrame[, attr(terms(modFrame), "response")] <- newresp 1381 } 1382 1383 rr <- if(isLMM(object)) 1384 mkRespMod(model.frame(object), REML = isREML(object)) 1385 else if(isGLMM(object)) { 1386 mkRespMod(model.frame(object), family = family(object)) 1387 } else 1388 stop("refit.merMod not working for nonlinear mixed models.\n", 1389 "try update.merMod instead.") 1390 1391 if(!is.null(newresp)) { 1392 if(family(object)$family == "binomial") { 1393 ## re-do conversion of two-column matrix and factor 1394 ## responses to proportion/weights format 1395 if (is.matrix(newresp) && ncol(newresp) == 2) { 1396 ntot <- rowSums(newresp) 1397 ## FIXME: test what happens for (0,0) rows 1398 newresp <- newresp[,1]/ntot 1399 rr$setWeights(ntot) 1400 } 1401 if (is.factor(newresp)) { 1402 ## FIXME: would be better to do this consistently with 1403 ## whatever machinery is used in glm/glm.fit/glmer ?? 1404 newresp <- as.numeric(newresp)-1 1405 } 1406 } 1407 ## if (isGLMM(object) && rr$family$family=="binomial") { 1408 1409 ## } 1410 stopifnot(length(newresp <- as.numeric(as.vector(newresp))) == 1411 length(rr$y)) 1412 1413 } 1414 1415 1416 if (isGLMM(object)) { 1417 GQmat <- GHrule(nAGQ) 1418 1419 if (nAGQ <= 1) { 1420 glmerPwrssUpdate(pp,rr, control$tolPwrss, GQmat, maxit=maxit) 1421 } else { 1422 glmerPwrssUpdate(pp,rr, control$tolPwrss, GQmat, maxit=maxit, 1423 grpFac = object@flist[[1]]) 1424 } 1425 } 1426 1427 devlist <- 1428 if (isGLMM(object)) { 1429 baseOffset <- forceCopy(object@resp$offset) 1430 1431 list(tolPwrss= dc$cmp [["tolPwrss"]], 1432 compDev = dc$dims[["compDev"]], 1433 nAGQ = unname(nAGQ), 1434 lp0 = pp$linPred(1), ## object@resp$eta - baseOffset, 1435 baseOffset = baseOffset, 1436 pwrssUpdate = glmerPwrssUpdate, 1437 ## save GQmat in the object and use that instead of nAGQ 1438 GQmat = GHrule(nAGQ), 1439 fac = object@flist[[1]], 1440 pp=pp, resp=rr, u0=pp$u0, verbose=verbose, dpars=seq_len(nth)) 1441 } else 1442 list(pp=pp, resp=rr, u0=pp$u0, verbose=verbose, dpars=seq_len(nth)) 1443 ff <- mkdevfun(list2env(devlist), nAGQ=nAGQ, maxit=maxit, verbose=verbose) 1444 ## rho <- environment(ff) == list2env(devlist) 1445 xst <- rep.int(0.1, nth) 1446 x0 <- pp$theta 1447 lower <- object@lower 1448 if (!is.na(nAGQ) && nAGQ > 0L) { 1449 xst <- c(xst, sqrt(diag(pp$unsc()))) 1450 x0 <- c(x0, unname(fixef(object))) 1451 lower <- c(lower, rep(-Inf,length(x0)-length(lower))) 1452 } 1453 ## control <- c(control,list(xst=0.2*xst, xt=xst*0.0001)) 1454 ## FIX ME: allow use.last.params to be passed through 1455 calc.derivs <- !is.null(object@optinfo$derivs) 1456 ## if(isGLMM(object)) { 1457 ## rho$resp$updateWts() 1458 ## rho$pp$updateDecomp() 1459 ## rho$lp0 <- rho$pp$linPred(1) 1460 ## } 1461 optimizer <- object@optinfo$optimizer 1462 if (!is.null(newopt <- ctrl.arg$optimizer)) { 1463 ## we might end up with a length-2 optimizer vector ... 1464 ## use the *last* element 1465 optimizer <- newopt[length(newopt)] 1466 } 1467 opt <- optwrap(optimizer, 1468 ff, x0, lower=lower, control=control$optCtrl, 1469 calc.derivs=calc.derivs) 1470 cc <- checkConv(attr(opt,"derivs"),opt$par, 1471 ## FIXME: was there a reason that ctrl was passed 1472 ## via the call slot? it was causing problems 1473 ## when optTheta called refit (github issue #173) 1474 # ctrl = eval(object@call$control)$checkConv, 1475 ctrl = control$checkConv, 1476 lbound=lower) 1477 if (isGLMM(object)) rr$setOffset(baseOffset) 1478 mkMerMod(environment(ff), opt, 1479 list(flist=object@flist, cnms=object@cnms, 1480 Gp=object@Gp, lower=object@lower), 1481 object@frame, getCall(object), cc) 1482} 1483 1484refitML.merMod <- function (x, optimizer="bobyqa", ...) { 1485 ## FIXME: optimizer is set to 'bobyqa' for back-compatibility, but that's not 1486 ## consistent with lmer (default NM). Should be based on internally stored 'optimizer' value 1487 if (!isREML(x)) return(x) 1488 stopifnot(is(rr <- x@resp, "lmerResp")) 1489 rho <- new.env(parent=parent.env(environment())) 1490 rho$resp <- new(class(rr), y=rr$y, offset=rr$offset, weights=rr$weights, REML=0L) 1491 xpp <- x@pp$copy() 1492 rho$pp <- new(class(xpp), X=xpp$X, Zt=xpp$Zt, Lambdat=xpp$Lambdat, 1493 Lind=xpp$Lind, theta=xpp$theta, n=nrow(xpp$X)) 1494 devfun <- mkdevfun(rho, 0L) # FIXME? also pass {verbose, maxit, control} 1495 opt <- ## "smart" calc.derivs rules 1496 if(optimizer == "bobyqa" && !any("calc.derivs" == ...names())) 1497 optwrap(optimizer, devfun, x@theta, lower=x@lower, calc.derivs=TRUE, ...) 1498 else 1499 optwrap(optimizer, devfun, x@theta, lower=x@lower, ...) 1500 ## FIXME: Should be able to call mkMerMod() here, and be done 1501 n <- length(rr$y) 1502 pp <- rho$pp 1503 p <- ncol(pp$X) 1504 dims <- c(N=n, n=n, p=p, nmp=n-p, q=nrow(pp$Zt), nth=length(pp$theta), 1505 useSc=1L, reTrms=length(x@cnms), 1506 spFe=0L, REML=0L, GLMM=0L, NLMM=0L)#, nAGQ=NA_integer_) 1507 wrss <- rho$resp$wrss() 1508 ussq <- pp$sqrL(1) 1509 pwrss <- wrss + ussq 1510 cmp <- c(ldL2=pp$ldL2(), ldRX2=pp$ldRX2(), wrss=wrss, ussq=ussq, 1511 pwrss=pwrss, drsum=NA, dev=opt$fval, REML=NA, 1512 sigmaML=sqrt(pwrss/n), sigmaREML=sqrt(pwrss/(n-p))) 1513 ## modify the call to have REML=FALSE. (without evaluating the call!) 1514 cl <- x@call 1515 cl[["REML"]] <- FALSE 1516 new("lmerMod", call = cl, frame=x@frame, flist=x@flist, 1517 cnms=x@cnms, theta=pp$theta, beta=pp$delb, u=pp$delu, 1518 optinfo = .optinfo(opt), 1519 lower=x@lower, devcomp=list(cmp=cmp, dims=dims), pp=pp, resp=rho$resp, 1520 Gp=x@Gp) 1521} 1522 1523##' residuals of merMod objects --> ../man/residuals.merMod.Rd 1524##' @param object a fitted [g]lmer (\code{merMod}) object 1525##' @param type type of residuals 1526##' @param scaled scale residuals by residual standard deviation (=scale parameter)? 1527##' @param \dots additional arguments (ignored: for method compatibility) 1528residuals.merMod <- function(object, 1529 type = if(isGLMM(object)) "deviance" else "response", 1530 scaled = FALSE, ...) { 1531 r <- residuals(object@resp, type,...) 1532 fr <- model.frame(object) 1533 if (is.null(nm <- rownames(fr))) nm <- seq_along(r) 1534 names(r) <- nm 1535 if (scaled) r <- r/sigma(object) 1536 if (!is.null(na.action <- attr(fr, "na.action"))) 1537 naresid(na.action, r) 1538 else 1539 r 1540} 1541 1542##' @rdname residuals.merMod 1543##' @S3method residuals lmResp 1544##' @method residuals lmResp 1545residuals.lmResp <- function(object, 1546 type = c("working", "response", "deviance", 1547 "pearson", "partial"), ...) { 1548 y <- object$y 1549 r <- object$wtres 1550 mu <- object$mu 1551 switch(match.arg(type), 1552 working =, 1553 response = y-mu, 1554 deviance =, 1555 pearson = r, 1556 partial = stop(gettextf("partial residuals are not implemented yet"), 1557 call. = FALSE) 1558 ) 1559} 1560 1561##' @rdname residuals.merMod 1562##' @S3method residuals glmResp 1563##' @method residuals glmResp 1564residuals.glmResp <- function(object, type = c("deviance", "pearson", 1565 "working", "response", "partial"), 1566 ...) { 1567 type <- match.arg(type) 1568 y <- object$y 1569 mu <- object$mu 1570 switch(type, 1571 deviance = { 1572 d.res <- sqrt(object$devResid()) 1573 ifelse(y > mu, d.res, -d.res) 1574 }, 1575 pearson = object$wtres, 1576 working = object$wrkResids(), 1577 response = y - mu, 1578 partial = stop(gettextf("partial residuals are not implemented yet"), 1579 call. = FALSE) 1580 ) 1581 1582} 1583 1584hatvalues.merMod <- function(model, fullHatMatrix = FALSE, ...) { 1585 if(isGLMM(model)) warning("the hat matrix may not make sense for GLMMs") 1586 ## FIXME: add restriction for NLMMs? 1587 1588 ## prior weights, W ^ {1/2} : 1589 sqrtW <- Diagonal(x = sqrt(weights(model, type = "prior"))) 1590 res <- with(getME(model, c("L", "Lambdat", "Zt", "RX", "X", "RZX")), { 1591 ## CL:= right factor of the random-effects component of the hat matrix (64) 1592 CL <- solve(L, solve(L, Lambdat %*% Zt %*% sqrtW, 1593 system = "P"), system = "L") 1594 ## CR:= right factor of the fixed-effects component of the hat matrix (65) 1595 ## {MM (FIXME Matrix): t(.) %*% here faster than crossprod()} 1596 CR <- solve(t(RX), t(X) %*% sqrtW - crossprod(RZX, CL)) 1597 if(fullHatMatrix) ## H = (C_L^T C_L + C_R^T C_R) (63) 1598 crossprod(CL) + crossprod(CR) 1599 else ## diagonal of the hat matrix, diag(H) : 1600 colSums(CR^2) + colSums(CL^2) 1601 }) 1602 napredict(attr(model.frame(model),"na.action"),res) 1603} 1604 1605 1606###----- Printing etc ---------------------------- 1607 1608## lme4.0, for GLMM had 1609## 'Generalized linear mixed model fit by the Laplace approximation' 1610## 'Generalized linear mixed model fit by the adaptive Gaussian Hermite approximation' 1611## so did *not* mention "maximum likelihood" at all in the GLMM case 1612 1613 1614methTitle <- function(dims) { # dims == object@devcomp$dims 1615 GLMM <- dims[["GLMM"]] 1616 kind <- switch(1L + GLMM * 2L + dims[["NLMM"]], 1617 "Linear", "Nonlinear", 1618 "Generalized linear", "Generalized nonlinear") 1619 paste(kind, "mixed model fit by", 1620 if(dims[["REML"]]) "REML" 1621 else paste("maximum likelihood", 1622 if(GLMM) { 1623 ## TODO? Use shorter wording here, for (new) 'long = FALSE' argument 1624 if((nAGQ <- dims[["nAGQ"]]) == 1) 1625 "(Laplace Approximation)" 1626 else 1627 sprintf("(Adaptive Gauss-Hermite Quadrature, nAGQ = %d)", 1628 nAGQ) 1629 })) 1630} 1631 1632 1633cat.f <- function(...) cat(..., fill = TRUE) 1634 1635famlink <- function(object, resp = object@resp) { 1636 if(is(resp, "glmResp")) 1637 resp$family[c("family", "link")] 1638 else list(family = NULL, link = NULL) 1639} 1640 1641##' @title print method title 1642##' @param mtit the result of methTitle(obj) 1643##' @param class typically class(obj) 1644.prt.methTit <- function(mtit, class) { 1645 if(nchar(mtit) + 5 + nchar(class) > (w <- getOption("width"))) { 1646 ## wrap around 1647 mtit <- strwrap(mtit, width = w - 2, exdent = 2) 1648 cat(mtit, " [",class,"]", sep = "", fill = TRUE) 1649 } else ## previous: simple one-liner 1650 cat(sprintf("%s ['%s']\n", mtit, class)) 1651} 1652 1653.prt.family <- function(famL) { 1654 if (!is.null(f <- famL$family)) { 1655 cat.f(" Family:", f, 1656 if(!is.null(ll <- famL$link)) paste(" (", ll, ")")) 1657 } 1658} 1659 1660.prt.resids <- function(resids, digits, title = "Scaled residuals:", ...) { 1661 cat(title,"\n") 1662 ## FIXME: need testing code 1663 rq <- setNames(zapsmall(quantile(resids, na.rm=TRUE), digits + 1L), 1664 c("Min", "1Q", "Median", "3Q", "Max")) 1665 print(rq, digits = digits, ...) 1666 cat("\n") 1667} 1668 1669.prt.call <- function(call, long = TRUE) { 1670 if (!is.null(cc <- call$formula)) 1671 cat.f("Formula:", deparse(cc)) 1672 if (!is.null(cc <- call$data)) 1673 cat.f(" Data:", deparse(cc)) 1674 if (!is.null(cc <- call$weights)) 1675 cat.f("Weights:", deparse(cc)) 1676 if (!is.null(cc <- call$offset)) 1677 cat.f(" Offset:", deparse(cc)) 1678 if (long && length(cc <- call$control) && 1679 !identical((dc <- deparse(cc)), "lmerControl()")) 1680 ## && !identical(eval(cc), lmerControl())) 1681 cat.f("Control:", dc) 1682 if (!is.null(cc <- call$subset)) 1683 cat.f(" Subset:", deparse(cc)) 1684} 1685 1686##' @title Extract Log Likelihood, AIC, and related statics from a Fitted LMM 1687##' @param object a LMM model fit 1688##' @param devianceFUN the function to be used for computing the deviance; should not be changed for \pkg{lme4} created objects. 1689##' @param chkREML optional logical indicating \code{object} maybe a REML fit. 1690##' @param devcomp 1691##' @return a list with components \item{logLik} and \code{AICtab} where the first is \code{\link{logLik(object)}} and \code{AICtab} is a "table" of AIC, BIC, logLik, deviance, df.residual() values. 1692llikAIC <- function(object, devianceFUN = devCrit, chkREML = TRUE, devcomp = object@devcomp) { 1693 llik <- logLik(object) # returns NA for a REML fit - maybe change? 1694 AICstats <- { 1695 if(chkREML && devcomp$dims[["REML"]]) 1696 devcomp$cmp["REML"] # *no* likelihood stats here 1697 else { 1698 c(AIC = AIC(llik), BIC = BIC(llik), logLik = c(llik), 1699 deviance = devianceFUN(object), 1700 df.resid = df.residual(object)) 1701 } 1702 } 1703 list(logLik = llik, AICtab = AICstats) 1704} 1705 1706.prt.aictab <- function(aictab, digits = 1) { 1707 t.4 <- round(aictab, digits) 1708 if (length(aictab) == 1 && names(aictab) == "REML") 1709 cat.f("REML criterion at convergence:", t.4) 1710 else { 1711 ## slight hack to get residual df formatted as an integer 1712 t.4F <- format(t.4) 1713 t.4F["df.resid"] <- format(t.4["df.resid"]) 1714 print(t.4F, quote = FALSE) 1715 } 1716} 1717 1718.prt.VC <- function(varcor, digits, comp, formatter = format, ...) { 1719 cat("Random effects:\n") 1720 fVC <- if(missing(comp)) 1721 formatVC(varcor, digits = digits, formatter = formatter) 1722 else 1723 formatVC(varcor, digits = digits, formatter = formatter, comp = comp) 1724 print(fVC, quote = FALSE, digits = digits, ...) 1725} 1726 1727.prt.grps <- function(ngrps, nobs) { 1728 cat(sprintf("Number of obs: %d, groups: ", nobs), 1729 paste(paste(names(ngrps), ngrps, sep = ", "), collapse = "; "), 1730 fill = TRUE) 1731} 1732 1733## FIXME: print header ("Warnings:\n") ? 1734## change position in output? comes at the very end, could get lost ... 1735.prt.warn <- function(optinfo, summary=FALSE, ...) { 1736 if(length(optinfo) == 0) return() # currently, e.g., from refitML() 1737 ## check all warning slots: print numbers of warnings (if any) 1738 cc <- optinfo$conv$opt 1739 msgs <- unlist(optinfo$conv$lme4$messages) 1740 ## can't put nmsgs/nwarnings compactly into || expression 1741 ## because of short-circuiting 1742 nmsgs <- length(msgs) 1743 warnings <- optinfo$warnings 1744 nwarnings <- length(warnings) 1745 if (cc > 0 || nmsgs > 0 || nwarnings > 0) { 1746 m <- if (cc==0) { 1747 "(OK)" 1748 } else if (!is.null(optinfo$message)) { 1749 sprintf("(%s)",optinfo$message) 1750 } else "" 1751 convmsg <- sprintf("optimizer (%s) convergence code: %d %s", 1752 optinfo$optimizer, cc, m) 1753 if (summary) { 1754 cat(convmsg,sprintf("; %d optimizer warnings; %d lme4 warnings", 1755 nwarnings,nmsgs),"\n") 1756 } else { 1757 cat(convmsg, 1758 msgs, 1759 unlist(warnings), 1760 sep="\n") 1761 cat("\n") 1762 } 1763 } 1764} 1765 1766## options(lme4.summary.cor.max = 20) --> ./hooks.R 1767## ~~~~~~~~ 1768## was .summary.cor.max <- 20 a lme4-namespace hidden global variable 1769 1770## This is modeled a bit after print.summary.lm : 1771## Prints *both* 'mer' and 'merenv' - as it uses summary(x) mainly 1772##' @S3method print summary.merMod 1773print.summary.merMod <- function(x, digits = max(3, getOption("digits") - 3), 1774 correlation = NULL, symbolic.cor = FALSE, 1775 signif.stars = getOption("show.signif.stars"), 1776 ranef.comp = c("Variance", "Std.Dev."), 1777 show.resids = TRUE, ...) 1778{ 1779 .prt.methTit(x$methTitle, x$objClass) 1780 .prt.family(x) 1781 .prt.call(x$call); cat("\n") 1782 .prt.aictab(x$AICtab); cat("\n") 1783 if (show.resids) 1784 ## need residuals.merMod() rather than residuals(): 1785 ## summary.merMod has no residuals method 1786 .prt.resids(x$residuals, digits = digits) 1787 .prt.VC(x$varcor, digits = digits, useScale = x$useScale, 1788 comp = ranef.comp, ...) 1789 .prt.grps(x$ngrps, nobs = x$devcomp$dims[["n"]]) 1790 1791 p <- nrow(x$coefficients) 1792 if (p > 0) { 1793 cat("\nFixed effects:\n") 1794 printCoefmat(x$coefficients, # too radical: zap.ind = 3, #, tst.ind = 4 1795 digits = digits, signif.stars = signif.stars) 1796 ## do not show correlation when summary(*, correlation=FALSE) was used: 1797 hasCor <- !is.null(VC <- x$vcov) && !is.null(VC@factors$correlation) 1798 if(is.null(correlation)) { # default 1799 cor.max <- getOption("lme4.summary.cor.max") 1800 correlation <- hasCor && p <= cor.max 1801 if(!correlation && p > cor.max) { 1802 nam <- deparse(substitute(x)) 1803 if(length(nam) > 1 || nchar(nam) >= 32) nam <- "...." 1804 message(sprintf(paste( 1805 "\nCorrelation matrix not shown by default, as p = %d > %d.", 1806 "Use print(%s, correlation=TRUE) or", 1807 " vcov(%s) if you need it\n", sep = "\n"), 1808 p, cor.max, nam, nam)) 1809 } 1810 } 1811 else if(!is.logical(correlation)) stop("'correlation' must be NULL or logical") 1812 if(correlation) { 1813 if(is.null(VC)) VC <- vcov(x, correlation = TRUE) 1814 corF <- VC@factors$correlation 1815 if (is.null(corF)) { # can this still happen? 1816 message("\nCorrelation of fixed effects could have been required in summary()") 1817 corF <- cov2cor(VC) 1818 } 1819 p <- ncol(corF) 1820 if (p > 1) { 1821 rn <- rownames(x$coefficients) 1822 rns <- abbreviate(rn, minlength = 11) 1823 cat("\nCorrelation of Fixed Effects:\n") 1824 if (is.logical(symbolic.cor) && symbolic.cor) { 1825 corf <- as(corF, "matrix") 1826 dimnames(corf) <- list(rns, 1827 abbreviate(rn, minlength = 1, strict = TRUE)) 1828 print(symnum(corf)) 1829 } else { 1830 corf <- matrix(format(round(corF@x, 3), nsmall = 3), 1831 ncol = p, 1832 dimnames = list(rns, abbreviate(rn, minlength = 6))) 1833 corf[!lower.tri(corf)] <- "" 1834 print(corf[-1, -p, drop = FALSE], quote = FALSE) 1835 } ## !symbolic.cor 1836 } ## if (p > 1) 1837 } ## if (correlation) 1838 } ## if (p>0) 1839 1840 if(length(x$fitMsgs) && any(nchar(x$fitMsgs) > 0)) { 1841 cat("fit warnings:\n"); writeLines(x$fitMsgs) 1842 } 1843 .prt.warn(x$optinfo,summary=FALSE) 1844 invisible(x) 1845}## print.summary.merMod 1846 1847 1848##' @S3method print merMod 1849print.merMod <- function(x, digits = max(3, getOption("digits") - 3), 1850 correlation = NULL, symbolic.cor = FALSE, 1851 signif.stars = getOption("show.signif.stars"), 1852 ranef.comp = "Std.Dev.", ...) 1853{ 1854 dims <- x@devcomp$dims 1855 .prt.methTit(methTitle(dims), class(x)) 1856 .prt.family(famlink(x, resp = x@resp)) 1857 .prt.call(x@call, long = FALSE) 1858 ## useScale <- as.logical(dims[["useSc"]]) 1859 1860 llAIC <- llikAIC(x) 1861 .prt.aictab(llAIC$AICtab, 4) 1862 varcor <- VarCorr(x) 1863 .prt.VC(varcor, digits = digits, comp = ranef.comp, ...) 1864 ngrps <- vapply(x@flist, nlevels, 0L) 1865 .prt.grps(ngrps, nobs = dims[["n"]]) 1866 if(length(cf <- fixef(x)) > 0) { 1867 cat("Fixed Effects:\n") 1868 print.default(format(cf, digits = digits), 1869 print.gap = 2L, quote = FALSE, ...) 1870 } else cat("No fixed effect coefficients\n") 1871 1872 fitMsgs <- .merMod.msgs(x) 1873 if(any(nchar(fitMsgs) > 0)) { 1874 cat("fit warnings:\n"); writeLines(fitMsgs) 1875 } 1876 .prt.warn(x@optinfo,summary=TRUE) 1877 invisible(x) 1878} 1879 1880##' @exportMethod show 1881setMethod("show", "merMod", function(object) print.merMod(object)) 1882 1883##' Return the deviance component list 1884devcomp <- function(x) { 1885 .Deprecated("getME(., \"devcomp\")") 1886 stopifnot(is(x, "merMod")) 1887 x@devcomp 1888} 1889 1890##' @exportMethod getL 1891setMethod("getL", "merMod", function(x) { 1892 .Deprecated("getME(., \"L\")") 1893 getME(x, "L") 1894}) 1895 1896##' used by tnames 1897mkPfun <- function(diag.only = FALSE, old = TRUE, prefix = NULL){ 1898 local({ 1899 function(g,e) { 1900 mm <- outer(e,e,paste,sep = ".") 1901 if(old) { 1902 diag(mm) <- e 1903 } else { 1904 mm[] <- paste(mm,g,sep = "|") 1905 if (!is.null(prefix)) mm[] <- paste(prefix[2],mm,sep = "_") 1906 diag(mm) <- paste(e,g,sep = "|") 1907 if (!is.null(prefix)) diag(mm) <- paste(prefix[1],diag(mm),sep = "_") 1908 } 1909 mm <- if (diag.only) diag(mm) else mm[lower.tri(mm,diag = TRUE)] 1910 if(old) paste(g,mm,sep = ".") else mm 1911 } 1912 }) 1913} 1914 1915##' Construct names of individual theta/sd:cor components 1916##' 1917##' @param object a fitted model 1918##' @param diag.only include only diagonal elements? 1919##' @param old (logical) give backward-compatible results? 1920##' @param prefix a character vector with two elements giving the prefix 1921##' for diagonal (e.g. "sd") and off-diagonal (e.g. "cor") elements 1922tnames <- function(object, diag.only = FALSE, old = TRUE, prefix = NULL) { 1923 pfun <- mkPfun(diag.only=diag.only, old=old, prefix=prefix) 1924 c(unlist(mapply(pfun, names(object@cnms), object@cnms))) 1925} 1926 1927## -> ../man/getME.Rd 1928getME <- function(object, name, ...) UseMethod("getME") 1929 1930 1931##' Extract or Get Generalized Components from a Fitted Mixed Effects Model 1932getME.merMod <- function(object, 1933 name = c("X", "Z","Zt", "Ztlist", "mmList", 1934 "y", "mu", "u", "b", 1935 "Gp", "Tp", 1936 "L", "Lambda", "Lambdat", "Lind", "Tlist", "A", 1937 "RX", "RZX", "sigma", 1938 "flist", 1939 "fixef", "beta", "theta", "ST", 1940 "REML", "is_REML", 1941 "n_rtrms", "n_rfacs", 1942 "N", "n", "p", "q", 1943 "p_i", "l_i", "q_i", "k", "m_i", "m", 1944 "cnms", 1945 "devcomp", "offset", "lower", 1946 "devfun", 1947 "glmer.nb.theta" 1948 ), ...) 1949{ 1950 if(missing(name)) stop("'name' must not be missing") 1951 ## Deal with multiple names -- "FIXME" is inefficiently redoing things 1952 if (length(name <- as.character(name)) > 1) { 1953 names(name) <- name 1954 return(lapply(name, getME, object = object)) 1955 } 1956 if(name == "ALL") ## recursively get all provided components 1957 return(sapply(eval(formals()$name), 1958 getME.merMod, object=object, simplify=FALSE)) 1959 1960 stopifnot(is(object,"merMod")) 1961 name <- match.arg(name) 1962 rsp <- object@resp 1963 PR <- object@pp 1964 dc <- object@devcomp 1965 th <- object@theta 1966 cnms <- object@cnms 1967 dims <- dc $ dims 1968 Tpfun <- function(cnms) { 1969 ltsize <- function(n) n*(n+1)/2 # lower triangle size 1970 cLen <- cumsum(ltsize(lengths(cnms))) 1971 setNames(c(0, cLen), 1972 c("beg__", names(cnms))) ## such that diff( Tp ) is well-named 1973 } 1974 ## mmList. <- mmList(object) 1975 if(any(name == c("p_i", "q_i", "m_i"))) 1976 p_i <- vapply(mmList(object), ncol, 1L) 1977 if(any(name == c("l_i", "q_i"))) 1978 l_i <- vapply(object@flist, nlevels, 1L) 1979 switch(name, 1980 "X" = PR$X, ## ok ? - check -- use model.matrix() method instead? 1981 "Z" = t(PR$Zt), 1982 "Zt" = PR$Zt, 1983 "Ztlist" = 1984 { 1985 getInds <- function(i) { 1986 n <- diff(object@Gp)[i] ## number of elements in this block 1987 nt <- length(cnms[[i]]) ## number of REs 1988 inds <- lapply(seq(nt), seq, to = n, by = nt) ## pull out individual RE indices 1989 inds <- lapply(inds,function(x) x + object@Gp[i]) ## add group offset 1990 } 1991 inds <- do.call(c,lapply(seq_along(cnms),getInds)) 1992 setNames(lapply(inds,function(i) PR$Zt[i,]), 1993 tnames(object,diag.only = TRUE)) 1994 }, 1995 "mmList" = mmList.merMod(object), 1996 "y" = rsp$y, 1997 "mu" = rsp$mu, 1998 "u" = object@u, 1999 "b" = crossprod(PR$Lambdat, object@u), # == Lambda %*% u 2000 "L" = PR$ L(), 2001 "Lambda" = t(PR$ Lambdat), 2002 "Lambdat" = PR$ Lambdat, 2003 "A" = PR$Lambdat %*% PR$Zt, 2004 "Lind" = PR$ Lind, 2005 "RX" = structure(PR$RX(), dimnames = list(colnames(PR$X), colnames(PR$X))), ## maybe add names elsewhere? 2006 "RZX" = structure(PR$RZX, dimnames = list(NULL, colnames(PR$X))), ## maybe add names elsewhere? 2007 "sigma" = sigma(object), 2008 "Gp" = object@Gp, 2009 "Tp" = Tpfun(cnms), # "term-wise theta pointer" 2010 "flist" = object@flist, 2011 "fixef" = fixef(object), 2012 "beta" = object@beta, 2013 "theta" = setNames(th, tnames(object)), 2014 "ST" = setNames(vec2STlist(object@theta, n = lengths(cnms)), 2015 names(cnms)), 2016 "Tlist" = { 2017 nc <- lengths(cnms) # number of columns per term 2018 nt <- length(nc) # number of random-effects terms 2019 ans <- vector("list",nt) 2020 names(ans) <- names(nc) 2021 pos <- 0L 2022 for (i in 1:nt) { # will need modification for more general Lambda 2023 nci <- nc[i] 2024 tt <- matrix(0., nci, nci) 2025 inds <- lower.tri(tt, diag=TRUE) 2026 nthi <- sum(inds) 2027 tt[inds] <- th[pos + seq_len(nthi)] 2028 pos <- pos + nthi 2029 ans[[i]] <- tt 2030 } 2031 ans 2032 }, 2033 "REML" = dims[["REML"]], 2034 "is_REML" = isREML(object), 2035 ## number of random-effects terms 2036 "n_rtrms" = length(cnms), 2037 ## number of random-effects grouping factors 2038 "n_rfacs" = length(object@flist), 2039 "N" = dims[["N"]], 2040 "n" = dims[["n"]], 2041 "p" = dims[["p"]], 2042 "q" = dims[["q"]], 2043 "p_i" = p_i, 2044 "l_i" = l_i, 2045 "q_i" = p_i * l_i, 2046 "k" = length(cnms), 2047 "m_i" = choose(p_i + 1, 2), 2048 "m" = dims[["nth"]], 2049 "cnms" = cnms, 2050 "devcomp" = dc, 2051 "offset" = rsp$offset, 2052 "lower" = setNames(object@lower, tnames(object)), 2053 "devfun" = { 2054 verbose <- getCall(object)$verbose; if (is.null(verbose)) verbose <- 0L 2055 if (isGLMM(object)) { 2056 reTrms <- getME(object,c("Zt","theta","Lambdat","Lind","flist","cnms")) 2057 d1 <- mkGlmerDevfun(object@frame, getME(object,"X"), 2058 reTrms=reTrms, family(object), 2059 verbose=verbose) 2060 nAGQ <- object@devcomp$dims[["nAGQ"]] 2061 updateGlmerDevfun(d1, reTrms, nAGQ=nAGQ) 2062 } else { 2063 ## copied from refit ... DRY ... 2064 devlist <- list(pp=PR, resp=rsp, u0=PR$u0, dpars=seq_along(PR$theta), verbose=verbose) 2065 mkdevfun(rho=list2env(devlist), 2066 ## FIXME: fragile ... // also pass 'maxit' ? 2067 verbose=verbose, control=object@optinfo$control) 2068 } 2069 }, 2070 ## FIXME: current version gives lower bounds for theta parameters only: 2071 ## -- these must be extended for [GN]LMMs -- give extended value including -Inf values for beta values? 2072 2073 "glmer.nb.theta" = if(isGLMM(object) && isNBfamily(rsp$family$family)) 2074 getNBdisp(object) else NA, 2075 2076 "..foo.." = # placeholder! 2077 stop(gettextf("'%s' is not implemented yet", 2078 sprintf("getME(*, \"%s\")", name))), 2079 ## otherwise 2080 stop(sprintf("Mixed-Effects extraction of '%s' is not available for class \"%s\"", 2081 name, class(object)))) 2082}## {getME} 2083 2084##' @importMethodsFrom Matrix t %*% crossprod diag tcrossprod 2085##' @importClassesFrom Matrix dgCMatrix dpoMatrix corMatrix 2086NULL 2087 2088## Extract the conditional variance-covariance matrix of the fixed-effects 2089## parameters 2090vcov.merMod <- function(object, correlation = TRUE, sigm = sigma(object), 2091 use.hessian = NULL, ...) 2092{ 2093 2094 hess.avail <- 2095 ## (1) numerical Hessian computed? 2096 (!is.null(h <- object@optinfo$derivs$Hessian) && 2097 ## (2) does Hessian include fixed-effect parameters? 2098 nrow(h) > (ntheta <- length(getME(object,"theta")))) 2099 if (is.null(use.hessian)) use.hessian <- hess.avail 2100 if (use.hessian && !hess.avail) stop(shQuote("use.hessian"), 2101 "=TRUE specified, ", 2102 "but Hessian is unavailable") 2103 calc.vcov.hess <- function(h) { 2104 ## invert 2*Hessian, catching errors and forcing symmetric result 2105 ## ~= forceSymmetric(solve(h/2)[i,i]) : solve(h/2) = 2*solve(h) 2106 h <- tryCatch(solve(h), 2107 error=function(e) matrix(NA,nrow=nrow(h),ncol=ncol(h))) 2108 i <- -seq_len(ntheta) ## drop var-cov parameters 2109 h <- h[i,i] 2110 forceSymmetric(h + t(h)) 2111 } 2112 2113 ## alternately: calculate var-cov from implicit (RX) information 2114 ## provided by fit (always the case for lmerMods) 2115 V <- sigm^2 * object@pp$unsc() 2116 2117 if (hess.avail) { 2118 V.hess <- calc.vcov.hess(h) 2119 bad.V.hess <- any(is.na(V.hess)) 2120 if (!bad.V.hess) { 2121 ## another 'bad var-cov' check: positive definite? 2122 e.hess <- eigen(V.hess,symmetric = TRUE,only.values = TRUE)$values 2123 if (min(e.hess) <= 0) bad.V.hess <- TRUE 2124 } 2125 } 2126 if (!use.hessian && hess.avail) { 2127 ## if hessian is available, go ahead and check 2128 ## for similarity with the RX-based estimate 2129 var.hess.tol <- 1e-4 # FIXME: should var.hess.tol be user controlled? 2130 if (!bad.V.hess && any(abs(V-V.hess) > var.hess.tol * V.hess)) 2131 warning("variance-covariance matrix computed ", 2132 "from finite-difference Hessian\nand ", 2133 "from RX differ by >",var.hess.tol,": ", 2134 "consider ",shQuote("use.hessian=TRUE")) 2135 } 2136 if (use.hessian) { 2137 if (!bad.V.hess) { 2138 V <- V.hess 2139 } else { 2140 warning("variance-covariance matrix computed ", 2141 "from finite-difference Hessian is\n", 2142 "not positive definite or contains NA values: falling back to ", 2143 "var-cov estimated from RX") 2144 } 2145 } 2146 2147 ## FIXME: try to catch non-PD matrices 2148 rr <- tryCatch(as(V, "dpoMatrix"), error = function(e)e) 2149 if (inherits(rr, "error")) { 2150 warning(gettextf("Computed variance-covariance matrix problem: %s;\nreturning NA matrix", 2151 rr$message), domain = NA) 2152 rr <- matrix(NA,nrow(V),ncol(V)) 2153 } 2154 2155 nmsX <- colnames(object@pp$X) 2156 dimnames(rr) <- list(nmsX,nmsX) 2157 2158 if(correlation) 2159 rr@factors$correlation <- 2160 if(!is.na(sigm)) as(rr, "corMatrix") else rr # (is NA anyway) 2161 rr 2162} 2163 2164##' @importFrom stats vcov 2165##' @S3method vcov summary.merMod 2166vcov.summary.merMod <- function(object, ...) { 2167 if(is.null(object$vcov)) stop("logic error in summary of merMod object") 2168 object$vcov 2169} 2170 2171##' Make variance and correlation matrices from \code{theta} 2172##' 2173##' @param sc scale factor (residual standard deviation) 2174##' @param cnms component names 2175##' @param nc numeric vector: number of terms in each RE component 2176##' @param theta theta vector (lower-triangle of Cholesky factors) 2177##' @param nms component names (FIXME: nms/cnms redundant: nms=names(cnms)?) 2178##' @seealso \code{\link{VarCorr}} 2179##' @return A matrix 2180##' @export 2181mkVarCorr <- function(sc, cnms, nc, theta, nms) { 2182 ncseq <- seq_along(nc) 2183 thl <- split(theta, rep.int(ncseq, (nc * (nc + 1))/2)) 2184 if(!all(nms == names(cnms))) ## the above FIXME 2185 warning("nms != names(cnms) -- whereas lme4-authors thought they were --\n", 2186 "Please report!", immediate. = TRUE) 2187 ans <- lapply(ncseq, function(i) 2188 { 2189 ## Li := \Lambda_i, the i-th block diagonal of \Lambda(\theta) 2190 Li <- diag(nrow = nc[i]) 2191 Li[lower.tri(Li, diag = TRUE)] <- thl[[i]] 2192 rownames(Li) <- cnms[[i]] 2193 ## val := \Sigma_i = \sigma^2 \Lambda_i \Lambda_i', the 2194 val <- tcrossprod(sc * Li) # variance-covariance 2195 stddev <- sqrt(diag(val)) 2196 corr <- t(val / stddev)/stddev 2197 diag(corr) <- 1 2198 structure(val, stddev = stddev, correlation = corr) 2199 }) 2200 if(is.character(nms)) { 2201 ## FIXME: do we want this? Maybe not. 2202 ## Potential problem: the names of the elements of the VarCorr() list 2203 ## are not necessarily unique (e.g. fm2 from example("lmer") has *two* 2204 ## Subject terms, so the names are "Subject", "Subject". The print method 2205 ## for VarCorrs handles this just fine, but it's a little awkward if we 2206 ## want to dig out elements of the VarCorr list ... ??? 2207 if (anyDuplicated(nms)) 2208 nms <- make.names(nms, unique = TRUE) 2209 names(ans) <- nms 2210 } 2211 structure(ans, sc = sc) 2212} 2213 2214##' Extract variance and correlation components 2215##' 2216VarCorr.merMod <- function(x, sigma = 1, ...) 2217{ 2218 ## TODO: now that we have '...', add type=c("varcov","sdcorr","logs" ? 2219 if (is.null(cnms <- x@cnms)) 2220 stop("VarCorr methods require reTrms, not just reModule") 2221 if(missing(sigma)) 2222 sigma <- sigma(x) 2223 nc <- lengths(cnms) # no. of columns per term 2224 structure(mkVarCorr(sigma, cnms = cnms, nc = nc, theta = x@theta, 2225 nms = { fl <- x@flist; names(fl)[attr(fl, "assign")]}), 2226 useSc = as.logical(x@devcomp$dims[["useSc"]]), 2227 class = "VarCorr.merMod") 2228} 2229 2230if(FALSE)## *NOWHERE* used _FIXME_ ?? 2231## Compute standard errors of fixed effects from an merMod object 2232## 2233## @title Standard errors of fixed effects 2234## @param object "merMod" object, 2235## @param ... additional, optional arguments. None are used at present. 2236## @return numeric vector of length length(fixef(.)) 2237unscaledVar <- function(object, ...) { 2238 stopifnot(is(object, "merMod")) 2239 sigma(object) * diag(object@pp$unsc()) 2240} 2241 2242##' @S3method print VarCorr.merMod 2243print.VarCorr.merMod <- function(x, digits = max(3, getOption("digits") - 2), 2244 comp = "Std.Dev.", formatter = format, ...) { 2245 print(formatVC(x, digits = digits, comp = comp, formatter = formatter), quote = FALSE, ...) 2246 invisible(x) 2247} 2248 2249##' __NOT YET EXPORTED__ 2250##' "format()" the 'VarCorr' matrix of the random effects -- for 2251##' print()ing and show()ing 2252##' 2253##' @title Format the 'VarCorr' Matrix of Random Effects 2254##' @param varc a \code{\link{VarCorr}} (-like) matrix with attributes. 2255##' @param digits the number of significant digits. 2256##' @param comp character vector of length one or two indicating which 2257##' columns out of "Variance" and "Std.Dev." should be shown in the 2258##' formatted output. 2259##' @param formatter the \code{\link{function}} to be used for 2260##' formatting the standard deviations and or variances (but 2261##' \emph{not} the correlations which (currently) are always formatted 2262##' as "0.nnn" 2263##' @param ... optional arguments for \code{formatter(*)} in addition 2264##' to the first (numeric vector) and \code{digits}. 2265##' @return a character matrix of formatted VarCorr entries from \code{varc}. 2266formatVC <- function(varcor, digits = max(3, getOption("digits") - 2), 2267 comp = "Std.Dev.", formatter = format, 2268 useScale = attr(varcor, "useSc"), 2269 ...) 2270{ 2271 c.nms <- c("Groups", "Name", "Variance", "Std.Dev.") 2272 avail.c <- c.nms[-(1:2)] 2273 if(anyNA(mcc <- pmatch(comp, avail.c))) 2274 stop("Illegal 'comp': ", comp[is.na(mcc)]) 2275 nc <- length(colnms <- c(c.nms[1:2], (use.c <- avail.c[mcc]))) 2276 if(length(use.c) == 0) 2277 stop("Must *either* show variances or standard deviations") 2278 reStdDev <- c(lapply(varcor, attr, "stddev"), 2279 if(useScale) list(Residual = unname(attr(varcor, "sc")))) 2280 reLens <- lengths(reStdDev) 2281 nr <- sum(reLens) 2282 reMat <- array('', c(nr, nc), list(rep.int('', nr), colnms)) 2283 reMat[1+cumsum(reLens)-reLens, "Groups"] <- names(reLens) 2284 reMat[,"Name"] <- c(unlist(lapply(varcor, colnames)), if(useScale) "") 2285 if(any("Variance" == use.c)) 2286 reMat[,"Variance"] <- formatter(unlist(reStdDev)^2, digits = digits, ...) 2287 if(any("Std.Dev." == use.c)) 2288 reMat[,"Std.Dev."] <- formatter(unlist(reStdDev), digits = digits, ...) 2289 if (any(reLens > 1)) { 2290 maxlen <- max(reLens) 2291 recorr <- lapply(varcor, attr, "correlation") 2292 corr <- 2293 do.call(rbind, 2294 lapply(recorr, 2295 function(x) { 2296 x <- as.matrix(x) 2297 dig <- max(2, digits - 2) # use 'digits' ! 2298 ## not using formatter() for correlations 2299 cc <- format(round(x, dig), nsmall = dig) 2300 cc[!lower.tri(cc)] <- "" 2301 nr <- nrow(cc) 2302 if (nr >= maxlen) return(cc) 2303 cbind(cc, matrix("", nr, maxlen-nr)) 2304 }))[, -maxlen, drop = FALSE] 2305 if (nrow(corr) < nrow(reMat)) 2306 corr <- rbind(corr, matrix("", nrow(reMat) - nrow(corr), ncol(corr))) 2307 colnames(corr) <- c("Corr", rep.int("", max(0L, ncol(corr)-1L))) 2308 cbind(reMat, corr) 2309 } else reMat 2310} 2311 2312##' @S3method summary merMod 2313summary.merMod <- function(object, 2314 correlation = (p <= getOption("lme4.summary.cor.max")), 2315 use.hessian = NULL, 2316 ...) 2317{ 2318 if (...length() > 0) { 2319 ## FIXME: need testing code 2320 warning("additional arguments ignored") 2321 } 2322 ## se.calc: 2323 hess.avail <- (!is.null(h <- object@optinfo$derivs$Hessian) && 2324 nrow(h) > length(getME(object,"theta"))) 2325 if (is.null(use.hessian)) use.hessian <- hess.avail 2326 if (use.hessian && !hess.avail) 2327 stop("'use.hessian=TRUE' specified, but Hessian is unavailable") 2328 resp <- object@resp 2329 devC <- object@devcomp 2330 dd <- devC$dims 2331 ## cmp <- devC$cmp 2332 useSc <- as.logical(dd[["useSc"]]) 2333 sig <- sigma(object) 2334 ## REML <- isREML(object) 2335 2336 famL <- famlink(resp = resp) 2337 p <- length(coefs <- fixef(object)) 2338 2339 vc <- vcov(object, use.hessian = use.hessian) 2340 stdError <- sqrt(diag(vc)) 2341 coefs <- cbind("Estimate" = coefs, 2342 "Std. Error" = stdError) 2343 if (p > 0) { 2344 coefs <- cbind(coefs, (cf3 <- coefs[,1]/coefs[,2]), deparse.level = 0) 2345 colnames(coefs)[3] <- paste(if(useSc) "t" else "z", "value") 2346 if (isGLMM(object)) # FIXME: if "t" above, cannot have "z" here 2347 coefs <- cbind(coefs, "Pr(>|z|)" = 2348 2*pnorm(abs(cf3), lower.tail = FALSE)) 2349 } 2350 2351 llAIC <- llikAIC(object) 2352 ## FIXME: You can't count on object@re@flist, 2353 ## nor compute VarCorr() unless is(re, "reTrms"): 2354 varcor <- VarCorr(object) 2355 # use S3 class for now 2356 structure(list(methTitle = methTitle(dd), 2357 objClass = class(object), 2358 devcomp = devC, 2359 isLmer = is(resp, "lmerResp"), useScale = useSc, 2360 logLik = llAIC[["logLik"]], 2361 family = famL$family, link = famL$link, 2362 ngrps = ngrps(object), 2363 coefficients = coefs, sigma = sig, 2364 vcov = vcov(object, correlation = correlation, sigm = sig), 2365 varcor = varcor, # and use formatVC(.) for printing. 2366 AICtab = llAIC[["AICtab"]], call = object@call, 2367 residuals = residuals(object,"pearson",scaled = TRUE), 2368 fitMsgs = .merMod.msgs(object), 2369 optinfo = object@optinfo 2370 ), class = "summary.merMod") 2371} 2372 2373## TODO: refactor? 2374##' @S3method summary summary.merMod 2375summary.summary.merMod <- function(object, varcov = TRUE, ...) { 2376 if(varcov && is.null(object$vcov)) 2377 object$vcov <- vcov(object, correlation = TRUE, sigm = object$sigma) 2378 object 2379} 2380 2381### Plots for the ranef.mer class ---------------------------------------- 2382 2383##' @importFrom lattice dotplot 2384##' @S3method dotplot ranef.mer 2385dotplot.ranef.mer <- function(x, data, main = TRUE, transf=I, ...) 2386{ 2387 prepanel.ci <- function(x, y, se, subscripts, ...) { 2388 if (is.null(se)) return(list()) 2389 x <- as.numeric(x) 2390 hw <- 1.96 * as.numeric(se[subscripts]) 2391 list(xlim = range(transf(x - hw), transf(x + hw), finite = TRUE)) 2392 } 2393 panel.ci <- function(x, y, se, subscripts, pch = 16, 2394 horizontal = TRUE, col = dot.symbol$col, 2395 lty = dot.line$lty, lwd = dot.line$lwd, 2396 col.line = dot.line$col, levels.fos = unique(y), 2397 groups = NULL, ...) 2398 { 2399 x <- as.numeric(x) 2400 y <- as.numeric(y) 2401 dot.line <- trellis.par.get("dot.line") 2402 dot.symbol <- trellis.par.get("dot.symbol") 2403 sup.symbol <- trellis.par.get("superpose.symbol") 2404 panel.abline(h = levels.fos, col = col.line, lty = lty, lwd = lwd) 2405 panel.abline(v = 0, col = col.line, lty = lty, lwd = lwd) 2406 if (!is.null(se)) { 2407 se <- as.numeric(se[subscripts]) 2408 panel.segments( transf(x - 1.96 * se), y, 2409 transf(x + 1.96 * se), y, col = 'black') 2410 } 2411 panel.xyplot(transf(x), y, pch = pch, ...) 2412 } 2413 f <- function(nx, ...) { 2414 ss <- asDf0(x,nx) 2415 mtit <- if(main) nx # else NULL 2416 dotplot(.nn ~ values | ind, ss, se = ss$se, 2417 prepanel = prepanel.ci, panel = panel.ci, 2418 xlab = NULL, main = mtit, ...) 2419 } 2420 setNames(lapply(names(x), f, ...), names(x)) 2421} 2422 2423##' @importFrom graphics plot 2424##' @S3method plot ranef.mer 2425plot.ranef.mer <- function(x, y, ...) 2426{ 2427 lapply(x, function(x) { 2428 cn <- lapply(colnames(x), as.name) 2429 switch(min(ncol(x), 3), 2430 qqmath(eval(substitute(~ x, list(x = cn[[1]]))), x, ...), 2431 xyplot(eval(substitute(y ~ x, 2432 list(y = cn[[1]], 2433 x = cn[[2]]))), x, ...), 2434 splom(~ x, ...)) 2435 }) 2436} 2437 2438##' @importFrom lattice qqmath 2439##' @S3method qqmath ranef.mer 2440qqmath.ranef.mer <- function(x, data, main = TRUE, ...) 2441{ 2442 prepanel.ci <- function(x, y, se, subscripts, ...) { 2443 x <- as.numeric(x) 2444 se <- as.numeric(se[subscripts]) 2445 hw <- 1.96 * se 2446 list(xlim = range(x - hw, x + hw, finite = TRUE)) 2447 } 2448 panel.ci <- function(x, y, se, subscripts, pch = 16, ...) { 2449 panel.grid(h = -1,v = -1) 2450 panel.abline(v = 0) 2451 x <- as.numeric(x) 2452 y <- as.numeric(y) 2453 se <- as.numeric(se[subscripts]) 2454 panel.segments(x - 1.96 * se, y, x + 1.96 * se, y, col = 'black') 2455 panel.xyplot(x, y, pch = pch, ...) 2456 } 2457 f <- function(nx) { 2458 xt <- x[[nx]] 2459 mtit <- if(main) nx # else NULL 2460 if (!is.null(pv <- attr(xt, "postVar"))) 2461 { 2462 d <- dim(pv) 2463 se <- vapply(seq_len(d[1]), function(i) sqrt(pv[i, i, ]), numeric(d[3])) 2464 nr <- nrow(xt) 2465 nc <- ncol(xt) 2466 ord <- unlist(lapply(xt, order)) + rep((0:(nc - 1)) * nr, each = nr) 2467 rr <- 1:nr 2468 ind <- gl(nc, nr, labels = names(xt)) 2469 xyplot(rep(qnorm((rr - 0.5)/nr), nc) ~ unlist(xt)[ord] | ind[ord], 2470 se = se[ord], prepanel = prepanel.ci, panel = panel.ci, 2471 scales = list(x = list(relation = "free")), 2472 ylab = "Standard normal quantiles", 2473 xlab = NULL, main = mtit, ...) 2474 } else { 2475 qqmath(~values|ind, stack(xt), 2476 scales = list(y = list(relation = "free")), 2477 xlab = "Standard normal quantiles", 2478 ylab = NULL, main = mtit, ...) 2479 } 2480 } 2481 sapply(names(x), f, simplify = FALSE) 2482} 2483 2484##' @importFrom graphics plot 2485##' @S3method plot coef.mer 2486plot.coef.mer <- function(x, y, ...) 2487{ 2488 ## remove non-varying columns from frames 2489 reduced <- lapply(x, function(el) 2490 el[, !vapply(el, function(cc) all(cc == cc[1L]), NA)]) 2491 plot.ranef.mer(reduced, ...) 2492} 2493 2494##' @importFrom lattice dotplot 2495##' @S3method dotplot coef.mer 2496dotplot.coef.mer <- function(x, data, ...) { 2497 mc <- match.call() 2498 mc[[1]] <- as.name("dotplot.ranef.mer") 2499 eval(mc) 2500} 2501 2502##' @importFrom stats weights 2503##' @S3method weights merMod 2504weights.merMod <- function(object, type = c("prior","working"), ...) { 2505 type <- match.arg(type) 2506 isPrior <- type == "prior" 2507 if(!isGLMM(object) && !isPrior) 2508 stop("working weights only available for GLMMs") 2509 res <- if(isPrior) object@resp$weights else object@pp$Xwts^2 2510 ## the working weights available through pp$Xwts should be 2511 ## equivalent to: 2512 ## object@resp$weights*(object@resp$muEta()^2)/object@resp$variance() 2513 ## however, the unit tests in tests/glmmWeights.R suggest that this 2514 ## equivalence is approximate. this may be fine, however, if the 2515 ## discrepancy is due to another instance of the general problem of 2516 ## reference class fields not being updated at the optimum, then this 2517 ## could cause real problems. see for example: 2518 ## https://github.com/lme4/lme4/issues/166 2519 2520 2521 ## FIXME: what to do about missing values (see stats:::weights.glm)? 2522 ## FIXME: add unit tests 2523 return(res) 2524} 2525 2526## utility function: x is a ranef.mer object, nx is the name of an element 2527asDf0 <- function(x,nx,id=FALSE) { 2528 xt <- x[[nx]] 2529 ss <- stack(xt) 2530 ss$ind <- factor(as.character(ss$ind), levels = colnames(xt)) 2531 ss$.nn <- rep.int(reorder(factor(rownames(xt)), xt[[1]], 2532 FUN = mean,sort = sort), ncol(xt)) 2533 ## allow 'postVar' *or* 'condVar' names 2534 pv <- attr(xt,"postVar") 2535 if (is.null(pv)) { 2536 pv <- attr(xt,"condVar") 2537 } 2538 if (!is.null(pv)) { 2539 tmpfun <- function(pvi) { 2540 unlist(lapply(1:nrow(pvi), function(i) sqrt(pvi[i, i, ]))) 2541 } 2542 if (!is.list(pv)) { 2543 ss$se <- tmpfun(pv) 2544 } else { 2545 ## rely on ordering when unpacking! 2546 ss$se <- unlist(lapply(pv,tmpfun)) 2547 } 2548 } 2549 if (id) ss$id <- nx 2550 return(ss) 2551} 2552 2553## convert ranef object to a long-format data frame, e.g. suitable 2554## for ggplot2 (or homemade lattice plots) 2555## FIXME: have some gymnastics to do if terms, levels are different 2556## for different grouping variables - want to maintain ordering 2557## but still allow rbind()ing 2558as.data.frame.ranef.mer <- function(x, 2559 ..., 2560 stringsAsFactors = default.stringsAsFactors()) { 2561 xL <- lapply(names(x),asDf0,x=x,id=TRUE) 2562 ## combine 2563 xD <- do.call(rbind,xL) 2564 ## rename ... 2565 oldnames <- c("values","ind",".nn","se","id") 2566 newnames <- c("condval","term","grp","condsd","grpvar") 2567 names(xD) <- newnames[match(names(xD),oldnames)] 2568 ## reorder ... 2569 neworder <- c("grpvar","term","grp","condval") 2570 if ("condsd" %in% names(xD)) neworder <- c(neworder,"condsd") 2571 return(xD[neworder]) 2572} 2573 2574dim.merMod <- function(x) { 2575 getME(x, c("n", "p", "q", "p_i", "l_i", "q_i", "k", "m_i", "m")) 2576} 2577 2578## Internal utility, only used in optwrap() : 2579##' @title Get the optimizer function and check it minimally 2580##' @param optimizer character string ( = function name) *or* function 2581getOptfun <- function(optimizer) { 2582 if (((is.character(optimizer) && optimizer == "optimx") || 2583 deparse(substitute(optimizer)) == "optimx")) { 2584 if (!requireNamespace("optimx")) { 2585 stop(shQuote("optimx")," package must be installed order to ", 2586 "use ",shQuote('optimizer="optimx"')) 2587 } 2588 optfun <- optimx::optimx 2589 } else if (is.character(optimizer)) { 2590 optfun <- tryCatch(get(optimizer), error = function(e) NULL) 2591 } else optfun <- optimizer 2592 if (is.null(optfun)) stop("couldn't find optimizer function ",optimizer) 2593 if (!is.function(optfun)) stop("non-function specified as optimizer") 2594 needArgs <- c("fn","par","lower","control") 2595 if (anyNA(match(needArgs, names(formals(optfun))))) 2596 stop("optimizer function must use (at least) formal parameters ", 2597 paste(sQuote(needArgs), collapse = ", ")) 2598 optfun 2599} 2600 2601optwrap <- function(optimizer, fn, par, lower = -Inf, upper = Inf, 2602 control = list(), adj = FALSE, calc.derivs = TRUE, 2603 use.last.params = FALSE, 2604 verbose = 0L) 2605{ 2606 ## control must be specified if adj==TRUE; 2607 ## otherwise this is a fairly simple wrapper 2608 optfun <- getOptfun(optimizer) 2609 optName <- if(is.character(optimizer)) optimizer 2610 else ## "good try": 2611 deparse(substitute(optimizer))[[1L]] 2612 2613 lower <- rep_len(lower, length(par)) 2614 upper <- rep_len(upper, length(par)) 2615 2616 if (adj) 2617 ## control parameter tweaks: only for second round in nlmer, glmer 2618 switch(optName, 2619 "bobyqa" = { 2620 if(!is.numeric(control$rhobeg)) control$rhobeg <- 0.0002 2621 if(!is.numeric(control$rhoend)) control$rhoend <- 2e-7 2622 }, 2623 "Nelder_Mead" = { 2624 if (is.null(control$xst)) { 2625 thetaStep <- 0.1 2626 nTheta <- length(environment(fn)$pp$theta) 2627 betaSD <- sqrt(diag(environment(fn)$pp$unsc())) 2628 control$xst <- 0.2* c(rep.int(thetaStep, nTheta), 2629 pmin(betaSD, 10)) 2630 } 2631 if (is.null(control$xt)) control$xt <- control$xst*5e-4 2632 }) 2633 switch(optName, 2634 "bobyqa" = { 2635 if(all(par == 0)) par[] <- 0.001 ## minor kludge 2636 if(!is.numeric(control$iprint)) control$iprint <- min(verbose, 3L) 2637 }, 2638 "Nelder_Mead" = control$verbose <- verbose, 2639 "nloptwrap" = control$print_level <- min(as.numeric(verbose),3L), 2640 ## otherwise: 2641 if(verbose) warning(gettextf( 2642 "'verbose' not yet passed to optimizer '%s'; consider fixing optwrap()", 2643 optName), domain = NA) 2644 ) 2645 arglist <- list(fn = fn, par = par, lower = lower, upper = upper, control = control) 2646 ## optimx: must pass method in control (?) because 'method' was previously 2647 ## used in lme4 to specify REML vs ML 2648 if (optName == "optimx") { 2649 if (is.null(method <- control$method)) 2650 stop("must specify 'method' explicitly for optimx") 2651 arglist$control$method <- NULL 2652 arglist <- c(arglist, list(method = method)) 2653 } 2654 ## FIXME: test! effects of multiple warnings?? 2655 ## may not need to catch warnings after all?? 2656 curWarnings <- list() 2657 opt <- withCallingHandlers(do.call(optfun, arglist), 2658 warning = function(w) { 2659 curWarnings <<- append(curWarnings,list(w$message)) 2660 }) 2661 ## cat("***",unlist(tail(curWarnings,1))) 2662 ## FIXME: set code to warn on convergence !=0 2663 ## post-fit tweaking 2664 if (optName == "bobyqa") { 2665 opt$convergence <- opt$ierr 2666 } else if (optName == "Nelder_Mead") { 2667 ## fix-up: Nelder_Mead treats running out of iterations as "convergence" (!?) 2668 if (opt$NM.result==4) opt$convergence <- 4 2669 } else if (optName == "optimx") { 2670 opt <- list(par = coef(opt)[1,], 2671 fvalues = opt$value[1], 2672 method = method, 2673 conv = opt$convcode[1], 2674 feval = opt$fevals + opt$gevals, 2675 message = attr(opt,"details")[,"message"][[1]]) 2676 } 2677 if ((optconv <- getConv(opt)) != 0) { 2678 wmsg <- paste("convergence code",optconv,"from",optName) 2679 if (!is.null(getMsg(opt))) wmsg <- paste0(wmsg,": ",getMsg(opt)) 2680 warning(wmsg) 2681 curWarnings <<- append(curWarnings,list(wmsg)) 2682 } 2683 ## pp_before <- environment(fn)$pp 2684 ## save(pp_before,file="pp_before.RData") 2685 2686 if (calc.derivs) { 2687 if (use.last.params) { 2688 ## +0 tricks R into doing a deep copy ... 2689 ## otherwise element of ref class changes! 2690 ## FIXME:: clunky!! 2691 orig_pars <- opt$par 2692 orig_theta <- environment(fn)$pp$theta+0 2693 orig_pars[seq_along(orig_theta)] <- orig_theta 2694 } 2695 if (verbose > 10) cat("computing derivatives\n") 2696 derivs <- deriv12(fn,opt$par,fx = opt$value) 2697 if (use.last.params) { 2698 ## run one more evaluation of the function at the optimized 2699 ## value, to reset the internal/environment variables in devfun ... 2700 fn(orig_pars) 2701 } 2702 } else derivs <- NULL 2703 2704 if (!use.last.params) { 2705 ## run one more evaluation of the function at the optimized 2706 ## value, to reset the internal/environment variables in devfun ... 2707 fn(opt$par) 2708 } 2709 structure(opt, ## store all auxiliary information 2710 optimizer = optimizer, 2711 control = control, 2712 warnings = curWarnings, 2713 derivs = derivs) 2714} 2715 2716as.data.frame.VarCorr.merMod <- function(x,row.names = NULL, 2717 optional = FALSE, 2718 order = c("cov.last", "lower.tri"), 2719 ...) { 2720 order <- match.arg(order) 2721 tmpf <- function(v,grp) { 2722 vcov <- c(diag(v), v[lt.v <- lower.tri(v, diag = FALSE)]) 2723 sdcor <- c(attr(v,"stddev"), 2724 attr(v,"correlation")[lt.v]) 2725 nm <- rownames(v) 2726 n <- nrow(v) 2727 dd <- data.frame(grp = grp, 2728 var1 = nm[c(seq(n), col(v)[lt.v])], 2729 var2 = c(rep(NA,n), nm[row(v)[lt.v]]), 2730 vcov, 2731 sdcor, 2732 stringsAsFactors = FALSE) 2733 if (order=="lower.tri") { 2734 ## reorder *back* to lower.tri order 2735 m <- matrix(NA,n,n) 2736 diag(m) <- seq(n) 2737 m[lower.tri(m)] <- (n+1):(n*(n+1)/2) 2738 dd <- dd[m[lower.tri(m, diag=TRUE)],] 2739 } 2740 dd 2741 } 2742 r <- do.call(rbind, 2743 c(mapply(tmpf, x,names(x), SIMPLIFY = FALSE), 2744 deparse.level = 0)) 2745 if (attr(x,"useSc")) { 2746 ss <- attr(x,"sc") 2747 r <- rbind(r,data.frame(grp = "Residual",var1 = NA,var2 = NA, 2748 vcov = ss^2, 2749 sdcor = ss), 2750 deparse.level = 0) 2751 } 2752 rownames(r) <- NULL 2753 r 2754} 2755 2756