[Lme4-commits] r1669 - in pkg/lme4: R man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Mar 16 18:49:24 CET 2012
Author: mmaechler
Date: 2012-03-16 18:49:24 +0100 (Fri, 16 Mar 2012)
New Revision: 1669
Modified:
pkg/lme4/R/bootMer.R
pkg/lme4/R/lmer.R
pkg/lme4/R/profile.R
pkg/lme4/R/utilities.R
pkg/lme4/man/Dyestuff.Rd
pkg/lme4/tests/nbinom.R
Log:
using paste0() and other cosmetics
Modified: pkg/lme4/R/bootMer.R
===================================================================
--- pkg/lme4/R/bootMer.R 2012-03-16 17:49:01 UTC (rev 1668)
+++ pkg/lme4/R/bootMer.R 2012-03-16 17:49:24 UTC (rev 1669)
@@ -1,6 +1,5 @@
.simpleCap <- function(x) {
- paste(toupper(substring(x, 1,1)), substring(x, 2),
- sep="", collapse=" ")
+ paste0(toupper(substr(x, 1,1)), substr(x, 2, 1000000L), collapse=" ")
}
### bootMer() --- <==> (TODO: semi-*)parametric bootstrap
@@ -93,8 +92,8 @@
{
stopifnot((nsim <- as.integer(nsim[1])) > 0)
if (.progress!="none") { ## progress bar
- pbfun <- get(paste(.progress,"ProgressBar",sep=""))
- setpbfun <- get(paste("set",.simpleCap(.progress),"ProgressBar",sep=""))
+ pbfun <- get(paste0(.progress,"ProgressBar"))
+ setpbfun <- get(paste0("set",.simpleCap(.progress),"ProgressBar"))
pb <- do.call(pbfun,PBargs)
}
FUN <- match.fun(FUN)
@@ -150,4 +149,4 @@
as.data.frame.boot <- function(x,...) {
as.data.frame(x$t)
}
-
+
Modified: pkg/lme4/R/lmer.R
===================================================================
--- pkg/lme4/R/lmer.R 2012-03-16 17:49:01 UTC (rev 1668)
+++ pkg/lme4/R/lmer.R 2012-03-16 17:49:24 UTC (rev 1669)
@@ -221,7 +221,7 @@
##' (including \code{L-BFGS-B} from base \code{optim} and
##' \code{nlminb}), pass the \code{method} argument to \code{optim}
##' in the \code{control} argument.
-##'
+##'
##' @param mustart optional starting values on the scale of the conditional mean,
##' as in \code{\link[stats]{glm}}; see there for details.
##' @param etastart optional starting values on the scale of the unbounded
@@ -351,7 +351,7 @@
opt <- optwrap(optimizer[[1]],devfun,rho$pp$theta, rho$lower,
control=control, rho=rho,
adj=FALSE, verbose=verbose)
-
+
rho$nAGQ <- nAGQ
if (nAGQ > 0L) {
rho$lower <- c(rho$lower, rep.int(-Inf, length(rho$beta0)))
@@ -452,8 +452,8 @@
opt <- optwrap(optimizer[[2]],devfun, par=c(rho$pp$theta, rho$beta0),
lower=rho$lower, control=control, rho=rho,
adj=TRUE, verbose=verbose)
-
-
+
+
}
mkMerMod(environment(devfun), opt, vals$reTrms, vals$frame, mc)
}## {nlmer}
@@ -720,9 +720,8 @@
coefMer <- function(object, ...)
{
if (length(list(...)))
- warning(paste('arguments named "',
- paste(names(list(...)), collapse = ", "),
- '" ignored', sep = ''))
+ warning('arguments named "', paste(names(list(...)), collapse = ", "),
+ '" ignored')
fef <- data.frame(rbind(fixef(object)), check.names = FALSE)
ref <- ranef(object)
val <- lapply(ref, function(x)
@@ -1022,7 +1021,7 @@
rr <- object at resp$copy()
if (!is.null(newresp)) {
-
+
if (!is.null(na.act <- attr(object at frame,"na.action"))) {
## will only get here if na.action is 'na.omit' or 'na.exclude'
if (is.matrix(newresp)) {
@@ -1031,7 +1030,7 @@
}
-
+
if (isGLMM(object) && rr$family$family=="binomial") {
if (is.matrix(newresp) && ncol(newresp)==2) {
ntot <- rowSums(newresp)
@@ -1050,7 +1049,7 @@
rr$setResp(newresp)
}
-
+
pp <- object at pp$copy()
dc <- object at devcomp
nAGQ <- dc$dims["nAGQ"]
@@ -1923,7 +1922,7 @@
lower <- rep(lower, length.out=length(par))
upper <- rep(upper, length.out=length(par))
-
+
if (adj && is.character(optimizer))
switch(optimizer,
bobyqa = {
@@ -1973,7 +1972,7 @@
}
if (opt$conv!=0) {
wmsg <- paste("convergence code",opt$conv,"from",optimizer)
- if (!is.null(opt$msg)) wmsg <- paste(wmsg,": ",opt$msg,sep="")
+ if (!is.null(opt$msg)) wmsg <- paste0(wmsg,": ",opt$msg)
warning(wmsg)
}
opt
Modified: pkg/lme4/R/profile.R
===================================================================
--- pkg/lme4/R/profile.R 2012-03-16 17:49:01 UTC (rev 1668)
+++ pkg/lme4/R/profile.R 2012-03-16 17:49:24 UTC (rev 1669)
@@ -1,9 +1,9 @@
##' Methods for profile() of [ng]lmer fitted models
-##'
+##'
##' Methods for function \code{\link{profile}} (package \pkg{stats}), here for
##' profiling (fitted) mixed effect models.
-##'
-##'
+##'
+##'
##' @name profile-methods
##' @aliases profile-methods profile.merMod
##' @docType methods
@@ -48,7 +48,7 @@
## be careful with scale parameter;
## profile all parameters at once rather than RE first and then fixed)
-
+
dd <- devfun2(fitted)
base <- attr(dd, "basedev")
thopt <- attr(dd, "thopt")
@@ -57,10 +57,10 @@
X.orig <- pp$X
n <- environment(dd)$n
p <- length(pp$beta0)
-
+
ans <- lapply(opt <- attr(dd, "optimum"), function(el) NULL)
bakspl <- forspl <- ans
-
+
nptot <- length(opt)
nvp <- nptot - p # number of variance-covariance pars
fe.orig <- opt[-seq_len(nvp)]
@@ -69,11 +69,11 @@
dimnames = list(NULL, names(res)), byrow = TRUE)
## FIXME: why is cutoff based on nptot (i.e. boundary of simultaneous LRT conf region for nptot values)
## when we are computing (at most) 2-parameter profiles here?
-
+
cutoff <- sqrt(qchisq(1 - alphamax, nptot))
-
+
## helper functions
-
+
## nextpar calculates the next value of the parameter being
## profiled based on the desired step in the profile zeta
## (absstep) and the values of zeta and column cc for rows
@@ -100,7 +100,7 @@
prev=mat[r,1+seqnvp][-pind],
extrap=stop("stub")) ## do something with mat[r-(1:0),1+seqnvp])[-pind] ...
}
-
+
## mkpar generates the parameter vector of theta and
## sigma from the values being profiled in position w
mkpar <- function(np, w, pw, pmw) {
@@ -109,7 +109,7 @@
par[-w] <- pmw
par
}
-
+
## fillmat fills the third and subsequent rows of the matrix
## using nextpar and zeta
### FIXME: add code to evaluate more rows near the minimum if that
@@ -129,14 +129,14 @@
}
mat
}
-
+
lower <- c(pmax(fitted at lower,-1.0), 0, rep.int(-Inf, p))
upper <- c(ifelse(fitted at lower==0,Inf,1.0), Inf, rep.int(Inf, p))
seqnvp <- seq_len(nvp)
lowvp <- lower[seqnvp]
upvp <- upper[seqnvp]
form <- .zeta ~ foo # pattern for interpSpline formula
-
+
for (w in seqnvp) {
if (verbose) cat("var-cov parameter ",w,":\n",sep="")
@@ -156,13 +156,14 @@
## NA/NaN detected
warning("NAs detected in profiling")
} else if (devdiff < (-devtol))
- stop("profiling detected new, lower deviance")
- if (devdiff<0) warning(paste("slightly lower deviances (diff=",devdiff,") detected",sep=""))
+ 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)
c(zz, mkpar(nvp, w, xx, ores$par), pp$beta(1))
}
-
+
### FIXME: The starting values for the conditional optimization should
### be determined from recent starting values, not always the global
### optimum values.
@@ -170,7 +171,7 @@
### Can do this most easily by taking the most recent by taking the change in the other parameter values at
### the two last points and extrapolating.
-
+
## intermediate storage for pos. and neg. increments
pres <- nres <- res
## assign one row, determined by inc. sign, from a small shift
@@ -182,14 +183,14 @@
bres$.par <- names(opt)[w]
ans[[w]] <- bres[order(bres[, wp1]), ]
form[[3]] <- as.name(names(opt)[w])
-
+
## FIXME: test for nearly-vertical slopes here ...
bakspl[[w]] <- try(backSpline(forspl[[w]] <- interpSpline(form, bres)),silent=TRUE)
if (inherits(bakspl[[w]],"try-error")) {
warning("non-monotonic profile")
}
- }
+ } ## for(w in ..)
offset.orig <- fitted at resp$offset
for (j in seq_len(p)) {
@@ -238,9 +239,9 @@
form[[3]] <- as.name(thisnm)
bakspl[[thisnm]] <-
try(backSpline(forspl[[thisnm]] <- interpSpline(form, bres)),silent=TRUE)
- if (inherits(bakspl[[thisnm]],"try-error")) warning("non-monotonic profile")
- }
-
+ if (inherits(bakspl[[thisnm]],"try-error")) warning("non-monotonic profile")
+ } ## for(j in 1..p)
+
ans <- do.call(rbind, ans)
ans$.par <- factor(ans$.par)
attr(ans, "forward") <- forspl
@@ -309,7 +310,7 @@
## Assumption: all parameters, including the residual SD on SD-scale
sigma <- pars[np]
## .Call(lmer_Deviance, pp$ptr(), resp$ptr(), pars[-np]/sigma)
- ## convert from sdcor vector back to 'unscaled theta'
+ ## convert from sdcor vector back to 'unscaled theta'
thpars <- Sv_to_Cv(pars,n=sapply(fm at cnms,length),s=sigma)
.Call(lmer_Deviance, pp$ptr(), resp$ptr(), thpars)
sigsq <- sigma^2
@@ -322,17 +323,17 @@
{
stopifnot(is.numeric(pars), length(pars) == np)
thpars <- pars[seq(np)]
-
+
## FIXME: allow useSc (i.e. NLMMs)
sigma <- pars[np]
## .Call(lmer_Deviance, pp$ptr(), resp$ptr(), pars[-np]/sigma)
- ## convert from sdcor vector back to 'unscaled theta'
+ ## convert from sdcor vector back to 'unscaled theta'
thpars <- Sv_to_Cv(pars,n=sapply(fm at cnms,length))
stop("STUB!")
sigsq <- sigma^2
pp$ldL2() + (resp$wrss() + pp$sqrL(1))/sigsq + n * log(2 * pi * sigsq)
}
- }
+ }
attr(ans, "optimum") <- opt # w/ names()
attr(ans, "basedev") <- basedev
attr(ans, "thopt") <- pp$theta
@@ -411,7 +412,7 @@
lsegments(x, y, x, 0, ...)
lims <- current.panel.limits()$xlim
myspl <- spl[[panel.number()]]
- if (inherits(myspl,"try-error")) browser()
+ if (inherits(myspl,"try-error")) browser() ## << FIXME
krange <- range(myspl$knots)
pr <- predict(myspl,
seq(max(lims[1], krange[1]),
@@ -594,7 +595,7 @@
axis.line.alpha = axis.line$alpha,
axis.line.lty = axis.line$lty,
axis.line.lwd = axis.line$lwd,
- i, j,
+ i, j,
...)
{
n.var <- eval.parent(expression(n.var))
Modified: pkg/lme4/R/utilities.R
===================================================================
--- pkg/lme4/R/utilities.R 2012-03-16 17:49:01 UTC (rev 1668)
+++ pkg/lme4/R/utilities.R 2012-03-16 17:49:24 UTC (rev 1669)
@@ -1,3 +1,6 @@
+if(getRversion() < "2.15")
+ paste0 <- function(...) paste(..., sep = '')
+
### Utilities for parsing and manipulating mixed-model formulas
##' From the result of \code{\link{findbars}} applied to a model formula and
@@ -3,5 +6,5 @@
##' and the evaluation frame, create the model matrix, etc. associated with
##' random-effects terms. See the description of the returned value for a
-##' detailed list.
+##' detailed list.
##'
##' @title Create Z, Lambda, Lind, etc.
@@ -388,7 +391,7 @@
## @return NULL. The function is executed for its side effect.
chck1 <- function(expr) {
if ((le <- length(expr)) == 1) {
- if (is.numeric(expr) && expr == 1)
+ if (is.numeric(expr) && expr == 1)
stop("1 is not meaningful in a nonlinear model formula")
return()
} else
@@ -405,7 +408,7 @@
##' parameter matrix. If the formula is to be used for optimizing
##' designs, the \code{"resp"} part can be omitted.
##'
-##'
+##'
##' @title Manipulate a nonlinear model formula.
##' @param mc matched call from the calling function. Should have arguments named
##' \describe{
@@ -454,7 +457,7 @@
nlenv <- list2env(fr, parent=parent.frame(2L))
lapply(pnames, function(nm) nlenv[[nm]] <- rep.int(nlpars[[nm]], n))
respMod <- mkRespMod(fr, nlenv=nlenv, nlmod=nlmod)
-
+
chck1(meform <- form[[3L]])
pnameexpr <- parse(text=paste(pnames, collapse='+'))[[1]]
nb <- nobars(meform)
@@ -465,7 +468,7 @@
frE[[nm]] <- as.numeric(rep(nm == pnames, each = n))
X <- model.matrix(fe, frE)
rownames(X) <- NULL
-
+
reTrms <- mkReTrms(lapply(findbars(meform),
function(expr) {
expr[[2]] <- substitute(0+foo, list(foo=expr[[2]]))
Modified: pkg/lme4/man/Dyestuff.Rd
===================================================================
--- pkg/lme4/man/Dyestuff.Rd 2012-03-16 17:49:01 UTC (rev 1668)
+++ pkg/lme4/man/Dyestuff.Rd 2012-03-16 17:49:24 UTC (rev 1669)
@@ -60,10 +60,11 @@
Sys.getenv("R_PROFILE")
cat("R_LIBS:\\n"); (RL <- strsplit(Sys.getenv("R_LIBS"), ":")[[1]])
nRL <- normalizePath(RL)
- cat("and extra .libPaths():\\n")
+ cat("and extra(:= not in R_LIBS) .libPaths():\\n")
.libPaths()[is.na(match(.libPaths(), nRL))]
sessionInfo()
+ searchpaths()
pkgI <- function(pkgname) {
pd <- packageDescription(pkgname)
cat(sprintf("\%s -- built: \%s\\n\%*s -- dir : \%s\\n",
@@ -72,6 +73,7 @@
}
pkgI("Matrix")
pkgI("Rcpp")
+ ## 2012-03-12{MM}: fails with --as-cran
pkgI("RcppEigen")
pkgI("minqa")
pkgI("lme4")
Modified: pkg/lme4/tests/nbinom.R
===================================================================
--- pkg/lme4/tests/nbinom.R 2012-03-16 17:49:01 UTC (rev 1668)
+++ pkg/lme4/tests/nbinom.R 2012-03-16 17:49:24 UTC (rev 1669)
@@ -3,7 +3,7 @@
getNBdisp <- lme4:::getNBdisp
refitNB <- lme4:::refitNB
-simfun <- function(sd.u=1,NBtheta=0.5,
+simfun <- function(sd.u=1, NBtheta=0.5,
nblock=25,
fform=~x,
beta=c(1,2),
@@ -117,20 +117,19 @@
if(!(Sys.getenv("USER") %in% c("maechler")))
q("no")
-## "too slow" for regular testing -- 49 (MM at lynne: 33) seconds
-
+## "too slow" for regular testing -- 49 (MM at lynne: 33, then 26) seconds:
(t4 <- system.time(g4 <- glmer.nb(y~ Base*trt + Age + Visit + (Visit|subject),
data=epil2, verbose=TRUE)))
g4
(Lg4 <- logLik(g4))
attributes(Lg4) <- attributes(Lg4)[c("class","df","nobs")]
stopifnot(
- all.equal(getNBdisp(g4), glmmADMB_epil_vals$ theta, tol = 0.002)
- ,
- all.equal(fixef (g4), glmmADMB_epil_vals$ fixef, tol = 0.004)
- ,
- all.equal(logLik.m (g4), - glmmADMB_epil_vals$ NLL, tol = 0.0002)
- )
+ all.equal(getNBdisp(g4), glmmADMB_epil_vals$ theta, tol= 0.0022)# was 0.002
+ ,
+ all.equal(fixef (g4), glmmADMB_epil_vals$ fixef, tol= 0.004)
+ ,
+ all.equal(logLik.m (g4), - glmmADMB_epil_vals$ NLL, tol= 0.0002)
+ )
cat('Time elapsed: ', proc.time(),'\n') # for ``statistical reasons''
More information about the Lme4-commits
mailing list