[Lme4-commits] r1759 - in pkg/lme4: R tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed May 30 10:39:23 CEST 2012
Author: bbolker
Date: 2012-05-30 10:39:23 +0200 (Wed, 30 May 2012)
New Revision: 1759
Modified:
pkg/lme4/R/profile.R
pkg/lme4/R/vcconv.R
pkg/lme4/tests/profile.R
Log:
tweaked safe_chol; only pivot if absolutely necessary
xyplot.thpr more robust, skips bad profiles
more profile examples, skip bad ones for now
Modified: pkg/lme4/R/profile.R
===================================================================
--- pkg/lme4/R/profile.R 2012-05-29 13:04:48 UTC (rev 1758)
+++ pkg/lme4/R/profile.R 2012-05-30 08:39:23 UTC (rev 1759)
@@ -411,7 +411,10 @@
}
## extract only the y component from a prediction
-predy <- function(sp, vv) predict(sp, vv)$y
+predy <- function(sp, vv) {
+ if (inherits(sp,"try-error")) rep(NA,length(vv))
+ else predict(sp, vv)$y
+}
stripExpr <- function(ll, nms) {
stopifnot(inherits(ll, "list"), is.character(nms))
@@ -456,14 +459,15 @@
bspl <- attr(x, "backward")
zeta <- c(-rev(levels), 0, levels)
fr <- data.frame(zeta = rep.int(zeta, length(spl)),
- pval = unlist(lapply(bspl, predy, zeta)),
+ pval = unlist(lapply(bspl,predy,zeta)),
pnm = gl(length(spl), length(zeta), labels = names(spl)))
if (length(ind <- which(is.na(fr$pval)))) {
fr[ind, "zeta"] <- 0
for (i in ind)
### FIXME: Should check which bound has been violated, although it
### will almost always be the minimum.
- fr[i, "pval"] <- min(spl[[fr[i, "pnm"] ]]$knots)
+ if (!is.null(curspl <- spl[[fr[i, "pnm"] ]]))
+ fr[i, "pval"] <- min(curspl$knots)
}
ylab <- expression(zeta)
if (absVal) {
@@ -477,22 +481,26 @@
panel = function(x, y, ...)
{
panel.grid(h = -1, v = -1)
- lsegments(x, y, x, 0, ...)
- lims <- current.panel.limits()$xlim
myspl <- spl[[panel.number()]]
- if (inherits(myspl,"try-error")) browser() ## << FIXME
- krange <- range(myspl$knots)
- pr <- predict(myspl,
- seq(max(lims[1], krange[1]),
- min(lims[2], krange[2]), len = 101))
- if (absVal) {
- pr$y <- abs(pr$y)
- y[y == 0] <- NA
- lsegments(x, y, rev(x), y)
+ if (inherits(myspl,"try-error") ||
+ is.null(myspl)) {
+ warning(sprintf("bad profile for variable %d: skipped",panel.number()))
} else {
- panel.abline(h = 0, ...)
+ lsegments(x, y, x, 0, ...)
+ lims <- current.panel.limits()$xlim
+ krange <- range(myspl$knots)
+ pr <- predict(myspl,
+ seq(max(lims[1], krange[1]),
+ min(lims[2], krange[2]), len = 101))
+ if (absVal) {
+ pr$y <- abs(pr$y)
+ y[y == 0] <- NA
+ lsegments(x, y, rev(x), y)
+ } else {
+ panel.abline(h = 0, ...)
+ }
+ panel.lines(pr$x, pr$y)
}
- panel.lines(pr$x, pr$y)
}))
do.call(xyplot, stripExpr(ll, names(spl)))
}
@@ -518,25 +526,27 @@
ci
}
+## FIXME: make bootMer more robust; make profiling more robust;
+## more warnings; documentation ...
##' @importFrom stats confint
##' @S3method confint merMod
confint.merMod <- function(object, parm, method=c("profile","quadratic","bootMer"), level = 0.95, zeta, nsim=500, boot.type="perc", ...)
{
method <- match.arg(method)
if (method=="profile") {
- if (missing(parm)) {
+ if (!missing(parm)) {
pp <- profile(object,which=parm)
} else {
pp <- profile(object)
}
- return(confint(pp,parm=parm,level=level,zeta=zeta,...))
+ return(confint(pp,level=level,zeta=zeta,...))
} else if (method=="quadratic") {
stop("stub")
## only give confidence intervals on fixed effects?
## or warn??
} else if (method=="bootMer") {
message("Computing bootstrap confidence intervals ...")
- bb <- bootMer(object, fixef, nsim=nsum)
+ bb <- bootMer(object, fixef, nsim=nsim)
boot::boot.ci(bb,type=boot.type,conf=level)
}
}
Modified: pkg/lme4/R/vcconv.R
===================================================================
--- pkg/lme4/R/vcconv.R 2012-05-29 13:04:48 UTC (rev 1758)
+++ pkg/lme4/R/vcconv.R 2012-05-30 08:39:23 UTC (rev 1759)
@@ -1,7 +1,12 @@
## convert matrix list to concatenated vector of lower triangles
mlist2vec <- function(L) {
n <- sapply(L,nrow)
- r <- unlist(lapply(L,function(x) x[lower.tri(x,diag=TRUE)]))
+ ## allow for EITHER upper- or lower-triangular input
+ ff <- function(x) {
+ if (all(x[upper.tri(x)]==0)) t(x[lower.tri(x,diag=TRUE)])
+ else x[upper.tri(x,diag=TRUE)]
+ }
+ r <- unlist(lapply(L,ff))
attr(r,"clen") <- n
r
}
@@ -56,12 +61,15 @@
if (all(dmult(m,0)==0)) { ## diagonal
return(diag(sqrt(diag(m))))
}
+ ## attempt regular Chol. decomp
+ if (!inherits(try(cc <- chol(m),silent=TRUE),"try-error"))
+ return(cc)
+ ## ... pivot if necessary ...
cc <- suppressWarnings(chol(m,pivot=TRUE))
- cc[,order(attr(cc,"pivot"))]
- ## FIXME: pivot is here to deal with semidefinite cases: do we have
- ## to be more careful? test!
- ## does pivoting carry a performance cost?
- ## FIXME: we
+ oo <- order(attr(cc,"pivot"))
+ cc[,oo]
+ ## FIXME: pivot is here to deal with semidefinite cases,
+ ## but results might be returned in a strange format: TEST
}
Vv_to_Cv <- function(v,n=NULL,s=1) {
@@ -73,6 +81,7 @@
attr(r,"clen") <- get_clen(v,n)
r
}
+
Sv_to_Cv <- function(v,n=NULL,s=1) {
if (!missing(s)) {
v <- v[-length(v)]
Modified: pkg/lme4/tests/profile.R
===================================================================
--- pkg/lme4/tests/profile.R 2012-05-29 13:04:48 UTC (rev 1758)
+++ pkg/lme4/tests/profile.R 2012-05-30 08:39:23 UTC (rev 1759)
@@ -51,7 +51,10 @@
profile(gm1,which=3)
xyplot(pr4,layout=c(5,1),as.table=TRUE)
-splom(pr4)
+if (FALSE) {
+ ## FIXME! fails because of NAs
+ splom(pr4)
+}
nm1 <- nlmer(circumference ~ SSlogis(age, Asym, xmid, scal) ~ Asym|Tree,
Orange, start = c(Asym = 200, xmid = 725, scal = 350))
@@ -71,3 +74,12 @@
print(splom(pr3))
}
+
+## NOT RUN
+if (FALSE) {
+ data("Contraception",package="mlmRev")
+ fm2 <- glmer(use ~ urban+age+livch+(urban|district), Contraception, binomial)
+ load("pr5.RData")
+ pr5 <- profile(fm2,verbose=10)
+ xyplot(pr5)
+}
More information about the Lme4-commits
mailing list