R version 4.4.0 beta (2024-04-14 r86421 ucrt) -- "Puppy Cup"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> #### Testing  refit()
> #### ----------------
> 
> library(lme4)
Loading required package: Matrix
> set.seed(101)
> testLevel <- if (nzchar(s <- Sys.getenv("LME4_TEST_LEVEL"))) as.numeric(s) else 1
> 
> ## for each type of model, should be able to
> ##  (1) refit with same data and get the same answer,
> ##     at least structurally (small numerical differences
> ##     are probably unavoidable)
> ##  (2) refit with simulate()d data
> 
> if (testLevel>1) {
+ getinfo <- function(x) {
+   c(fixef(x), logLik(x), unlist(ranef(x)), unlist(VarCorr(x)))
+ }
+ 
+ dropterms <- function(x) {
+     attr(x@frame,"terms") <- NULL
+     x
+ }
+ 
+ if (getRversion() >= "3.0.0") {
+     attach(system.file("testdata", "lme-tst-fits.rda", package="lme4"))
+ } else {
+     ## saved fits are not safe with old R versions; just re-compute ("cheat"!):
+ 
+     fit_sleepstudy_2 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
+ 
+     cbpp$obs <- factor(seq(nrow(cbpp)))
+     ## intercept-only fixed effect
+     fit_cbpp_0 <- glmer(cbind(incidence, size-incidence) ~ 1 + (1|herd),
+                         cbpp, family=binomial)
+     ## include fixed effect of period
+     fit_cbpp_1 <- update(fit_cbpp_0, . ~ . + period)
+     if(FALSE) ## include observation-level RE
+         fit_cbpp_2 <- update(fit_cbpp_1, . ~ . + (1|obs))
+     ## specify formula by proportion/weights instead
+     fit_cbpp_3 <- update(fit_cbpp_1, incidence/size ~ period + (1 | herd), weights = size)
+ }
+ 
+ 
+ ## LMM
+ fm1 <- fit_sleepstudy_2
+ fm1R <- refit(fm1, sleepstudy$Reaction)
+ fm1S <- refit(fm1, simulate(fm1)[[1]])
+ 
+ stopifnot(all.equal(getinfo(fm1 ),
+                     getinfo(fm1R), tolerance = 6e-3),
+           all.equal(getinfo(fm1 ),
+                     getinfo(fm1S), tolerance = 0.5) # <- simulate()
+           )
+ 
+ if(FALSE) { ## show all differences
+     sapply(slotNames(fm1), function(.)
+         all.equal( slot(fm1,.), slot(fm1R,.), tolerance=0))
+ }
+ 
+ if (getRversion() >= "3.4.0") {
+     ## differences: FALSE for resp, theta, u, devcomp, pp, optinfo?
+     ## FIXME: this isn't actually tested in any way ...
+     sapply(slotNames(fm1),
+            function(.) isTRUE(all.equal( slot(fm1,.), slot(fm1R,.), tolerance= 1.5e-5)))
+     str(fm1 @ optinfo)
+     str(fm1R@ optinfo)
+ }
+ 
+ fm1ML <- refitML(fm1)
+ stopifnot(
+     all.equal(getinfo(fm1), getinfo(fm1ML), tolerance=0.05)# 0.029998
+ )
+ 
+ 
+ ## binomial GLMM (two-column)
+ gm1 <- fit_cbpp_1
+ gm1R <- refit(gm1, with(cbpp, cbind(incidence,size-incidence)))
+ 
+  sim1Z <- simulate(gm1)[[1]]
+  sim1Z[4,] <- c(0,0)
+  (gm1. <- refit(gm1, sim1Z)) # earlier gave Error:  ... PIRLS ... failed ...
+ 
+ all.equal(getinfo(gm1), getinfo(gm1R), tolerance=0) # to see it --> 5.52e-4
+ # because glmer() uses Laplace approx. (? -- still, have *same* y !)
+ stopifnot(all.equal(getinfo(gm1), getinfo(gm1R), tolerance = 1e-4))
+ 
+ gm1S <- refit(gm1, simulate(gm1)[[1]])
+           all.equal(getinfo(gm1), getinfo(gm1S), tolerance=0) # to see:
+ stopifnot(all.equal(getinfo(gm1), getinfo(gm1S), tolerance = 0.4))
+ 
+ ## binomial GLMM (prob/weights)
+ formula(gm2 <- fit_cbpp_3)
+ ## glmer(incidence/size ~ period + (1 | herd), cbpp, binomial, weights=size)
+ gm2R <- refit(gm2,  with(cbpp, incidence/size))
+           all.equal(getinfo(gm2), getinfo(gm2R), tolerance= 0)
+ stopifnot(all.equal(getinfo(gm2), getinfo(gm2R), tolerance= 6e-4))
+ 
+ ## FIXME: check on Windows == 2015-06: be brave
+ gm2S <- refit(gm2, simulate(gm2)[[1]])
+           all.equal(getinfo(gm2), getinfo(gm2S), tolerance=0)# 0.17 .. upto 0.28
+ stopifnot(all.equal(getinfo(gm2), getinfo(gm2S), tolerance=0.40))
+ 
+ ## from Alexandra Kuznetsova
+ set.seed(101)
+ Y <- matrix(rnorm(1000),ncol=2)
+ d <- data.frame(y1=Y[,1],  x=rnorm(100), f=rep(1:10,10))
+ fit1 <- lmer(y1 ~ x+(1|f),data=d)
+ fit2 <- refit(fit1, newresp = Y[,2], rename.response=TRUE)
+ ## check, but ignore terms attribute of model frame ...
+ tools::assertWarning(refit(fit1, newresp = Y[,2], junk=TRUE))
+ if (isTRUE(all.equal(fit1,fit2))) stop("fit1 and fit2 should not be equal")
+ ## hack number of function evaluations
+ u2 <- update(fit2)
+ fit2@optinfo$feval <- u2@optinfo$feval <-  NA
+ 
+ d1 <- dropterms(fit2)
+ d2 <- dropterms( u2 )
+ ## They are not "all equal", but mostly :
+ for (i in slotNames(d1)) {
+     ae <- all.equal(slot(d1,i), slot(d2,i))
+     cat(sprintf("%10s: %s\n", i, if(isTRUE(ae)) "all.equal"
+ 		else paste(ae, collapse="\n  ")))
+ }
+           all.equal(getinfo(d1), getinfo(d2),  tolerance = 0)# -> 0.00126
+ stopifnot(all.equal(getinfo(d1), getinfo(d2),  tolerance = 0.005))
+ 
+ 
+ ## Bernoulli GLMM (specified as factor)
+ if (requireNamespace("mlmRev")) {
+     data(Contraception, package="mlmRev")
+     gm3 <- glmer(use ~ urban + age + livch + (1|district),
+                  Contraception, binomial)
+     gm3R <- refit(gm3, Contraception$use)
+     gm3S <- refit(gm3, simulate(gm3)[[1]])
+     stopifnot(all.equal(getinfo(gm3 ),
+                         getinfo(gm3R), tolerance = 1e-5),# 64b_Lx: 7.99e-7
+               all.equal(getinfo(gm3 ),
+                         getinfo(gm3S), tolerance = 0.05) # <- simulated data
+               )
+     cat("gm3: glmer(..):\n"        ); print(getinfo(gm3))
+     cat("gm3R: refit(*, y):\n"     ); print(getinfo(gm3R))
+     cat("gm3S: refit(*, sim.()):\n"); print(getinfo(gm3S))
+ 
+     data(Mmmec, package="mlmRev")
+     if (lme4:::testLevel() > 1) {
+         gm4 <- glmer(deaths ~ uvb + (1|region), data=Mmmec,
+                      family = poisson,
+                      offset = log(expected))
+         ## FIXME: Fails to converge (with larger maxit: "downdate .. not pos.def..")
+         try( gm4R <- refit(gm4, Mmmec $ deaths) )
+         try( gm4S <- refit(gm4, simulate(gm4)[[1]]) )
+         if(FALSE) { ## FIXME (above)
+             cat("gm4R: refit(*,y):\n" ); print( getinfo(gm4R) )
+             cat("gm4S: refit(*,y):\n" ); print( getinfo(gm4S) )
+             stopifnot(all.equal(getinfo(gm4),getinfo(gm4R),tolerance=6e-5))
+         }
+     }
+ }
+ 
+ ## ----------------------------------------------------------------------
+ ## issue: #231, http://ms.mcmaster.ca/~bolker/misc/boot_reset.html
+ ## commits: 1a34cd0, e33d698, 53ce966, 7dbfff1, 73aa1bb, a693ba9, 8dc8cf0
+ ## ----------------------------------------------------------------------
+ 
+ formGrouse <- TICKS ~ YEAR + scale(HEIGHT) + (1 | BROOD) + (1 | INDEX) + (1 | LOCATION)
+ gmGrouse <- glmer(formGrouse, family = "poisson", data = grouseticks)
+ set.seed(105)
+ simTICKS <- simulate(gmGrouse)[[1]]
+ newdata <- transform(grouseticks, TICKS = simTICKS)
+ gmGrouseUpdate <-  update(gmGrouse, data = newdata)
+ gmGrouseRefit  <-  refit(gmGrouse, newresp = simTICKS)
+ 
+ ## compute and print tolerances
+ all.equal(bet.U <- fixef(gmGrouseUpdate),
+           bet.R <- fixef(gmGrouseRefit), tolerance = 0)
+ all.equal(th.U <- getME(gmGrouseUpdate, "theta"),
+           th.R <- getME(gmGrouseRefit,  "theta"), tolerance = 0)
+ all.equal(dev.U <- deviance(gmGrouseUpdate),
+           dev.R <- deviance(gmGrouseRefit),	tolerance = 0)
+ stopifnot(
+     all.equal(bet.U, bet.R, tolerance = 6e-5), # saw 1.0e-5
+     all.equal( th.U,  th.R, tolerance = 4e-5), # saw 1.2e-5
+     all.equal(dev.U, dev.R, tolerance = 2e-5)) # saw 4.6e-6
+ } ## testLevel>1
> 
> proc.time()
   user  system elapsed 
   1.09    0.17    1.25