1## --> ../man/profile-methods.Rd 2 3profnames <- function(object, signames=TRUE, 4 useSc=isLMM(object), prefix=c("sd","cor")) { 5 ntp <- length(getME(object,"theta")) 6 ## return 7 c(if(signames) sprintf(".sig%02d", seq(ntp)) 8 else tnames(object, old=FALSE, prefix=prefix), 9 if(useSc) if (signames) ".sigma" else "sigma") 10} 11 12 13##' @importFrom splines backSpline interpSpline periodicSpline 14##' @importFrom stats profile 15##' @method profile merMod 16##' @export 17profile.merMod <- function(fitted, 18 which=NULL, 19 alphamax = 0.01, 20 maxpts = 100, 21 delta = NULL, 22 delta.cutoff = 1/8, 23 verbose=0, devtol=1e-9, 24 maxmult = 10, 25 startmethod = "prev", 26 optimizer = NULL, 27 control = NULL, 28 signames = TRUE, 29 parallel = c("no", "multicore", "snow"), 30 ncpus = getOption("profile.ncpus", 1L), 31 cl = NULL, 32 prof.scale = c("sdcor","varcov"), 33 ...) 34{ 35 36 ## FIXME: allow choice of nextstep/nextstart algorithm? 37 ## FIXME: allow selection of individual variables to profile by name? 38 ## FIXME: allow for failure of bounds (non-pos-definite correlation matrices) when >1 cor parameter 39 40 prof.scale <- match.arg(prof.scale) 41 parallel <- match.arg(parallel) 42 do_parallel <- have_mc <- have_snow <- NULL # "-Wall" are set here: 43 eval(initialize.parallel)# (parallel, ncpus) 44 45 if (is.null(optimizer)) optimizer <- fitted@optinfo$optimizer 46 ## hack: doesn't work to set bobyqa parameters to *ending* values stored 47 ## in @optinfo$control 48 ignore.pars <- c("xst", "xt") 49 control.internal <- fitted@optinfo$control 50 if (length(ign <- which(names(control.internal) %in% ignore.pars)) > 0) 51 control.internal <- control.internal[-ign] 52 if (!is.null(control)) { 53 i <- names(control) 54 control.internal[[i]] <- control[[i]] 55 } 56 control <- control.internal 57 useSc <- isLMM(fitted) || isNLMM(fitted) 58 prof.prefix <- 59 switch(prof.scale, 60 "sdcor" = { 61 transfuns <- list(from.chol= Cv_to_Sv, 62 to.chol = Sv_to_Cv, 63 to.sd = identity) 64 c("sd", "cor") 65 }, 66 "varcov" = { 67 transfuns <- list(from.chol= Cv_to_Vv, 68 to.chol = Vv_to_Cv, 69 to.sd = sqrt) 70 c("var", "cov") 71 }, 72 stop("internal error, prof.scale=", prof.scale)) 73 dd <- devfun2(fitted, useSc=useSc, signames=signames, 74 transfuns=transfuns, prefix=prof.prefix, ...) 75 ## FIXME: figure out to what do here ... 76 if (isGLMM(fitted) && fitted@devcomp$dims[["useSc"]]) 77 stop("can't (yet) profile GLMMs with non-fixed scale parameters") 78 stopifnot(devtol >= 0) 79 base <- attr(dd, "basedev") 80 ## protect against accidental tampering by later devfun calls 81 thopt <- forceCopy(attr(dd, "thopt")) 82 stderr <- attr(dd, "stderr") 83 pp <- environment(dd)$pp 84 X.orig <- pp$X 85 n <- environment(dd)$n 86 p <- length(pp$beta0) 87 88 opt <- attr(dd, "optimum") 89 nptot <- length(opt) 90 nvp <- nptot - p # number of variance-covariance pars 91 wi.vp <- seq_len(nvp) 92 if(nvp > 0) fe.orig <- opt[- wi.vp] 93 94 which <- get.which(which, nvp, nptot, names(opt), verbose) 95 96 res <- c(.zeta = 0, opt) 97 res <- matrix(res, nrow = maxpts, ncol = length(res), 98 dimnames = list(NULL, names(res)), byrow = TRUE) 99 100 ## FIXME: why is cutoff based on nptot 101 ## (i.e. boundary of simultaneous LRT conf region for nptot values) 102 ## when we are computing (at most) 2-parameter profiles here? 103 cutoff <- sqrt(qchisq(1 - alphamax, nptot)) 104 105 if (is.null(delta)) 106 delta <- cutoff*delta.cutoff 107 108 ## helper functions 109 110 ## nextpar calculates the next value of the parameter being 111 ## profiled based on the desired step in the profile zeta 112 ## (absstep) and the values of zeta and column cc for rows 113 ## r-1 and r. The parameter may not be below lower (or above upper) 114 nextpar <- function(mat, cc, r, absstep, 115 lower = -Inf, upper = Inf, minstep=1e-6) { 116 rows <- r - (1:0) # previous two row numbers 117 pvals <- mat[rows, cc] 118 zeta <- mat[rows, ".zeta"] 119 num <- diff(pvals) 120 if (is.na(denom <- diff(zeta)) || denom==0) { 121 warning("Last two rows have identical or NA .zeta values: using minstep") 122 step <- minstep 123 } else { 124 step <- absstep*num/denom 125 if (step<0) { 126 warning("unexpected decrease in profile: using minstep") 127 step <- minstep 128 } else { 129 if (r>1) { 130 if (abs(step) > (maxstep <- abs(maxmult*num))) { 131 maxstep <- sign(step)*maxstep 132 if (verbose) cat(sprintf("capped step at %1.2f (multiplier=%1.2f > %1.2f)\n", 133 maxstep,abs(step/num),maxmult)) 134 step <- maxstep 135 } 136 } 137 } 138 } 139 min(upper, max(lower, pvals[2] + sign(num) * step)) 140 } 141 142 nextstart <- function(mat, pind, r, method) { 143 ## FIXME: indexing may need to be checked (e.g. for fixed-effect parameters) 144 switch(method, 145 global= opt[seqpar1][-pind], ## address opt, no zeta column 146 prev = mat[r,1+seqpar1][-pind], 147 extrap = stop("not yet implemented"),## do something with mat[r-(1:0),1+seqnvp])[-pind] 148 stop("invalid nextstart method")) 149 } 150 151 ## mkpar generates the parameter vector of theta and 152 ## sigma from the values being profiled in position w 153 mkpar <- function(np, w, pw, pmw) { 154 par <- numeric(np) 155 par[w] <- pw 156 par[-w] <- pmw 157 par 158 } 159 160 ## fillmat fills the third and subsequent rows of the matrix 161 ## using nextpar and zeta 162 ## FIXME: add code to evaluate more rows near the minimum if that 163 ## constraint was active. 164 fillmat <- function(mat, lowcut, upcut, zetafun, cc) { 165 nr <- nrow(mat) 166 i <- 2L 167 while (i < nr && mat[i, cc] > lowcut && mat[i,cc] < upcut && 168 (is.na(curzeta <- abs(mat[i, ".zeta"])) || curzeta <= cutoff)) { 169 np <- nextpar(mat, cc, i, delta, lowcut, upcut) 170 ns <- nextstart(mat, pind = cc-1, r=i, method=startmethod) 171 mat[i + 1L, ] <- zetafun(np,ns) 172 if (verbose>0) { 173 cat(i,cc,mat[i+1L,],"\n") 174 } 175 i <- i + 1L 176 } 177 if (mat[i-1,cc]==lowcut) { 178 ## fill in more values near the minimum 179 } 180 if (mat[i-1,cc]==upcut) { 181 ## fill in more values near the maximum 182 } 183 mat 184 } 185 186 ## bounds on Cholesky (== fitted@lower): [0,Inf) for diag, (-Inf,Inf) for off-diag 187 ## bounds on sd-corr: [0,Inf) for diag, (-1.0,1.0) for off-diag 188 ## bounds on var-corr: [0,Inf) for diag, (-Inf,Inf) for off-diag 189 if (prof.scale=="sdcor") { 190 lower <- pmax(fitted@lower, -1.) 191 upper <- ifelse(fitted@lower==0, Inf, 1.0) 192 } else { 193 lower <- fitted@lower 194 upper <- rep(Inf,length.out=length(fitted@lower)) 195 } 196 if (useSc) { # bounds for sigma 197 lower <- c(lower,0) 198 upper <- c(upper,Inf) 199 } 200 ## bounds on fixed parameters (TODO: allow user-specified bounds, e.g. for NLMMs) 201 lower <- c(lower,rep.int(-Inf, p)) 202 upper <- c(upper, rep.int(Inf, p)) 203 npar1 <- if (isLMM(fitted)) nvp else nptot 204 ## check that devfun2() computation for the base parameters is (approx.) the 205 ## same as the original devfun() computation 206 if(!isTRUE(all.equal(unname(dd(opt[seq(npar1)])), base, tolerance=1e-5))){ 207 stop("Profiling over both the residual variance and\n", 208 "fixed effects is not numerically consistent with\n", 209 "profiling over the fixed effects only")} 210 211 ## sequence of variance parameters to profile 212 seqnvp <- intersect(seq_len(npar1), which) 213 ## sequence of 'all' parameters 214 seqpar1 <- seq_len(npar1) 215 lowvp <- lower[seqpar1] 216 upvp <- upper[seqpar1] 217 form <- .zeta ~ foo # pattern for interpSpline formula 218 219 ## as in bootMer.R, define FUN as a 220 ## closure containing the referenced variables 221 ## in its scope to avoid explicit clusterExport statement 222 ## in the PSOCKcluster case 223 FUN <- local({ 224 function(w) { 225 if (verbose) cat(if(isLMM(fitted)) "var-cov " else "", 226 "parameter ",w,":\n",sep="") 227 wp1 <- w + 1L 228 start <- opt[seqpar1][-w] 229 pw <- opt[w] 230 lowcut <- lower[w] 231 upcut <- upper[w] 232 zeta <- function(xx,start) { 233 ores <- tryCatch(optwrap(optimizer, par=start, 234 fn=function(x) dd(mkpar(npar1, w, xx, x)), 235 lower = lowvp[-w], 236 upper = upvp [-w], 237 control = control), 238 error=function(e) NULL) 239 if (is.null(ores)) { 240 devdiff <- NA 241 pars <- NA 242 } else { 243 devdiff <- ores$fval - base 244 pars <- ores$par 245 } 246 if (is.na(devdiff)) { 247 warning("NAs detected in profiling") 248 } else { 249 if(verbose && devdiff < 0) 250 cat("old deviance ",base,",\n", 251 "new deviance ",ores$fval,",\n", 252 "new params ", 253 paste(mkpar(npar1,w,xx,ores$par), 254 collapse=","),"\n") 255 if (devdiff < (-devtol)) 256 stop("profiling detected new, lower deviance") 257 if(devdiff < 0) 258 warning(gettextf("slightly lower deviances (diff=%g) detected", 259 devdiff), domain=NA) 260 } 261 devdiff <- max(0,devdiff) 262 zz <- sign(xx - pw) * sqrt(devdiff) 263 r <- c(zz, mkpar(npar1, w, xx, pars)) 264 if (isLMM(fitted)) c(r, pp$beta(1)) else r 265 }## {zeta} 266 267 ## intermediate storage for pos. and neg. increments 268 pres <- nres <- res 269 ## assign one row, determined by inc. sign, from a small shift 270 ## FIXME:: do something if pw==0 ??? 271 shiftpar <- if (pw==0) 1e-3 else pw*1.01 272 ## Since both the pos- and neg-increment matrices are already 273 ## filled with the opt. par. results, this sets the first 274 ## two rows of the positive-increment matrix 275 ## to (opt. par, shiftpar) and the first two rows of 276 ## the negative-increment matrix to (shiftpar, opt. par), 277 ## which sets up two points going in the right direction 278 ## for each matrix (since the profiling algorithm uses increments 279 ## between rows to choose the next parameter increment) 280 nres[1, ] <- pres[2, ] <- zeta(shiftpar, start=opt[seqpar1][-w]) 281 ## fill in the rest of the arrays and collapse them 282 upperf <- fillmat(pres, lowcut, upcut, zeta, wp1) 283 lowerf <- 284 if (pw > lowcut) 285 fillmat(nres, lowcut, upcut, zeta, wp1) 286 else ## don't bother to fill in 'nres' matrix 287 nres 288 ## this will throw away the extra 'opt. par' and 'shiftpar' 289 ## rows introduced above: 290 bres <- as.data.frame(unique(rbind2(upperf,lowerf))) 291 pname <- names(opt)[w] 292 bres$.par <- pname 293 bres <- bres[order(bres[, wp1]), ] 294 295 ## FIXME: test for bad things here?? 296 form[[3]] <- as.name(pname) 297 forspl <- NULL # (in case of error) 298 ## bakspl 299 bakspl <- 300 tryCatch(backSpline( 301 forspl <- interpSpline(form, bres, na.action=na.omit)), 302 error=function(e)e) 303 if (inherits(bakspl, "error")) 304 warning("non-monotonic profile for ",pname) 305 ## return: 306 namedList(bres,bakspl,forspl) # namedList() -> lmerControl.R 307 308 }}) ## FUN() 309 310 ## copied from bootMer: DRY! 311 L <- if (do_parallel) { 312 if (have_mc) { 313 parallel::mclapply(seqnvp, FUN, mc.cores = ncpus) 314 } else if (have_snow) { 315 if (is.null(cl)) { 316 cl <- parallel::makePSOCKcluster(rep("localhost", ncpus)) 317 ## explicit export of the lme4 namespace since most FUNs will probably 318 ## use some of them 319 parallel::clusterExport(cl, varlist=getNamespaceExports("lme4")) 320 if(RNGkind()[1L] == "L'Ecuyer-CMRG") 321 parallel::clusterSetRNGStream(cl) 322 pres <- parallel::parLapply(cl, seqnvp, FUN) 323 parallel::stopCluster(cl) 324 pres 325 } else parallel::parLapply(cl, seqnvp, FUN) 326 } 327 } else lapply(seqnvp, FUN) 328 nn <- names(opt[seqnvp]) 329 ans <- setNames(lapply(L, `[[`, "bres"), nn) 330 bakspl <- setNames(lapply(L, `[[`,"bakspl"), nn) 331 forspl <- setNames(lapply(L, `[[`,"forspl"), nn) 332 333 ## profile fixed effects separately (for LMMs) 334 if (isLMM(fitted)) { 335 offset.orig <- fitted@resp$offset 336 fp <- seq_len(p) 337 fp <- fp[(fp+nvp) %in% which] 338 ## FIXME: parallelize this too ... 339 for (j in fp) { 340 if (verbose) cat("fixed-effect parameter ",j,":\n",sep="") 341 pres <- # intermediate results for pos. incr. 342 nres <- res # and negative increments 343 est <- opt[nvp + j] 344 std <- stderr[j] 345 Xw <- X.orig[, j, drop=TRUE] 346 Xdrop <- .modelMatrixDrop(X.orig, j) 347 pp1 <- new(class(pp), 348 X = Xdrop, 349 Zt = pp$Zt, 350 Lambdat = pp$Lambdat, 351 Lind = pp$Lind, 352 theta = pp$theta, 353 n = nrow(Xdrop)) 354### FIXME Change this to use the deep copy and setWeights, setOffset, etc. 355 rr <- new(Class=class(fitted@resp), y=fitted@resp$y) 356 rr$setWeights(fitted@resp$weights) 357 fe.zeta <- function(fw, start) { 358 ## (start parameter ignored) 359 rr$setOffset(Xw * fw + offset.orig) 360 rho <- list2env(list(pp=pp1, resp=rr), parent = parent.frame()) 361 ores <- optwrap(optimizer, par = thopt, fn = mkdevfun(rho, 0L), 362 lower = fitted@lower) 363 ## this optimization is done on the ORIGINAL 364 ## theta scale (i.e. not the sigma/corr scale) 365 ## upper=Inf for all cases 366 ## lower = pmax(fitted@lower, -1.0), 367 ## upper = 1/(fitted@lower != 0))## = ifelse(fitted@lower==0, Inf, 1.0) 368 fv <- ores$fval 369 sig <- sqrt((rr$wrss() + pp1$sqrL(1))/n) 370 c(sign(fw - est) * sqrt(fv - base), 371 Cv_to_Sv(ores$par, lengths(fitted@cnms), s=sig), 372 ## ores$par * sig, sig, 373 mkpar(p, j, fw, pp1$beta(1))) 374 } 375 nres[1, ] <- pres[2, ] <- fe.zeta(est + delta * std) 376 poff <- nvp + 1L + j 377 ## Workaround R bug [rbind2() is S4 generic; cannot catch warnings in its arg] 378 ## see lme4 GH issue #304 379 upperf <- fillmat(pres, -Inf, Inf, fe.zeta, poff) 380 lowerf <- fillmat(nres, -Inf, Inf, fe.zeta, poff) 381 bres <- as.data.frame(unique(rbind2(upperf, lowerf))) 382 bres$.par <- n.j <- names(fe.orig)[j] 383 ans[[n.j]] <- bres[order(bres[, poff]), ] 384 form[[3]] <- as.name(n.j) 385 bakspl[[n.j]] <- 386 tryCatch(backSpline(forspl[[n.j]] <- interpSpline(form, bres)), 387 error=function(e)e) 388 if (inherits(bakspl[[n.j]], "error")) 389 warning("non-monotonic profile for ", n.j) 390 } ## for(j in 1..p) 391 } ## if isLMM 392 393 ans <- do.call(rbind, ans) 394 row.names(ans) <- NULL ## TODO: rbind(*, make.row.names=FALSE) 395 ans$.par <- factor(ans$.par) 396 structure(ans, 397 forward = forspl, 398 backward = bakspl, 399 lower = lower[seqnvp], 400 upper = upper[seqnvp], 401 class = c("thpr", "data.frame"))# 'thpr': see ../man/profile-methods.Rd 402} ## profile.merMod 403 404##' Transform 'which' \in {parnames | integer | "beta_" | "theta_"} 405##' into integer indices 406##' @param which numeric or character vector 407##' @param nvp number of variance-covariance parameters 408##' @param nptot total number of parameters 409##' @param parnames vector of parameter names 410##' @param verbose print messages? 411##' @examples 412##' fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy) 413##' tn <- names(getME(fm1,"theta")) 414##' nn <- c(tn,names(fixef(fm1))) 415##' get.which("theta_",length(tn),length(nn),nn, verbose=TRUE) 416##' 417get.which <- function(which, nvp, nptot, parnames, verbose=FALSE) { 418 if (is.null(which)) 419 seq_len(nptot) 420 else if (is.character(which)) { 421 wi <- integer(); wh <- which 422 if(any(j <- wh == "theta_")) { 423 wi <- seq_len(nvp); wh <- wh[!j] } 424 if(any(j <- wh == "beta_") && nptot > nvp) { 425 wi <- c(wi, seq(nvp+1, nptot)); wh <- wh[!j] } 426 if(any(j <- parnames %in% wh)) { ## which containing param.names 427 wi <- sort(unique(c(wi, seq_len(nptot)[j]))) 428 } 429 if(verbose) message(gettextf("From original which = %s: new which <- %s", 430 deparse(which, nlines=1), deparse(wi, nlines=1)), 431 domain=NA) 432 if(length(wi) == 0) 433 warning(gettextf("Nothing selected by 'which=%s'", deparse(which)), 434 domain=NA) 435 wi 436 } else # stopifnot( .. numeric ..) 437 which 438} 439 440## This is a hack. The preferred approach is to write a 441## subset method for the ddenseModelMatrix and dsparseModelMatrix 442## classes 443.modelMatrixDrop <- function(mm, w) { 444 if (isS4(mm)) { 445 nX <- slotNames(X <- mm[, -w, drop = FALSE]) 446 do.call(new, 447 c(list(Class = class(mm), 448 assign = attr(mm,"assign")[-w], 449 contrasts = NULL 450 ## FIXME: where did the contrasts information go?? 451 ## mm@contrasts 452 ), 453 lapply(structure(nX, .Names=nX), 454 function(nm) slot(X, nm)))) 455 } else { 456 structure(mm[, -w, drop=FALSE], 457 assign = attr(mm, "assign")[-w]) 458 } 459} 460 461## The deviance is profiled with respect to the fixed-effects 462## parameters but not with respect to sigma. The other parameters 463## are on the standard deviation scale, not the theta scale. 464## 465## @title Return a function for evaluation of the deviance. 466## @param fm a fitted model of class merMod 467## @param useSc (logical) whether a scale parameter is used 468## @param \dots arguments passed to profnames (\code{signames=TRUE} 469## to use old-style .sigxx names, FALSE uses (sd_cor|xx); 470## also \code{prefix=c("sd","cor")}) 471## @return a function for evaluating the deviance in the extended 472## parameterization. This is profiled with respect to the 473## variance-covariance parameters (fixed-effects done separately). 474devfun2 <- function(fm, useSc = if(isLMM(fm)) TRUE else NA, 475 transfuns = list(from.chol = Cv_to_Sv, 476 to.chol = Sv_to_Cv, to.sd = identity), 477 ...) 478{ 479 ## FIXME: have to distinguish between 480 ## 'useSc' (GLMM: report profiled scale parameter) and 481 ## 'useSc' (NLMM/LMM: scale theta by sigma) 482 ## hasSc := GLMMuseSc <- fm@devcomp$dims["useSc"] 483 stopifnot(is(fm, "merMod")) 484 fm <- refitML(fm) 485 basedev <- -2*c(logLik(fm)) ## no longer deviance() 486 vlist <- lengths(fm@cnms) 487 ## "has scale" := isLMM or GLMM with scale parameter 488 hasSc <- as.logical(fm@devcomp$dims[["useSc"]]) 489 stdErr <- unname(coef(summary(fm))[,2]) 490 pp <- fm@pp$copy() 491 if (useSc) { 492 sig <- sigma(fm) ## only if hasSc is TRUE? 493 ## opt <- c(pp$theta*sig, sig) 494 opt <- transfuns$from.chol(pp$theta, n=vlist, s=sig) 495 } else { 496 opt <- transfuns$from.chol(pp$theta, n=vlist) 497 } 498 names(opt) <- profnames(fm, useSc=useSc, ...) 499 opt <- c(opt, fixef(fm)) 500 resp <- fm@resp$copy() 501 np <- length(pp$theta) 502 nf <- length(fixef(fm)) 503 if(hasSc) np <- np + 1L # was if(!isGLMM(fm)) np <- np + 1L 504 n <- nrow(pp$V) # use V, not X so it works with nlmer 505 if (isLMM(fm)) { # ==> hasSc 506 ans <- function(pars) 507 { 508 stopifnot(is.numeric(pars), length(pars) == np) 509 ## Assumption: we can translate the *last* parameter back 510 ## to sigma (SD) scale ... 511 sigma <- transfuns$to.sd(pars[np]) 512 ## .Call(lmer_Deviance, pp$ptr(), resp$ptr(), pars[-np]/sigma) 513 ## convert from sdcor vector back to 'unscaled theta' 514 thpars <- transfuns$to.chol(pars, n=vlist, s=sigma) 515 .Call(lmer_Deviance, pp$ptr(), resp$ptr(), thpars) 516 sigsq <- sigma^2 517 pp$ldL2() - ldW + (resp$wrss() + pp$sqrL(1))/sigsq + n * log(2 * pi * sigsq) 518 } 519 ldW <- sum(log(environment(ans)$resp$weights)) 520 assign("ldW", ldW, envir = environment(ans)) 521 } else { # GLMM *and* NLMMs 522 d0 <- getME(fm, "devfun") 523 ## from glmer: 524 ## rho <- new.env(parent=parent.env(environment())) 525 ## rho$pp <- do.call(merPredD$new, c(reTrms[c("Zt","theta","Lambdat","Lind")], n=nrow(X), list(X=X))) 526 ## rho$resp <- mkRespMod(fr, if(REML) p else 0L) 527 ans <- function(pars) 528 { 529 stopifnot(is.numeric(pars), length(pars) == np+nf) 530 thpars <- if(!useSc) 531 transfuns$to.chol(pars[seq(np)], n=vlist) 532 else 533 transfuns$to.chol(pars[seq(np)], n=vlist, s=pars[np]) 534 fixpars <- pars[-seq(np)] 535 d0(c(thpars,fixpars)) 536 } 537 } 538 attr(ans, "optimum") <- opt # w/ names() 539 attr(ans, "basedev") <- basedev 540 attr(ans, "thopt") <- pp$theta 541 attr(ans, "stderr") <- stdErr 542 class(ans) <- "devfun" 543 ans 544} 545 546## extract only the y component from a prediction 547predy <- function(sp, vv) { 548 if(inherits(sp,"error")) rep(NA_real_, length(vv)) 549 else predict(sp, vv)$y 550} 551 552stripExpr <- function(ll, nms) { 553 stopifnot(is.list(ll), is.character(nms)) 554 fLevs <- as.expression(nms) # variable names; now replacing log(sigma[.]) etc: 555 fLevs[nms == ".sigma"] <- expression(sigma) 556 fLevs[nms == ".lsigma"] <- expression(log(sigma)) 557 fLevs[nms == ".sigmasq"] <- expression(sigma^2) 558 sigNms <- grep("^\\.sig[0-9]+", nms) 559 lsigNms <- grep("^\\.lsig[0-9]+", nms) 560 sig2Nms <- grep("^\\.sigsq[0-9]+", nms) 561 ## the <n> in ".sig0<n>" and then in ".lsig0<n>": 562 sigsub <- as.integer(substring(nms[ sigNms], 5)) 563 lsigsub <- as.integer(substring(nms[lsigNms], 6)) 564 sig2sub <- as.integer(substring(nms[sig2Nms], 7)) 565 fLevs[ sigNms] <- lapply( sigsub, function(i) bquote( sigma[.(i)])) 566 fLevs[lsigNms] <- lapply(lsigsub, function(i) bquote(log(sigma[.(i)]))) 567 fLevs[sig2Nms] <- lapply(sig2sub, function(i) bquote( {sigma[.(i)]}^2)) 568 ## result of using { .. }^2 is easier to understand == == 569 levsExpr <- substitute(strip.custom(factor.levels=foo), list(foo=fLevs)) 570 snames <- c("strip", "strip.left") 571 if(!any(.in <- snames %in% names(ll))) { 572 ll$strip <- levsExpr 573 } else { 574 for(nm in snames[.in]) 575 if(is.logical(v <- ll[[nm]]) && v) 576 ll[[nm]] <- levsExpr 577 } 578 ll 579} 580 581panel.thpr <- function(x, y, spl, absVal, ...) 582{ 583 panel.grid(h = -1, v = -1) 584 myspl <- spl[[panel.number()]] 585 lsegments(x, y, x, 0, ...) 586 if (absVal) { 587 y[y == 0] <- NA 588 lsegments(x, y, rev(x), y) 589 } else { 590 panel.abline(h = 0, ...) 591 } 592 if (!is(myspl,"spline")) { 593 ## 'difficult' data 594 if (absVal) myspl$y <- abs(myspl$y) 595 panel.lines (myspl$x, myspl$y) 596 panel.points(myspl$x, myspl$y, pch="+") 597 warning(gettextf("bad profile for variable %d: using linear interpolation", 598 panel.number()), domain=NA) 599 } else { 600 lims <- current.panel.limits()$xlim 601 krange <- range(myspl$knots) 602 pr <- predict(myspl, 603 seq(max(lims[1], krange[1]), 604 min(lims[2], krange[2]), len = 101)) 605 if (absVal) pr$y <- abs(pr$y) 606 panel.lines(pr$x, pr$y) 607 } 608} 609 610## A lattice-based plot method for profile objects 611##' @importFrom lattice xyplot 612##' @S3method xyplot thpr 613xyplot.thpr <- 614 function (x, data = NULL, 615 levels = sqrt(qchisq(pmax.int(0, pmin.int(1, conf)), df = 1)), 616 conf = c(50, 80, 90, 95, 99)/100, 617 absVal = FALSE, scales = NULL, 618 which = 1:nptot, ...) 619{ 620 if(any(!is.finite(conf) | conf <= 0 | conf >= 1)) 621 stop("values of 'conf' must be strictly between 0 and 1") 622 stopifnot(1 <= (nptot <- length(nms <- levels(x[[".par"]])))) 623 ## FIXME: is this sufficiently reliable? 624 ## (include "sigma" in 'theta' parameters) 625 nvp <- length(grep("^(\\.sig[0-9]+|.sigma|sd_|cor_)", nms)) 626 which <- get.which(which, nvp, nptot, nms) 627 levels <- sort(levels[is.finite(levels) & levels > 0]) 628 spl <- attr(x, "forward") [which] 629 bspl <- attr(x, "backward")[which] 630 ## for parameters for which spline couldn't be computed, 631 ## replace the 'spl' element with the raw profile data 632 if(any(badSpl <- vapply(spl, is.null, NA))) { 633 spl[badSpl] <- lapply(which(badSpl), function(i) { 634 n <- names(badSpl)[i] 635 r <- x[x[[".par"]] == n, ] 636 data.frame(y = r[[".zeta"]], x = r[[n]]) 637 }) 638 bspl[badSpl] <- lapply(spl[badSpl], function(d) data.frame(x=d$y,y=d$x)) 639 ## FIXME: more efficient, not yet ok ? 640 ## ibad <- which(badSpl) 641 ## spl[ibad] <- lapply(names(ibad), function(n) { 642 ## r <- x[x[[".par"]]==n,] 643 ## data.frame(y = r[[".zeta"]], x = r[[n]]) 644 ## }) 645 ## bspl[ibad] <- lapply(spl[ibad], function(d) data.frame(x=d$y,y=d$x)) 646 } 647 zeta <- c(-rev(levels), 0, levels) 648 mypred <- function(bs, zeta) { ## use linear approximation if backspline doesn't work 649 if (inherits(bs,"spline")) 650 predy(bs, zeta) 651 else if(is.numeric(x <- bs$x) && is.numeric(y <- bs$y) && length(x) == length(y)) 652 approx(x, y, xout = zeta)$y 653 else 654 rep_len(NA, length(zeta)) 655 } 656 fr <- data.frame(zeta = rep.int(zeta, length(spl)), 657 pval = unlist(lapply(bspl, mypred, zeta)), 658 pnm = gl(length(spl), length(zeta), labels = names(spl))) 659 if (length(ind <- which(is.na(fr$pval)))) { 660 fr[ind, "zeta"] <- 0 661 for (i in ind) 662### FIXME: Should check which bound has been violated, although it 663### will almost always be the minimum. 664 if (inherits(curspl <- spl[[fr[i, "pnm"] ]], "spline")) { 665 fr[i, "pval"] <- min(curspl$knots) 666 } 667 } 668 ylab <- if (absVal) { 669 fr$zeta <- abs(fr$zeta) 670 expression("|" * zeta * "|") 671 } else 672 expression(zeta) 673 intscales <- list(x = list(relation = 'free')) 674 ## FIXME: is there something less clunky we can do here 675 ## that allows for all possible user inputs 676 ## (may want to (1) override x$relation (2) add elements to $x 677 ## (3) add elements to scales) 678 if (!is.null(scales)) { 679 if (!is.null(scales[["x"]])) { 680 if (!is.null(scales[["x"]]$relation)) { 681 intscales[["x"]]$relation <- scales[["x"]]$relation 682 scales[["x"]]$relation <- NULL 683 } 684 intscales[["x"]] <- c(intscales[["x"]],scales[["x"]]) 685 scales[["x"]] <- NULL 686 } 687 intscales <- c(intscales,scales) 688 } 689 ll <- c(list(...), 690 list(x = zeta ~ pval | pnm, data=fr, 691 scales = intscales, 692 ylab = ylab, xlab = NULL, panel=panel.thpr, 693 spl = spl, absVal = absVal)) 694 695 do.call(xyplot, stripExpr(ll, names(spl))) 696} 697 698## copy of stats:::format.perc (not exported, and ":::" being forbidden nowadays): 699format.perc <- function (probs, digits) { 700 paste(format(100 * probs, trim = TRUE, 701 scientific = FALSE, digits = digits), 702 "%") 703} 704 705##' confint() method for our profile() results 'thpr' 706##' @importFrom stats confint 707confint.thpr <- function(object, parm, level = 0.95, zeta, 708 ## tolerance for non-monotonic profiles 709 ## (raw values, not splines) 710 non.mono.tol=1e-2, 711 ...) 712{ 713 bak <- attr(object, "backward") 714 ## fallback strategy for old profiles that don't have a lower/upper 715 ## attribute saved ... 716 if (is.null(lower <- attr(object,"lower"))) 717 lower <- rep(NA,length(parm)) 718 if (is.null(upper <- attr(object,"upper"))) 719 upper <- rep(NA,length(parm)) 720 ## FIXME: work a little harder to add -Inf/Inf for fixed effect 721 ## parameters? (Should only matter for really messed-up profiles) 722 bnms <- names(bak) 723 parm <- if (missing(parm)) 724 bnms 725 else if(is.numeric(parm)) # e.g., when called from confint.merMod() 726 bnms[parm] 727 else if (length(parm <- as.character(parm)) == 1) { 728 if (parm == "theta_") 729 grep("^(sd_|cor_|.sig|sigma$)", bnms, value=TRUE) 730 else if (parm == "beta_") 731 grep("^(sd_|cor_|.sig|sigma$)", bnms, value=TRUE, invert=TRUE) 732 } else 733 intersect(parm, bnms) 734 cn <- 735 if (missing(zeta)) { 736 a <- (1 - level)/2 737 a <- c(a, 1 - a) 738 zeta <- qnorm(a) 739 format.perc(a, 3) 740 } ## else NULL 741 ci <- matrix(NA_real_, nrow=length(parm), ncol=2L, dimnames = list(parm,cn)) 742 for (i in seq_along(parm)) { 743 ## would like to build this machinery into predy, but 744 ## predy is used in many places and it's much harder to 745 ## tell in general whether an NA indicates a lower or an 746 ## upper bound ... 747 badprof <- FALSE 748 p <- rep(NA,2) 749 if (!inherits(b <- bak[[parm[i]]], "error")) { 750 p <- predy(b, zeta) 751 } else { 752 obj1 <- object[object$.par==parm[[i]],c(parm[[i]],".zeta")] 753 if (all(is.na(obj1[,2]))) { 754 badprof <- TRUE 755 warning("bad profile for ",parm[i]) 756 } else if (min(diff(obj1[,2])<(-non.mono.tol),na.rm=TRUE)) { 757 badprof <- TRUE 758 warning("non-monotonic profile for ",parm[i]) 759 } else { 760 warning("bad spline fit for ",parm[i],": falling back to linear interpolation") 761 p <- approxfun(obj1[,2],obj1[,1])(zeta) 762 } 763 } 764 if (!badprof) { 765 if (is.na(p[1])) p[1] <- lower[i] 766 if (is.na(p[2])) p[2] <- upper[i] 767 } 768 ci[i,] <- p 769 } 770 ci 771} 772 773## FIXME: make bootMer more robust; make profiling more robust; 774## more warnings; documentation ... 775 776##' Compute confidence intervals on the parameters of an lme4 fit 777##' @param object a fitted [ng]lmer model 778##' @param parm parameters (specified by integer position) 779##' @param level confidence level 780##' @param method for computing confidence intervals 781##' @param zeta likelihood cutoff 782##' (if not specified, computed from \code{level}: "profile" only) 783##' @param nsim number of simulations for parametric bootstrap intervals 784##' @param boot.type bootstrap confidence interval type 785##' @param quiet (logical) suppress messages about computationally intensive profiling? 786##' @param oldNames (logical) use old-style names for \code{method="profile"}? (See \code{signames} argument to \code{\link{profile}} 787##' @param \dots additional parameters to be passed to \code{\link{profile.merMod}} or \code{\link{bootMer}} 788##' @return a numeric table of confidence intervals 789confint.merMod <- function(object, parm, level = 0.95, 790 method = c("profile","Wald","boot"), 791 zeta, nsim=500, boot.type = c("perc","basic","norm"), 792 FUN = NULL, quiet=FALSE, oldNames=TRUE, ...) 793{ 794 method <- match.arg(method) 795 boot.type <- match.arg(boot.type) 796 ## 'parm' corresponds to 'which' in other contexts 797 if (method=="boot" && !is.null(FUN)) { 798 ## custom boot function, don't expand parameter names 799 } else { 800 ## "use scale" = GLMM with scale parameter *or* LMM .. 801 useSc <- as.logical(object@devcomp$dims[["useSc"]]) 802 vn <- profnames(object, oldNames, useSc=useSc) 803 an <- c(vn,names(fixef(object))) 804 parm <- if(missing(parm)) seq_along(an) 805 else 806 get.which(parm, nvp=length(vn), nptot=length(an), parnames=an) 807 if (!quiet && method %in% c("profile","boot")) { 808 mtype <- switch(method, profile="profile", boot="bootstrap") 809 message("Computing ",mtype," confidence intervals ...") 810 flush.console() 811 } 812 } 813 switch(method, 814 "profile" = { 815 pp <- profile(object, which=parm, signames=oldNames, ...) 816 ## confint.thpr() with missing(parm) using all names: 817 confint(pp, level=level, zeta=zeta) 818 }, 819 "Wald" = { 820 a <- (1 - level)/2 821 a <- c(a, 1 - a) 822 ci.vcov <- array(NA,dim = c(length(vn), 2L), 823 dimnames = list(vn, format.perc(a,3))) 824 ## copied with small changes from confint.default 825 cf <- fixef(object) ## coef() -> fixef() 826 pnames <- names(cf) 827 ## N.B. can't use sqrt(...)[parm] (diag() loses names) 828 ses <- sqrt(diag(vcov(object))) 829 ci.fixed <- array(cf + ses %o% qnorm(a), 830 dim = c(length(pnames), 2L), 831 dimnames = list(pnames, format.perc(a, 3))) 832 vnames <- tnames(object) 833 ci.all <- rbind(ci.vcov,ci.fixed) 834 ci.all[parm,,drop=FALSE] 835 }, 836 "boot" = { 837 bootFun <- function(x) { 838 th <- x@theta 839 nvec <- lengths(x@cnms) 840 scaleTh <- (isLMM(x) || isNLMM(x)) 841 ## FIXME: still ugly. Best cleanup via Cv_to_Sv ... 842 ss <- if (scaleTh) { ## scale variances by sigma and include it 843 Cv_to_Sv(th, n=nvec, s=sigma(x)) 844 } else if (useSc) { ## don't scale variances but do include sigma 845 c(Cv_to_Sv(th, n=nvec), sigma(x)) 846 } else { ## no scaling, no sigma 847 Cv_to_Sv(th, n=nvec) 848 } 849 c(setNames(ss, vn), fixef(x)) 850 } 851 if (is.null(FUN)) FUN <- bootFun 852 bb <- bootMer(object, FUN=FUN, nsim=nsim,...) 853 if (all(is.na(bb$t))) stop("*all* bootstrap runs failed!") 854 print.bootWarnings(bb, verbose=FALSE) 855 citab <- confint(bb,level=level,type=boot.type) 856 if (missing(parm)) { 857 ## only happens if we have custom boot method 858 if (is.null(parm <- rownames(citab))) { 859 parm <- seq(nrow(citab)) 860 } 861 } 862 citab[parm, , drop=FALSE] 863 }, 864 ## should never get here ... 865 stop("unknown confidence interval method")) 866} 867 868##' Convert x-cosine and y-cosine to average and difference. 869##' 870##' Convert the x-cosine and the y-cosine to an average and difference 871##' ensuring that the difference is positive by flipping signs if 872##' necessary 873##' @param xc x-cosine 874##' @param yc y-cosine 875ad <- function(xc, yc) 876{ 877 a <- (xc + yc)/2 878 d <- (xc - yc) 879 cbind(sign(d)* a, abs(d)) 880} 881 882##' convert d versus a (as an xyVector) and level to a matrix of taui and tauj 883##' @param xy an xyVector 884##' @param lev the level of the contour 885tauij <- function(xy, lev) lev * cos(xy$x + outer(xy$y/2, c(-1, 1))) 886 887##' @title safe arc-cosine 888##' @param x numeric vector argument 889##' @return acos(x) being careful of boundary conditions 890sacos <- function(x) acos(pmax.int(-0.999, pmin.int(0.999, x))) 891 892##' Generate a contour 893##' 894##' @title Generate a contour 895##' @param sij the arc-cosines of i on j 896##' @param sji the arc-cosines of j on i 897##' @param levels numeric vector of levels at which to interpolate 898##' @param nseg number of segments in the interpolated contour 899##' @return a list with components 900##' \item{tki}{the tau-scale predictions of i on j at the contour levels} 901##' \item{tkj}{the tau-scale predictions of j on i at the contour levels} 902##' \item{pts}{an array of dimension (length(levels), nseg, 2) containing the points on the contours} 903cont <- function(sij, sji, levels, nseg = 101) 904{ 905 ada <- array(0, c(length(levels), 2L, 4L)) 906 ada[, , 1] <- ad(0, sacos(predy(sij, levels)/levels)) 907 ada[, , 2] <- ad( sacos(predy(sji, levels)/levels), 0) 908 ada[, , 3] <- ad(pi, sacos(predy(sij, -levels)/levels)) 909 ada[, , 4] <- ad(sacos(predy(sji, -levels)/levels), pi) 910 pts <- array(0, c(length(levels), nseg + 1, 2)) 911 for (i in seq_along(levels)) 912 pts[i, ,] <- tauij(predict(periodicSpline(ada[i, 1, ], ada[i, 2, ]), 913 nseg = nseg), levels[i]) 914 levs <- c(-rev(levels), 0, levels) 915 list(tki = predict(sij, levs), tkj = predict(sji, levs), pts = pts) 916} 917 918## copied from lattice:::chooseFace 919chooseFace <- function (fontface = NULL, font = 1) 920{ 921 if (is.null(fontface)) 922 font 923 else fontface 924} 925 926 927##' Draws profile pairs plots. Contours are for the marginal 928##' two-dimensional regions (i.e. using df = 2). 929##' 930##' @title Profile pairs plot 931##' @param x the result of \code{\link{profile}} (or very similar structure) 932##' @param data unused - only for compatibility with generic 933##' @param levels the contour levels to be shown; usually derived from \code{conf} 934##' @param conf numeric vector of confidence levels to be shown as contours 935##' @param ... further arguments passed to \code{\link{splom}} 936##' @importFrom grid gpar viewport 937##' @importFrom lattice splom 938##' @method splom thpr 939##' @export 940splom.thpr <- function (x, data, 941 levels = sqrt(qchisq(pmax.int(0, pmin.int(1, conf)), 2)), 942 conf = c(50, 80, 90, 95, 99)/100, which = 1:nptot, 943 draw.lower = TRUE, draw.upper = TRUE, ...) 944{ 945 stopifnot(1 <= (nptot <- length(nms <- names(attr(x, "forward"))))) 946 singfit <- FALSE 947 for (i in grep("^(\\.sig[0-9]+|sd_)", names(x))) 948 singfit <- singfit || any(x[,".zeta"] == 0 & x[,i] == 0) 949 if (singfit) warning("splom is unreliable for singular fits") 950 951 nvp <- length(grep("^(\\.sig[0-9]+|.sigma|sd_|cor_)", nms)) 952 which <- get.which(which, nvp, nptot, nms) 953 if (length(which) == 1) 954 stop("can't draw a scatterplot matrix for a single variable") 955 956 mlev <- max(levels) 957 spl <- attr(x, "forward")[which] 958 frange <- sapply(spl, function(x) range(x$knots)) 959 bsp <- attr(x, "backward")[which] 960 x <- x[x[[".par"]] %in% nms[which],c(".zeta",nms[which],".par")] 961 ## brange <- sapply(bsp, function(x) range(x$knots)) 962 pfr <- do.call(cbind, lapply(bsp, predy, c(-mlev, mlev))) 963 pfr[1, ] <- pmax.int(pfr[1, ], frange[1, ], na.rm = TRUE) 964 pfr[2, ] <- pmin.int(pfr[2, ], frange[2, ], na.rm = TRUE) 965 nms <- names(spl) 966 ## Create data frame fr of par. vals in zeta coordinates 967 fr <- x[, -1] 968 for (nm in nms) fr[[nm]] <- predy(spl[[nm]], na.omit(fr[[nm]])) 969 fr1 <- fr[1, nms] 970 ## create a list of lists with the names of the parameters 971 traces <- lapply(fr1, function(el) lapply(fr1, function(el1) list())) 972 .par <- NULL ## => no R CMD check warning 973 for (j in seq_along(nms)[-1]) { 974 frj <- subset(fr, .par == nms[j]) 975 for (i in seq_len(j - 1L)) { 976 fri <- subset(fr, .par == nms[i]) 977 sij <- interpSpline(fri[ , i], fri[ , j]) 978 sji <- interpSpline(frj[ , j], frj[ , i]) 979 ll <- cont(sij, sji, levels) 980 traces[[j]][[i]] <- list(sij = sij, sji = sji, ll = ll) 981 } 982 } 983 if(draw.lower) ## panel function for lower triangle 984 lp <- function(x, y, groups, subscripts, i, j, ...) { 985 tr <- traces[[j]][[i]] 986 grid::pushViewport(viewport(xscale = c(-1.07, 1.07) * mlev, 987 yscale = c(-1.07, 1.07) * mlev)) 988 dd <- sapply(current.panel.limits(), diff)/50 989 psij <- predict(tr$sij) 990 ll <- tr$ll 991 ## now do the actual plotting 992 panel.grid(h = -1, v = -1) 993 llines(psij$y, psij$x, ...) 994 llines(predict(tr$sji), ...) 995 with(ll$tki, lsegments(y - dd[1], x, y + dd[1], x, ...)) 996 with(ll$tkj, lsegments(x, y - dd[2], x, y + dd[2], ...)) 997 for (k in seq_along(levels)) llines(ll$pts[k, , ], ...) 998 grid::popViewport(1) 999 } 1000 if(draw.upper) ## panel function for upper triangle 1001 up <- function(x, y, groups, subscripts, i, j, ...) { 1002 ## panels are transposed so reverse i and j 1003 jj <- i 1004 ii <- j 1005 tr <- traces[[jj]][[ii]] 1006 ll <- tr$ll 1007 pts <- ll$pts 1008 ## limits <- current.panel.limits() 1009 psij <- predict(tr$sij) 1010 psji <- predict(tr$sji) 1011 ## do the actual plotting 1012 panel.grid(h = -1, v = -1) 1013 llines(predy(bsp[[ii]], psij$x), predy(bsp[[jj]], psij$y), ...) 1014 llines(predy(bsp[[ii]], psji$y), predy(bsp[[jj]], psji$x), ...) 1015 for (k in seq_along(levels)) 1016 llines(predy(bsp[[ii]], pts[k, , 2]), 1017 predy(bsp[[jj]], pts[k, , 1]), ...) 1018 } 1019 dp <- function(x = NULL, # diagonal panel 1020 varname = NULL, limits, at = NULL, lab = NULL, 1021 draw = TRUE, 1022 1023 varname.col = add.text$col, 1024 varname.cex = add.text$cex, 1025 varname.lineheight = add.text$lineheight, 1026 varname.font = add.text$font, 1027 varname.fontfamily = add.text$fontfamily, 1028 varname.fontface = add.text$fontface, 1029 1030 axis.text.col = axis.text$col, 1031 axis.text.alpha = axis.text$alpha, 1032 axis.text.cex = axis.text$cex, 1033 axis.text.font = axis.text$font, 1034 axis.text.fontfamily = axis.text$fontfamily, 1035 axis.text.fontface = axis.text$fontface, 1036 1037 axis.line.col = axis.line$col, 1038 axis.line.alpha = axis.line$alpha, 1039 axis.line.lty = axis.line$lty, 1040 axis.line.lwd = axis.line$lwd, 1041 i, j, 1042 ...) 1043 { 1044 n.var <- eval.parent(expression(n.var)) 1045 add.text <- trellis.par.get("add.text") 1046 axis.line <- trellis.par.get("axis.line") 1047 axis.text <- trellis.par.get("axis.text") 1048 if (!is.null(varname)) 1049 grid::grid.text(varname, 1050 gp = 1051 gpar(col = varname.col, 1052 cex = varname.cex, 1053 lineheight = varname.lineheight, 1054 fontface = chooseFace(varname.fontface, varname.font), 1055 fontfamily = varname.fontfamily)) 1056 if (draw) 1057 { 1058 at <- pretty(limits) 1059 sides <- c("left", "top") 1060 if (j == 1) sides <- "top" 1061 if (j == n.var) sides <- "left" 1062 for (side in sides) 1063 panel.axis(side = side, 1064 at = at, 1065 labels = format(at, trim = TRUE), 1066 ticks = TRUE, 1067 check.overlap = TRUE, 1068 half = side == "top" && j > 1, 1069 1070 tck = 1, rot = 0, 1071 1072 text.col = axis.text.col, 1073 text.alpha = axis.text.alpha, 1074 text.cex = axis.text.cex, 1075 text.font = axis.text.font, 1076 text.fontfamily = axis.text.fontfamily, 1077 text.fontface = axis.text.fontface, 1078 1079 line.col = axis.line.col, 1080 line.alpha = axis.line.alpha, 1081 line.lty = axis.line.lty, 1082 line.lwd = axis.line.lwd) 1083 lims <- c(-1.07, 1.07) * mlev 1084 grid::pushViewport(viewport(xscale = lims, yscale = lims)) 1085 side <- if(j == 1) "right" else "bottom" 1086 which.half <- if(j == 1) "lower" else "upper" 1087 at <- pretty(lims) 1088 panel.axis(side = side, at = at, labels = format(at, trim = TRUE), 1089 ticks = TRUE, half = TRUE, which.half = which.half, 1090 tck = 1, rot = 0, 1091 1092 text.col = axis.text.col, 1093 text.alpha = axis.text.alpha, 1094 text.cex = axis.text.cex, 1095 text.font = axis.text.font, 1096 text.fontfamily = axis.text.fontfamily, 1097 text.fontface = axis.text.fontface, 1098 1099 line.col = axis.line.col, 1100 line.alpha = axis.line.alpha, 1101 line.lty = axis.line.lty, 1102 line.lwd = axis.line.lwd) 1103 grid::popViewport(1) 1104 } 1105 } 1106 1107 panel.blank <- function(...) {} 1108 splom(~ pfr, 1109 lower.panel = if(draw.lower) lp else panel.blank, 1110 upper.panel = if(draw.upper) up else panel.blank, 1111 diag.panel = dp, ...) 1112} 1113 1114## return an lmer profile like x with all the .sigNN parameters 1115## replaced by .lsigNN. The forward and backward splines for 1116## these parameters are recalculated. -> ../man/profile-methods.Rd 1117logProf <- function (x, base = exp(1), ranef=TRUE, 1118 sigIni = if(ranef) "sig" else "sigma") 1119{ 1120 stopifnot(inherits(x, "thpr")) 1121 cn <- colnames(x) 1122 sigP <- paste0("^\\.", sigIni) 1123 if (length(sigs <- grep(sigP, cn))) { 1124 repP <- sub("sig", ".lsig", sigIni) 1125 colnames(x) <- cn <- sub(sigP, repP, cn) 1126 levels(x[[".par"]]) <- sub(sigP, repP, levels(x[[".par"]])) 1127 names(attr(x, "backward")) <- 1128 names(attr(x, "forward")) <- 1129 sub(sigP, repP, names(attr(x, "forward"))) 1130 for (nm in cn[sigs]) { 1131 x[[nm]] <- log(x[[nm]], base = base) 1132 fr <- x[x[[".par"]] == nm & is.finite(x[[nm]]), TRUE, drop=FALSE] 1133 form <- eval(substitute(.zeta ~ nm, list(nm = as.name(nm)))) 1134 attr(x, "forward")[[nm]] <- isp <- interpSpline(form, fr) 1135 attr(x, "backward")[[nm]] <- backSpline(isp) 1136 } 1137 ## eliminate rows that produced non-finite logs 1138 x <- x[apply(is.finite(as.matrix(x[, sigs])), 1, all), , drop=FALSE] 1139 } 1140 x 1141} 1142## the log() method must have (x, base); no other arguments 1143log.thpr <- function (x, base = exp(1)) logProf(x, base=base) 1144 1145 1146 1147 1148##' Create an approximating density from a profile object -- called only from densityplot.thpr() 1149##' 1150##' @title Approximate densities from profiles 1151##' @param pr a profile object 1152##' @param npts number of points at which to evaluate the density 1153##' @param upper upper bound on cumulative for a cutoff 1154##' @return a data frame 1155dens <- function(pr, npts=201, upper=0.999) { 1156 npts <- as.integer(npts) 1157 spl <- attr(pr, "forward") 1158 bspl <- attr(pr, "backward") 1159 stopifnot(inherits(pr, "thpr"), npts > 0, 1160 is.numeric(upper), 0.5 < upper, upper < 1, 1161 identical((vNms <- names(spl)), names(bspl))) 1162 bad_vars <- character(0) 1163 zeta <- c(-1,1) * qnorm(upper) 1164 rng <- vector(mode="list", length=length(bspl)) 1165 names(rng) <- vNms 1166 dd <- rng # empty named list to be filled 1167 for (i in seq_along(bspl)) { 1168 if (inherits(bspl[[i]], "error")) { 1169 rng[[i]] <- rep(NA,npts) 1170 bad_vars <- c(bad_vars, vNms[i]) 1171 next 1172 } 1173 rng0 <- predy(bspl[[i]], zeta) 1174 if (is.na(rng0[1])) rng0[1] <- 0 1175 if (is.na(rng0[2])) { 1176 ## try harder to pick an upper bound 1177 for(upp in 1 - 10^seq(-4,-1, length=21)) { # <<-- TODO: should be related to 'upper' 1178 if(!is.na(rng0[2L] <- predy(bspl[[i]], qnorm(upp)))) 1179 break 1180 } 1181 if (is.na(rng0[2])) { 1182 warning("can't find an upper bound for the profile") 1183 bad_vars <- c(bad_vars, vNms[i]) 1184 next 1185 } 1186 } 1187 rng[[i]] <- seq(rng0[1], rng0[2], len=npts) 1188 } 1189 fr <- data.frame(pval = unlist(rng), 1190 pnm = gl(length(rng), npts, labels=vNms)) 1191 for (nm in vNms) { 1192 dd[[nm]] <- 1193 if (!inherits(spl[[nm]], "spline")) { 1194 rep(NA_real_, npts) 1195 } else { 1196 zz <- predy(spl[[nm]], rng[[nm]]) 1197 dnorm(zz) * predict(spl[[nm]], rng[[nm]], deriv = 1L)$y 1198 } 1199 } 1200 fr$density <- unlist(dd) 1201 if (length(bad_vars)) 1202 warning("unreliable profiles for some variables skipped: ", 1203 paste(bad_vars, collapse=", ")) 1204 fr 1205} 1206 1207##' Densityplot method for a mixed-effects model profile 1208## 1209##' @title densityplot from a mixed-effects profile 1210##' @param x a mixed-effects profile 1211##' @param data not used - for compatibility with generic 1212##' @param ... optional arguments to \code{\link[lattice]{densityplot}()} 1213##' from package \pkg{lattice}. 1214##' @return a density plot 1215##' @examples ## see example("profile.merMod") 1216##' @importFrom lattice densityplot 1217##' @method densityplot thpr 1218##' @export 1219densityplot.thpr <- function(x, data, npts=201, upper=0.999, ...) { 1220 dd <- dens(x, npts=npts, upper=upper) 1221 ll <- c(list(...), 1222 list(x=density ~ pval|pnm, 1223 data=dd, 1224 type=c("l","g"), 1225 scales=list(relation="free"), 1226 xlab=NULL)) 1227 do.call(xyplot, stripExpr(ll, names(attr(x, "forward")))) 1228} 1229 1230##' Transform a mixed-effects profile to the variance scale 1231varianceProf <- function(x, ranef=TRUE) { 1232 ## "parallel" to logProf() 1233 stopifnot(inherits(x, "thpr")) 1234 cn <- colnames(x) 1235 if(length(sigs <- grep(paste0("^\\.", if(ranef)"sig" else "sigma"), cn))) { 1236 ## s/sigma/sigmasq/ ; s/sig01/sig01sq/ etc 1237 sigP <- paste0("^(\\.sig", if(ranef) "(ma)?" else "ma", ")") 1238 repP <- "\\1sq" 1239 colnames(x) <- cn <- sub(sigP, repP, cn) 1240 levels(x[[".par"]]) <- sub(sigP, repP, levels(x[[".par"]])) 1241 names(attr(x, "backward")) <- 1242 names(attr(x, "forward")) <- 1243 sub(sigP, repP, names(attr(x, "forward"))) 1244 for (nm in cn[sigs]) { 1245 x[[nm]] <- x[[nm]]^2 1246 ## select rows (and currently drop extra attributes) 1247 fr <- x[x[[".par"]] == nm, TRUE, drop=FALSE] 1248 form <- eval(substitute(.zeta ~ nm, list(nm = as.name(nm)))) 1249 attr(x, "forward")[[nm]] <- isp <- interpSpline(form, fr) 1250 attr(x, "backward")[[nm]] <- backSpline(isp) 1251 } 1252 } 1253 x 1254} 1255 1256## convert profile to data frame, adding a .focal parameter to simplify lattice/ggplot plotting 1257##' @method as.data.frame thpr 1258##' @param x the result of \code{\link{profile}} (or very similar structure) 1259##' @export 1260##' @rdname profile-methods 1261as.data.frame.thpr <- function(x,...) { 1262 class(x) <- "data.frame" 1263 m <- as.matrix(x[,seq(ncol(x))-1]) ## omit .par 1264 x.p <- x[[".par"]] 1265 x[[".focal"]] <- m[cbind(seq(nrow(x)),match(x.p,names(x)))] 1266 x[[".par"]] <- factor(x.p, levels=unique(as.character(x.p))) ## restore order 1267 x 1268} 1269 1270 1271