[Lme4-commits] r1757 - in pkg/lme4: R src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue May 29 15:03:27 CEST 2012
Author: bbolker
Date: 2012-05-29 15:03:26 +0200 (Tue, 29 May 2012)
New Revision: 1757
Modified:
pkg/lme4/R/lmer.R
pkg/lme4/R/profile.R
pkg/lme4/src/external.cpp
Log:
(see previous commit message)
Modified: pkg/lme4/R/lmer.R
===================================================================
--- pkg/lme4/R/lmer.R 2012-05-29 13:02:46 UTC (rev 1756)
+++ pkg/lme4/R/lmer.R 2012-05-29 13:03:26 UTC (rev 1757)
@@ -500,6 +500,7 @@
##' points. A value of 0 indicates that both the fixed-effects parameters
##' and the random effects are optimized by the iteratively reweighted least
##' squares algorithm.
+##' @param verbose Logical: print verbose output?
##' @return A function of one numeric argument.
##' @seealso \code{\link{lmer}}, \code{\link{glmer}} and \code{\link{nlmer}}
##' @keywords models
@@ -508,7 +509,7 @@
##' (dd <- lmer(Yield ~ 1|Batch, Dyestuff, devFunOnly=TRUE))
##' dd(0.8)
##' minqa::bobyqa(1, dd, 0)
-mkdevfun <- function(rho, nAGQ=1L) {
+mkdevfun <- function(rho, nAGQ=1L, verbose=0) {
## FIXME: should nAGQ be automatically embedded in rho?
stopifnot(is.environment(rho), is(rho$resp, "lmResp"))
ff <- NULL
@@ -520,14 +521,14 @@
function(theta) {
resp$updateMu(lp0)
pp$setTheta(theta)
- pwrssUpdate(pp, resp, 1e-7, GHrule(0L), compDev)
+ pwrssUpdate(pp, resp, 1e-7, GHrule(0L), compDev, verbose)
}
else
function(pars) {
resp$updateMu(lp0)
pp$setTheta(as.double(pars[dpars])) # theta is first part of pars
resp$setOffset(baseOffset + pp$X %*% as.numeric(pars[-dpars]))
- pwrssUpdate(pp, resp, tolPwrss, GQmat, compDev, fac)
+ pwrssUpdate(pp, resp, tolPwrss, GQmat, compDev, fac, verbose)
}
} else if (is(rho$resp, "nlsResp")) {
if (nAGQ < 2L) {
@@ -592,8 +593,7 @@
resp$resDev() + pp$sqrL(1)
}
-glmerPwrssUpdate <- function(pp, resp, tol, GQmat, compDev=TRUE, grpFac=NULL) {
- verbose <- TRUE ## hack
+glmerPwrssUpdate <- function(pp, resp, tol, GQmat, compDev=TRUE, grpFac=NULL, verbose=0) {
nAGQ <- nrow(GQmat)
if (compDev) {
if (nAGQ < 2L)
@@ -602,19 +602,31 @@
}
oldpdev <- .Machine$double.xmax
uOnly <- nAGQ == 0L
+ i <- 0
repeat {
- oldu <- pp$delu
- if (uOnly) oldbeta <- pp$delb
+ ## oldu <- pp$delu
+ ## olddelb <- pp$delb
pdev <- RglmerWrkIter(pp, resp, uOnly=uOnly)
+ if (verbose>2) cat(i,": ",pdev,"\n",sep="")
## check convergence first so small increases don't trigger errors
if (abs((oldpdev - pdev) / pdev) < tol)
break
- if (pdev > oldpdev) {
- stop("PIRLS update failed")
- }
+ ## if (pdev > oldpdev) {
+ ## ## try step-halving
+ ## ## browser()
+ ## k <- 0
+ ## while (k < 10 && pdev > oldpdev) {
+ ## pp$setDelu((oldu + pp$delu)/2.)
+ ## if (!uOnly) pp$setDelb((olddelb + pp$delb)/2.)
+ ## pdev <- RglmerWrkIter(pp, resp, uOnly=uOnly)
+ ## k <- k+1
+ ## }
+ ## }
+ if (pdev>oldpdev) stop("PIRLS update failed")
oldpdev <- pdev
+ i <- i+1
}
- resp$Laplace(pp$ldL2(), 0., pp$sqrL(1))
+ resp$Laplace(pp$ldL2(), 0., pp$sqrL(1)) ## FIXME: should 0. be pp$ldRX2 ?
}
## create a deviance evaluation function that uses the sigma parameters
Modified: pkg/lme4/R/profile.R
===================================================================
--- pkg/lme4/R/profile.R 2012-05-29 13:02:46 UTC (rev 1756)
+++ pkg/lme4/R/profile.R 2012-05-29 13:03:26 UTC (rev 1757)
@@ -110,13 +110,13 @@
step <- minstep
} else {
step <- absstep*num/denom
- }
- if (r>1) {
- if (abs(step) > (maxstep <- abs(maxmult*num))) {
- maxstep <- sign(step)*maxstep
- if (verbose) cat(sprintf("capped step at %1.2f (multiplier=%1.2f > %1.2f)\n",
- maxstep,abs(step/num),maxmult))
- step <- maxstep
+ if (r>1) {
+ if (abs(step) > (maxstep <- abs(maxmult*num))) {
+ maxstep <- sign(step)*maxstep
+ if (verbose) cat(sprintf("capped step at %1.2f (multiplier=%1.2f > %1.2f)\n",
+ maxstep,abs(step/num),maxmult))
+ step <- maxstep
+ }
}
}
min(upper, max(lower, pvals[2] + sign(num) * step))
@@ -148,7 +148,7 @@
i <- 2L
while (i < nr && mat[i, cc] > lowcut && mat[i,cc] < upcut &&
(is.na(curzeta <- abs(mat[i, ".zeta"])) || curzeta <= cutoff)) {
- np <- nextpar(mat, cc, i, delta, lowcut, upcut, step)
+ np <- nextpar(mat, cc, i, delta, lowcut, upcut)
ns <- nextstart(mat, cc-1, i, startmethod)
mat[i + 1L, ] <- zetafun(np,ns)
if (verbose>0) {
@@ -195,20 +195,29 @@
lowcut <- lower[w]
upcut <- upper[w]
zeta <- function(xx,start) {
- ores <- optwrap(optimizer, par=start,
- fn=function(x) dd(mkpar(npar1, w, xx, x)),
- lower = lowvp[-w],
- upper = upvp[-w])
- devdiff <- ores$fval-base
- if (is.na(ores$fval)) {
+ ores <- try(optwrap(optimizer, par=start,
+ fn=function(x) dd(mkpar(npar1, w, xx, x)),
+ lower = lowvp[-w],
+ upper = upvp[-w]),silent=TRUE)
+ if (inherits(ores,"try-error"))
+ {
+ devdiff <- NA
+ pars <- NA
+ } else {
+ devdiff <- ores$fval-base
+ pars <- ores$par
+ }
+ if (is.na(devdiff)) {
warning("NAs detected in profiling")
- } else if (devdiff < (-devtol))
- stop("profiling detected new, lower deviance")
- if(devdiff<0)
- warning("slightly lower deviances (diff=",devdiff,") detected")
+ } else {
+ if (devdiff < (-devtol))
+ stop("profiling detected new, lower deviance")
+ if(devdiff<0)
+ warning("slightly lower deviances (diff=",devdiff,") detected")
+ }
devdiff <- max(0,devdiff)
zz <- sign(xx - pw) * sqrt(devdiff)
- r <- c(zz, mkpar(npar1, w, xx, ores$par))
+ r <- c(zz, mkpar(npar1, w, xx, pars))
if (isLMM(fitted)) r <- c(r,pp$beta(1))
r
}
@@ -227,16 +236,16 @@
## FIXME:: do something if pw==0 ???
nres[1, ] <- pres[2, ] <- zeta(pw * 1.01, start=opt[seqpar1][-w])
## fill in the rest of the arrays and collapse them
- bres <-
- as.data.frame(unique(rbind2(fillmat(pres,lowcut, upcut, zeta, wp1),
- fillmat(nres,lowcut, upcut, zeta, wp1))))
+ upperf <- fillmat(pres,lowcut, upcut, zeta, wp1)
+ lowerf <- fillmat(nres,lowcut, upcut, zeta, wp1)
+ bres <- as.data.frame(unique(rbind2(upperf,lowerf)))
pname <- names(opt)[w]
bres$.par <- pname
ans[[pname]] <- bres[order(bres[, wp1]), ]
form[[3]] <- as.name(pname)
## FIXME: test for bad things here??
- bakspl[[pname]] <- try(backSpline(forspl[[pname]] <- interpSpline(form, bres)),silent=TRUE)
+ bakspl[[pname]] <- try(backSpline(forspl[[pname]] <- interpSpline(form, bres,na.action=na.omit)),silent=TRUE)
if (inherits(bakspl[[pname]],"try-error")) {
warning("non-monotonic profile")
}
@@ -509,6 +518,29 @@
ci
}
+##' @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)) {
+ pp <- profile(object,which=parm)
+ } else {
+ pp <- profile(object)
+ }
+ return(confint(pp,parm=parm,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)
+ boot::boot.ci(bb,type=boot.type,conf=level)
+ }
+}
+
## Convert x-cosine and y-cosine to average and difference.
## Convert the x-cosine and the y-cosine to an average and difference
@@ -582,7 +614,7 @@
nms <- names(spl)
## Create data frame fr of par. vals in zeta coordinates
fr <- x[, -1]
- for (nm in nms) fr[[nm]] <- predy(spl[[nm]], fr[[nm]])
+ for (nm in nms) fr[[nm]] <- predy(spl[[nm]], na.omit(fr[[nm]]))
fr1 <- fr[1, nms]
## create a list of lists with the names of the parameters
traces <- lapply(fr1, function(el) lapply(fr1, function(el1) list()))
Modified: pkg/lme4/src/external.cpp
===================================================================
--- pkg/lme4/src/external.cpp 2012-05-29 13:02:46 UTC (rev 1756)
+++ pkg/lme4/src/external.cpp 2012-05-29 13:03:26 UTC (rev 1757)
@@ -561,6 +561,20 @@
XPtr<merPredD>(ptr)->setBeta0(as<MVec>(beta0));
END_RCPP;
}
+
+
+ SEXP merPredDsetDelu(SEXP ptr, SEXP delu) {
+ BEGIN_RCPP;
+ XPtr<merPredD>(ptr)->setDelu(as<MVec>(delu));
+ END_RCPP;
+ }
+
+
+ SEXP merPredDsetDelb(SEXP ptr, SEXP delb) {
+ BEGIN_RCPP;
+ XPtr<merPredD>(ptr)->setDelb(as<MVec>(delb));
+ END_RCPP;
+ }
// getters
SEXP merPredDCcNumer(SEXP ptr) {
BEGIN_RCPP;
@@ -891,6 +905,9 @@
CALLDEF(merPredDsetTheta, 2), // setters
CALLDEF(merPredDsetBeta0, 2),
+ CALLDEF(merPredDsetDelu, 2), // setters
+ CALLDEF(merPredDsetDelb, 2),
+
CALLDEF(merPredDCcNumer, 1), // getters
CALLDEF(merPredDL, 1),
CALLDEF(merPredDPvec, 1),
More information about the Lme4-commits
mailing list