diff -Nru lme4-1.1-8/debian/changelog lme4-1.1-10/debian/changelog --- lme4-1.1-8/debian/changelog 2015-10-26 17:11:23.000000000 +0000 +++ lme4-1.1-10/debian/changelog 2015-10-26 17:11:23.000000000 +0000 @@ -1,3 +1,17 @@ +lme4 (1.1-10-1) unstable; urgency=low + + * New upstream release + + -- Dirk Eddelbuettel Tue, 06 Oct 2015 22:33:37 -0500 + +lme4 (1.1-9-1) unstable; urgency=low + + * New upstream release + + * debian/control: Set Build-Depends: to current R version + + -- Dirk Eddelbuettel Thu, 20 Aug 2015 18:48:03 -0500 + lme4 (1.1-8-1) unstable; urgency=low * New upstream release diff -Nru lme4-1.1-8/debian/control lme4-1.1-10/debian/control --- lme4-1.1-8/debian/control 2015-10-26 17:11:23.000000000 +0000 +++ lme4-1.1-10/debian/control 2015-10-26 17:11:23.000000000 +0000 @@ -2,7 +2,7 @@ Section: gnu-r Priority: optional Maintainer: Dirk Eddelbuettel -Build-Depends: debhelper (>= 7.0.0), cdbs, r-base-dev (>= 3.2.1), r-cran-matrix (>= 1.1-0), r-cran-lattice, r-cran-nlme, r-cran-mass, r-cran-rcppeigen (>= 0.3.2.0.2-2), r-cran-minqa (>= 1.2.2-2), r-cran-nloptr +Build-Depends: debhelper (>= 7.0.0), cdbs, r-base-dev (>= 3.2.2), r-cran-matrix (>= 1.1-0), r-cran-lattice, r-cran-nlme, r-cran-mass, r-cran-rcppeigen (>= 0.3.2.0.2-2), r-cran-minqa (>= 1.2.2-2), r-cran-nloptr Standards-Version: 3.9.6 Package: r-cran-lme4 diff -Nru lme4-1.1-8/DESCRIPTION lme4-1.1-10/DESCRIPTION --- lme4-1.1-8/DESCRIPTION 2015-06-22 07:09:58.000000000 +0000 +++ lme4-1.1-10/DESCRIPTION 2015-10-06 06:08:29.000000000 +0000 @@ -1,5 +1,5 @@ Package: lme4 -Version: 1.1-8 +Version: 1.1-10 Title: Linear Mixed-Effects Models using 'Eigen' and S4 Authors@R: c(person("Douglas","Bates", role="aut"), person("Martin","Maechler", role="aut"), @@ -23,7 +23,7 @@ Depends: R (>= 3.0.2), Matrix (>= 1.1.1), methods, stats LinkingTo: Rcpp (>= 0.10.5), RcppEigen Imports: graphics, grid, splines, utils, parallel, MASS, nlme, lattice, - minqa (>= 1.1.15), nloptr + minqa (>= 1.1.15), nloptr (>= 1.0.4) Suggests: knitr, boot, PKPDmodels, MEMSS, testthat (>= 0.8.1), ggplot2, mlmRev, optimx (>= 2013.8.6), gamm4, pbkrtest, HSAUR2, numDeriv VignetteBuilder: knitr @@ -32,6 +32,6 @@ URL: https://github.com/lme4/lme4/ http://lme4.r-forge.r-project.org/ BugReports: https://github.com/lme4/lme4/issues NeedsCompilation: yes -Packaged: 2015-06-20 21:23:09 UTC; bolker +Packaged: 2015-10-05 17:56:37 UTC; bolker Repository: CRAN -Date/Publication: 2015-06-22 09:09:58 +Date/Publication: 2015-10-06 08:08:29 diff -Nru lme4-1.1-8/inst/CITATION lme4-1.1-10/inst/CITATION --- lme4-1.1-8/inst/CITATION 2015-01-20 12:47:19.000000000 +0000 +++ lme4-1.1-10/inst/CITATION 2015-10-05 14:17:31.000000000 +0000 @@ -1,23 +1,25 @@ -year <- sub("-.*", "", meta$Date) -note <- sprintf("R package version %s", meta$Version) - -bibentry(bibtype = "Manual", - title = "{lme4}: Linear mixed-effects models using {Eigen} and {S4}", - author = c(person("Douglas", "Bates"), - person("Martin", "Maechler"), - person("Ben", "Bolker"), - person("Steven", "Walker")), - year = year, - note = note, - url = "http://CRAN.R-project.org/package=lme4") - -bibentry(bibtype = "Misc", - title = "Fitting Linear Mixed-Effects Models using {lme4}", - year = "2015", - author = c(person("Douglas", "Bates"), - person("Martin", "Maechler"), - person(c("Benjamin", "M."), "Bolker"), - person("Steven", "Walker")), - note = "ArXiv e-print; in press, \\emph{Journal of Statistical Software}", - url = "http://arxiv.org/abs/1406.5823") - +bibentry(bibtype = "Article", + title = "Fitting Linear Mixed-Effects Models Using {lme4}", + author = c(person(given = "Douglas", + family = "Bates"), + person(given = "Martin", + family = "M{\\\"a}chler"), + person(given = "Ben", + family = "Bolker"), + person(given = "Steve", + family = "Walker")), + journal = "Journal of Statistical Software", + year = "2015", + volume = "67", + number = "1", + pages = "1--48", + doi = "10.18637/jss.v067.i01", + + header = "To cite lme4 in publications use:", + textVersion = + paste("Douglas Bates, Martin Maechler, Ben Bolker, Steve Walker (2015).", + "Fitting Linear Mixed-Effects Models Using lme4.", + "Journal of Statistical Software, 67(1), 1-48.", + "doi:10.18637/jss.v067.i01.") +) + diff -Nru lme4-1.1-8/inst/NEWS.Rd lme4-1.1-10/inst/NEWS.Rd --- lme4-1.1-8/inst/NEWS.Rd 2015-06-13 23:08:09.000000000 +0000 +++ lme4-1.1-10/inst/NEWS.Rd 2015-10-05 14:19:14.000000000 +0000 @@ -3,7 +3,76 @@ \name{NEWS} \title{lme4 News} \encoding{UTF-8} -\section{CHANGES IN VERSION 1.1-8}{ + +\section{CHANGES IN VERSION 1.1-10 (2015-10-05)}{ + This release is primarily a version bump for the release of + the paper in J. Stat. Software. + \subsection{USER-VISIBLE CHANGES}{ + \itemize{ + \item updated CITATION file. + } + } + \subsection{NEW FEATURES}{ + \itemize{ + \item We export set of about a dozen printing utility functions + which are used in our \code{print} methods. + \item \code{bootMer} now allows the use of \code{re.form}. + } + } + \subsection{BUG FIXES}{ + \itemize{ + \item fixed reordering bug in names of \code{getME(.,"Ztlist")} + (terms are reordered in decreasing order of the number of + levels of the grouping variable, but names were not being + reordered) + \item fixed issue with simulation when complex forms (such + as nested random effects terms) are included in the model + (Github #335) + } + } +} +\section{CHANGES IN VERSION 1.1-9 (2015-08-20)}{ + \subsection{USER-VISIBLE CHANGES}{ + \itemize{ + \item explicit \code{maxit} arguments for various functions + (\code{refit}, \code{mkGlmerDevfun}, ...) + } + } + \subsection{NEW FEATURES}{ + \itemize{ + \item \code{terms} and \code{formula} methods now have + \code{random.only} options + + \item \code{getME} gains a \code{glmer.nb.theta} option. + It is now (an S3) generic with an \code{"merMod"} method in + \pkg{lme4} and potentially other methods in dependent packages. + + \item \code{simulate} now works for \code{glmer.nb} models + (Github #284: idea from @aosmith16) + } + } + \subsection{BUG FIXES}{ + \itemize{ + \item prediction and simulation now work when random-effects + terms have data-dependent bases (e.g., \code{poly(.)} or + \code{ns(.)} terms) (Github #313, Edgar Gonzalez) + \item \code{logLik} for \code{glmer.nb} models now + includes the overdispersion parameter in the + parameter count (\code{df} attribute) + \item \code{lmList} handles offsets and weights better + \item lots of fixes to \code{glmer.nb} (Github #176, #266, #287, + #318). \strong{Please note that glmer.nb is still somewhat + unstable/under construction.} + } + } + \subsection{CRAN-COMPATIBILITY UPDATES}{ + \itemize{ + \item import functions from base packages to pass CRAN checks + \item tweak to failing tests on Windows + } + } +} +\section{CHANGES IN VERSION 1.1-8 (2015-06-22)}{ \subsection{NEW FEATURES}{ \itemize{ \item \code{getME} gains a \code{"Tlist"} option @@ -71,13 +140,12 @@ \item As with \code{lm()}, users now get an error for non-finite (\code{Inf}, \code{NA}, or \code{NaN}) values in the response unless \code{na.action} is set to - exclude or omit them. Github #310 - + exclude or omit them (Github #310) } } } -\section{CHANGES IN VERSION 1.1-7}{ +\section{CHANGES IN VERSION 1.1-7 (2014-07-19)}{ \subsection{NEW FEATURES}{ \itemize{ \item the \code{nloptr} package is now imported; @@ -98,11 +166,11 @@ \item the use of \code{deviance} to return the REML criterion is now deprecated; users should use \code{REMLcrit()} instead (Github #211) - \item changed the default value of \code{check.nobs.vs.rankZ} to - \code{"ignore"} (Github #214) + \item changed the default value of \code{check.nobs.vs.rankZ} to + \code{"ignore"} (Github #214) } } - \subsection{BUG FIXES}{ + \subsection{BUG FIXES}{ \itemize{ \item change gradient testing from absolute to relative \item fix \code{confint(.,method="boot")} to allow/work @@ -119,7 +187,8 @@ } } } -\section{CHANGES IN VERSION 1.1-6}{ + +\section{CHANGES IN VERSION 1.1-6 (2014-04-13)}{ This version incorporates no changes in functionality, just modifications to testing and dependencies for CRAN/backward compatibility. \subsection{BUG FIXES}{ @@ -133,7 +202,8 @@ } } } -\section{CHANGES IN VERSION 1.1-5}{ + +\section{CHANGES IN VERSION 1.1-5 (2014-03-14)}{ \subsection{BUG FIXES}{ \itemize{ \item improved NA handling in \code{simulate} and \code{refit} @@ -155,6 +225,7 @@ } } } + \section{CHANGES IN VERSION 1.1-4}{ \subsection{BUG FIXES}{ \itemize{ @@ -296,14 +367,15 @@ } } -\section{CHANGES IN VERSION 1.0-6 (2013-10-27)}{ +\section{CHANGES IN VERSION 1.0-6 (2014-02-02)}{ \subsection{BUG FIXES}{ \itemize{ \item prediction now works when new data have fewer factor levels than are present in the original data (Github issue #143, reported by Rune Haubo) \item the existence of a variable "new" in the global environment - would mess \code{lme4} up: reported at http://stackoverflow.com/questions/19801070/error-message-glmer-using-r-what-must-be-a-character-string-or-a-function + would mess \code{lme4} up: reported at + \url{http://stackoverflow.com/questions/19801070/error-message-glmer-using-r-what-must-be-a-character-string-or-a-function} } } } @@ -554,23 +626,23 @@ \section{CHANGES IN VERSION 0.999375-16 (2008-06-23)}{ \subsection{MAJOR USER-VISIBLE CHANGES}{ \itemize{ - \item The underlying algorithms and representations for all the - mixed-effects models fit by this package have changed - for - the better, we hope. The class "mer" is a common - mixed-effects model representation for linear, generalized - linear, nonlinear and generalized nonlinear mixed-effects - models. - \item ECME iterations are no longer used at all, nor are analytic - gradients. Components named 'niterEM', 'EMverbose', or - 'gradient' can be included in the 'control' argument to - lmer(), glmer() or nlmer() but have no effect. - \item PQL iterations are no longer used in glmer() and nlmer(). - Only the Laplace approximation is currently available. AGQ, - for certain classes of GLMMs or NLMMs, is being added. - \item The 'method' argument to lmer(), glmer() or nlmer() is - deprecated. Use the 'REML = FALSE' in lmer() to obtain ML - estimates. Selection of AGQ in glmer() and nlmer() will be - controlled by the argument 'nAGQ', when completed. + \item The underlying algorithms and representations for all the + mixed-effects models fit by this package have changed - for + the better, we hope. The class "mer" is a common + mixed-effects model representation for linear, generalized + linear, nonlinear and generalized nonlinear mixed-effects + models. + \item ECME iterations are no longer used at all, nor are analytic + gradients. Components named 'niterEM', 'EMverbose', or + 'gradient' can be included in the 'control' argument to + lmer(), glmer() or nlmer() but have no effect. + \item PQL iterations are no longer used in glmer() and nlmer(). + Only the Laplace approximation is currently available. AGQ, + for certain classes of GLMMs or NLMMs, is being added. + \item The 'method' argument to lmer(), glmer() or nlmer() is + deprecated. Use the 'REML = FALSE' in lmer() to obtain ML + estimates. Selection of AGQ in glmer() and nlmer() will be + controlled by the argument 'nAGQ', when completed. } } \subsection{NEW FEATURES}{ Binary files /tmp/YPVG8_S46B/lme4-1.1-8/inst/testdata/badprof.rds and /tmp/wzFoJ4Js8o/lme4-1.1-10/inst/testdata/badprof.rds differ diff -Nru lme4-1.1-8/inst/tests/test-formulaEval.R lme4-1.1-10/inst/tests/test-formulaEval.R --- lme4-1.1-8/inst/tests/test-formulaEval.R 2015-06-06 17:12:42.000000000 +0000 +++ lme4-1.1-10/inst/tests/test-formulaEval.R 2015-08-19 12:47:37.000000000 +0000 @@ -18,10 +18,19 @@ modStr <- (paste("x ~", "y +", F, "+", rF)) modForm <- as.formula(modStr) + ## WARNING: these drop/environment tests are extremely sensitive to environment + ## they may fail/not fail, or fail differently, within a "testthat" environment vs. + ## when run interactively + ## AICvec <- c(77.0516381151634, 75.0819116367084, 75.1915023640827) expect_that(m_data.3 <- glmer( modStr , data=d, family="binomial"), is_a("glmerMod")) - expect_error(drop1(m_data.3),"'data' not found") expect_that(m_data.4 <- glmer( "x ~ y + z + (1|r)" , data=d, family="binomial"), is_a("glmerMod")) + ## interactively: + ## expect_equal(drop1(m_data.3)$AIC,AICvec) + ## expect_equal(drop1(m_data.4)$AIC,AICvec) + ## in test environment: + expect_error(drop1(m_data.3),"'data' not found") expect_error(drop1(m_data.4),"'data' not found") + }) test_that("glmerForm", { @@ -41,6 +50,7 @@ ## formulas have environments associated, but character vectors don't ## data argument not specified: ## should work, but documentation warns against it + expect_that(m_nodata.0 <- glmer( x ~ y + z + (1|r) , family="binomial"), is_a("glmerMod")) expect_that(m_nodata.1 <- glmer( as.formula(modStr) , family="binomial"), is_a("glmerMod")) expect_that(m_nodata.2 <- glmer( modForm , family="binomial"), is_a("glmerMod")) @@ -48,11 +58,12 @@ expect_that(m_nodata.4 <- glmer( "x ~ y + z + (1|r)" , family="binomial"), is_a("glmerMod")) ## apply drop1 to all of these ... - m_nodata_List <- list(m_nodata.0,m_nodata.1,m_nodata.2,m_nodata.3,m_nodata.4) + m_nodata_List <- list(m_nodata.0, + m_nodata.1,m_nodata.2,m_nodata.3,m_nodata.4) d_nodata_List <- lapply(m_nodata_List,drop1) rm(list=c("x","y","z","r")) - + ## data argument specified expect_that(m_data.0 <- glmer( x ~ y + z + (1|r) , data=d, family="binomial"), is_a("glmerMod")) expect_that(m_data.1 <- glmer( as.formula(modStr) , data=d, family="binomial"), is_a("glmerMod")) diff -Nru lme4-1.1-8/inst/tests/test-glmernb.R lme4-1.1-10/inst/tests/test-glmernb.R --- lme4-1.1-8/inst/tests/test-glmernb.R 1970-01-01 00:00:00.000000000 +0000 +++ lme4-1.1-10/inst/tests/test-glmernb.R 2015-08-19 12:47:37.000000000 +0000 @@ -0,0 +1,22 @@ +library("testthat") +library("lme4") + +context("glmer.nb") +test_that("basic", { + set.seed(101) + dd <- expand.grid(f1 = factor(1:3), + f2 = LETTERS[1:2], g=1:9, rep=1:15, + KEEP.OUT.ATTRS=FALSE) + mu <- 5*(-4 + with(dd, as.integer(f1) + 4*as.numeric(f2))) + dd$y <- rnbinom(nrow(dd), mu = mu, size = 0.5) + require("MASS") + m.glm <- glm.nb(y ~ f1*f2, data=dd) + m.nb <- glmer.nb(y ~ f1*f2 + (1|g), data=dd, verbose=TRUE) + expect_null(m.nb@call$verbose) ## check: GH #321 + if (FALSE) { + ## FIXME: still trying to work this one out (GH #319) + expect_equal(fixef(m.nb), + coef (m.glm), tol=1e-5) + } +} +) diff -Nru lme4-1.1-8/inst/tests/test-glmer.R lme4-1.1-10/inst/tests/test-glmer.R --- lme4-1.1-8/inst/tests/test-glmer.R 2015-06-06 17:12:42.000000000 +0000 +++ lme4-1.1-10/inst/tests/test-glmer.R 2015-08-19 12:47:37.000000000 +0000 @@ -101,9 +101,10 @@ expect_is(glmer(prop ~ period + (1 | herd), data = cbppX, family = binomial, weights=size, na.action="na.exclude"), "glmerMod") + expect_is(glmer(prop ~ period + (1 | herd), - data = cbppX, family = binomial, weights=size, offset=rep(0,nrow(cbppX))), - "glmerMod") + data = cbppX, family = binomial, weights=size, offset=rep(0,nrow(cbppX))), + "glmerMod") expect_is(glmer(prop ~ period + (1 | herd), data = cbppX, family = binomial, weights=size, contrasts=NULL), "glmerMod") @@ -123,7 +124,7 @@ "extra argument.*disregarded") if(FALSE) { ## Hadley broke this expect_warning(glmer(cbind(incidence, size - incidence) ~ period + (1 | herd), - data = cbpp, family = binomial, + data = cbpp, family = binomial, control=list()), "instead of passing a list of class") expect_warning(glmer(cbind(incidence, size - incidence) ~ period + (1 | herd), @@ -136,8 +137,8 @@ mod <- glmer(presabs~predictor+(1|species),family=binomial, radinger_dat) expect_is(mod,"merMod") - ## TODO: is this reliable across platforms or do we have to loosen? - expect_equal(unname(fixef(mod)),c(0.5425528,6.4289962)) + ## tolerance: 32-bit Windows (CRAN) reported ave.diff of 5.33e-8 + expect_equal(unname(fixef(mod)), c(0.5425528,6.4289962), tolerance = 4e-7) set.seed(101) ## complete separation case d <- data.frame(y=rbinom(1000,size=1,p=0.5), diff -Nru lme4-1.1-8/inst/tests/test-lmer.R lme4-1.1-10/inst/tests/test-lmer.R --- lme4-1.1-8/inst/tests/test-lmer.R 2015-06-06 17:12:42.000000000 +0000 +++ lme4-1.1-10/inst/tests/test-lmer.R 2015-08-23 16:44:48.000000000 +0000 @@ -205,6 +205,24 @@ coefs=c(1,1), ctrl=lmerControl()$checkConv,lbound=0), "Problem with Hessian check") + + ## test ordering of Ztlist names + ## this is a silly model, just using it for a case + ## where nlevs(RE term 1) < nlevs(RE term 2)x + data(cbpp) + cbpp <- transform(cbpp,obs=factor(1:nrow(cbpp))) + fm0 <- lmer(incidence~1+(1|herd)+(1|obs),cbpp, + control=lmerControl(check.nobs.vs.nlev="ignore", + check.nobs.vs.rankZ="ignore", + check.nobs.vs.nRE="ignore", + check.conv.grad="ignore", + check.conv.singular="ignore", + check.conv.hess="ignore")) + fm0B <- update(fm0, .~1+(1|obs)+(1|herd)) + expect_equal(names(getME(fm0,"Ztlist")), + c("obs.(Intercept)", "herd.(Intercept)")) + ## stable regardless of order in formula + expect_equal(getME(fm0,"Ztlist"),getME(fm0B,"Ztlist")) }) ## test_that(..) diff -Nru lme4-1.1-8/inst/tests/test-methods.R lme4-1.1-10/inst/tests/test-methods.R --- lme4-1.1-8/inst/tests/test-methods.R 2015-06-20 01:00:03.000000000 +0000 +++ lme4-1.1-10/inst/tests/test-methods.R 2015-10-05 17:56:04.000000000 +0000 @@ -3,9 +3,8 @@ L <- load(system.file("testdata", "lme-tst-fits.rda", package="lme4", mustWork=TRUE)) -if (getRversion() > "3.0.0") { - ## saved fits are not safe with old R versions - +## FIXME: should test for old R versions, skip reloading test data in that +## case? fm0 <- fit_sleepstudy_0 fm1 <- fit_sleepstudy_1 fm2 <- fit_sleepstudy_2 @@ -41,9 +40,27 @@ context("summary") test_that("summary", { - ## test for multiple-correlation-warning bug - cc <- capture.output(print(summary(fit_agridat_archbold),correlation=FALSE)) - expect_true(length(g <- grep("not shown by default",cc))==0 || g<=1) + ## test for multiple-correlation-warning bug and other 'correlation = *' variants + ## Have 2 summary() versions, each with 3 print(.) ==> 6 x capture.output(.) + sf.aa <- summary(fit_agridat_archbold) + msg1 <- "Correlation.* not shown by default" + ## message => *not* part of capture.*(.) + expect_message(c1 <- capture.output(sf.aa), msg1) + # correlation = NULL - default + cF <- capture.output(print(sf.aa, correlation=FALSE)) + ## TODO? ensure the above gives *no* message/warning/error + expect_identical(c1, cF) + expect_message( + cT <- capture.output(print(sf.aa, correlation=TRUE)) + , "Correlation.* could have been required in summary()") + expect_identical(cF, cT[seq_along(cF)]) + sfT.aa <- summary(fit_agridat_archbold, correlation=TRUE) + expect_message(cT2 <- capture.output(sfT.aa), msg1) + expect_identical(cF, cT2) + cT3 <- capture.output(print(sfT.aa, correlation=TRUE)) + expect_identical(cT, cT3) + cF2 <- capture.output(print(sfT.aa, correlation=FALSE)) + expect_identical(cF, cF2) }) context("anova") @@ -162,6 +179,18 @@ suppressWarnings(confint(fm1, method="boot", FUN=testFun, nsim=10, quiet=TRUE)))), c(243.7551,256.9104),tol=1e-3) + + ## passing re.form to bootMer + FUN <- function(.){ + predict(., type="response") + } + fm2 <- lmer(strength ~ (1|batch/cask), Pastes) + expect_is(bootMer(fm2, predict, nsim=3),"boot") + expect_is(bootMer(fm2, predict, re.form=NULL, nsim=3),"boot") + expect_is(bootMer(fm2, predict, re.form=~(1|batch)+(1|cask:batch), nsim=3), + "boot") + expect_is(b3 <- bootMer(fm2, predict, re.form=~(1|batch), nsim=3), + "boot") }) context("confint_other") test_that("confint", { @@ -198,9 +227,48 @@ expect_equal(dimnames(ci1.p.n),dimnames(ci1.w.n)) expect_equal(dimnames(ci1.p.n),dimnames(ci1.b.n)) + ## test case of slightly wonky (spline fit fails) but monotonic profiles: + ## + simfun <- function(J,n_j,g00,g10,g01,g11,sig2_0,sig01,sig2_1){ + N <- sum(rep(n_j,J)) + x <- rnorm(N) + z <- rnorm(J) + mu <- c(0,0) + sig <- matrix(c(sig2_0,sig01,sig01,sig2_1),ncol=2) + u <- MASS::mvrnorm(J,mu=mu,Sigma=sig) + b_0j <- g00 + g01*z + u[,1] + b_1j <- g10 + g11*z + u[,2] + y <- rep(b_0j,each=n_j)+rep(b_1j,each=n_j)*x + rnorm(N,0,sqrt(0.5)) + sim_data <- data.frame(Y=y,X=x,Z=rep(z,each=n_j), + group=rep(1:J,each=n_j)) + } + set.seed(102) + dat <- simfun(10,5,1,.3,.3,.3,(1/18),0,(1/18)) + fit <- lmer(Y~X+Z+X:Z+(X||group),data=dat) + if (Sys.info()["sysname"] != "SunOS") { + ## doesn't produce warnings on Solaris ... + expect_warning(pp <- profile(fit,"theta_",quiet=TRUE), + "non-monotonic profile") + expect_warning(cc <- confint(pp),"falling back to linear interpolation") + ## very small/unstable problem, needs large tolerance + expect_equal(unname(cc[2,]),c(0,0.5427609),tolerance=1e-2) + } + badprof <- readRDS(system.file("testdata","badprof.rds", + package="lme4")) + expect_warning(cc <- confint(badprof), "falling back to linear") + expect_equal(cc, + structure(c(0, -1, 2.50856219044636, 48.8305727797906, NA, NA, + 33.1204478717389, 1, 7.33374326592662, 68.7254711217912, + -6.90462047196017, + NA), .Dim = c(6L, 2L), + .Dimnames = list(c(".sig01", ".sig02", + ".sig03", ".sigma", "(Intercept)", "cYear"), + c("2.5 %", "97.5 %"))), + tol=1e-3) }) + context("refit") test_that("refit", { s1 <- simulate(fm1) @@ -240,20 +308,21 @@ data("Orthodont", package = "MEMSS") # (differently "coded" from the 'default' "nlme" one) silly <- glmer(Sex ~ distance + (1|Subject), data = Orthodont, family = binomial) - sillypred <- data.frame(distance = c(20, 25)) - options(warn = 2) # no warnings! - ps <- predict(silly, sillypred, re.form=NA, type = "response") - expect_is(ps, "numeric") - expect_equal(unname(ps), c(0.999989632, 0.999997201)) - ## a case with interactions (failed in one temporary version): - ## fails differently on Windows and on other platforms? - expect_warning(fmPixS <<- update(fmPix, .~. + Side), - "(nearly unidentifiable|not uniquely determined)") } + sillypred <- data.frame(distance = c(20, 25)) + op <- options(warn = 2) # no warnings! + ps <- predict(silly, sillypred, re.form=NA, type = "response") + expect_is(ps, "numeric") + expect_equal(unname(ps), c(0.999989632, 0.999997201)) + ## a case with interactions (failed in one temporary version): + expect_warning(fmPixS <<- update(fmPix, .~. + Side), + "nearly unidentifiable|unable to evaluate scaled gradient|failed to converge") + ## (1|2|3); 2 and 3 seen (as Error??) on CRAN's Windows 32bit + options(op) + set.seed(1); ii <- sample(nrow(Pixel), 16) expect_equal(predict(fmPix, newdata = Pixel[ii,]), fitted(fmPix )[ii]) expect_equal(predict(fmPixS, newdata = Pixel[ii,]), fitted(fmPixS)[ii]) - options(warn = 0) set.seed(7); n <- 100; y <- rnorm(n) dd <- data.frame(id = factor(sample(10, n, replace = TRUE)), @@ -294,12 +363,14 @@ expect_equal(predict(fm3, newdata = model.frame(fm3)[2:3, ])[2], predict(fm3, newdata = model.frame(fm3)[3, ])) + ## complex-basis functions in RANDOM effect: (currently) + fm5 <- lmer(Reaction~Days+(poly(Days,2)|Subject),sleepstudy) + expect_equal(predict(fm5,sleepstudy[1,]),fitted(fm5)[1]) + ## complex-basis functions in FIXED effect are fine + fm6 <- lmer(Reaction~poly(Days,2)+(1|Subject),sleepstudy) + expect_equal(predict(fm6,sleepstudy[1,]),fitted(fm6)[1]) }) - - - - context("simulate") test_that("simulate", { expect_is(simulate(gm2), "data.frame") @@ -316,27 +387,38 @@ p9 <- simulate(gm2, re.form = NA, seed = 101) p10 <- simulate(gm2,use.u = FALSE, seed = 101) p11 <- simulate(gm2,use.u = TRUE, seed = 101) + ## minimal check of content: + expect_identical(colSums(p1[,1]), c(incidence = 95, 747)) + expect_identical(colSums(p2[,1]), c(incidence = 109, 733)) + ## equivalences: + ## group ~0: expect_equal(p2,p3) expect_equal(p2,p5) expect_equal(p2,p6) expect_equal(p2,p8) expect_equal(p2,p9) expect_equal(p2,p10) + ## group 1: expect_equal(p1,p4) expect_equal(p1,p7) expect_equal(p1,p11) expect_error(simulate(gm2,use.u = TRUE, re.form = NA), "should specify only one") + ## ## hack: test with three REs p1 <- lmer(diameter ~ (1|plate) + (1|plate) + (1|sample), Penicillin, - control = lmerControl(check.conv.hess = "ignore")) - expect_is(sp1 <- simulate(p1), "data.frame") - expect_true(all(dim(sp1) == c(nrow(Penicillin), 1))) + control = lmerControl(check.conv.hess = "ignore", + check.conv.grad = "ignore")) + expect_is(sp1 <- simulate(p1, seed=123), "data.frame") + expect_identical(dim(sp1), c(nrow(Penicillin), 1L)) + expect_equal(fivenum(sp1[,1]), + c(20.9412, 22.5805, 23.5575, 24.6095, 27.6997), tol=1e-5) ## Pixel example - expect_true(all(dim(simulate(fmPixS)) == c(nPix,1))) - expect_true(all(dim(simulate(fmPix )) == c(nPix,1))) + + expect_identical(dim(simulate(fmPixS)), c(nPix, 1L)) + expect_identical(dim(simulate(fmPix )), c(nPix, 1L)) ## simulation with newdata smaller/larger different from original - fm <- lmer(diameter ~ 1 + (1|plate) + (1|sample),Penicillin) + fm <- lmer(diameter ~ 1 + (1|plate) + (1|sample), Penicillin) expect_is(simulate(fm,newdata=Penicillin[1:10,],allow.new.levels=TRUE),"data.frame") expect_is(simulate(fm,newdata=do.call(rbind,replicate(4,Penicillin,simplify=FALSE))),"data.frame") @@ -345,19 +427,26 @@ dd <- data.frame(f=factor(rep(1:10,each=20)), x=runif(200), y=rnbinom(200,size=2,mu=2)) - g1 <- glmer.nb(y~x+(1|f),data=dd) - s1 <- simulate(g1) - expect_equal(mean(s1[[1]]),2.35) - - d <- sleepstudy - d$Subject <- factor(rep(1:18, each=10)) + g1 <- glmer.nb(y ~ x + (1|f), data=dd) + th.g1 <- getME(g1, "glmer.nb.theta") + ts1 <- table(s1 <- simulate(g1)[,1]) + expect_equal(fixef(g1), c("(Intercept)" = 0.630067, x = -0.0167248), tol = 1e-5) + expect_equal(th.g1, 2.013, tol = 1e-4) + expect_equal(th.g1, g1@call$family[["theta"]])# <- important for pkg{effects} eval() + expect_identical(sum(s1), 403) + expect_identical(as.vector(ts1[as.character(0:5)]), + c(51L, 54L, 36L, 21L, 14L, 9L)) + + ## Simulate with newdata with *new* RE levels: + d <- sleepstudy[-1] # droping the response ("Reaction") + ## d$Subject <- factor(rep(1:18, each=10)) ## Add 18 new subjects: - d <- rbind(sleepstudy, sleepstudy) + d <- rbind(d, d) d$Subject <- factor(rep(1:36, each=10)) - d$simulated <- simulate(fm1, seed=1, newdata=d[-1], + d$simulated <- simulate(fm1, seed=1, newdata = d, re.form=NULL, - allow.new.levels=TRUE)[[1]] - expect_equal(mean(d$simulated),299.9384608) + allow.new.levels = TRUE)[,1] + expect_equal(mean(d$simulated), 299.9384608) }) context("misc") @@ -369,7 +458,6 @@ } expect_is(as.data.frame(VarCorr(fm1)), "data.frame") }) -}# R >= 3.0.0 context("plot") test_that("plot", { @@ -406,3 +494,4 @@ data=grouseticks) expect_is(family(gnb),"family") }) + diff -Nru lme4-1.1-8/man/bootMer.Rd lme4-1.1-10/man/bootMer.Rd --- lme4-1.1-8/man/bootMer.Rd 2014-11-11 00:31:19.000000000 +0000 +++ lme4-1.1-10/man/bootMer.Rd 2015-09-20 00:07:32.000000000 +0000 @@ -2,7 +2,7 @@ \alias{bootMer} \title{Model-based (Semi-)Parametric Bootstrap for Mixed Models} \usage{ -bootMer(x, FUN, nsim = 1, seed = NULL, use.u = FALSE, +bootMer(x, FUN, nsim = 1, seed = NULL, use.u = FALSE, re.form=NA, type = c("parametric", "semiparametric"), verbose = FALSE, .progress = "none", PBargs = list(), parallel = c("no", "multicore", "snow"), @@ -24,6 +24,11 @@ inference is conditional on these values. If \code{FALSE}, new normal deviates are drawn (see Details).} + \item{re.form}{formula, \code{NA} (equivalent to \code{use.u=FALSE}), + or \code{NULL} (equivalent to \code{use.u=TRUE}): + alternative to \code{use.u} for + specifying which random effects to incorporate. + See \code{\link{simulate.merMod}} for details.} \item{type}{character string specifying the type of bootstrap, \code{"parametric"} or \code{"semiparametric"}; partial matching is allowed.} diff -Nru lme4-1.1-8/man/devcomp.Rd lme4-1.1-10/man/devcomp.Rd --- lme4-1.1-8/man/devcomp.Rd 2012-11-17 18:14:27.000000000 +0000 +++ lme4-1.1-10/man/devcomp.Rd 2015-08-23 16:36:36.000000000 +0000 @@ -21,7 +21,7 @@ \code{devcomp} slot as described in the value section. } \note{ - This function is deprecated, use \code{getME(., - "devcomp")} + This function is deprecated, use \code{\link{getME}(., "devcomp")} + %% ---------- } diff -Nru lme4-1.1-8/man/getME.Rd lme4-1.1-10/man/getME.Rd --- lme4-1.1-8/man/getME.Rd 2014-11-11 00:31:19.000000000 +0000 +++ lme4-1.1-10/man/getME.Rd 2015-08-20 12:25:28.000000000 +0000 @@ -1,17 +1,21 @@ \name{getME} +\title{Extract or Get Generalized Components from a Fitted Mixed Effects Model} \alias{getL} \alias{getL,merMod-method} \alias{getME} -\title{Extract or Get Generalized Components from a Fitted Mixed Effects Model} +\alias{getME.merMod} \usage{ -getME(object, +getME(object, name, ...) + +\S3method{getME}{merMod}(object, name = c("X", "Z", "Zt", "Ztlist", "mmList", "y", "mu", "u", "b", "Gp", "Tp", "L", "Lambda", "Lambdat", "Lind", "Tlist", "A", "RX", "RZX", "sigma", "flist", "fixef", "beta", "theta", "ST", "REML", "is_REML", "n_rtrms", "n_rfacs", "N", "n", "p", "q", "p_i", "l_i", "q_i", "k", "m_i", "m", - "cnms", "devcomp", "offset", "lower", "devfun")) + "cnms", "devcomp", "offset", "lower", "devfun", "glmer.nb.theta"), + \dots) } \arguments{ \item{object}{a fitted mixed-effects model of class @@ -19,83 +23,82 @@ \code{\link{lmer}()}, \code{\link{glmer}()} or \code{\link{nlmer}()}.} \item{name}{a character vector specifying the name(s) of - the \dQuote{component}. If \code{length(name)}>1, a named - list of components will be returned. Possible values are:\cr + the \dQuote{component}. If \code{length(name) > 1} or if \code{name + = "ALL"}, a named \code{\link{list}} of components will be returned. Possible values are:\cr \describe{ - \item{X}{fixed-effects model matrix} - \item{Z}{random-effects model matrix} - \item{Zt}{transpose + \item{\code{"X"}:}{fixed-effects model matrix} + \item{\code{"Z"}:}{random-effects model matrix} + \item{\code{"Zt"}:}{transpose of random-effects model matrix. Note that the structure of \code{Zt} has changed since \code{lme4.0}; to get a backward-compatible structure, use \code{do.call(Matrix::rBind,getME(.,"Ztlist"))}} - \item{Ztlist}{list of components of the transpose of the + \item{\code{"Ztlist"}:}{list of components of the transpose of the random-effects model matrix, separated by individual variance component} - \item{mmList}{list of raw model matrices associated with random + \item{\code{"mmList"}:}{list of raw model matrices associated with random effects terms} - \item{y}{response vector} - \item{mu}{conditional mean of the response} - \item{u}{conditional mode of the \dQuote{spherical} + \item{\code{"y"}:}{response vector} + \item{\code{"mu"}:}{conditional mean of the response} + \item{\code{"u"}:}{conditional mode of the \dQuote{spherical} random effects variable} - \item{b}{conditional mode of the + \item{\code{"b"}:}{conditional mode of the random effects variable} - \item{Gp}{groups pointer vector. + \item{\code{"Gp"}:}{groups pointer vector. A pointer to the beginning of each group of random effects corresponding to the random-effects terms, beginning with 0 and including a final element giving the total number of random effects} - \item{Tp}{theta pointer vector. A pointer to the beginning of the theta + \item{\code{"Tp"}:}{theta pointer vector. A pointer to the beginning of the theta sub-vectors corresponding to the random-effects terms, beginning with 0 and including a final element giving the number of thetas.} - \item{L}{sparse Cholesky factor of the penalized random-effects model.} - \item{Lambda}{relative covariance factor \eqn{\Lambda}{Lambda} of the random effects.} - \item{Lambdat}{transpose \eqn{\Lambda'}{Lambda'} of \eqn{\Lambda}{Lambda} above.} - \item{Lind}{index vector for inserting elements of + \item{\code{"L"}:}{sparse Cholesky factor of the penalized random-effects model.} + \item{\code{"Lambda"}:}{relative covariance factor \eqn{\Lambda}{Lambda} of the random effects.} + \item{\code{"Lambdat"}:}{transpose \eqn{\Lambda'}{Lambda'} of \eqn{\Lambda}{Lambda} above.} + \item{\code{"Lind"}:}{index vector for inserting elements of \eqn{\theta}{theta} into the nonzeros of \eqn{\Lambda}{Lambda}.} - \item{Tlist}{vector of template matrices from which the blocks of + \item{\code{"Tlist"}:}{vector of template matrices from which the blocks of \eqn{\Lambda}{Lambda} are generated.} - \item{A}{Scaled sparse model matrix (class + \item{\code{"A"}:}{Scaled sparse model matrix (class \code{"\link[Matrix:dgCMatrix-class]{dgCMatrix}"}) for the unit, orthogonal random effects, \eqn{U}, equal to \code{getME(.,"Zt") \%*\% getME(.,"Lambdat")}} - \item{RX}{Cholesky factor for the fixed-effects parameters} - \item{RZX}{cross-term in the full Cholesky factor} - \item{sigma}{residual standard error; note that \code{sigma(object)} is preferred.} - \item{flist}{a list of the grouping variables (factors) + \item{\code{"RX"}:}{Cholesky factor for the fixed-effects parameters} + \item{\code{"RZX"}:}{cross-term in the full Cholesky factor} + \item{\code{"sigma"}:}{residual standard error; note that \code{sigma(object)} is preferred.} + \item{\code{"flist"}:}{a list of the grouping variables (factors) involved in the random effect terms} - \item{fixef}{fixed-effects parameter estimates} - \item{beta}{fixed-effects parameter estimates (identical + \item{\code{"fixef"}:}{fixed-effects parameter estimates} + \item{\code{"beta"}:}{fixed-effects parameter estimates (identical to the result of \code{\link{fixef}}, but without names)} - \item{theta}{random-effects parameter estimates: these + \item{\code{"theta"}:}{random-effects parameter estimates: these are parameterized as the relative Cholesky factors of each random effect term} - \item{ST}{A list of S and T factors in the TSST' Cholesky + \item{\code{"ST"}:}{A list of S and T factors in the TSST' Cholesky factorization of the relative variance matrices of the random effects associated with each random-effects term. The unit lower triangular matrix, \eqn{T}, and the diagonal matrix, \eqn{S}, for each term are stored as a single matrix with diagonal elements from \eqn{S} and off-diagonal elements from \eqn{T}.} - \item{n_rtrms}{number of random-effects terms} - \item{n_rfacs}{number of distinct random-effects grouping factors} - \item{N}{number of rows of \code{X}} - \item{n}{length of the response vector, \code{y}} - \item{p}{number of columns of the fixed effects model matrix, \code{X}} - \item{q}{number of columns of the random effects model matrix, \code{Z}} - \item{p_i}{numbers of columns of the raw model matrices, - \code{mmList}} - \item{l_i}{numbers of levels of the grouping factors} - \item{q_i}{numbers of columns of the term-wise model matrices, \code{ZtList}} - \item{k}{number of random effects terms} - \item{m_i}{numbers of covariance parameters in each term} - \item{m}{total number of covariance parameters} - \item{cnms}{the \dQuote{component names}, a \code{\link{list}}.} - \item{REML}{\code{0} indicates the model was fitted by maximum - likelihood, any other positive integer indicates fitting by + \item{\code{"n_rtrms"}:}{number of random-effects terms} + \item{\code{"n_rfacs"}:}{number of distinct random-effects grouping factors} + \item{\code{"N"}:}{number of rows of \code{X}} + \item{\code{"n"}:}{length of the response vector, \code{y}} + \item{\code{"p"}:}{number of columns of the fixed effects model matrix, \code{X}} + \item{\code{"q"}:}{number of columns of the random effects model matrix, \code{Z}} + \item{\code{"p_i"}:}{numbers of columns of the raw model matrices, \code{mmList}} + \item{\code{"l_i"}:}{numbers of levels of the grouping factors} + \item{\code{"q_i"}:}{numbers of columns of the term-wise model matrices, \code{ZtList}} + \item{\code{"k"}:}{number of random effects terms} + \item{\code{"m_i"}:}{numbers of covariance parameters in each term} + \item{\code{"m"}:}{total number of covariance parameters} + \item{\code{"cnms"}:}{the \dQuote{component names}, a \code{\link{list}}.} + \item{\code{"REML"}:}{\code{0} indicates the model was fitted by maximum + likelihood, any other positive integer indicates fitting by restricted maximum likelihood} - \item{is_REML}{same as the result of \code{\link{isREML}(.)}} - \item{devcomp}{a list consisting of a named numeric vector, + \item{\code{"is_REML"}:}{same as the result of \code{\link{isREML}(.)}} + \item{\code{"devcomp"}:}{a list consisting of a named numeric vector, \code{cmp}, and a named integer vector, \code{dims}, describing the fitted model. The elements of \code{cmp} are:\cr \describe{ @@ -127,18 +130,25 @@ \item{useSc}{\code{TRUE} if model has a scale parameter} \item{reTrms}{number of random effects terms} \item{REML}{\code{0} indicates the model was fitted by maximum - likelihood, any other positive integer indicates fitting by + likelihood, any other positive integer indicates fitting by restricted maximum likelihood} \item{GLMM}{\code{TRUE} if a GLMM} \item{NLMM}{\code{TRUE} if an NLMM} } } - \item{offset}{model offset} - \item{lower}{lower bounds on model parameters (random effects + \item{\code{"offset"}:}{model offset} + \item{\code{"lower"}:}{lower bounds on model parameters (random effects parameters only).} - \item{devfun}{deviance function (so far only available for LMMs)} + \item{\code{"devfun"}:}{deviance function (so far only available for LMMs)} + \item{\code{"glmer.nb.theta"}:}{negative binomial \eqn{\theta} parameter, + only for \code{\link{glmer.nb}}.} + + %% -- keep at the very end: + \item{\code{"ALL"}:}{get all of the above as a \code{\link{list}}.} } } + \item{\dots}{currently unused in \pkg{lme4}, potentially further + arguments in methods.} } \value{ Unspecified, as very much depending on the \code{\link{name}}. @@ -169,16 +179,15 @@ Z <- getME(fm1, "Z") stopifnot(is(Z, "CsparseMatrix"), c(180,36) == dim(Z), - all.equal(fixef(fm1), getME(fm1, "beta"), + all.equal(fixef(fm1), b1 <- getME(fm1, "beta"), check.attributes=FALSE, tolerance = 0)) ## A way to get *all* getME()s : -getME.all <- function(obj) { - sapply(eval(formals(getME)$name), getME, object=obj, simplify=FALSE) -} ## internal consistency check ensuring that all work: -parts <- getME.all(fm1) +parts <- getME(fm1, "ALL") str(parts, max=2) +stopifnot(identical(Z, parts $ Z), + identical(b1, parts $ beta)) } \keyword{utilities} diff -Nru lme4-1.1-8/man/glmer.nb.Rd lme4-1.1-10/man/glmer.nb.Rd --- lme4-1.1-8/man/glmer.nb.Rd 2012-11-17 18:14:27.000000000 +0000 +++ lme4-1.1-10/man/glmer.nb.Rd 2015-08-19 12:47:37.000000000 +0000 @@ -1,21 +1,81 @@ \name{glmer.nb} +\title{Fitting GLMM's for Negative Binomial} \alias{glmer.nb} -\title{glmer() for Negative Binomial} +\alias{negative.binomial}% re-exported, needed e.g. in update() +\concept{GLMM} +\description{ + Fits a generalized linear mixed-effects model (GLMM) for the negative + binomial family, building on \code{\link{glmer}}, and initializing via + \code{\link[MASS]{theta.ml}} from \pkg{MASS}. +} \usage{ - glmer.nb(..., interval = log(th) + c(-3, 3), - verbose = FALSE) +glmer.nb(\dots, interval = log(th) + c(-3, 3), + tol = 5e-5, verbose = FALSE, nb.control = NULL, + initCtrl = list(limit = 20, eps = 2*tol, trace = verbose)) } \arguments{ - \item{...}{formula, data, etc: the arguments for - \code{\link{glmer}(..)} (apart from \code{family}!).} + \item{\dots}{arguments as for \code{glmer(.)} such as \code{formula}, + \code{data}, \code{control}, etc, but \emph{not} \code{family}!} + \item{interval}{interval in which to start the optimization. The + default is symmetric on log scale around the initially estimated theta.} + \item{tol}{tolerance for the optimization via \code{\link{optimize}}.} + \item{verbose}{\code{\link{logical}} indicating how much + progress information should be printed during the optimization. Use + \code{verbose = 2} (or larger) to enable \code{verbose=TRUE} in the + \code{\link{glmer}()} calls.} + \item{nb.control}{optional \code{\link{list}}, like \code{\link{glmerControl}()}, + used in \code{\link{refit}(*, control = control.nb)} during the + optimization.} + \item{initCtrl}{(\emph{\bold{experimental}, do not rely on this}:) a + \code{\link{list}} with named components as in the default, passed to + \code{\link[MASS]{theta.ml}} (package \pkg{MASS}) for the initial + value of the negative binomial parameter \code{theta}.} +} +\value{ + An object of class \code{glmerMod}, for which many + methods are available (e.g. \code{methods(class="glmerMod")}), see + \code{\link{glmer}}. +} +\note{For historical reasons, the shape parameter of the negative + binomial and the random effects parameters in our (G)LMM models are + both called \code{theta} (\eqn{\theta}), but are unrelated here. - \item{interval}{interval in which to start the - optimization} + The negative binomial \eqn{\theta} can be extracted from a fit + \code{g <- glmer.nb()} by \code{\link{getME}(g, "glmer.nb.theta")}. - \item{verbose}{logical indicating how much progress - information should be printed.} -} -\description{ - glmer() for Negative Binomial + + Parts of \code{glmer.nb()} are still experimental and methods are + still missing or suboptimal. In particular, there is no inference + available for the dispersion parameter \eqn{\theta}, yet. } +\seealso{ + \code{\link{glmer}}; from package \pkg{MASS}, + \code{\link[MASS]{negative.binomial}} (which we re-export currently) and + \code{\link[MASS]{theta.ml}}, the latter for initialization of + optimization. + The \sQuote{Details} of \code{\link{pnbinom}} for the definition of + the negative binomial distribution. +} +\examples{ +set.seed(101) +dd <- expand.grid(f1 = factor(1:3), + f2 = LETTERS[1:2], g=1:9, rep=1:15, + KEEP.OUT.ATTRS=FALSE) +summary(mu <- 5*(-4 + with(dd, as.integer(f1) + 4*as.numeric(f2)))) +dd$y <- rnbinom(nrow(dd), mu = mu, size = 0.5) +str(dd) +require("MASS")## and use its glm.nb() - as indeed we have zero random effect: +m.glm <- glm.nb(y ~ f1*f2, data=dd, trace=TRUE) +summary(m.glm) +m.nb <- glmer.nb(y ~ f1*f2 + (1|g), data=dd, verbose=TRUE) +m.nb +## The neg.binomial theta parameter: +getME(m.nb, "glmer.nb.theta") +LL <- logLik(m.nb) +## mixed model has 1 additional parameter (RE variance) +stopifnot(attr(LL,"df")==attr(logLik(m.glm),"df")+1) + +plot(m.nb, resid(.) ~ g)# works, as long as data 'dd' is found +} +\keyword{models} diff -Nru lme4-1.1-8/man/glmer.Rd lme4-1.1-10/man/glmer.Rd --- lme4-1.1-8/man/glmer.Rd 2015-06-18 22:19:15.000000000 +0000 +++ lme4-1.1-10/man/glmer.Rd 2015-08-19 12:47:37.000000000 +0000 @@ -4,7 +4,7 @@ \concept{ GLMM } \description{ Fit a generalized linear mixed-effects model (GLMM). Both fixed - effects and random effects are specified via the model \code{formula} + effects and random effects are specified via the model \code{formula}. } \usage{ glmer(formula, data = NULL, family = gaussian, control = glmerControl(), @@ -116,7 +116,7 @@ integral over the random effects space. For a linear mixed-effects model (LMM), as fit by \code{\link{lmer}}, this integral can be evaluated exactly. For a GLMM the integral must be approximated. The - most reliable approximation for GLMMs + most reliable approximation for GLMMs is adaptive Gauss-Hermite quadrature, at present implemented only for models with a single scalar random effect. The @@ -131,7 +131,7 @@ % deviance. For 20-dimensional evaluations of the GLM deviance per % evaluation of the approximate GLMM deviance. % The default approximation is the Laplace approximation, -% corresponding to \code{nAGQ=1}. +% corresponding to \code{nAGQ=1}. } \seealso{ @@ -139,6 +139,8 @@ parameterization); \code{\link[stats]{glm}} for Generalized Linear Models (\emph{without} random effects). \code{\link{nlmer}} for nonlinear mixed-effects models. + + \code{\link{glmer.nb}} to fit negative binomial GLMMs. } \examples{ ## generalized linear mixed model diff -Nru lme4-1.1-8/man/golden-class.Rd lme4-1.1-10/man/golden-class.Rd --- lme4-1.1-8/man/golden-class.Rd 2012-11-17 18:14:27.000000000 +0000 +++ lme4-1.1-10/man/golden-class.Rd 2015-08-23 16:36:36.000000000 +0000 @@ -1,10 +1,23 @@ -\docType{class} \name{golden-class} +\docType{class} +\title{Class \code{"golden"} and Generator for Golden Search Optimizer Class} \alias{golden-class} -\title{Class \code{"golden"}} +\alias{golden} \description{ - A reference class for a golden search scalar optimizer - using reverse communication. + \code{"golden"} is a reference class for a golden search scalar optimizer, + for a parameter within an interval. + + \code{golden()} is the generator for the \code{"golden"} + class. The optimizer uses reverse communications. +} +\usage{ +golden(...) +} +\arguments{ + \item{\dots}{(partly optional) arguments passed to + \code{\link{new}()} must be named arguments. \code{lower} and + \code{upper} are the bounds for the scalar parameter; they must be + finite.} } \section{Extends}{ All reference classes extend and inherit methods from @@ -12,6 +25,8 @@ } \examples{ showClass("golden") + +golden(lower= -100, upper= 1e100) } \keyword{classes} diff -Nru lme4-1.1-8/man/golden.Rd lme4-1.1-10/man/golden.Rd --- lme4-1.1-8/man/golden.Rd 2013-05-10 01:10:06.000000000 +0000 +++ lme4-1.1-10/man/golden.Rd 1970-01-01 00:00:00.000000000 +0000 @@ -1,30 +0,0 @@ -\name{golden} -\alias{golden} -\title{Generator object for the golden search optimizer class.} -\usage{ - golden(...) -} -\arguments{ - \item{\dots}{additional, optional arguments. None are - used at present.} -} -\description{ - The generator objects for the \code{\linkS4class{golden}} - class of a scalar optimizer for a parameter within an - interval. The optimizer uses reverse communications. -} -\note{ - Arguments to the \code{new} methods must be named - arguments. \code{lower} and \code{upper} are the bounds - for the scalar parameter; they must be finite. -} -\section{Methods}{ - \describe{ \item{\code{new(lower=lower, - upper=upper)}}{Create a new \code{\linkS4class{golden}} - object.} } -} -\seealso{ - \code{\linkS4class{golden}} -} -\keyword{classes} - diff -Nru lme4-1.1-8/man/lmerControl.Rd lme4-1.1-10/man/lmerControl.Rd --- lme4-1.1-8/man/lmerControl.Rd 2015-04-29 21:20:03.000000000 +0000 +++ lme4-1.1-10/man/lmerControl.Rd 2015-09-04 18:45:01.000000000 +0000 @@ -26,7 +26,7 @@ use.last.params=FALSE, sparseX = FALSE, ## input checking options - check.nobs.vs.rankZ = "ignore", + check.nobs.vs.rankZ = "ignore", check.nobs.vs.nlev = "stop", check.nlev.gtreq.5 = "ignore", check.nlev.gtr.1 = "stop", @@ -225,7 +225,7 @@ \code{.makeCC} returns a list containing the check specification (action, tolerance, and optionally relative tolerance). - + } \details{ Note that (only!) the pre-fitting \dQuote{checking options} @@ -246,16 +246,17 @@ str(lmerControl()) str(glmerControl()) \dontrun{ - ## fit with default Nelder-Mead algorithm ... - fm0 <- lmer(Reaction ~ Days + (1 | Subject), sleepstudy) + ## fit with default "bobyqa" algorithm ... + fm0 <- lmer(Reaction ~ Days + ( 1 | Subject), sleepstudy) fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy) - ## or with minqa::bobyqa ... - fm1_bobyqa <- update(fm1,control=lmerControl(optimizer="bobyqa")) + ## or with "Nelder_Mead" (the previous default) ... + fm1_bobyqa <- update(fm1, control = lmerControl(optimizer="Nelder_Mead")) ## or with the nlminb function used in older (<1.0) versions of lme4; ## this will usually replicate older results require(optimx) - fm1_nlminb <- update(fm1,control=lmerControl(optimizer="optimx", - optCtrl=list(method="nlminb"))) + fm1_nlminb <- update(fm1, + control = lmerControl(optimizer= "optimx", + optCtrl = list(method="nlminb"))) ## The other option here is method="L-BFGS-B". ## Or we can wrap base::optim(): optimwrap <- function(fn,par,lower,upper,control=list(), @@ -266,16 +267,17 @@ ## "Brent" requires finite upper values (lower bound will always ## be zero in this case) if (method=="Brent") upper <- pmin(1e4,upper) - res <- optim(par=par,fn=fn,lower=lower,upper=upper, + res <- optim(par=par, fn=fn, lower=lower,upper=upper, control=control,method=method,...) - with(res,list(par=par, - fval=value, - feval=counts[1], - conv=convergence, - message=message)) + with(res, list(par = par, + fval = value, + feval= counts[1], + conv = convergence, + message = message)) } - fm0_brent <- update(fm0,control=lmerControl(optimizer="optimwrap", - optCtrl=list(method="Brent"))) + fm0_brent <- update(fm0, + control = lmerControl(optimizer = "optimwrap", + optCtrl = list(method="Brent"))) ## You can also use functions from the nloptr package. if (require(nloptr)) { defaultControl <- list(algorithm="NLOPT_LN_BOBYQA", @@ -284,15 +286,15 @@ for (n in names(defaultControl)) if (is.null(control[[n]])) control[[n]] <- defaultControl[[n]] res <- nloptr(x0=par,eval_f=fn,lb=lower,ub=upper,opts=control,...) - with(res,list(par=solution, - fval=objective, - feval=iterations, - conv=if (status>0) 0 else status, - message=message)) + with(res, list(par = solution, + fval = objective, + feval= iterations, + conv = if (status>0) 0 else status, + message = message)) } - fm1_nloptr <- update(fm1,control=lmerControl(optimizer="nloptwrap")) - fm1_nloptr_NM <- update(fm1,control=lmerControl(optimizer="nloptwrap", - optCtrl=list(algorithm="NLOPT_LN_NELDERMEAD"))) + fm1_nloptr <- update(fm1, control=lmerControl(optimizer="nloptwrap")) + fm1_nloptr_NM <- update(fm1, control=lmerControl(optimizer="nloptwrap", + optCtrl=list(algorithm="NLOPT_LN_NELDERMEAD"))) } ## other algorithm options include NLOPT_LN_COBYLA, NLOPT_LN_SBPLX } diff -Nru lme4-1.1-8/man/lmer.Rd lme4-1.1-10/man/lmer.Rd --- lme4-1.1-8/man/lmer.Rd 2014-11-13 20:34:53.000000000 +0000 +++ lme4-1.1-10/man/lmer.Rd 2015-08-19 12:47:37.000000000 +0000 @@ -140,7 +140,7 @@ \examples{ ## linear mixed models - reference values from older code (fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)) -summary(fm1)# (with its own print method) +summary(fm1)# (with its own print method; see class?merMod % ./merMod-class.Rd str(terms(fm1)) stopifnot(identical(terms(fm1, fixed.only=FALSE), diff -Nru lme4-1.1-8/man/lmResp-class.Rd lme4-1.1-10/man/lmResp-class.Rd --- lme4-1.1-8/man/lmResp-class.Rd 2013-02-15 03:30:53.000000000 +0000 +++ lme4-1.1-10/man/lmResp-class.Rd 2015-08-19 12:47:37.000000000 +0000 @@ -1,10 +1,9 @@ \name{lmResp-class} +\title{Reference Classes for Response Modules, \code{"(lm|glm|nls|lmer)Resp"}} \alias{glmResp-class} \alias{lmerResp-class} \alias{lmResp-class} \alias{nlsResp-class} -\title{Classes \code{"lmResp"}, \code{"glmResp"}, \code{"nlsResp"} and -\code{"lmerResp"}} \description{ Reference classes for response modules, including linear models, \code{"lmResp"}, generalized linear models, @@ -14,6 +13,12 @@ same name. As is customary, the generator object for each class has the same name as the class. } +\section{Extends}{ + All reference classes extend and inherit methods from + \code{"\linkS4class{envRefClass}"}. Furthermore, + \code{"glmResp"}, \code{"nlsResp"} and \code{"lmerResp"} + all extend the \code{"lmResp"} class. +} \note{ Objects from these reference classes correspond to objects in C++ classes. Methods are invoked on the C++ @@ -28,12 +33,6 @@ method. Access to the external pointer should be through the method, not through the field. } -\section{Extends}{ - All reference classes extend and inherit methods from - \code{"\linkS4class{envRefClass}"}. Furthermore, - \code{"glmResp"}, \code{"nlsResp"} and \code{"lmerResp"} - all extend the \code{"lmResp"} class. -} \examples{ showClass("lmResp") str(lmResp$new(y=1:4)) diff -Nru lme4-1.1-8/man/merMod-class.Rd lme4-1.1-10/man/merMod-class.Rd --- lme4-1.1-8/man/merMod-class.Rd 2015-06-18 21:22:07.000000000 +0000 +++ lme4-1.1-10/man/merMod-class.Rd 2015-08-19 12:47:37.000000000 +0000 @@ -50,12 +50,12 @@ REMLcrit(object) \S3method{extractAIC}{merMod}(fit, scale = 0, k = 2, ...) \S3method{family}{merMod}(object, ...) -\S3method{formula}{merMod}(x, fixed.only = FALSE, ...) +\S3method{formula}{merMod}(x, fixed.only = FALSE, random.only = FALSE, ...) \S3method{fitted}{merMod}(object, ...) \S3method{logLik}{merMod}(object, REML = NULL, ...) \S3method{nobs}{merMod}(object, ...) \S3method{ngrps}{merMod}(object, ...) -\S3method{terms}{merMod}(x, fixed.only = TRUE, \dots) +\S3method{terms}{merMod}(x, fixed.only = TRUE, random.only = FALSE, \dots) \S3method{vcov}{merMod}(object, correlation = TRUE, sigm = sigma(object), use.hessian = NULL, \dots) \S3method{model.frame}{merMod}(formula, fixed.only = FALSE, ...) @@ -90,17 +90,26 @@ anova table.} \item{scale}{Not currently used (see \code{\link{extractAIC}}).} \item{k}{see \code{\link{extractAIC}}.} - \item{REML}{Logical. If \code{TRUE}, return the restricted log-likelihood + \item{REML}{Logical. If \code{TRUE}, return the restricted log-likelihood rather than the log-likelihood. If \code{NULL} (the default), set \code{REML} to \code{isREML(object)} (see \code{\link{isREML}}).} \item{fixed.only}{logical indicating if only the fixed effects - \code{\link{terms}} are sought, defaults to true. If false, all - terms, including random ones are returned.} - \item{correlation}{(logical) for \code{vcov}, indicates whether the - correlation matrix as well as the variance-covariance matrix is - desired; for \code{print.summary.merMod}, indicates whether the - correlation matrix of the fixed-effects parameters should be - printed.} + components (terms or formula elements) are sought. If false, all + components, including random ones, are returned.} + \item{random.only}{complement of \code{fixed.only}; indicates + whether random components only are sought. (Trying to specify + \code{fixed.only} and \code{random.only} at the same time + will produce an error.)} + \item{correlation}{(logical) + for \code{vcov}, indicates whether the correlation matrix as well as + the variance-covariance matrix is desired; + for \code{summary.merMod}, indicates whether the correlation matrix + should be computed and stored along with the covariance; + for \code{print.summary.merMod}, indicates whether the correlation + matrix of the fixed-effects parameters should be printed. In the + latter case, when \code{NULL} (the default), the correlation matrix + is printed when it has been computed by \code{summary(.)}, and when + \eqn{p <= 20}.}% and '20' can be changed by options(lme4.summary.cor.max = ) \item{use.hessian}{(logical) indicates whether to use the finite-difference Hessian of the deviance function to compute standard errors of the fixed effects, rather estimating @@ -184,7 +193,7 @@ } } \section{Deviance and log-likelihood of GLMMs}{ - + One must be careful when defining the deviance of a GLM. For example, should the deviance be defined as minus twice the log-likelihood or does it involve subtracting the deviance for a saturated model? To @@ -227,19 +236,19 @@ %% We define the conditional likelihood as the density of the response %% variable, \eqn{y}, conditional on the spherical random effects, %% \eqn{u}: - + %% \deqn{f_{\theta,\beta}(y|u)}{f(y|u)} - + %% which depends on the fixed effect and covariance parameters, %% \eqn{\beta} and \eqn{\theta}. We define the likelihood itself as: - + %% \deqn{f_{\theta,\beta}(y) = \int f_{\theta,\beta}(y|u) f(u)}{f(y) = integral(f(y|u)f(u))} - + %% where \eqn{f(u)} is an independent multivariate normal. It is not %% typically possible to evaluate this integral for GLMMs. In %% \code{lme4} we approximate it using the Laplace approximation and %% adaptive Gauss-Hermite quadrature. - + %% For canonical links, the Laplace approximation to the log-likelihood %% is: @@ -264,7 +273,7 @@ %% quadrature node (\code{\link{GHrule}}), and \eqn{u_{jk}} is the value %% of the spherical random effect coefficient associated with the %% \eqn{j}th level of the grouping factor and \eqn{k}th quadrature node. - + } \section{Slots}{ \describe{ diff -Nru lme4-1.1-8/man/mkMerMod.Rd lme4-1.1-10/man/mkMerMod.Rd --- lme4-1.1-8/man/mkMerMod.Rd 2015-06-05 17:51:40.000000000 +0000 +++ lme4-1.1-10/man/mkMerMod.Rd 2015-08-19 12:47:37.000000000 +0000 @@ -1,13 +1,18 @@ \name{mkMerMod} +\title{Create a 'merMod' Object} \alias{mkMerMod} -\title{Create a merMod object} +\description{ + Create an object of (a subclass of) class \code{\linkS4class{merMod}} + from the environment of the objective function and the value returned + by the optimizer. +} \usage{ - mkMerMod(rho, opt, reTrms, fr, mc, lme4conv=NULL) + mkMerMod(rho, opt, reTrms, fr, mc, lme4conv = NULL) } \arguments{ \item{rho}{the environment of the objective function} \item{opt}{the optimization result returned by the optimizer - (a list: see \code{\link{lmerControl}} for required elements)} + (a \code{\link{list}}: see \code{\link{lmerControl}} for required elements)} \item{reTrms}{random effects structure from the calling function (see \code{\link{mkReTrms}} for required elements)} \item{fr}{model frame (see \code{\link{model.frame}})} @@ -16,13 +21,10 @@ (results of \code{checkConv})} } \value{ - an object from a class that inherits from - \code{\linkS4class{merMod}} -} -\description{ - Create an object in a subclass of - \code{\linkS4class{merMod}} from the environment of the - objective function and the value returned by the - optimizer. + an object from a class that inherits from \code{\linkS4class{merMod}}. } - +%% FIXME: TODO +%% \examples{ +%% ## 1) An "empty" merMod object : +%% ## 2) A "lmer()-like" merMod object, using our "modular" approach instead of lmer() +%% } diff -Nru lme4-1.1-8/man/mkReTrms.Rd lme4-1.1-10/man/mkReTrms.Rd --- lme4-1.1-8/man/mkReTrms.Rd 2014-11-25 04:22:28.000000000 +0000 +++ lme4-1.1-10/man/mkReTrms.Rd 2015-08-19 12:47:37.000000000 +0000 @@ -17,28 +17,45 @@ \value{ a \code{\link{list}} with components \item{Zt}{transpose of the sparse model matrix for the random effects} - \item{Lambdat}{transpose of the sparse relative covariance factor} + \item{theta}{initial values of the covariance parameters} \item{Lind}{an integer vector of indices determining the mapping of the elements of the \code{theta} vector to the \code{"x"} slot of \code{Lambdat}} - \item{theta}{initial values of the covariance parameters} + \item{Gp}{} \item{lower}{lower bounds on the covariance parameters} + \item{Lambdat}{transpose of the sparse relative covariance factor} \item{flist}{list of grouping factors used in the random-effects terms} \item{cnms}{a list of column names of the random effects according to the grouping factors} + \item{Ztlist}{list of components of the transpose of the + random-effects model matrix, separated by random-effects term} } \seealso{ Other utilities: \code{\link{findbars}}, \code{\link{mkRespMod}}, \code{\link{nlformula}}, - \code{\link{nobars}}, \code{\link{subbars}} + \code{\link{nobars}}, \code{\link{subbars}}. + \code{\link{getME}} can retrieve these components + from a fitted model, although their values and/or forms + may be slightly different in the final fitted model from + their original values as returned from \code{mkReTrms}. } \examples{ data("Pixel", package="nlme") mform <- pixel ~ day + I(day^2) + (day | Dog) + (1 | Side/Dog) (bar.f <- findbars(mform)) # list with 3 terms -%% >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> FIXME -## Unfinished _FIXME_ -## Pfrm <- model.frame(bar.f, data = Pixel) +mf <- model.frame(subbars(mform),data=Pixel) +rt <- mkReTrms(bar.f,mf) +names(rt) } +%fm1 <- lmer(mform,Pixel) +%rt2 <- getME(fm1,names(rt)) +%for (i in seq_along(rt)) +% cat(names(rt)[[i]], +% isTRUE(all.equal(rt[[i]],rt2[[i]])), +% "\n") +% ## theta and Lambda components have the same structure, but have +% ## been updated in the fitting; flist is a data frame in reTrms +% ## but converted to a list; Ztlist is quite different (has been +% ## decomposed in getME() \keyword{utilities} \ No newline at end of file diff -Nru lme4-1.1-8/man/utilities.Rd lme4-1.1-10/man/utilities.Rd --- lme4-1.1-8/man/utilities.Rd 1970-01-01 00:00:00.000000000 +0000 +++ lme4-1.1-10/man/utilities.Rd 2015-09-04 18:39:15.000000000 +0000 @@ -0,0 +1,168 @@ +\name{prt-utilities} +\title{Print and Summary Method Utilities for Mixed Effects} +\alias{.prt.methTit} +\alias{.prt.VC} +\alias{.prt.aictab} +\alias{.prt.call} +\alias{.prt.family} +\alias{.prt.grps} +\alias{.prt.methTit} +\alias{.prt.resids} +\alias{.prt.warn} +\alias{formatVC} +\alias{llikAIC} +\alias{methTitle} +\description{ + The \code{\link{print}}, \code{\link{summary}} methods (including the + \code{print} for the \code{summary()} result) in \pkg{lme4} are + modular, using about ten small utility functions. Other packages, + building on \pkg{lme4} can use the same utilities for ease of + programming and consistency of output. + + Notably see the Examples. + + \code{llikAIC()} extracts the log likelihood, AIC, and related + statics from a Fitted LMM. + + \code{.prt.*()} all use \code{\link{cat}} and \code{\link{print}} to + produce output. +} +\usage{ +llikAIC(object, devianceFUN = devCrit, chkREML = TRUE, + devcomp = object@devcomp) + +methTitle(dims) + +.prt.methTit(mtit, class) +.prt.family (famL) +.prt.resids (resids, digits, title = "Scaled residuals:", \dots) +.prt.call (call, long = TRUE) +.prt.aictab (aictab, digits = 1) +.prt.VC (varcor, digits, comp, formatter = format, \dots) +.prt.grps (ngrps, nobs) +.prt.warn (optinfo, summary = FALSE, \dots) +formatVC (varcor, digits = max(3, getOption("digits") - 2), + comp = "Std.Dev.", formatter = format, + useScale = attr(varcor, "useSc"), + ...) + +} +\arguments{ + %% llikAIC() : + \item{object}{a LMM model fit} + \item{devianceFUN}{the function to be used for computing the deviance; + should not be changed for \pkg{lme4} created objects.} + \item{chkREML}{optional logical indicating if \code{object} maybe a REML + fit.}% use TRUE for \pkg{lme4} fits + \item{devcomp}{for \pkg{lme4} always the equivalent of + \code{object@devcomp}; here a \code{\link{list}}}%... FIXME + + %% methTitle(): + \item{dims}{for \pkg{lme4} always the equivalent of + \code{object@devcomp$dims}, a named vector or list with components + \code{"GLMM"}, \code{"NLMM"}, \code{"REML"}, and \code{"nAGQ"} of + which the first two are \code{\link{logical}} scalars, and the latter + two typically are \code{FALSE} or \code{\link{numeric}}.} + + %% .prt.methTit + \item{mtit}{the result of \code{methTitle(object)}} + \item{class}{typically \code{\link{class}(object)}.} + + %% .prt.family (famL) + \item{famL}{a \code{\link{list}} with components \code{family} and + \code{link}, each a \code{\link{character}} string; note that standard + \R \code{\link{family}} objects can be used directly, as well.} + + %% .prt.resids (resids, digits, title = "Scaled residuals:", \dots) + \item{resids}{numeric vector of model \code{\link{residuals}}.} + \item{digits}{non-negative integer of (significant) digits to print minimally.} + \item{title}{\code{\link{character}} string.} + \item{\dots}{optional arguments passed on, e.g., to \code{\link{residuals}()}.} + + %% .prt.call (call, long = TRUE) + \item{call}{the \code{\link{call}} of the model fit; e.g., available + via (generic) function \code{\link{getCall}()}.} + \item{long}{logical indicating if the output may be long, e.g., + printing the \code{control} part of the call if there is one.} + + %% .prt.aictab (aictab, digits = 1) + \item{aictab}{typically the \code{AICtab} component of the result of + \code{llikAIC()}.} + + %% .prt.VC (varcor, digits, comp, formatter = format, \dots) + \item{varcor}{typically the result of \code{\link{VarCorr}()}.} + \item{comp}{optional ...}% ... FIXME + \item{formatter}{a \code{\link{function}} used for formatting the numbers.} + + %% .prt.grps (ngrps, nobs) + \item{ngrps}{integer (vector), typically the result of + \code{\link{ngrps}(object)}.} + \item{nobs}{integer; the number of observations, e.g., the result + of \code{\link{nobs}}.} + + %% .prt.warn (optinfo, summary = FALSE, \dots) + \item{optinfo}{typically \code{object @ optinfo}, the optimization + infos, including warnings if there were.} +\item{summary}{logical} + +%% formatVC +\item{useScale}{(logical) whether the parent model estimates a scale + parameter} + +} +\value{ + \code{llikAIC()} returns a \code{\link{list}} with components + + \item{logLik}{which is \code{\link{logLik}(object)}, and} + \item{AICtab}{ a \dQuote{table} of \code{\link{AIC}}, \code{\link{BIC}}, + \code{\link{logLik}}, deviance and \code{\link{df.residual}()} values.} +} +\examples{ +## Create a few "lme4 standard" models ------------------------------ +fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy) +fmM <- update(fm1, REML=FALSE) # -> Maximum Likelihood + +gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd), + data = cbpp, family = binomial) +gmA <- update(gm1, nAGQ = 5) + + +(lA1 <- llikAIC(fm1)) +(lAM <- llikAIC(fmM)) +(lAg <- llikAIC(gmA)) + +(m1 <- methTitle(fm1 @ devcomp $ dims)) +(mM <- methTitle(fmM @ devcomp $ dims)) +(mG <- methTitle(gm1 @ devcomp $ dims)) +(mA <- methTitle(gmA @ devcomp $ dims)) + +.prt.methTit(m1, class(fm1)) +.prt.methTit(mA, class(gmA)) + +.prt.family(gaussian()) +.prt.family(binomial()) +.prt.family( poisson()) + +.prt.resids(residuals(fm1), digits = 4) +.prt.resids(residuals(fmM), digits = 2) + +.prt.call(getCall(fm1)) +.prt.call(getCall(gm1)) + +.prt.aictab ( lA1 $ AICtab ) # REML +.prt.aictab ( lAM $ AICtab ) # ML --> AIC, BIC, ... + +.prt.VC( VarCorr(fm1), digits = 4, formatter = format) +## Random effects: +## Groups Name Std.Dev. Corr +## Subject (Intercept) 24.740 +## Days 5.922 0.07 +## Residual 25.592 + +.prt.grps(ngrps = ngrps(fm1), + nobs = nobs (fm1)) +## --> Number of obs: 180, groups: Subject, 18 + +.prt.warn(fm1 @ optinfo) # nothing .. had no warnings +} +\keyword{utilities} diff -Nru lme4-1.1-8/MD5 lme4-1.1-10/MD5 --- lme4-1.1-8/MD5 2015-06-22 07:09:59.000000000 +0000 +++ lme4-1.1-10/MD5 2015-10-06 06:08:29.000000000 +0000 @@ -1,28 +1,28 @@ 2e5b726f280690ac96aa8279ec72c848 *ChangeLog -3373cc2aa6ad07468cba6e702ee74b85 *DESCRIPTION -bd3b98991254194f6d50684f8a57480c *NAMESPACE -69ed5127bd7ce619c475d5728035444e *R/AllClass.R +752cc39087f1824d78b98dde873e901a *DESCRIPTION +3d88311f9b329890d7282bf221b209a0 *NAMESPACE +9af71f88a54ea40f543360462351ac17 *R/AllClass.R c1eda7d4300d934fccc1271aba17f3a1 *R/AllGeneric.R 1e7b92f8de2a1188c3fa118c94fbc548 *R/GHrule.R -d87e18ddea0fd76232462531b0689adc *R/bootMer.R +45722688b9f97e33d9bb05ed68c06455 *R/bootMer.R 22fe5d7328a1dac86b29130e99d75e41 *R/checkConv.R 8f986a71400c586e9bf7c02dd1093a35 *R/deriv.R -d612a69940e4d68a2fc53917a606ae33 *R/hooks.R -919f774c793186520e6e794e2159d717 *R/lmList.R -4e51d65b8cf5274a97bfa48b862dc022 *R/lmer.R +7990e82a2970913d774fda874dc4dd90 *R/hooks.R +12302ab656c7c74b5b906a1f578fe7ed *R/lmList.R +4df68dc47919b258d40bf83d14f6ad36 *R/lmer.R 775978845374cd9f493154ddf9b7db48 *R/lmerControl.R 4bde5438bd9ddba5d3cd26e5a43a494b *R/mcmcsamp.R -37d473df2de188555c32e5a14c70bd01 *R/modular.R -77d0b243966f2b47d735acf51e5d18bd *R/nbinom.R +6b877b08cd1936284b18341cd81ed531 *R/modular.R +b8ef3b6a8d0030d626d0106a2d0b2108 *R/nbinom.R 89be54f6700b0d9dfbd5b5c28dbc6282 *R/optimizer.R -798a6f3ea8995550459ea3dddf024218 *R/plot.R -a684ee05ea762259224f79afc56b7f6e *R/predict.R -1a76b4a817ee13117cf3c71aa944241d *R/profile.R +8b14979bea44c5a5b589d4d17fd60e1f *R/plot.R +07bbb1242ffcaa82b9109a02fcffd8b2 *R/predict.R +d0bfef2199c25ed277896ca151455553 *R/profile.R 02ea75a63f3728709a68a702801ac6d3 *R/sparsegrid.R e057bf248efa620aa2c5c7ad122cd98d *R/sysdata.rda -99662168cf71db53187d9b55c5fdd168 *R/utilities.R +d1e158a17d5d70cbaff8d7198e74b30b *R/utilities.R ca217187200eb85c6e7ff5bac37a383a *R/vcconv.R -54da34dcf54f0ea6715afa1ebbf71c96 *README.md +04f4a9822cda60a932b43cf6de7756a7 *README.md d5c218fd9ffcb288bb236f782b63c58a *ToDo b4a18df55812c43d0230a7faf2abf1a0 *build/vignette.rds 99c921d83f9ba51e1653a5275fa8760c *configure.ac @@ -37,8 +37,8 @@ 09bbd2170531a7256a2cf743f4bec399 *data/cbpp.rda 28d1ffd0a450108d274431f7335e859d *data/grouseticks.rda 3eb17ae1551f7a264c4214460041908c *data/sleepstudy.rda -a6619e4aba637a90556c6c9bddca244a *inst/CITATION -f2afb4ad0beec44dabf1c7d8e0952d28 *inst/NEWS.Rd +8c81a4cdd4cd634e4102999545422b7e *inst/CITATION +af9e44c261b244b6e77678d54eace919 *inst/NEWS.Rd a1f824867c6c98ebe87a82f8ce95e3f6 *inst/doc/Doxyfile 7ea2caa8013f7ffc6e4c333fa1db1534 *inst/doc/PLSvGLS.R f84eeec527edf425f54bb4754e23da61 *inst/doc/PLSvGLS.Rnw @@ -58,6 +58,7 @@ 0ad4852aad9639f10ad13606ecf09908 *inst/testdata/SO_sep25.RData 56173e2cb895c04358bf1bf81bab1495 *inst/testdata/Thailand.rda 6dd6554b65cb2bfeb420364b1641a6f3 *inst/testdata/badhess.RData +aedcb80c37b6d01f4db5bd076ebc3298 *inst/testdata/badprof.rds da186eafec8322a55ec78451e19fd01a *inst/testdata/boo01L.RData b4b99dec3880c755f42a8dc9e77a067b *inst/testdata/colonizer_rand.rda 46c36252a311ba97a3413f84a3633d41 *inst/testdata/confint_ex.rda @@ -93,14 +94,15 @@ 60481bb722fbe5d948cd37ebdc42744c *inst/tests/test-catch.R f65087832e44b34225efdad862a90200 *inst/tests/test-doubleVertNotation.R e35e77e6b29924f0c5fb59e9a85225ed *inst/tests/test-factors.R -83cd3fc546575a219fab41c5182e8bc4 *inst/tests/test-formulaEval.R +a5b409d5411b21d07169d523c9fbbc62 *inst/tests/test-formulaEval.R 4fc4d941f9e8e52fcd53e30e687a9382 *inst/tests/test-glmFamily.R -c11a47d52fc1d0fdecefee78e16cd7e3 *inst/tests/test-glmer.R +a52dc2f4fd8e06c5b1d8ed87eec2ec83 *inst/tests/test-glmer.R 969dab36e89b1174d56c6fc84cd42217 *inst/tests/test-glmer.Rout +d2e5d3affd0c9bc9fbeffeb9f926a5fc *inst/tests/test-glmernb.R 77035e2c69747009e655f43ffdad341f *inst/tests/test-glmmFail.R -8cd345dc7efb1e689072252a93c011e9 *inst/tests/test-lmer.R +a8dd1c64fd3a3754625e518d2e2b5ecd *inst/tests/test-lmer.R 1ef8bc825d617df86ba8128cef02eff4 *inst/tests/test-lmerResp.R -fd4778c9c8fc81d459c9923ee008192f *inst/tests/test-methods.R +71c01379c8947d21d600929025597022 *inst/tests/test-methods.R c5a2a10463306115ca7aa399dbcc0fb6 *inst/tests/test-oldRZXfailure.R 09d4ffd86392d2e6851ac968ab6dc195 *inst/tests/test-rank.R 96da1f4dfe7466fe1320ed435199c984 *inst/tests/test-resids.R @@ -124,12 +126,12 @@ 2dfbaad44da43caaf910a65a19ecb2f7 *man/Penicillin.Rd 8cf2a0e6edc0891b30c7e15e06e60e1f *man/VarCorr.Rd 041c90508c36866544fec23b39ba1c25 *man/VerbAgg.Rd -c137b485e3b00950017f7caaa3937b78 *man/bootMer.Rd +2c8f37a212c8e08be26b17db8bd5c6a8 *man/bootMer.Rd f6191a9e9dbe98f609eeb90f8ca19133 *man/cake.Rd fd93f2971f3ebf948961188570e94b82 *man/cbpp.Rd 8c86f88fe902f884dbd1e15b6c5518ca *man/confint.merMod.Rd a354550440181ace499e4bc09f16bb67 *man/convergence.Rd -ca5d619637f766872600f92b39531c1f *man/devcomp.Rd +08522dd3a52eb723bb0520cc3ad66879 *man/devcomp.Rd 4387831c1db6d08a62c36ccf6b058732 *man/drop1.merMod.Rd 0ffadf84fce25e56463e86ccf165e04b *man/dummy.Rd 00b8ea9f50b7e2bef6e42dacf48348fb *man/expandDoubleVerts.Rd @@ -137,30 +139,29 @@ f699fae826b328940860943c133ac02e *man/findbars.Rd b6f832b80984b12fe181c79428184895 *man/fixef.Rd b6277de71fa73cdaa9b9259389971d18 *man/fortify.Rd -f5729a5008ce8ad5450c42a914f96440 *man/getME.Rd +7027a0d2904e83bc59511d3e104b4028 *man/getME.Rd c0c66bc34fc7cd0b242d62901a9d4cb8 *man/glmFamily-class.Rd b595fec321f7ca73ed88b843126baeef *man/glmFamily.Rd -5cba1347b6c11c5e430c3f39e0b663a8 *man/glmer.Rd -4ee5ed791ede41adfd03842a7fe1243f *man/glmer.nb.Rd +2a51f5cfc199cd1c81385b69e52291fa *man/glmer.Rd +8bad43f4cb7fc9046110384ffaa73eb9 *man/glmer.nb.Rd 36383a0393e3f8f285f7c1d3fe04160b *man/glmerLaplaceHandle.Rd -1262b2a7fa1c626eeddd08aff11b66b6 *man/golden-class.Rd -acbd4ee20f34021d218f551626852afc *man/golden.Rd +dc26eb9f86b5056b90e2e90d531490c1 *man/golden-class.Rd 44bb2d1bc6568facec7c4f3d0c0f2cf8 *man/grouseticks.Rd c4236aab077a04046402a23b86635a7b *man/hatvalues.merMod.Rd 3af905a45dfae46c1eda3ac1be6dd36c *man/isNested.Rd e10407c28ced55696cb726daac2d9c05 *man/isREML.Rd dad5bd62532208cd257f8390284f2f40 *man/lmList.Rd e4faeace55fd64773ae425a8d0f72713 *man/lmList4-class.Rd -227a86299cbe6c9fae0d7624b3ed09f7 *man/lmResp-class.Rd +7d8ebc023d6f58fd43e3a35047240342 *man/lmResp-class.Rd 83c485d70631b2f161fbd85e32b47f0b *man/lmResp.Rd 5be68d6000e53966c5c7fbb470609929 *man/lme4-package.Rd -a428ce63cb64e739f521740ef7721cdb *man/lmer.Rd -ae8aae9ff686fe1aa98221afecb50898 *man/lmerControl.Rd -39524e8e130013c2307a0b523b8adb22 *man/merMod-class.Rd +4d2ab197f565800b94f5bfd0d437d725 *man/lmer.Rd +66fcbb4186a5e10c7244513aff809659 *man/lmerControl.Rd +a932c248c59f641d36e0c31d13eb7b5c *man/merMod-class.Rd 67daa8376c15b9f3eebf77f83cd5b6c7 *man/merPredD-class.Rd 2c4a82fb8f0e4d1d3a5f754c242ccee5 *man/merPredD.Rd -a804d1cb7c32b7c48cefe02bad4e2257 *man/mkMerMod.Rd -1caf1dfffd3559b68b9fad70d7e3a860 *man/mkReTrms.Rd +09352c1fb6121b13186be99c445751a0 *man/mkMerMod.Rd +79017852aa08c419552fbb049c62abea *man/mkReTrms.Rd 0b8a7c8eb451e46d3c2f589fe257faab *man/mkRespMod.Rd 97cd5db2899f0c8587e1205df415cd21 *man/mkSimulationTemplate.Rd c061c134315aa205c8f412e3114e6617 *man/mkVarCorr.Rd @@ -188,6 +189,7 @@ 07b9b7279df8bb4ceeab7b390b4213c4 *man/sleepstudy.Rd 74597be45b8bc7612dfd345902340bfb *man/subbars.Rd a000d83fa06dd8e6c7cb84bb85f6547b *man/troubleshooting.Rd +f153e7778a5d4adcc8ad132656946e60 *man/utilities.Rd 71d3aefeda044a276a668f032a3a6b66 *man/varianceProf.Rd 3565765b67b88fddd3df9f82fd82738e *man/vcconv.Rd d36ee08626d65f230c54f0e26fe19284 *src/Makevars @@ -215,7 +217,7 @@ 6161bbbe3647bab0d31ef43e38a46818 *tests/agridat_gotway.R 75f09bf23e2ae1a3f09aa6a01631eeb0 *tests/bootMer.R c7f6dd652db43941cbe4fba46b4965ba *tests/bootMer.Rout -6fae7d7b84ff26885b137932064906cd *tests/boundary.R +d715f699ff662fe2b0455682fed3b5e9 *tests/boundary.R 827076d9e188c9cd0bfce43a42493081 *tests/condVar.R c6edf910eef9140bd26c1b213d39959e *tests/confint.R d4812b1e5a87481da2f1ffef165f4e77 *tests/confint.Rout @@ -229,7 +231,7 @@ 65c300f2156a21d8492b0a050800057f *tests/extras.R 39ef0a3821bd138f881e1edb6c38a0ba *tests/falsezero_dorie.R 099dc03c7490ef92e6343fa4892fd289 *tests/fewlevels.R -fce6689bab7ec8240432b46da15bc5c6 *tests/getME.R +897cfd778c8d948f842d6fbca0ec9624 *tests/getME.R cec2ff7851b2e2c430c5d9c643be0396 *tests/glmer-1.R 3a32381a5397cb2df8130d18665c3e65 *tests/glmerControlPass.R 03e7594ca74c8566534bc412e05b1ac5 *tests/glmerWarn.R @@ -237,7 +239,7 @@ 043f79ace204c4a59e19352dfcdd156a *tests/glmmWeights.R 26ec7a0086bdb879cad3a9f49b639fea *tests/hatvalues.R 51d9a3dbdcece6163e76fb867fce6ec4 *tests/is.R -deba44e35b937c0875124e3fd2749b33 *tests/lmList-tst.R +1d56eef1e6c4e909ee841122fddf4bb6 *tests/lmList-tst.R 0294c12565b3fc2f4cccbcf2f9940e63 *tests/lme-tst-fits.rda c322d71ff1cb5aae7f7752259339781d *tests/lme4_nlme.R 0796e3c062c28b5aa160f36e9de767ab *tests/lmer-0.R @@ -249,7 +251,7 @@ 4e30d737875ce21593cc4ddfb24a4311 *tests/minval.R 96ca4f4b592712ec3510bc4028a51bbe *tests/mkRout 5aeaf5b07676f61dba91683b4ea0a1cf *tests/modFormula.R -a93cd55731860f9d5ba4e76df0a9e212 *tests/nbinom.R +eb9d0da40437ec407001e1455d31c816 *tests/nbinom.R 3af38d4ca473ebcd29306eee0e5c5ec7 *tests/nlmer-conv.R 8f8452a6f2f995a1cb2f3d98ce6495c4 *tests/nlmer.R e6993cb36c085cea5d4bd324c96e912e *tests/nlmer.Rout.save @@ -272,7 +274,7 @@ f41d72b2048efd5d4aee555f78e1292f *tests/resids.R 71936dde95ed77155e38f49177119801 *tests/respiratory.R 3e319657e6f69219b1e18cdb9d2603a0 *tests/simulate.R -f9f746963a905227459b6d86598ee78d *tests/testOptControl.R +8375fd1d85d82c024a5a2bbab8b01a14 *tests/testOptControl.R 304b26b3344af0ec8c65b59449dd0f8f *tests/test_times.R c01f9cfa35d54301d5ddf216832170a0 *tests/test_timings.RData 5e20fb18025336c8eb692d017103506d *tests/testcolonizer.R diff -Nru lme4-1.1-8/NAMESPACE lme4-1.1-10/NAMESPACE --- lme4-1.1-8/NAMESPACE 2015-06-05 17:51:40.000000000 +0000 +++ lme4-1.1-10/NAMESPACE 2015-09-08 23:41:51.000000000 +0000 @@ -1,19 +1,30 @@ -useDynLib(lme4,.registration=TRUE) +useDynLib(lme4, .registration=TRUE) ## base packages -importFrom("graphics", plot) -importFrom("grid", gpar,viewport) -importFrom("splines", backSpline,interpSpline,periodicSpline) -importFrom("methods", show) +importFrom("graphics", par, plot) +importFrom("grid", gpar, viewport) +importFrom("splines", backSpline, interpSpline, periodicSpline) +importFrom("methods", as, getDataPart, is, new, rbind2, show, slot, slotNames) importFrom("stats", - anova, coef, confint, deviance, df.residual, drop1, - extractAIC, family, fitted, formula, hatvalues, logLik, - model.matrix, nobs, profile, residuals, simulate, terms, - update, vcov, weights) -importFrom("utils", flush.console) + AIC, BIC, anova, approx, approxfun, + as.formula, asOneSidedFormula, ave, + coef, confint, delete.response, + deviance, df.residual, dnorm, drop.scope, drop1, extractAIC, + family, fitted, formula, gaussian, getCall, glm, + hatvalues, lm, logLik, model.extract, model.frame, + model.matrix, model.offset, model.response, model.weights, + na.exclude, na.omit, na.pass, napredict, naresid, nlminb, + nobs, optimize, pchisq, pnorm, poisson, ppoints, + predict, printCoefmat, profile, pt, qchisq, qnorm, + qqnorm, qt, quantile, rbinom, reformulate, reorder, + resid, residuals, rgamma, rnbinom, rnorm, rpois, + runif, sd, setNames, simulate, symnum, terms, terms.formula, + update, update.formula, var, vcov, weights) + +importFrom("utils", flush.console, packageVersion, stack, str) + ## Recommended packages -importClassesFrom("Matrix", corMatrix, dgCMatrix, dpoMatrix, dCHMsimpl) importFrom("lattice", bwplot, current.panel.limits, @@ -27,17 +38,25 @@ splom, strip.custom, strip.default, trellis.par.get, xyplot) importFrom("MASS", negative.binomial, theta.ml) + ## generics we provide methods for and re-export: importFrom("nlme", fixef,ranef, VarCorr, getGroups) importFrom("Matrix", drop0, rankMatrix, sparseMatrix, cBind,rBind, forceSymmetric, fac2sparse, KhatriRao, Diagonal) -importMethodsFrom("Matrix", coerce, "%*%",crossprod,tcrossprod, - t,diag, chol2inv,solve, colSums,rowSums) +importClassesFrom("Matrix", corMatrix, dgCMatrix, dpoMatrix, dCHMsimpl) +## methods incl. S4 generics: +importMethodsFrom("Matrix", coerce, cov2cor, "%*%", crossprod,tcrossprod, + t, diag, chol2inv, solve, colSums,rowSums) ## other CRAN packages: importFrom("minqa", bobyqa) importFrom("nloptr", nloptr) +## Re-Exports : +export(negative.binomial) +## currently needed for some eval()ing methods for glmer.nb() objects + +## Our Exports: export(bootMer) export(devcomp) export(dummy) @@ -106,17 +125,28 @@ } export(.simulateFun) export(subbars) -## export(tnames) export(updateGlmerDevfun) export(VarCorr) export(varianceProf) -exportClasses(glmerMod) -exportClasses(lmerMod) + +## print() and print.summary() utilities -- other [NG]LMM packages can import +export(.prt.VC, .prt.aictab, .prt.call, .prt.family, + .prt.grps, .prt.methTit, .prt.resids, .prt.warn, + formatVC, + llikAIC, methTitle) +## export(tnames) + + +##------ Our S4 Classes ------------------------ +exportClasses(glmerMod, lmerMod, nlmerMod, merMod) exportClasses(lmList4) -exportClasses(merMod) -exportClasses(nlmerMod) + +##------ Our S4 Generics / Methods -------------- exportMethods(getL) exportMethods(show) + + +##------ Our S3 Methods ------------------------- S3method(anova,merMod) S3method(as.data.frame,boot) S3method(as.data.frame,thpr) @@ -146,6 +176,7 @@ S3method(formula,merMod) ## S3method(fortify,merMod)# but don't hide: export(fortify.merMod) +S3method(getME, merMod) S3method(hatvalues,merMod) S3method(isGLMM,merMod) S3method(isLMM,merMod) @@ -192,6 +223,7 @@ S3method(weights,merMod) S3method(xyplot,thpr) + S3method(sigma,lmList4) ## Re-using S3 methods from nlme (for 'lmList') as methods for our 'lmList4': @@ -207,3 +239,4 @@ S3method(summary,lmList4) ## Auxiliaries: S3method(getGroups,lmList4) + diff -Nru lme4-1.1-8/R/AllClass.R lme4-1.1-10/R/AllClass.R --- lme4-1.1-8/R/AllClass.R 2015-06-06 17:12:42.000000000 +0000 +++ lme4-1.1-10/R/AllClass.R 2015-08-23 16:36:36.000000000 +0000 @@ -24,24 +24,6 @@ ## ## MM: _Or_ have "merPred" with X of class "mMatrix" := classUnion {matrix, Matrix} -##' Class \code{"merPredD"} - a dense predictor reference class -##' -##' A reference class for a mixed-effects model predictor module with a dense -##' model matrix for the fixed-effects parameters. The reference class is -##' associated with a C++ class of the same name. As is customary, the -##' generator object, \code{\link{merPredD}}, for the class has the same name as -##' the class. -##' @name merPredD-class -##' @note Objects from this reference class correspond to objects in a C++ -##' class. Methods are invoked on the C++ class object using the external -##' pointer in the \code{Ptr} field. When saving such an object the external -##' pointer is converted to a null pointer, which is why there are redundant -##' fields containing enough information as R objects to be able to regenerate -##' the C++ object. The convention is that a field whose name begins with an -##' upper-case letter is an R object and the corresponding field, whose name -##' begins with the lower-case letter is a method. References to the -##' external pointer should be through the method, not directly through the -##' \code{Ptr} field. merPredD <- setRefClass("merPredD", # Predictor class for mixed-effects models with dense X fields = @@ -262,38 +244,8 @@ "X", "Xwts", "Zt", "beta0", "delb", "delu", "theta", "u0") -##' Generator objects for the response classes -##' -##' The generator objects for the \code{\linkS4class{lmResp}}, -##' \code{\linkS4class{lmerResp}}, \code{\linkS4class{glmResp}} and -##' \code{\linkS4class{nlsResp}} reference classes. Such objects are -##' primarily used through their \code{new} methods. -##' -##' @aliases lmResp lmerResp glmResp nlsResp -##' @param ... List of arguments (see Note). -##' @note Arguments to the \code{new} methods must be named arguments. -##' \itemize{ -##' \item{y}{ the numeric response vector} -##' \item{family}{ a \code{\link{family}} object} -##' \item{nlmod}{ the nonlinear model function} -##' \item{nlenv}{ an environment holding data objects for evaluation of \code{nlmod}} -##' \item{pnames}{ a character vector of parameter names} -##' \item{gam}{ a numeric vector - the initial linear predictor} -##'} -##' @section Methods: -##' \describe{ -##' \item{\code{new(y=y)}:}{Create a new -##' \code{\linkS4class{lmResp}} or \code{\linkS4class{lmerResp}} object.} -##' \item{\code{new(family=family, y=y)}:}{Create a new -##' \code{\linkS4class{glmResp}} object.} -##' \item{\code{new(y=y, nlmod=nlmod, nlenv=nlenv, pnames=pnames, -##' gam=gam)}:}{Create a new -##' \code{\linkS4class{nlsResp}} object.} -##' } -##' @seealso \code{\linkS4class{lmResp}}, \code{\linkS4class{lmerResp}}, -##' \code{\linkS4class{glmResp}}, \code{\linkS4class{nlsResp}} -##' @keywords classes -##' @export +## -> ../man/lmResp-class.Rd +## ~~~~~~~~~~~~~~~~~~~~~~ lmResp <- # base class for response modules setRefClass("lmResp", fields = @@ -345,7 +297,7 @@ assign(field, forceCopy(current), envir = vEnv) } } - do.call("new", c(as.list(vEnv), Class=def)) + do.call(new, c(as.list(vEnv), Class=def)) }, initializePtr = function() { Ptr <<- .Call(lm_Create, y, weights, offset, mu, sqrtXwt, @@ -383,49 +335,9 @@ lmResp$lock("mu", "offset", "sqrtXwt", "sqrtrwt", "weights", "wtres")#, "y") -##' Classes \code{"lmResp"}, \code{"glmResp"}, \code{"nlsResp"} and -##' \code{"lmerResp"} -##' -##' Reference classes for response modules, including linear models, -##' \code{"lmResp"}, generalized linear models, \code{"glmResp"}, nonlinear -##' models, \code{"nlsResp"} and linear mixed-effects models, \code{"lmerResp"}. -##' Each reference class is associated with a C++ class of the same name. As is -##' customary, the generator object for each class has the same name as the -##' class. -##' @name lmResp-class -##' @aliases lmResp-class glmResp-class lmerResp-class nlsResp-class -##' @note Objects from these reference classes correspond to objects in C++ -##' classes. Methods are invoked on the C++ classes using the external pointer -##' in the \code{ptr} field. When saving such an object the external pointer is -##' converted to a null pointer, which is why there are redundant fields -##' containing enough information as R objects to be able to regenerate the C++ -##' object. The convention is that a field whose name begins with an upper-case -##' letter is an R object and the corresponding field whose name begins with the -##' lower-case letter is a method. Access to the external pointer should be -##' through the method, not through the field. -##' @section Extends: All reference classes extend and inherit methods from -##' \code{"\linkS4class{envRefClass}"}. Furthermore, \code{"glmResp"}, -##' \code{"nlsResp"} and \code{"lmerResp"} all extend the \code{"lmResp"} class. -##' @seealso \code{\link{lmer}}, \code{\link{glmer}}, \code{\link{nlmer}}, -##' \code{\linkS4class{merMod}}. -##' @examples -##' -##' showClass("lmResp") -##' str(lmResp$new(y=1:4)) -##' showClass("glmResp") -##' str(glmResp$new(family=poisson(), y=1:4)) -##' showClass("nlsResp") -##' showClass("lmerResp") -##' str(lmerResp$new(y=1:4)) -##' @keywords classes -NULL - -##' @export lmerResp <- - setRefClass("lmerResp", - fields= - list(REML="integer"), - contains="lmResp", + setRefClass("lmerResp", contains = "lmResp", + fields = list(REML = "integer"), methods= list(initialize = function(...) { REML <<- as.integer(list(...)$REML) @@ -454,10 +366,8 @@ ##' @export glmResp <- - setRefClass("glmResp", - fields= - list(eta="numeric", family="family", n="numeric"), - contains="lmResp", + setRefClass("glmResp", contains = "lmResp", + fields = list(eta = "numeric", family = "family", n = "numeric"), methods= list(initialize = function(...) { callSuper(...) @@ -551,19 +461,16 @@ ) ) - glmResp$lock("family", "n", "eta") + ##' @export nlsResp <- - setRefClass("nlsResp", - fields= - list(gam="numeric", - nlmod="formula", - nlenv="environment", - pnames="character" - ), - contains="lmResp", + setRefClass("nlsResp", contains = "lmResp", + fields= list(gam = "numeric", + nlmod = "formula", + nlenv = "environment", + pnames= "character"), methods= list(initialize = function(...) { callSuper(...) @@ -601,6 +508,7 @@ nlsResp$lock("nlmod", "nlenv", "pnames") + ##' Generator object for the \code{\linkS4class{glmFamily}} class ##' ##' The generator object for the \code{\linkS4class{glmFamily}} reference class. @@ -919,7 +827,7 @@ length(asgn <- as.integer(attr(Flist, "assign"))) == ntrms) cnms <<- Cnms flist <<- Flist - ncols <<- unname(vapply(cnms, length, 0L)) + ncols <<- unname(lengths(cnms)) nctot <<- unname(as.vector(tapply(ncols, asgn, sum))) nlevs <<- unname(vapply(flist, function(el) length(levels(el)), 0L)) # why not replace the sapply with ncols*nlevs[asgn] ?? @@ -979,7 +887,7 @@ Gp <<- getME(mer, "Gp") cnms <<- Cnms flist <<- Flist - ncols <<- unname(vapply(cnms, length, 0L)) + ncols <<- unname(lengths(cnms)) nctot <<- unname(as.vector(tapply(ncols, asgn, sum))) nlevs <<- unname(vapply(flist, function(el) length(levels(el)), 0L)) offsets <<- c(0L, cumsum(sapply(seq_along(asgn), @@ -1029,7 +937,7 @@ all(names(CV) == names(covar)), all(sapply(CV, isSymmetric)), all(sapply(CV, ncol) == covsiz)) - if (!all(vapply(cnms, length, 1L) == covsiz)) + if (!all(lengths(cnms) == covsiz)) error("setRECovar currently requires distinct grouping factors") theta <<- sapply(CV, function(mm) { diff -Nru lme4-1.1-8/R/bootMer.R lme4-1.1-10/R/bootMer.R --- lme4-1.1-8/R/bootMer.R 2015-04-15 19:40:55.000000000 +0000 +++ lme4-1.1-10/R/bootMer.R 2015-09-19 18:33:54.000000000 +0000 @@ -28,7 +28,8 @@ ##' ##' The semi-parametric variant is not yet implemented, and we only ##' provide a method for \code{\link{lmer}} and \code{\link{glmer}} results. -bootMer <- function(x, FUN, nsim = 1, seed = NULL, use.u = FALSE, +bootMer <- function(x, FUN, nsim = 1, seed = NULL, + use.u = FALSE, re.form = NA, type = c("parametric","semiparametric"), verbose = FALSE, .progress = "none", PBargs=list(), @@ -71,8 +72,17 @@ ## FIXME: remove prefix when incorporated in package if (type=="parametric") { - ss <- simulate(x, nsim=nsim, use.u=use.u, na.action=na.exclude) + argList <- list(x, nsim=nsim, na.action=na.exclude) + if (!missing(re.form)) { + argList <- c(argList,list(re.form=re.form)) + } else { + argList <- c(argList,list(use.u=use.u)) + } + ss <- do.call(simulate,argList) } else { + if (!missing(re.form)) + stop(paste(sQuote("re.form")), + "cannot be used with semiparametric bootstrapping") if (use.u) { if (isGLMM(x)) warning("semiparametric bootstrapping is questionable for GLMMs") ss <- replicate(nsim,fitted(x)+sample(residuals(x,"response"), diff -Nru lme4-1.1-8/R/hooks.R lme4-1.1-10/R/hooks.R --- lme4-1.1-8/R/hooks.R 2014-11-11 00:31:19.000000000 +0000 +++ lme4-1.1-10/R/hooks.R 2015-08-19 12:47:37.000000000 +0000 @@ -1,5 +1,9 @@ +.onLoad <- function(libname, pkgname) { + options(lme4.summary.cor.max = 20) +} + .onUnload <- function(libpath) { gc() - if (is.loaded("lmer_Deviance",PACKAGE="lme4")) + if (is.loaded("lmer_Deviance", PACKAGE="lme4")) library.dynam.unload("lme4", libpath) } diff -Nru lme4-1.1-8/R/lmer.R lme4-1.1-10/R/lmer.R --- lme4-1.1-8/R/lmer.R 2015-06-10 01:28:29.000000000 +0000 +++ lme4-1.1-10/R/lmer.R 2015-10-05 14:24:51.000000000 +0000 @@ -502,9 +502,9 @@ class(val) <- c("anova", class(val)) forms <- lapply(lapply(calls, `[[`, "formula"), deparse) structure(val, - heading = c(header, "Models:", - paste(rep(names(mods), times = vapply(forms, length, 1)), - unlist(forms), sep = ": "))) + heading = c(header, "Models:", + paste(rep(names(mods), times = lengths(forms)), + unlist(forms), sep = ": "))) } else { ## ------ single model --------------------- dc <- getME(object, "devcomp") @@ -518,7 +518,7 @@ nmeffects <- c("(Intercept)", nmeffects) ss <- unlist(lapply(split(ss, asgn), sum)) stopifnot(length(ss) == length(nmeffects)) - df <- vapply(split(asgn, asgn), length, 1L) + df <- lengths(split(asgn, asgn)) ## dfr <- unlist(lapply(split(dfr, asgn), function(x) x[1])) ms <- ss/df f <- ms/(sigma(object)^2) @@ -592,18 +592,55 @@ ##' @importFrom stats deviance ##' @S3method deviance merMod deviance.merMod <- function(object, REML = NULL, ...) { + ## type = c("conditional", "unconditional", "penalized"), + ## relative = TRUE, ...) { if (isGLMM(object)) { return(sum(residuals(object,type="deviance")^2)) + ## ------------------------------------------------------------ + ## proposed change to deviance function for GLMMs + ## ------------------------------------------------------------ + ## @param type Type of deviance (can be unconditional, + ## penalized, conditional) + ## @param relative Should deviance be shifted relative to a + ## saturated model? (only available with type == penalized or + ## conditional) + ## ------------------------------------------------------------ + ## ans <- switch(type[1], + ## unconditional = { + ## if (relative) { + ## stop("unconditional and relative deviance is undefined") + ## } + ## c(-2 * logLik(object)) + ## }, + ## penalized = { + ## sqrL <- object@pp$sqrL(1) + ## if (relative) { + ## object@resp$resDev() + sqrL + ## } else { + ## useSc <- unname(getME(gm1, "devcomp")$dims["useSc"]) + ## qLog2Pi <- unname(getME(object, "q")) * log(2 * pi) + ## object@resp$aic() - (2 * useSc) + sqrL + qLog2Pi + ## } + ## }, + ## conditional = { + ## if (relative) { + ## object@resp$resDev() + ## } else { + ## useSc <- unname(getME(gm1, "devcomp")$dims["useSc"]) + ## object@resp$aic() - (2 * useSc) + ## } + ## }) + ## return(ans) } if (isREML(object) && is.null(REML)) { warning("deviance() is deprecated for REML fits; use REMLcrit for the REML criterion or deviance(.,REML=FALSE) for deviance calculated at the REML fit") - return(devCritFun(object, REML=TRUE)) + return(devCrit(object, REML=TRUE)) } - devCritFun(object, REML=FALSE) + devCrit(object, REML=FALSE) } REMLcrit <- function(object) { - devCritFun(object, REML=TRUE) + devCrit(object, REML=TRUE) } ## original deviance.merMod -- now wrapped by REMLcrit @@ -617,7 +654,7 @@ ## REML=FALSE: ## if ML fit, return deviance ## if REML fit, compute and return deviance -devCritFun <- function(object, REML = NULL) { +devCrit <- function(object, REML = NULL) { ## cf. (1) lmerResp::Laplace in respModule.cpp ## (2) section 5.6 of lMMwR, listing lines 34-42 if (isTRUE(REML) && !isLMM(object)) @@ -852,7 +889,9 @@ ##' @importFrom stats formula ##' @S3method formula merMod -formula.merMod <- function(x, fixed.only=FALSE, ...) { +formula.merMod <- function(x, fixed.only=FALSE, random.only=FALSE, ...) { + if (missing(fixed.only) && random.only) fixed.only <- FALSE + if (fixed.only && random.only) stop("can't specify 'only fixed' and 'only random' terms") if (is.null(form <- attr(x@frame,"formula"))) { if (!grepl("lmer$",deparse(getCall(x)[[1]]))) stop("can't find formula stored in model frame or call") @@ -861,6 +900,10 @@ if (fixed.only) { form <- getFixedFormula(form) } + if (random.only) { + ## from predict.R + form <- reOnly(form,response=TRUE) + } form } @@ -886,8 +929,19 @@ } npar.merMod <- function(object) { - length(object@beta) + length(object@theta) + + n <- length(object@beta) + length(object@theta) + object@devcomp[["dims"]][["useSc"]] + ## FIXME: this is a bit of a hack: a user *might* have specified + ## negative binomial family with a known theta, in which case we + ## shouldn't count it as extra. Either glmer.nb needs to set a + ## flag somewhere, or we need class 'nbglmerMod' to extend 'glmerMod' ... + ## We do *not* want to use the 'useSc' slot (as above), because + ## although theta is in some sense a scale parameter, it's not + ## one in the formal sense (and isn't stored in the 'sigma' slot) + if (grepl("Negative Binomial",family(object)$family)) { + n <- n+1 + } + return(n) ## TODO: how do we feel about counting the scale parameter ??? } @@ -896,7 +950,7 @@ logLik.merMod <- function(object, REML = NULL, ...) { if (is.null(REML) || is.na(REML[1])) REML <- isREML(object) - val <- -devCritFun(object, REML = REML)/2 + val <- -devCrit(object, REML = REML)/2 ## dc <- object@devcomp nobs <- nobs.merMod(object) structure(val, @@ -1062,14 +1116,14 @@ sQuote("condVar")," instead") condVar <- postVar } - ans <- object@pp$b(1.) + ans <- object@pp$b(1) ## not always == c(matrix(unlist(getME(object,"b")))) if (!is.null(object@flist)) { ## evaluate the list of matrices levs <- lapply(fl <- object@flist, levels) asgn <- attr(fl, "assign") cnms <- object@cnms - nc <- vapply(cnms, length, 1L) - nb <- nc * vapply(levs, length, 1L)[asgn] + nc <- lengths(cnms) + nb <- nc * lengths(levs)[asgn] nbseq <- rep.int(seq_along(nb), nb) ml <- split(ans, nbseq) for (i in seq_along(ml)) @@ -1090,7 +1144,7 @@ if (condVar) { sigsqr <- sigma(object)^2 rp <- rePos$new(object) - if(any(sapply(rp$terms, length) > 1)){ + if(any(lengths(rp$terms) > 1L)) { # TODO: actually use condVar and then convert back to array format warning("conditional variances not currently available via ", "ranef when there are multiple terms per factor") @@ -1128,10 +1182,10 @@ refit.merMod <- function(object, newresp=NULL, rename.response=FALSE, maxit = 100L, ...) { - newControl <- NULL + ctrl.arg <- NULL if (ll <- length(l... <- list(...)) > 0) { if ((ll == 1L) && (names(l...)[1] == "control")) { - newControl <- l...$control + ctrl.arg <- l...$control } else { warning("additional arguments to refit.merMod ignored") @@ -1162,17 +1216,21 @@ ## somewhat repeated from profile.merMod, but sufficiently ## different that refactoring is slightly non-trivial ## "three minutes' thought would suffice ..." - ignore.pars <- c("xst", "xt") - control.internal <- object@optinfo$control - if (length(ign <- which(names(control.internal) %in% ignore.pars)) > 0) - control.internal <- control.internal[-ign] - if (!is.null(newControl)) { - control <- newControl - if (length(control$optCtrl)==0) - control$optCtrl <- control.internal - } else { - control <- if (isGLMM(object)) glmerControl() else lmerControl() - } + control <- + if (!is.null(ctrl.arg)) { + if (length(ctrl.arg$optCtrl) == 0) { ## use object's version: + obj.control <- object@optinfo$control + ignore.pars <- c("xst", "xt") + if (any(ign <- names(obj.control) %in% ignore.pars)) + obj.control <- obj.control[!ign] + ctrl.arg$optCtrl <- obj.control + } + ctrl.arg + } + else if (isGLMM(object)) + glmerControl() + else + lmerControl() ## we need this stuff defined before we call .glmerLaplace below ... pp <- object@pp$copy() @@ -1234,12 +1292,7 @@ length(rr$y)) } - ## hacking around to try to get internals properly set up - ## for refitting. This helps, but not all the way ... - ## oldresp <- rr$y # set this above from before copy - ## rr$setResp(newresp) - ## rr$setResp(oldresp) - ## rr$setResp(newresp) + if (isGLMM(object)) { GQmat <- GHrule(nAGQ) if (nAGQ <= 1) { @@ -1249,12 +1302,6 @@ grpFac = object@flist[[1]]) } } - ## .Call(glmerLaplace, pp$ptr(), rr$ptr(), nAGQ, - ## control$tolPwrss, as.integer(30), verbose) - ## nAGQ, - ## control$tolPwrss, as.integer(30), # maxit = 30 - ## verbose) - ## lp0 <- pp$linPred(1) # each pwrss opt begins at this eta devlist <- if (isGLMM(object)) { @@ -1319,26 +1366,30 @@ rho$pp <- new(class(xpp), X=xpp$X, Zt=xpp$Zt, Lambdat=xpp$Lambdat, Lind=xpp$Lind, theta=xpp$theta, n=nrow(xpp$X)) devfun <- mkdevfun(rho, 0L) # also pass ?? (verbose, maxit, control) - opt <- optwrap(optimizer, devfun, x@theta, lower=x@lower, - calc.derivs=TRUE) - ## FIXME: smarter calc.derivs rules - ## opt <- bobyqa(x@theta, devfun, x@lower) + opt <- ## "smart" calc.derivs rules + if(optimizer == "bobyqa" && !any("calc.derivs" == names(list(...)))) + optwrap(optimizer, devfun, x@theta, lower=x@lower, calc.derivs=TRUE, ...) + else + optwrap(optimizer, devfun, x@theta, lower=x@lower, ...) + ## FIXME: Should be able to call mkMerMod() here, and be done n <- length(rr$y) pp <- rho$pp p <- ncol(pp$X) - dims <- c(N=n, n=n, nmp=n-p, nth=length(pp$theta), p=p, q=nrow(pp$Zt), - nAGQ=NA_integer_, useSc=1L, reTrms=length(x@cnms), - spFe=0L, REML=0L, GLMM=0L, NLMM=0L) + dims <- c(N=n, n=n, p=p, nmp=n-p, q=nrow(pp$Zt), nth=length(pp$theta), + useSc=1L, reTrms=length(x@cnms), + spFe=0L, REML=0L, GLMM=0L, NLMM=0L)#, nAGQ=NA_integer_) wrss <- rho$resp$wrss() ussq <- pp$sqrL(1) pwrss <- wrss + ussq cmp <- c(ldL2=pp$ldL2(), ldRX2=pp$ldRX2(), wrss=wrss, ussq=ussq, pwrss=pwrss, drsum=NA, dev=opt$fval, REML=NA, sigmaML=sqrt(pwrss/n), sigmaREML=sqrt(pwrss/(n-p))) -### FIXME: Should modify the call slot to set REML=FALSE. It is -### tricky to do so without causing the call to be evaluated - new("lmerMod", call=x@call, frame=x@frame, flist=x@flist, + ## modify the call to have REML=FALSE. (without evaluating the call!) + cl <- x@call + cl[["REML"]] <- FALSE + new("lmerMod", call = cl, frame=x@frame, flist=x@flist, cnms=x@cnms, theta=pp$theta, beta=pp$delb, u=pp$delu, + optinfo = .optinfo(opt), lower=x@lower, devcomp=list(cmp=cmp, dims=dims), pp=pp, resp=rho$resp) } @@ -1434,13 +1485,20 @@ ##' @importFrom stats terms ##' @S3method terms merMod -terms.merMod <- function(x, fixed.only=TRUE, ...) { - if (fixed.only) { - tt <- terms.formula(formula(x,fixed.only=TRUE)) - attr(tt,"predvars") <- attr(attr(x@frame,"terms"),"predvars.fixed") - tt - } - else attr(x@frame,"terms") +terms.merMod <- function(x, fixed.only=TRUE, random.only=FALSE, ...) { + if (missing(fixed.only) && random.only) fixed.only <- FALSE + if (fixed.only && random.only) stop("can't specify 'only fixed' and 'only random' terms") + tt <- attr(x@frame,"terms") + if (fixed.only) { + tt <- terms.formula(formula(x,fixed.only=TRUE)) + attr(tt,"predvars") <- attr(terms(x@frame),"predvars.fixed") + } + if (random.only) { + tt <- terms.formula(subbars(formula(x,random.only=TRUE))) + ## FIXME: predvars should be random-only + attr(tt,"predvars") <- attr(terms(x@frame),"predvars.random") + } + tt } ##' @importFrom stats update @@ -1480,13 +1538,14 @@ ## 'Generalized linear mixed model fit by the adaptive Gaussian Hermite approximation' ## so did *not* mention "maximum likelihood" at all in the GLMM case -methTitle <- function(object, dims = object@devcomp$dims) { + +methTitle <- function(dims) { # dims == object@devcomp$dims GLMM <- dims[["GLMM"]] kind <- switch(1L + GLMM * 2L + dims[["NLMM"]], "Linear", "Nonlinear", "Generalized linear", "Generalized nonlinear") paste(kind, "mixed model fit by", - if(isREML(object)) "REML" + if(dims[["REML"]]) "REML" else paste("maximum likelihood", if(GLMM) { ## TODO? Use shorter wording here, for (new) 'long = FALSE' argument @@ -1550,15 +1609,22 @@ cat.f("Control:", dc) if (!is.null(cc <- call$subset)) cat.f(" Subset:", deparse(cc)) - } +} -getLlikAIC <- function(object, cmp = object@devcomp$cmp) { +##' @title Extract Log Likelihood, AIC, and related statics from a Fitted LMM +##' @param object a LMM model fit +##' @param devianceFUN the function to be used for computing the deviance; should not be changed for \pkg{lme4} created objects. +##' @param chkREML optional logical indicating \code{object} maybe a REML fit. +##' @param devcomp +##' @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. +llikAIC <- function(object, devianceFUN = devCrit, chkREML = TRUE, devcomp = object@devcomp) { llik <- logLik(object) # returns NA for a REML fit - maybe change? AICstats <- { - if(isREML(object)) cmp["REML"] # *no* likelihood stats here + if(chkREML && devcomp$dims[["REML"]]) + devcomp$cmp["REML"] # *no* likelihood stats here else { c(AIC = AIC(llik), BIC = BIC(llik), logLik = c(llik), - deviance = devCritFun(object), + deviance = devianceFUN(object), df.resid = df.residual(object)) } } @@ -1595,6 +1661,7 @@ ## FIXME: print header ("Warnings:\n") ? ## change position in output? comes at the very end, could get lost ... .prt.warn <- function(optinfo, summary=FALSE, ...) { + if(length(optinfo) == 0) return() # currently, e.g., from refitML() ## check all warning slots: print numbers of warnings (if any) cc <- optinfo$conv$opt msgs <- unlist(optinfo$conv$lme4$messages) @@ -1603,12 +1670,12 @@ nmsgs <- length(msgs) warnings <- optinfo$warnings nwarnings <- length(warnings) - if (cc>0 || nmsgs>0 || nwarnings>0) { + if (cc > 0 || nmsgs > 0 || nwarnings > 0) { if (summary) { cat(sprintf("convergence code %d; %d optimizer warnings; %d lme4 warnings", cc,nmsgs,nwarnings),"\n") } else { - cat(sprintf("convergence code: %d",cc), + cat(sprintf("convergence code: %d", cc), msgs, unlist(warnings), sep="\n") @@ -1617,7 +1684,9 @@ } } -.summary.cor.max <- 20 +## options(lme4.summary.cor.max = 20) --> ./hooks.R +## ~~~~~~~~ +## was .summary.cor.max <- 20 a lme4-namespace hidden global variable ## This is modeled a bit after print.summary.lm : ## Prints *both* 'mer' and 'merenv' - as it uses summary(x) mainly @@ -1645,26 +1714,29 @@ cat("\nFixed effects:\n") printCoefmat(x$coefficients, zap.ind = 3, #, tst.ind = 4 digits = digits, signif.stars = signif.stars) + ## do not show correlation when summary(*, correlation=FALSE) was used: + hasCor <- !is.null(VC <- x$vcov) && !is.null(VC@factors$correlation) if(is.null(correlation)) { # default - correlation <- p <= .summary.cor.max - if(!correlation) { + cor.max <- getOption("lme4.summary.cor.max") + correlation <- hasCor && p <= cor.max + if(!correlation && p > cor.max) { nam <- deparse(substitute(x)) if(length(nam) > 1 || nchar(nam) >= 32) nam <- "...." message(sprintf(paste( "\nCorrelation matrix not shown by default, as p = %d > %d.", "Use print(%s, correlation=TRUE) or", " vcov(%s) if you need it\n", sep = "\n"), - p, .summary.cor.max, nam, nam)) + p, cor.max, nam, nam)) } } else if(!is.logical(correlation)) stop("'correlation' must be NULL or logical") if(correlation) { - if(is.null(VC <- x$vcov)) VC <- vcov(x, correlation = TRUE) + if(is.null(VC)) VC <- vcov(x, correlation = TRUE) corF <- VC@factors$correlation - if (is.null(corF)) { + if (is.null(corF)) { # can this still happen? message("\nCorrelation of fixed effects could have been required in summary()") corF <- cov2cor(VC) - } ## else { + } p <- ncol(corF) if (p > 1) { rn <- rownames(x$coefficients) @@ -1701,12 +1773,12 @@ ranef.comp = "Std.Dev.", ...) { dims <- x@devcomp$dims - .prt.methTit(methTitle(x, dims = dims), class(x)) + .prt.methTit(methTitle(dims), class(x)) .prt.family(famlink(x, resp = x@resp)) .prt.call(x@call, long = FALSE) ## useScale <- as.logical(dims[["useSc"]]) - llAIC <- getLlikAIC(x) + llAIC <- llikAIC(x) .prt.aictab(llAIC$AICtab, 4) varcor <- VarCorr(x) .prt.VC(varcor, digits = digits, comp = ranef.comp, ...) @@ -1774,8 +1846,11 @@ } ## -> ../man/getME.Rd +getME <- function(object, name, ...) UseMethod("getME") + + ##' Extract or Get Generalize Components from a Fitted Mixed Effects Model -getME <- function(object, +getME.merMod <- function(object, name = c("X", "Z","Zt", "Ztlist", "mmList", "y", "mu", "u", "b", "Gp", "Tp", @@ -1789,15 +1864,21 @@ "p_i", "l_i", "q_i", "k", "m_i", "m", "cnms", "devcomp", "offset", "lower", - "devfun")) + "devfun", + "glmer.nb.theta" + ), ...) { if(missing(name)) stop("'name' must not be missing") - stopifnot(is(object,"merMod")) ## Deal with multiple names -- "FIXME" is inefficiently redoing things if (length(name <- as.character(name)) > 1) { names(name) <- name return(lapply(name, getME, object = object)) } + if(name == "ALL") ## recursively get all provided components + return(sapply(eval(formals()$name), + getME.merMod, object=object, simplify=FALSE)) + + stopifnot(is(object,"merMod")) name <- match.arg(name) rsp <- object@resp PR <- object@pp @@ -1808,7 +1889,7 @@ dims <- dc $ dims Tpfun <- function(cnms) { ltsize <- function(n) n*(n+1)/2 # lower triangle size - cLen <- cumsum(ltsize(vapply(cnms, length, 1L))) + cLen <- cumsum(ltsize(lengths(cnms))) setNames(c(0, cLen), c("beg__", names(cnms))) ## such that diff( Tp ) is well-named } @@ -1821,17 +1902,17 @@ "Z" = t(PR$Zt), "Zt" = PR$Zt, "Ztlist" = - { - getInds <- function(i) { - n <- diff(object@Gp)[i] ## number of elements in this block - nt <- length(cnms[[i]]) ## number of REs - inds <- lapply(seq(nt),seq,to = n,by = nt) ## pull out individual RE indices - inds <- lapply(inds,function(x) x + object@Gp[i]) ## add group offset - } - inds <- do.call(c,lapply(seq_along(cnms),getInds)) - setNames(lapply(inds,function(i) PR$Zt[i,]), - tnames(object,diag.only = TRUE)) - }, + { + getInds <- function(i) { + n <- diff(object@Gp)[i] ## number of elements in this block + nt <- length(cnms[[i]]) ## number of REs + inds <- lapply(seq(nt), seq, to = n, by = nt) ## pull out individual RE indices + inds <- lapply(inds,function(x) x + object@Gp[i]) ## add group offset + } + inds <- do.call(c,lapply(seq_along(cnms),getInds)) + setNames(lapply(inds,function(i) PR$Zt[i,]), + tnames(object,diag.only = TRUE)) + }, "mmList" = mmList.merMod(object), "y" = rsp$y, "mu" = rsp$mu, @@ -1851,11 +1932,11 @@ "fixef" = fixef(object), "beta" = object@beta, "theta" = setNames(th, tnames(object)), - "ST" = setNames(vec2STlist(object@theta, n = vapply(cnms, length, 0L)), + "ST" = setNames(vec2STlist(object@theta, n = lengths(cnms)), names(cnms)), "Tlist" = { - nc <- vapply(cnms, length, 1L) # number of columns per term - nt <- length(nc) # number of random-effects terms + nc <- lengths(cnms) # number of columns per term + nt <- length(nc) # number of random-effects terms ans <- vector("list",nt) names(ans) <- names(nc) pos <- 0L @@ -1891,33 +1972,37 @@ "devcomp" = dc, "offset" = rsp$offset, "lower" = object@lower, - "devfun" = - { - ## copied from refit ... DRY ... - verbose <- getCall(object)$verbose; if (is.null(verbose)) verbose <- 0L - devlist <- - if (isGLMM(object)) { - stop("getME('devfun') not yet working for GLMMs")## FIXME - baseOffset <- rsp$offset - nAGQ <- unname(dc$dims["nAGQ"]) - list(tolPwrss= dc$cmp ["tolPwrss"], - compDev = dc$dims["compDev"], - nAGQ = nAGQ, - lp0 = rsp$eta - baseOffset, - baseOffset = baseOffset, - pwrssUpdate = glmerPwrssUpdate, - GQmat = GHrule(nAGQ), - fac = object@flist[[1]], - pp=PR, resp=rsp, u0=PR$u0, dpars=seq_along(PR$theta), - verbose=verbose) - } - else - list(pp=PR, resp=rsp, u0=PR$u0, dpars=seq_along(PR$theta), verbose=verbose) - mkdevfun(rho=list2env(devlist), - ## FIXME: fragile ... - verbose=verbose, control=object@optinfo$control) - }, - ## FIXME: current version gives lower bounds for theta parameters only -- these must be extended for [GN]LMMs -- give extended value including -Inf values for beta values? + "devfun" = { + ## copied from refit ... DRY ... + verbose <- getCall(object)$verbose; if (is.null(verbose)) verbose <- 0L + devlist <- + if (isGLMM(object)) { + stop("getME('devfun') not yet working for GLMMs")## FIXME + baseOffset <- rsp$offset + nAGQ <- unname(dc$dims["nAGQ"]) + list(tolPwrss= dc$cmp ["tolPwrss"], + compDev = dc$dims["compDev"], + nAGQ = nAGQ, + lp0 = rsp$eta - baseOffset, + baseOffset = baseOffset, + pwrssUpdate = glmerPwrssUpdate, + GQmat = GHrule(nAGQ), + fac = object@flist[[1]], + pp=PR, resp=rsp, u0=PR$u0, dpars=seq_along(PR$theta), + verbose=verbose) + } + else + list(pp=PR, resp=rsp, u0=PR$u0, dpars=seq_along(PR$theta), verbose=verbose) + mkdevfun(rho=list2env(devlist), + ## FIXME: fragile ... + verbose=verbose, control=object@optinfo$control) + }, + ## FIXME: current version gives lower bounds for theta parameters only: + ## -- these must be extended for [GN]LMMs -- give extended value including -Inf values for beta values? + + "glmer.nb.theta" = if(isGLMM(object) && isNBfamily(rsp$family$family)) + getNBdisp(object) else NA, + "..foo.." = # placeholder! stop(gettextf("'%s' is not implemented yet", sprintf("getME(*, \"%s\")", name))), @@ -2075,7 +2160,7 @@ stop("VarCorr methods require reTrms, not just reModule") if(missing(sigma)) sigma <- sigma(x) - nc <- vapply(cnms, length, 1L) # no. of columns per term + nc <- lengths(cnms) # no. of columns per term structure(mkVarCorr(sigma, cnms = cnms, nc = nc, theta = x@theta, nms = { fl <- x@flist; names(fl)[attr(fl, "assign")]}), useSc = as.logical(x@devcomp$dims[["useSc"]]), @@ -2118,8 +2203,10 @@ ##' @param ... optional arguments for \code{formatter(*)} in addition ##' to the first (numeric vector) and \code{digits}. ##' @return a character matrix of formatted VarCorr entries from \code{varc}. -formatVC <- function(varc, digits = max(3, getOption("digits") - 2), - comp = "Std.Dev.", formatter = format, ...) +formatVC <- function(varcor, digits = max(3, getOption("digits") - 2), + comp = "Std.Dev.", formatter = format, + useScale = attr(varcor, "useSc"), + ...) { c.nms <- c("Groups", "Name", "Variance", "Std.Dev.") avail.c <- c.nms[-(1:2)] @@ -2128,23 +2215,22 @@ nc <- length(colnms <- c(c.nms[1:2], (use.c <- avail.c[mcc]))) if(length(use.c) == 0) stop("Must *either* show variances or standard deviations") - useScale <- attr(varc, "useSc") - reStdDev <- c(lapply(varc, attr, "stddev"), - if(useScale) list(Residual = unname(attr(varc, "sc")))) - reLens <- vapply(reStdDev, length, 1L) + reStdDev <- c(lapply(varcor, attr, "stddev"), + if(useScale) list(Residual = unname(attr(varcor, "sc")))) + reLens <- lengths(reStdDev) nr <- sum(reLens) reMat <- array('', c(nr, nc), list(rep.int('', nr), colnms)) reMat[1+cumsum(reLens)-reLens, "Groups"] <- names(reLens) - reMat[,"Name"] <- c(unlist(lapply(varc, colnames)), if(useScale) "") + reMat[,"Name"] <- c(unlist(lapply(varcor, colnames)), if(useScale) "") if(any("Variance" == use.c)) reMat[,"Variance"] <- formatter(unlist(reStdDev)^2, digits = digits, ...) if(any("Std.Dev." == use.c)) reMat[,"Std.Dev."] <- formatter(unlist(reStdDev), digits = digits, ...) if (any(reLens > 1)) { maxlen <- max(reLens) - recorr <- lapply(varc, attr, "correlation") + recorr <- lapply(varcor, attr, "correlation") corr <- - do.call("rBind", + do.call(rBind, lapply(recorr, function(x) { x <- as(x, "matrix") @@ -2165,7 +2251,7 @@ ##' @S3method summary merMod summary.merMod <- function(object, - correlation = (p <= .summary.cor.max), + correlation = (p <= getOption("lme4.summary.cor.max")), use.hessian = NULL, ...) { @@ -2177,10 +2263,8 @@ hess.avail <- (!is.null(h <- object@optinfo$derivs$Hessian) && nrow(h) > length(getME(object,"theta"))) if (is.null(use.hessian)) use.hessian <- hess.avail - if (use.hessian && !hess.avail) stop(shQuote("use.hessian"), - "=TRUE specified, ", - "but Hessian is unavailable") - + if (use.hessian && !hess.avail) + stop("'use.hessian=TRUE' specified, but Hessian is unavailable") resp <- object@resp devC <- object@devcomp dd <- devC$dims @@ -2203,12 +2287,12 @@ 2*pnorm(abs(cf3), lower.tail = FALSE)) } - llAIC <- getLlikAIC(object) + llAIC <- llikAIC(object) ## FIXME: You can't count on object@re@flist, ## nor compute VarCorr() unless is(re, "reTrms"): varcor <- VarCorr(object) # use S3 class for now - structure(list(methTitle = methTitle(object, dims = dd), + structure(list(methTitle = methTitle(dd), objClass = class(object), devcomp = devC, isLmer = is(resp, "lmerResp"), useScale = useSc, @@ -2278,7 +2362,7 @@ prepanel = prepanel.ci, panel = panel.ci, xlab = NULL, main = mtit, ...) } - setNames(lapply(names(x), f, ...),names(x)) + setNames(lapply(names(x), f, ...), names(x)) } ##' @importFrom graphics plot @@ -2543,17 +2627,19 @@ m <- matrix(NA,n,n) diag(m) <- seq(n) m[lower.tri(m)] <- (n+1):(n*(n+1)/2) - dd <- dd[m[lower.tri(m,diag=TRUE)],] + dd <- dd[m[lower.tri(m, diag=TRUE)],] } dd } r <- do.call(rbind, - mapply(tmpf,x,names(x),SIMPLIFY = FALSE)) + c(mapply(tmpf, x,names(x), SIMPLIFY = FALSE), + deparse.level = 0)) if (attr(x,"useSc")) { ss <- attr(x,"sc") r <- rbind(r,data.frame(grp = "Residual",var1 = NA,var2 = NA, vcov = ss^2, - sdcor = ss)) + sdcor = ss), + deparse.level = 0) } rownames(r) <- NULL r diff -Nru lme4-1.1-8/R/lmList.R lme4-1.1-10/R/lmList.R --- lme4-1.1-8/R/lmList.R 2015-04-15 15:40:54.000000000 +0000 +++ lme4-1.1-10/R/lmList.R 2015-08-19 12:47:37.000000000 +0000 @@ -53,25 +53,35 @@ ## FIXME: catch errors and pass them on as warnings? ## (simply passing them along silently gives confusing output) groups <- eval(mform$groups, frm) - val <- lapply(split(frm, groups), - if (!isGLM) ## lm(.) - function(dat, formula) { - tryCatch({ - data <- as.data.frame(dat) - lm(formula, data) - }, error = errorH) - } - else ## family ==> glm(.) - function(dat, formula) { - tryCatch({ - data <- as.data.frame(dat) - glm(formula, family=family, data) - }, error = errorH) - }, - formula = as.formula(mform$model)) + if (!is.factor(groups)) groups <- factor(groups) + ## FIXME: this splitting of data, weights, offset is really + ## ugly/brute force. I feel like there ought to be a way + ## to leverage the fact that 'weights' and 'offset' have + ## automatically been incorporated into the model frame ... + fit <- if (isGLM) glm else lm + mf2 <- if (missing(family)) NULL else list(family=family) + fitfun <- function(dat,formula) { + tryCatch({ + data <- as.data.frame(dat) + do.call(fit,c(list(formula, data, + weights=dat[["(weights)"]], + offset=dat[["(offset)"]]), + mf2)) + }, error=errorH) + } + frm.split <- split(frm, groups) + ## NB: levels() is only OK if grouping variable is a factor + nms <- names(frm.split) + ## null.split <- replicate(length(nms),NULL) + ## weights.split <- if (missing(weights)) null.split else split(weights, groups) + ## offset.split <- if (missing(offset)) null.split else split(offset, groups) + val <- ## mapply(fitfun, + lapply( + frm.split,fitfun, + ## weights.split, + ## offset.split, + formula = as.formula(mform$model)) - ## NB: identical(levels(groups), names(split(frm, groups))): - nms <- levels(groups) ## Contrary to nlme, we keep the erronous ones as well pool <- !isGLM || .hasScale(family2char(family)) new("lmList4", setNames(val, nms), diff -Nru lme4-1.1-8/R/modular.R lme4-1.1-10/R/modular.R --- lme4-1.1-8/R/modular.R 2015-06-06 17:12:42.000000000 +0000 +++ lme4-1.1-10/R/modular.R 2015-08-20 12:25:28.000000000 +0000 @@ -7,8 +7,9 @@ is.character(x) && !any(x == "ignore") } -RHSForm <- function(formula) { - formula[[length(formula)]] +RHSForm <- function(form,as.form=FALSE) { + rhsf <- form[[length(form)]] + if (as.form) reformulate(deparse(rhsf)) else rhsf } `RHSForm<-` <- function(formula,value) { @@ -330,7 +331,7 @@ ignoreArgs <- c("start","verbose","devFunOnly","control") l... <- list(...) l... <- l...[!names(l...) %in% ignoreArgs] - do.call("checkArgs", c(list("lmer"),l...)) + do.call(checkArgs, c(list("lmer"), l...)) if (!is.null(list(...)[["family"]])) { ## lmer(...,family=...); warning issued within checkArgs mc[[1]] <- quote(lme4::glFormula) @@ -394,6 +395,17 @@ fixedfr <- eval(mf, parent.frame()) attr(attr(fr,"terms"), "predvars.fixed") <- attr(attr(fixedfr,"terms"), "predvars") + + ## ran-effects model frame (for predvars) + ## important to COPY formula (and its environment)? + ranform <- formula + RHSForm(ranform) <- subbars(RHSForm(reOnly(formula))) + mf$formula <- ranform + ranfr <- eval(mf, parent.frame()) + attr(attr(fr,"terms"), "predvars.random") <- + attr(terms(ranfr), "predvars") + + ## FIXME: shouldn't we have this already in the full-frame predvars? X <- model.matrix(fixedform, fr, contrasts)#, sparse = FALSE, row.names = FALSE) ## sparseX not yet ## backward compatibility (keep no longer than ~2015): if(is.null(rankX.chk <- control[["check.rankX"]])) @@ -575,10 +587,10 @@ } } } - if (boundary.tol > 0) { - opt <- check.boundary(rho,opt,devfun,boundary.tol) - } - return(opt) + if (boundary.tol > 0) + check.boundary(rho, opt, devfun, boundary.tol) + else + opt } ## TODO: remove any arguments that aren't actually used by glFormula (same for lFormula) @@ -609,7 +621,7 @@ ignoreArgs <- c("start","verbose","devFunOnly","optimizer", "control", "nAGQ") l... <- list(...) l... <- l...[!names(l...) %in% ignoreArgs] - do.call("checkArgs", c(list("glmer"), l...)) + do.call(checkArgs, c(list("glmer"), l...)) cstr <- "check.formula.LHS" checkCtrlLevels(cstr, control[[cstr]]) @@ -662,6 +674,16 @@ fixedfr <- eval(mf, parent.frame()) attr(attr(fr,"terms"),"predvars.fixed") <- attr(attr(fixedfr,"terms"),"predvars") + + ## ran-effects model frame (for predvars) + ## important to COPY formula (and its environment)? + ranform <- formula + RHSForm(ranform) <- subbars(RHSForm(reOnly(formula))) + mf$formula <- ranform + ranfr <- eval(mf, parent.frame()) + attr(attr(fr,"terms"), "predvars.random") <- + attr(terms(ranfr), "predvars") + X <- model.matrix(fixedform, fr, contrasts)#, sparse = FALSE, row.names = FALSE) ## sparseX not yet ## backward compatibility (keep no longer than ~2015): if(is.null(rankX.chk <- control[["check.rankX"]])) @@ -691,9 +713,9 @@ c(reTrms[c("Zt","theta","Lambdat","Lind")], n=nrow(X), list(X=X))) rho$resp <- if (missing(fr)) - mkRespMod(family=family, ...) + mkRespMod(family=family, ...) else - mkRespMod(fr, family=family) + mkRespMod(fr, family=family) nAGQinit <- if(control$nAGQ0initStep) 0L else 1L ## allow trivial y if (length(y <- rho$resp$y) > 0) { diff -Nru lme4-1.1-8/R/nbinom.R lme4-1.1-10/R/nbinom.R --- lme4-1.1-8/R/nbinom.R 2015-04-03 15:21:36.000000000 +0000 +++ lme4-1.1-10/R/nbinom.R 2015-08-19 12:47:37.000000000 +0000 @@ -1,22 +1,27 @@ ##' @importFrom MASS negative.binomial ##' @importFrom MASS theta.ml -## should be getME(object,"NBdisp") ? -## MM: should the *user* use it? if yes, consider method sigma() or AIC() ? -getNBdisp <- function(object) { - get(".Theta",envir=environment(object@resp$family$aic)) -} +## ==> user should use getME(object, "glmer.nb.theta") +getNBdisp <- function(object) environment(object@resp$family$aic)[[".Theta"]] +## Hidden, originally used at least once, well tested (!), but +if(FALSE) # not needed anymore currently +getNBdisp.fam <- function(familyString) + as.numeric(sub(".*([-+]?\\d+(\\.\\d*)?([Ee][-+]?\\d+)?){1}.*", + "\\1", familyString)) + +## Package "constants" {only on depending the glmResp definition in ./AllClass.R}: +glmResp.f.nms <- names(glmResp$fields()) +glmNB.to.change <- setdiff(glmResp.f.nms, c("Ptr", "family")) -## should be setME(object,"NBdisp") ? +##' setNBdisp(object,theta) := +##' NB-object with changed [DISP]ersion parameter 'theta' (and all that entails) setNBdisp <- function(object,theta) { ## assign(".Theta",theta,envir=environment(object@resp$family$aic)) - ff <- setdiff(names(getRefClass("glmResp")$fields()),c("Ptr","family")) rr <- object@resp - arg1 <- lapply(ff,rr$field) - names(arg1) <- ff newresp <- do.call(glmResp$new, - c(arg1, list(family=negative.binomial(theta=theta)))) + c(lapply(setNames(nm=glmNB.to.change), rr$field), + list(family = negative.binomial(theta=theta)))) newresp$setOffset(rr$offset) newresp$updateMu(rr$eta - rr$offset) object@resp <- newresp @@ -27,59 +32,82 @@ refit(setNBdisp(object, theta), control = control) } +##' @title Optimize over the Negative Binomial Parameter Theta +##' @param object a "glmerMod" object, updated from poisson to negative.binomial() +##' @param interval and +##' @param tol are both passed to \code{\link{optimize}()}. +##' @param verbose ## show our own progress +##' @param control passed to \code{\link{refit}()} +##' @return the last fit, an object like 'object' optTheta <- function(object, interval = c(-5,5), - maxit = 20, + tol = .Machine$double.eps^0.25, verbose = FALSE, - control = NULL) { + control = NULL) +{ lastfit <- object - evalcnt <- 0 - optval <- optimize(function(t) { - ## FIXME: kluge to retain last value and evaluation count - ## Perhaps use a reference class object to keep track of this - ## auxilliary information? DB - dev <- deviance(lastfit <<- refitNB(lastfit, - theta = exp(t), - control = control)) - evalcnt <<- evalcnt+1 - if (verbose) cat(evalcnt,exp(t),dev,"\n") + it <- 0L + NBfun <- function(t) { + ## Kluge to retain last value and evaluation count {good enough for ..} + dev <- -2*logLik(lastfit <<- refitNB(lastfit, + theta = exp(t), + control = control)) + it <<- it+1L + if (verbose) + cat(sprintf("%2d: th=%#15.10g, dev=%#14.8f, beta[1]=%#14.8f\n", + it, exp(t), dev, lastfit@beta[1])) dev - }, interval=interval) - stopifnot(all.equal(optval$minimum,log(getNBdisp(lastfit)))) + } + optval <- optimize(NBfun, interval=interval, tol=tol) + stopifnot(all.equal(optval$minimum, log(getNBdisp(lastfit)))) ## FIXME: return eval count info somewhere else? MM: new slot there, why not? - attr(lastfit,"nevals") <- evalcnt + attr(lastfit,"nevals") <- it + ## fix up the 'th' expression, replacing it by the real number, + ## so effects:::mer.to.glm() can eval() it: + lastfit@call$family[["theta"]] <- exp(optval$minimum) lastfit } ## use MASS machinery to estimate theta from residuals -est_theta <- function(object) { +est_theta <- function(object, limit = 20, + eps = .Machine$double.eps^0.25, trace = 0) +{ Y <- model.response(model.frame(object)) mu <- fitted(object) - w <- object@resp$weights - control <- list(maxit=20,trace=0) - theta.ml(Y, mu, sum(w), w, limit = control$maxit, - trace = control$trace > 2) + theta.ml(Y, mu, weights = object@resp$weights, + limit = limit, eps = eps, trace = trace) } -## FIXME: really should document glmer.nb() on the same help page as glmer() -## I (MM) don't know how to use roxygen for that.. - +## -------> ../man/glmer.Rd ##' glmer() for Negative Binomial ##' @param ... formula, data, etc: the arguments for ##' \code{\link{glmer}(..)} (apart from \code{family}!). -##' @param interval interval in which to start the optimization -##' @param verbose logical indicating how much progress information -##' should be printed. -##' @export -glmer.nb <- function(..., interval = log(th)+c(-3,3), verbose=FALSE) +glmer.nb <- function(..., interval = log(th) + c(-3,3), + tol = 5e-5, verbose = FALSE, nb.control = NULL, + initCtrl = list(limit = 20, eps = 2*tol, trace = verbose)) { - control <- list(...)$control - g0 <- glmer(..., family = poisson) - th <- est_theta(g0) - if(verbose) cat("th := est_theta(glmer(..)) =", format(th),"\n") + dotE <- as.list(substitute(E(...))[-1]) + ## nE <- names(dotE <- as.list(substitute(E(...))[-1])) + ## i <- match(c("formula",""), nE) + ## i.frml <- i[!is.na(i)][[1]] # the index of "formula" in '...' + ## dots <- list(...) + g0 <- glmer(..., family = poisson, verbose = verbose >= 2) + th <- est_theta(g0, limit = initCtrl$limit, + eps = initCtrl$eps, trace = initCtrl$trace) + ## using format() on purpose, influenced by options(digits = *) : + if(verbose) cat("th := est_theta(glmer(..)) =", format(th)) g1 <- update(g0, family = negative.binomial(theta=th)) - ## if (is.null(interval)) interval <- log(th)+c(-3,3) - optTheta(g1, interval = interval, verbose = verbose, control = control) + if(verbose) cat(" --> dev.= -2*logLik(.) =", format(-2*logLik(g1)),"\n") + ## fix the 'data' part (only now!) + if("data" %in% names(g1@call)) + g1@call[["data"]] <- dotE[["data"]] + else + warning("no 'data = *' in glmer.nb() call.. Not much is guaranteed") + if ("verbose" %in% names(g1@call)) { + g1@call[["verbose"]] <- dotE[["verbose"]] + } + optTheta(g1, interval=interval, tol=tol, verbose=verbose, + control = nb.control) } ## do we want to facilitate profiling on theta?? diff -Nru lme4-1.1-8/R/plot.R lme4-1.1-10/R/plot.R --- lme4-1.1-8/R/plot.R 2015-01-29 22:47:13.000000000 +0000 +++ lme4-1.1-10/R/plot.R 2015-08-19 12:47:37.000000000 +0000 @@ -1,16 +1,19 @@ ## copied/modified from nlme -splitFormula <- - ## split, on the nm call, the rhs of a formula into a list of subformulas - function(form, sep = "/") + +##' split, on the nm call, the rhs of a formula into a list of subformulas +splitFormula <- function(form, sep = "/") { if (inherits(form, "formula") || - mode(form) == "call" && form[[1]] == as.name("~")) - return(splitFormula(form[[length(form)]], sep = sep)) - if (mode(form) == "call" && form[[1]] == as.name(sep)) - return(do.call("c", lapply(as.list(form[-1]), splitFormula, sep = sep))) - if (mode(form) == "(") return(splitFormula(form[[2]], sep = sep)) - if (length(form) < 1) return(NULL) - list(asOneSidedFormula(form)) + mode(form) == "call" && form[[1]] == as.name("~")) + splitFormula(form[[length(form)]], sep = sep) + else if (mode(form) == "call" && form[[1]] == as.name(sep)) + do.call(c, lapply(as.list(form[-1]), splitFormula, sep = sep)) + else if (mode(form) == "(") + splitFormula(form[[2]], sep = sep) + else if (length(form)) + list(asOneSidedFormula(form)) + ## else + ## NULL } ## Recursive version of all.vars @@ -23,27 +26,12 @@ } } -## crippled version of getData.gnls from nlme +## simple version of getData.gnls from nlme +## but we *should* and *can* work with environment(formula(.)) getData <- function(object) { mCall <- object@call - data <- eval(mCall$data,environment(formula(object))) - ## if (is.null(data)) return(data) - ## FIXME: deal with NAs, subset appropriately - ## naPat <- eval(mCall$naPattern) - ## if (!is.null(naPat)) { - ## data <- data[eval(naPat[[2]], data), , drop = FALSE] - ## } - ## naAct <- eval(mCall$na.action) - ## if (!is.null(naAct)) { - ## data <- naAct(data) - ## } - ## subset <- mCall@subset - ## if (!is.null(subset)) { - ## subset <- eval(asOneSidedFormula(subset)[[2]], data) - ## data <- data[subset, ] - ## } - return(data) + eval(mCall$data, environment(formula(object))) } asOneFormula <- @@ -220,10 +208,10 @@ ##' } ##' @S3method plot merMod ##' @method plot merMod -##' @export +##' @export plot.merMod <- function(x, form = resid(., type = "pearson") ~ fitted(.), abline, - id = NULL, idLabels = NULL, + id = NULL, idLabels = NULL, grid, ...) ## Diagnostic plots based on residuals and/or fitted values { diff -Nru lme4-1.1-8/R/predict.R lme4-1.1-10/R/predict.R --- lme4-1.1-8/R/predict.R 2015-06-20 01:01:18.000000000 +0000 +++ lme4-1.1-10/R/predict.R 2015-09-20 00:02:48.000000000 +0000 @@ -6,8 +6,10 @@ } ##' Random Effects formula only -reOnly <- function(f) { - reformulate(paste0("(", vapply(findbars(f), safeDeparse, ""), ")")) +reOnly <- function(f,response=FALSE) { + response <- if (response && length(f)==3) f[[2]] else NULL + reformulate(paste0("(", vapply(findbars(f), safeDeparse, ""), ")"), + response=response) } reFormHack <- function(re.form,ReForm,REForm,REform) { @@ -135,20 +137,52 @@ ## need rfd to inherit appropriate na.action; need grouping ## variables as well as any covariates that are included ## in RE terms + ## FIXME: mfnew is new data frame, rfd is processed new data + ## why do we need both/what is each doing/how do they differ? + ## rfd is *only* used in mkReTrms + ## mfnew is *only* used for its na.action attribute (!) [fixed only] + ## using model.frame would mess up matrix-valued predictors (GH #201) if (is.null(newdata)) { rfd <- mfnew <- model.frame(object) } else { - mfnew <- model.frame(delete.response(terms(object,fixed.only=TRUE)), - newdata, na.action=na.action) - ## make sure we pass na.action with new data - ## it would be nice to do something more principled like - ## rfd <- model.frame(~.,newdata,na.action=na.action) - ## but this adds complexities (stored terms, formula, etc.) - ## that mess things up later on ... - rfd <- na.action(newdata) - if (is.null(attr(rfd,"na.action"))) - attr(rfd,"na.action") <- na.action - } + mfnew <- model.frame(delete.response(terms(object,fixed.only=TRUE)), + newdata, na.action=na.action) + ## make sure we pass na.action with new data + ## it would be nice to do something more principled like + ## rfd <- model.frame(~.,newdata,na.action=na.action) + ## but this adds complexities (stored terms, formula, etc.) + ## that mess things up later on ... + ## rfd <- na.action(get_all_vars(delete.response(terms(object,fixed.only=FALSE)), newdata)) + old <- FALSE + if (old) { + rfd <- na.action(newdata) + if (is.null(attr(rfd,"na.action"))) + attr(rfd,"na.action") <- na.action + } else { + newdata.NA <- newdata + if (!is.null(fixed.na.action <- attr(mfnew,"na.action"))) { + newdata.NA <- newdata.NA[-fixed.na.action,] + } + tt <- delete.response(terms(object,random.only=TRUE)) + ## need to let NAs in RE components go through -- they're handled downstream + rfd <- model.frame(tt,newdata.NA,na.action=na.pass) + if (!is.null(fixed.na.action)) + attr(rfd,"na.action") <- fixed.na.action + } + ## + ## ## need terms to preserve info about spline/orthog polynomial bases + ## attr(rfd,"terms") <- terms(object) + ## ## ... but variables list messes things up; can we fix it? + ## vlist <- lapply(all.vars(terms(object)), as.name) + ## attr(attr(rfd,"terms"),"variables") <- as.call(c(quote(list), vlist)) + ## + ## take out variables that appear *only* in fixed effects + ## all.v <- all.vars(delete.response(terms(object,fixed.only=FALSE))) + ## ran.v <- vapply(findbars(formula(object)),all.vars,"") + ## fix.v <- all.vars(delete.response(terms(object,fixed.only=TRUE))) + ## rfd <- model.frame(delete.response(terms(object,fixed.only=FALSE)), + ## newdata,na.action=na.action) + } if (inherits(re.form, "formula")) { ## DROP values with NAs in fixed effects if (length(fit.na.action <- attr(mfnew,"na.action")) > 0) { @@ -325,6 +359,9 @@ xlev=orig_levs) X <- model.matrix(RHS, data=mfnew, contrasts.arg=attr(X,"contrasts")) + ## hack to remove unused interaction levels? + ## X <- X[,colnames(X0)] + offset <- 0 # rep(0, nrow(X)) tt <- terms(object) if (!is.null(off.num <- attr(tt, "offset"))) { @@ -485,13 +522,22 @@ ## compre.form <- noLHSform(formula(object)) ## construct RE formula ONLY: leave out fixed terms, ## which might have loose terms like offsets in them ... - fb <- findbars(formula(object)) - compReForm <- reformulate(vapply(fb, function(x) - paste("(",safeDeparse(x),")"), "")) + + ##' combine unary or binary operator + arguments (sugar for 'substitute') + makeOp <- function(x,y,op=NULL) { + if (is.null(op)) { ## unary + substitute(OP(X),list(X=x,OP=y)) + } else substitute(OP(X,Y), list(X=x,OP=op,Y=y)) + } + + ## this fails for complex RE terms such as (1|f/g): + ## findbars() is longer than the number of RE forms + compReForm <- reOnly(formula(object)) if (!noReForm(re.form)) { - rr <- re.form[[length(re.form)]] ## RHS of formula + rr <- reOnly(re.form)[[2]] ## expand RE and strip ~ ftemplate <- substitute(.~.-XX, list(XX=rr)) - compReForm <- update.formula(compReForm,ftemplate) + compReForm <- update.formula(compReForm,ftemplate)[-2] + ## update, then delete LHS } ## (1) random effect(s) @@ -499,9 +545,12 @@ newRE <- mkNewReTrms(object, newdata, compReForm, na.action=na.action, allow.new.levels=allow.new.levels) - ## paranoia ... - stopifnot(!is.null(newdata) || - isTRUE(all.equal(newRE$Lambdat,getME(object,"Lambdat")))) + ## this *can* justifiably happen, if we are using mkNewReTrms + ## in the context of predicting/simulating with a non-trivial + ## re.form ... + ## paranoia ... + ## stopifnot(!is.null(newdata) || + ## isTRUE(all.equal(newRE$Lambdat,getME(object,"Lambdat")))) U <- t(newRE$Lambdat %*% newRE$Zt) # == Z Lambda u <- rnorm(ncol(U)*nsim) ## UNSCALED random-effects contribution: @@ -517,12 +566,8 @@ ## GLMM ## n.b. DON'T scale random-effects (???) etasim <- etapred+sim.reff - family <- object@resp$family - ## hack (NB families have weird names) from @aosmith16 - if(gsub("[^[:alpha:]]", "", family$family) == "NegativeBinomial") { - family$family <- "negative.binomial" - } - musim <- family$linkinv(etasim) + family <- normalizeFamilyName(object@resp$family) + musim <- family$linkinv(etasim) #-> family$family == "negative.binomial" if(NB) ## ntot <- length(musim) ## FIXME: or could be dims["n"]? ## FIXME: is it possible to leverage family$simulate ... ??? ## @@ -707,15 +752,14 @@ ## in the original MASS version, .Theta is assigned into the environment ## (triggers a NOTE in R CMD check) ## modified from @aosmith16 GH contribution -negative.binomial_simfun <- function (object, nsim, ftd=fitted(object)) + +negative.binomial_simfun <- function (object, nsim, ftd = fitted(object)) { wts <- weights(object) if (any(wts != 1)) warning("ignoring prior weights") - theta <- as.numeric(gsub("[[:alpha:][:blank:]+?&/\\()]", "", - object@resp$family$family)) - rnbinom(nsim * length(ftd), mu = ftd, - size = theta) + theta <- getNBdisp(object) + rnbinom(nsim * length(ftd), mu = ftd, size = theta) } diff -Nru lme4-1.1-8/R/profile.R lme4-1.1-10/R/profile.R --- lme4-1.1-8/R/profile.R 2015-06-13 22:58:20.000000000 +0000 +++ lme4-1.1-10/R/profile.R 2015-08-23 16:36:36.000000000 +0000 @@ -294,7 +294,7 @@ forspl <- interpSpline(form, bres, na.action=na.omit)), error=function(e)e) if (inherits(bakspl, "error")) - warning("non-monotonic profile") + warning("non-monotonic profile for ",pname) ## return: namedList(bres,bakspl,forspl) # namedList() -> lmerControl.R @@ -361,7 +361,7 @@ fv <- ores$fval sig <- sqrt((rr$wrss() + pp1$sqrL(1))/n) c(sign(fw - est) * sqrt(fv - base), - Cv_to_Sv(ores$par, vapply(fitted@cnms,length, 1), s=sig), + Cv_to_Sv(ores$par, lengths(fitted@cnms), s=sig), ## ores$par * sig, sig, mkpar(p, j, fw, pp1$beta(1))) } @@ -406,7 +406,7 @@ ##' tn <- names(getME(fm1,"theta")) ##' nn <- c(tn,names(fixef(fm1))) ##' get.which("theta_",length(tn),length(nn),nn, verbose=TRUE) -##' +##' get.which <- function(which, nvp, nptot, parnames, verbose=FALSE) { if (is.null(which)) seq_len(nptot) @@ -435,19 +435,20 @@ ## classes .modelMatrixDrop <- function(mm, w) { if (isS4(mm)) { - ll <- list(Class = class(mm), - assign = attr(mm,"assign")[-w], - contrasts = NULL) - ## FIXME: where did the contrasts information go?? - ## mm@contrasts) - X <- mm[, -w, drop = FALSE] - ll <- c(ll, lapply(structure(slotNames(X), .Names=slotNames(X)), - function(nm) slot(X, nm))) - return(do.call("new", ll)) + nX <- slotNames(X <- mm[, -w, drop = FALSE]) + do.call(new, + c(list(Class = class(mm), + assign = attr(mm,"assign")[-w], + contrasts = NULL + ## FIXME: where did the contrasts information go?? + ## mm@contrasts + ), + lapply(structure(nX, .Names=nX), + function(nm) slot(X, nm)))) + } else { + structure(mm[, -w, drop=FALSE], + assign = attr(mm, "assign")[-w]) } - ans <- mm[, -w, drop=FALSE] - attr(ans, "assign") <- attr(mm, "assign")[-w] - ans } ## The deviance is profiled with respect to the fixed-effects @@ -468,7 +469,7 @@ stopifnot(is(fm, "merMod")) fm <- refitML(fm) basedev <- -2*c(logLik(fm)) ## no longer deviance() - vlist <- sapply(fm@cnms,length) + vlist <- lengths(fm@cnms) sig <- sigma(fm) ## only if useSc=TRUE? stdErr <- unname(coef(summary(fm))[,2]) pp <- fm@pp$copy() @@ -692,7 +693,11 @@ ##' confint() method for our profile() results 'thpr' ##' @importFrom stats confint -confint.thpr <- function(object, parm, level = 0.95, zeta, ...) +confint.thpr <- function(object, parm, level = 0.95, zeta, + ## tolerance for non-monotonic profiles + ## (raw values, not splines) + non.mono.tol=1e-2, + ...) { bak <- attr(object, "backward") ## fallback strategy for old profiles that don't have a lower/upper @@ -701,6 +706,8 @@ lower <- rep(NA,length(parm)) if (is.null(upper <- attr(object,"upper"))) upper <- rep(NA,length(parm)) + ## FIXME: work a little harder to add -Inf/Inf for fixed effect + ## parameters? (Should only matter for really messed-up profiles) bnms <- names(bak) if (missing(parm)) parm <- bnms else if (is.numeric(parm)) parm <- bnms[parm] @@ -719,12 +726,28 @@ ## predy is used in many places and it's much harder to ## tell in general whether an NA indicates a lower or an ## upper bound ... + badprof <- FALSE + p <- rep(NA,2) if (!inherits(b <- bak[[parm[i]]], "error")) { p <- predy(b, zeta) + } else { + obj1 <- object[object$.par==parm[[i]],c(parm[[i]],".zeta")] + if (all(is.na(obj1[,2]))) { + badprof <- TRUE + warning("bad profile for ",parm[i]) + } else if (min(diff(obj1[,2])<(-non.mono.tol),na.rm=TRUE)) { + badprof <- TRUE + warning("non-monotonic profile for ",parm[i]) + } else { + warning("bad spline fit for ",parm[i],": falling back to linear interpolation") + p <- approxfun(obj1[,2],obj1[,1])(zeta) + } + } + if (!badprof) { if (is.na(p[1])) p[1] <- lower[i] if (is.na(p[2])) p[2] <- upper[i] - ci[i,] <- p } + ci[i,] <- p } ci } @@ -798,7 +821,7 @@ { bootFun <- function(x) { th <- x@theta - nvec <- vapply(x@cnms, length, 1L) + nvec <- lengths(x@cnms) scaleTh <- (isLMM(x) || isNLMM(x)) useSc <- as.logical(x@devcomp$dims[["useSc"]]) ## FIXME: still ugly. Best cleanup via Cv_to_Sv ... diff -Nru lme4-1.1-8/R/utilities.R lme4-1.1-10/R/utilities.R --- lme4-1.1-8/R/utilities.R 2015-06-06 17:12:42.000000000 +0000 +++ lme4-1.1-10/R/utilities.R 2015-09-08 23:41:51.000000000 +0000 @@ -1,11 +1,15 @@ -if((Rv <- getRversion()) < "3.1.0") { - anyNA <- function(x) any(is.na(x)) - if(Rv < "3.0.0") { - rep_len <- function(x, length.out) rep(x, length.out=length.out) - if(Rv < "2.15") - paste0 <- function(...) paste(..., sep = '') - } -}; rm(Rv) +if((Rv <- getRversion()) < "3.2.0") { + lengths <- function (x, use.names = TRUE) vapply(x, length, 1L, USE.NAMES = use.names) + if(Rv < "3.1.0") { + anyNA <- function(x) any(is.na(x)) + if(Rv < "3.0.0") { + rep_len <- function(x, length.out) rep(x, length.out=length.out) + if(Rv < "2.15") + paste0 <- function(...) paste(..., sep = '') + } + } +} ## R < 3.2.0 +rm(Rv) ## From Matrix package isDiagonal(.) : all0 <- function(x) !anyNA(x) && all(!x) @@ -122,8 +126,9 @@ ## order terms stably by decreasing number of levels in the factor if (any(diff(nl) > 0)) { ord <- rev(order(nl)) - blist <- blist[ord] - nl <- nl[ord] + blist <- blist [ord] + nl <- nl [ord] + term.names <- term.names[ord] } Ztlist <- lapply(blist, `[[`, "sm") Zt <- do.call(rBind, Ztlist) ## eq. 7, JSS lmer paper @@ -134,7 +139,7 @@ ## any potential reordering of the terms. cnms <- lapply(blist, `[[`, "cnms") # list of column names of the # model matrix per term - nc <- vapply(cnms, length, 0L) # no. of columns per term + nc <- lengths(cnms) # no. of columns per term # (in lmer jss: p_i) nth <- as.integer((nc * (nc+1))/2) # no. of parameters per term # (in lmer jss: ??) @@ -161,7 +166,7 @@ ltri <- lower.tri(dd, diag = TRUE) ii <- row(dd)[ltri] jj <- col(dd)[ltri] - dd[cbind(ii, jj)] <- seq_along(ii) # FIXME: this line unnecessary? + ## unused: dd[cbind(ii, jj)] <- seq_along(ii) data.frame(i = as.vector(mm[, ii]) + boff[i], j = as.vector(mm[, jj]) + boff[i], x = as.double(rep.int(seq_along(ii), @@ -169,8 +174,8 @@ thoff[i])) })))) thet <- numeric(sum(nth)) - ll <- list(Zt=Matrix::drop0(Zt), theta=thet, Lind=as.integer(Lambdat@x), - Gp=unname(c(0L, cumsum(nb)))) + ll <- list(Zt = drop0(Zt), theta = thet, Lind = as.integer(Lambdat@x), + Gp = unname(c(0L, cumsum(nb)))) ## lower bounds on theta elements are 0 if on diagonal, else -Inf ll$lower <- -Inf * (thet + 1) ll$lower[unique(diag(Lambdat))] <- 0 @@ -576,16 +581,16 @@ sbnm(form) } -## Check for a constant term (a literal 1) in an expression +##' Check for a constant term (a literal 1) in an expression ## -## In the mixed-effects part of a nonlinear model formula, a constant -## term is not meaningful because every term must be relative to a -## nonlinear model parameter. This function recursively checks the -## expressions in the formula for a a constant, calling stop() if -## such a term is encountered. -## @title Check for constant terms. -## @param expr an expression -## @return NULL. The function is executed for its side effect. +##' In the mixed-effects part of a nonlinear model formula, a constant +##' term is not meaningful because every term must be relative to a +##' nonlinear model parameter. This function recursively checks the +##' expressions in the formula for a a constant, calling stop() if +##' such a term is encountered. +##' @title Check for constant terms. +##' @param expr an expression +##' @return NULL. The function is executed for its side effect. chck1 <- function(expr) { if ((le <- length(expr)) == 1) { if (is.numeric(expr) && expr == 1) @@ -602,7 +607,7 @@ start <- eval(mc$start, parent.frame(2L)) if (is.numeric(start)) start <- list(nlpars = start) stopifnot(is.numeric(nlpars <- start$nlpars), - vapply(nlpars, length, 0L) == 1L, + lengths(nlpars) == 1L, length(pnames <- names(nlpars)) == length(nlpars), length(form <- as.formula(mc$formula)) == 3L, is(nlform <- eval(form[[2]]), "formula"), @@ -696,11 +701,35 @@ ## } ################################################################################ -.minimalOptinfo <- function() { - list(conv = list( - opt = 0L, - lme4 = list(messages = character(0)))) -} +.minimalOptinfo <- function() + list(conv = list(opt = 0L, + lme4 = list(messages = character(0)))) + +.optinfo <- function(opt, lme4conv=NULL) + list(optimizer = attr(opt, "optimizer"), + control = attr(opt, "control"), + derivs = attr(opt, "derivs"), + conv = list(opt = opt$conv, lme4 = lme4conv), + feval = if (is.null(opt$feval)) NA else opt$feval, + warnings = attr(opt, "warnings"), + val = opt$par) + +##' Potentically needed in more than one place, be sure to keep consistency! +##' hack (NB families have weird names) from @aosmith16; then corrected +isNBfamily <- function(familyString) + grepl("^Negative ?Binomial", familyString, ignore.case=TRUE) +normalizeFamilyName <- function(family) { # such as object@resp$family + if(isNBfamily(family$family)) + family$family <- "negative.binomial" + family +} + +##' Is it a family with no scale parameter +hasNoScale <- function(family) + any(substr(family$family, 1L, 12L) + == c("poisson", "binomial", "negative.bin", "Negative Bin")) + + ##--> ../man/mkMerMod.Rd ---Create a merMod object ##' @param rho the environment of the objective function @@ -712,24 +741,24 @@ is(pp <- rho$pp, "merPredD"), is(resp <- rho$resp, "lmResp"), is.list(opt), "par" %in% names(opt), - c("conv","fval") %in% substr(names(opt),1,4), ## "conv[ergence]", "fval[ues]" + c("conv", "fval") %in% substr(names(opt),1,4), ## "conv[ergence]", "fval[ues]" is.list(reTrms), c("flist", "cnms", "Gp", "lower") %in% names(reTrms), length(rcl <- class(resp)) == 1) n <- nrow(pp$V) p <- ncol(pp$V) - dims <- c(N=nrow(pp$X), n=n, p=p, nmp=n-p, - nth=length(pp$theta), q=nrow(pp$Zt), - nAGQ=rho$nAGQ, + isGLMM <- (rcl == "glmResp") + dims <- c(N = nrow(pp$X), n=n, p=p, nmp = n-p, q = nrow(pp$Zt), + nth = length(pp$theta), + nAGQ= rho$nAGQ, compDev=rho$compDev, ## 'use scale' in the sense of whether dispersion parameter should ## be reported/used (*not* whether theta should be scaled by sigma) - useSc=(rcl != "glmResp" || - !resp$family$family %in% c("poisson","binomial")), + useSc = !(isGLMM && hasNoScale(resp$family)), reTrms=length(reTrms$cnms), - spFe=0L, - REML=if (rcl=="lmerResp") resp$REML else 0L, - GLMM=(rcl=="glmResp"), - NLMM=(rcl=="nlsResp")) + spFe= 0L, + REML = if (rcl=="lmerResp") resp$REML else 0L, + GLMM= isGLMM, + NLMM= (rcl=="nlsResp")) storage.mode(dims) <- "integer" fac <- as.numeric(rcl != "nlsResp") if (trivial.y <- (length(resp$y)==0)) { @@ -743,7 +772,7 @@ ## weights <- resp$weights beta <- pp$beta(fac) ## rescale - if (!is.null(sc <- attr(pp$X,"scaled:scale"))) { + if (!is.null(sc <- attr(pp$X, "scaled:scale"))) { warning("auto(un)scaling not yet finished/tested") ## FIXME: test/handle no-intercept models ## (only need to worry if we do centering as well as scaling) @@ -754,7 +783,7 @@ beta2[names(sc)] <- sc*beta2[names(sc)] beta <- beta2 } - if (!is.null(attr(pp$X,"scaled:center"))) { + if (!is.null(attr(pp$X, "scaled:center"))) { warning("auto(un)centering not yet implemented") } #sigmaML <- pwrss/sum(weights) @@ -776,19 +805,13 @@ tolPwrss=rho$tolPwrss) ## TODO: improve this hack to get something in frame slot (maybe need weights, etc...) if(missing(fr)) fr <- data.frame(resp$y) - new(switch(rcl, lmerResp="lmerMod", glmResp="glmerMod", nlsResp="nlmerMod"), + new(switch(rcl, lmerResp = "lmerMod", glmResp = "glmerMod", nlsResp = "nlmerMod"), call=mc, frame=fr, flist=reTrms$flist, cnms=reTrms$cnms, Gp=reTrms$Gp, theta=pp$theta, beta=beta, u=if (trivial.y) rep(NA_real_,nrow(pp$Zt)) else pp$u(fac), lower=reTrms$lower, devcomp=list(cmp=cmp, dims=dims), pp=pp, resp=resp, - optinfo = list (optimizer= attr(opt,"optimizer"), - control = attr(opt,"control"), - derivs = attr(opt,"derivs"), - conv = list(opt=opt$conv, lme4=lme4conv), - feval = if (is.null(opt$feval)) NA else opt$feval, - warnings = attr(opt,"warnings"), val = opt$par) - ) + optinfo = .optinfo(opt, lme4conv)) }## {mkMerMod} ## generic argument checking @@ -808,9 +831,11 @@ ## (different meanings/hints depending on glmer vs lmer) if (!is.null(l...[["method"]])) { msg <- paste("Argument", sQuote("method"), "is deprecated.") - if (type=="lmer") msg <- paste(msg,"Use the REML argument to specify ML or REML estimation.") - if (type=="glmer") msg <- paste(msg,"Use the nAGQ argument to specify Laplace (nAGQ=1) or adaptive", - "Gauss-Hermite quadrature (nAGQ>1). PQL is no longer available.") + if (type == "lmer") + msg <- paste(msg, "Use the REML argument to specify ML or REML estimation.") + else if (type == "glmer") + msg <- paste(msg, "Use the nAGQ argument to specify Laplace (nAGQ=1) or adaptive", + "Gauss-Hermite quadrature (nAGQ>1). PQL is no longer available.") warning(msg,call.=FALSE) l... <- l...[names(l...) != "method"] } @@ -892,7 +917,7 @@ OK <- allvarex(parent.frame(i)) cat("vars exist in parent frame ", i) if (i == glEnv) cat(" (global)") - cat(" ",OK,"\n") + cat(" ",OK, "\n") } cat("vars exist in env of formula ", allvarex(denv), "\n") } ## if (debug) @@ -941,8 +966,8 @@ ##' @return Sparse covariance matrix condVar <- function(object) { s2 <- sigma(object)^2 - Lamt <- getME(object,"Lambdat") - L <- getME(object,"L") + Lamt <- getME(object, "Lambdat") + L <- getME(object, "L") ## never do it this way! fortune("SOOOO") #V <- solve(L, system = "A") @@ -1069,38 +1094,6 @@ defaultControl <- list(algorithm="NLOPT_LN_BOBYQA", xtol_abs=1e-6, ftol_abs=1e-6, maxeval=1e5) function(par,fn,lower,upper,control=list(),...) { - if(packageVersion("nloptr") < "1.0.2") { - ## kluge: we would like to import nloptr.default.options from nloptr - ## up front (it's required by nloptr and isn't accessible if - ## nloptr::nloptr is only imported without the package being loaded) - ## (1) adding export(nloptr.default.options) to the NAMESPACE - ## of nloptr doesn't work (package install fails) - ## (2) the following code puts a copy of nloptr.default.options into the - ## environment of nloptwrap(), but nloptr can't see it - ## from there when called from within nloptwrap ... - ## if (!exists("nloptr.default.options")) { - ## data("nloptr.default.options", - ## package="nloptr", - ## envir=environment()) } - ## (3) solution used here is to load it *into the global environment* - ## every time nloptwrap() is - ## called (if it doesn't exist) -- ugly but works. - ## ... but it provokes a complaint from R CMD check - ## ... loading it into the environment of nloptwrap, or - ## the environment within nloptwrap (i.e. envir=environment() - ## or envir=parent.env(environment()) doesn't work because - ## nloptr doesn't look there for it. - ## ... loading it into the environment of nloptr doesn't work - ## because the environment is locked (and I don't know if/how - ## to unlock it: ?unlockBinding - ## http://stackoverflow.com/questions/19132492/how-to-unlock-environment-in-r - ## https://gist.github.com/wch/3280369#file-unlockenvironment-r - ## this seems like a lot of effort to go to - ## try to deceive R CMD check - ee <- .GlobalEnv ## environment(nloptr) - if (!exists("nloptr.default.options", envir=ee)) - data("nloptr.default.options", package="nloptr", envir=ee) - } for (n in names(defaultControl)) if (is.null(control[[n]])) control[[n]] <- defaultControl[[n]] res <- nloptr(x0=par, eval_f=fn, lb=lower,ub=upper, opts=control, ...) diff -Nru lme4-1.1-8/README.md lme4-1.1-10/README.md --- lme4-1.1-8/README.md 2015-06-19 21:51:50.000000000 +0000 +++ lme4-1.1-10/README.md 2015-08-23 16:57:14.000000000 +0000 @@ -7,11 +7,10 @@ ## Recent/release notes -* We have submitted release 1.1-8 to CRAN. There are no major user-visible changes. - * We have fixed some bugs in `predict`, `simulate`, and `refit`. - * Convergence and positive-definite-Hessian warnings are still overly conservative for large (>10^4 rows) data sets, but we are holding off on changing anything until we really understand the problem; see `help("convergence")`. - * The deviance computation has changed for GLMMs; see "Deviance and log-likelihood of GLMMs" in [merMod-class.Rd](https://github.com/lme4/lme4/blob/6203f71f4f6aa75e3a69f08c40e5d2fc176610d6/man/merMod-class.Rd) -* Otherwise, see the [NEWS file](https://github.com/lme4/lme4/blob/master/inst/NEWS.Rd) (or `news(Version=="1.1.8",package="lme4")`). +* Release 1.1-9 is on CRAN now. There are no major user-visible changes. + * Some improvements to \code{glmer.nb}, including the ability to simulate; note that \code{glmer.nb} is still somewhat unstable, e.g. for models with additional arguments such as \code{offset} or \code{contrasts}. + * Bug fixes to \code{predict} for models with orthogonal polynomial/spline bases, and to \code{lmList} models with offsets and weights +* Otherwise, see the [NEWS file](https://github.com/lme4/lme4/blob/master/inst/NEWS.Rd) (or `news(Version=="1.1.9",package="lme4")`). ## Features diff -Nru lme4-1.1-8/tests/boundary.R lme4-1.1-10/tests/boundary.R --- lme4-1.1-8/tests/boundary.R 2014-11-13 20:34:53.000000000 +0000 +++ lme4-1.1-10/tests/boundary.R 2015-08-19 12:47:37.000000000 +0000 @@ -9,55 +9,58 @@ if(!dev.interactive(orNone=TRUE)) pdf("boundary_plots.pdf") ## Stephane Laurent: -dat <- read.csv(system.file("testdata","dat20101314.csv",package="lme4")) +dat <- read.csv(system.file("testdata","dat20101314.csv", package="lme4")) -fit <- lmer(y ~ (1|Operator)+(1|Part)+(1|Part:Operator), data=dat, - control=lmerControl(optimizer="Nelder_Mead")) +fit <- lmer(y ~ (1|Operator)+(1|Part)+(1|Part:Operator), data=dat, + control= lmerControl(optimizer="Nelder_Mead")) fit_b <- lmer(y ~ (1|Operator)+(1|Part)+(1|Part:Operator), data=dat, - control=lmerControl(optimizer="bobyqa",restart_edge=FALSE)) + control= lmerControl(optimizer="bobyqa", restart_edge=FALSE)) fit_c <- lmer(y ~ (1|Operator)+(1|Part)+(1|Part:Operator), data=dat, - control=lmerControl(optimizer="Nelder_Mead", restart_edge=FALSE, - check.conv.hess="ignore")) + control= lmerControl(optimizer="Nelder_Mead", restart_edge=FALSE, + check.conv.hess="ignore")) ## final fit gives degenerate-Hessian warning ## FIXME: use fit_c with expect_warning() as a check on convergence tests ## tolerance=1e-5 seems OK in interactive use but not in R CMD check ... ?? -stopifnot(all.equal(getME(fit,"theta"),getME(fit_b,"theta"), tolerance=2e-5)) -stopifnot(all(getME(fit,"theta")>0)) +stopifnot(all.equal(getME(fit, "theta") -> th.f, + getME(fit_b,"theta"), tolerance=2e-5), + all(th.f > 0)) ## Manuel Koller -source(system.file("testdata","koller-data.R",package="lme4")) +source(system.file("testdata", "koller-data.R", package="lme4")) ldata <- getData(13) ## old (backward compatible/buggy) fm4 <- lmer(y ~ (1|Var2), ldata, control=lmerControl(optimizer="Nelder_Mead", use.last.params=TRUE), start=list(theta=1)) -fm4b <- lmer(y ~ (1|Var2), ldata, control=lmerControl(optimizer="Nelder_Mead", - use.last.params=TRUE, restart_edge=FALSE, - check.conv.hess="ignore", - check.conv.grad="ignore"), - start=list(theta=1)) +fm4b <- lmer(y ~ (1|Var2), ldata, + control = lmerControl(optimizer="Nelder_Mead", use.last.params=TRUE, + restart_edge = FALSE, + check.conv.hess="ignore", check.conv.grad="ignore"), + start = list(theta=1)) ## FIXME: use as convergence test check -stopifnot(getME(fm4b,"theta")==0) +stopifnot(getME(fm4b,"theta") == 0) fm4c <- lmer(y ~ (1|Var2), ldata, control=lmerControl(optimizer="bobyqa", - use.last.params=TRUE), + use.last.params=TRUE), start=list(theta=1)) -stopifnot(all.equal(getME(fm4,"theta"),getME(fm4c,"theta"),tolerance=1e-4)) -stopifnot(getME(fm4,"theta") > 0) +stopifnot(all.equal(getME(fm4, "theta") -> th4, + getME(fm4c,"theta"), tolerance=1e-4), + th4 > 0) + ## new: doesn't get stuck at edge any more, but gets stuck somewhere else ... fm5 <- lmer(y ~ (1|Var2), ldata, control=lmerControl(optimizer="Nelder_Mead", - check.conv.hess="ignore", - check.conv.grad="ignore"), - start=list(theta=1)) + check.conv.hess="ignore", + check.conv.grad="ignore"), + start=list(theta=1)) fm5b <- lmer(y ~ (1|Var2), ldata, control=lmerControl(optimizer="Nelder_Mead", - restart_edge=FALSE, - check.conv.hess="ignore", - check.conv.grad="ignore"), - start=list(theta=1)) + restart_edge=FALSE, + check.conv.hess="ignore", + check.conv.grad="ignore"), + start = list(theta = 1)) fm5c <- lmer(y ~ (1|Var2), ldata, control=lmerControl(optimizer="bobyqa"), - start=list(theta=1)) + start = list(theta = 1)) stopifnot(all.equal(unname(getME(fm5c,"theta")), 0.21067645, tolerance = 1e-7)) # 0.21067644264 [64-bit, lynne] diff -Nru lme4-1.1-8/tests/getME.R lme4-1.1-10/tests/getME.R --- lme4-1.1-8/tests/getME.R 2014-08-19 04:22:07.000000000 +0000 +++ lme4-1.1-10/tests/getME.R 2015-08-20 12:25:28.000000000 +0000 @@ -29,7 +29,7 @@ "theta") ## internal consistency check ensuring that all allowed 'name's work (and are not empty): -(nmME <- eval(formals(getME)$name)) +(nmME <- eval(formals(lme4:::getME.merMod)$name)) chkMEs <- function(fm, nms) { stopifnot(is.character(nms)) str(parts <- sapply(nms, getME, object = fm, simplify=FALSE)) diff -Nru lme4-1.1-8/tests/lmList-tst.R lme4-1.1-10/tests/lmList-tst.R --- lme4-1.1-8/tests/lmList-tst.R 2015-04-13 00:44:21.000000000 +0000 +++ lme4-1.1-10/tests/lmList-tst.R 2015-10-05 16:12:45.000000000 +0000 @@ -139,3 +139,18 @@ ## $ logLik : chr "log-likelihood not available with NULL fits" ## $ summary: chr "subscript out of bounds" stopifnot(length(errs4) <= 2) + +## from GH #320 +x <- c(1,2,3,4,5,6,7,8) +y <- c(2,2,5,4,3,1,2,1) +g <- c(1,1,1,2,2,3,3,3) +dat <- data.frame(x=x, y=y, g=g) +m1 <- lmList(y ~ x | g, data=dat) +stopifnot(!any(is.na(coef(m1)))) +m2 <- lmList(Reaction ~ Days | Subject, + weights=runif(nrow(sleepstudy)), sleepstudy) +m3 <- lmList(Reaction ~ Days | Subject, sleepstudy) +m4 <- lmList(Reaction ~ Days | Subject, + offset=runif(nrow(sleepstudy)), sleepstudy) +stopifnot(!identical(m2,m3)) +stopifnot(!identical(m4,m3)) diff -Nru lme4-1.1-8/tests/nbinom.R lme4-1.1-10/tests/nbinom.R --- lme4-1.1-8/tests/nbinom.R 2014-11-13 20:34:53.000000000 +0000 +++ lme4-1.1-10/tests/nbinom.R 2015-08-19 12:47:37.000000000 +0000 @@ -2,8 +2,8 @@ cat("lme4 testing level: ", testLevel <- lme4:::testLevel(), "\n") +getNBdisp <- function(x) getME(x,"glmer.nb.theta") ## for now, use hidden functions [MM: this is a sign, we should *export* them] -getNBdisp <- lme4:::getNBdisp refitNB <- lme4:::refitNB simfun <- function(sd.u=1, NBtheta=0.5, @@ -24,17 +24,15 @@ set.seed(102) d.1 <- simfun() t1 <- system.time(g1 <- glmer.nb(z ~ x + (1|f), data=d.1, verbose=TRUE)) - ## no longer: Model failed to converge with max|grad| = 0.00378242 (tol = 0.001) - g1 - ## ^^ FIXME: 'Data:' shows '..2' ![eval.parent() etc.?] - + d1 <- getNBdisp(g1) (g1B <- refitNB(g1, theta = d1)) (ddev <- deviance(g1) - deviance(g1B)) (reld <- (fixef(g1) - fixef(g1B)) / fixef(g1)) stopifnot(abs(ddev) < 1e-4, # was 6.18e-7 now 1.045e-6, now -6.367e-5 (!) abs(reld) < 1e-4)# 0, then 4.63e-6 + ## 2 Aug 2015: ddev==reld==0 on 32-bit Ubuntu 12.04 ## library(glmmADMB) @@ -60,9 +58,12 @@ ## no more at all ??! ## no all.equal( d1, glmmADMB_vals$theta, tolerance=0.0016) ## , - all.equal(fixef(g1B), glmmADMB_vals$ fixef, tolerance=0.1)# was 0.01 ! + all.equal(fixef(g1B), glmmADMB_vals$fixef, tolerance=0.1)# was 0.01 ! + ## Ubuntu 12.04/32-bit: 0.0094 , - all.equal(logLik.m(g1B), -glmmADMB_vals$ NLL, tolerance=0.4)# was 0.001 (!!) + all.equal(c(logLik.m(g1B)), c(-glmmADMB_vals$NLL), tolerance=0.4)# was 0.001 (!!) + ## Ubuntu 12.04/32-bit: 1.61e-5 + ## except that df=3 vs 4 ) }## end if( testLevel > 1 ) @@ -134,6 +135,7 @@ ## , all.equal(fixef (g4), glmmADMB_epil_vals$ fixef, tolerance= 0.03)#was 0.004 ## , +### still 0.00374 on Ubuntu 12.04 ## FIXME: even df differ ! ## all.equal(logLik.m (g4), - glmmADMB_epil_vals$ NLL, tolerance= 0.0) ## was 0.0002 ) diff -Nru lme4-1.1-8/tests/testOptControl.R lme4-1.1-10/tests/testOptControl.R --- lme4-1.1-8/tests/testOptControl.R 2014-09-17 18:24:43.000000000 +0000 +++ lme4-1.1-10/tests/testOptControl.R 2015-08-19 12:47:37.000000000 +0000 @@ -7,6 +7,15 @@ control= lmerControl("NMcopy", optCtrl= list(iprint=20)))) + ## check that printing goes through step 140 twice and up to 240 once -findStep <- function(str,n) sum(grepl(paste0("^\\(NM\\) ",n,": "),cc)) -stopifnot(findStep(cc,140)==2 && findStep(cc,240)==1) +## findStep <- function(str,n) sum(grepl(paste0("\\(NM\\) ",n,": "),str)) +cc <- paste(cc,collapse="") +countStep <- function(str,n) { + length(gregexpr(paste0("\\(NM\\) ",n,": "),str)[[1]]) +} + +stopifnot(countStep(cc,140)==2 && countStep(cc,240)==1) + +## testStr <- +## "(NM) 20: f = -53.3709 at 0.706667 0.813333 1.46444(NM) 40: f = -147.132 at 0 0 19.18(NM) 60: f = -147.159 at 0 0 17.4275(NM) 80: f = -147.159 at 0 0 17.5615(NM) 100: f = -147.159 at 0 0 17.5754(NM) 120: f = -147.159 at 0 0 17.5769(NM) 140: f = -147.159 at 0 0 17.5768(NM) 20: f = -165.55 at 0.0933333 0.573333 17.3168(NM) 40: f = -173.704 at 0.23799 1.4697 16.9728(NM) 60: f = -173.849 at 0.449634 1.39998 16.9452(NM) 80: f = -174.421 at 0.52329 1.69123 18.1534(NM) 100: f = -176.747 at 0.762043 1.88271 32.8993(NM) 120: f = -176.839 at 0.751206 1.75371 37.2128(NM) 140: f = -176.853 at 0.706425 1.7307 35.7528(NM) 160: f = -176.853 at 0.710803 1.73476 35.7032(NM) 180: f = -176.853 at 0.710159 1.73449 35.6699(NM) 200: f = -176.853 at 0.710271 1.73461 35.6689(NM) 220: f = -176.853 at 0.710259 1.7346 35.6684(NM) 240: f = -176.853 at 0.710257 1.73459 35.6685Linear mixed model fit by REML ['lmerMod']"