R Under development (unstable) (2023-11-02 r85465 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2023 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.

> if (.Platform$OS.type != "windows") {
+ 
+ library(lme4)
+ library(testthat)
+ 
+ source(system.file("testdata/lme-tst-funs.R", package="lme4", mustWork=TRUE))
+ ##-> gSim(), a general simulation function ...
+ 
+ ## hand-coded Pearson residuals {for  sumFun() }
+ mypresid <- function(x) {
+     mu <- fitted(x)
+     (getME(x,"y") - mu) * sqrt(weights(x)) / sqrt(x@resp$family$variance(mu))
+ }
+ 
+ ## should be equal (up to numerical error) to weights(.,type="working")
+ workingWeights <- function(mod) mod@resp$weights*(mod@resp$muEta()^2)/mod@resp$variance()
+ 
+ ##' Sum of weighted residuals, 4 ways; the last three are identical
+ sumFun <- function(m) {
+     wrss1 <- m@devcomp$cmp["wrss"]
+     wrss2 <- sum(residuals(m,type="pearson")^2)
+     wrss3 <- sum(m@resp$wtres^2)
+     ## compare to hand-fitted Pearson resids ...
+     wrss4 <- sum(mypresid(m)^2)
+     c(wrss1,wrss2,wrss3,wrss4)
+ }
+ ## The relative "error"/differences of the weights w[] entries
+ rel.diff <- function(w) abs(1 - w[-1]/w[1])
+ 
+ set.seed(101)
+ 
+ ## GAMMA
+ g0 <-  glmer(y~x+(1|block),data=gSim(),family=Gamma)
+ expect_true(all(rel.diff(sumFun(g0)) < 1e-13))
+ expect_equal(weights(g0, type = "working"), workingWeights(g0),
+              tolerance = 1e-4)  ## FIXME:  why is such a high tolerance required?
+ 
+ ## BERNOULLI
+ g1 <-  glmer(y~x+(1|block),data=gSim(family=binomial(),nbinom=1),
+              family=binomial)
+ expect_true(all(rel.diff(sumFun(g1)) < 1e-13))
+ expect_equal(weights(g1, type = "working"), workingWeights(g1),
+              tolerance = 1e-5)  ## FIXME:  why is such a high tolerance required?
+ 
+ 
+ ## POISSON
+ (n <- nrow(d.P <- gSim(family=poisson())))
+ g2 <-  glmer(y ~ x + (1|block), data = d.P, family=poisson)
+ g2W <- glmer(y ~ x + (1|block), data = d.P, family=poisson, weights = rep(2,n))
+ expect_true(all(rel.diff(sumFun(g2 )) < 1e-13))
+ expect_true(all(rel.diff(sumFun(g2W)) < 1e-13))
+ ## correct
+ expect_equal(weights(g2, type = "working"), workingWeights(g2),
+              tolerance = 1e-5)  ## FIXME:  why is such a high tolerance required?
+ expect_equal(weights(g2W, type = "working"), workingWeights(g2W),
+              tolerance = 1e-5)  ## FIXME:  why is such a high tolerance required?
+ 
+ 
+ ## non-Bernoulli BINOMIAL
+ g3 <- glmer(y ~ x + (1|block), data= gSim(family=binomial(), nbinom=10),
+             family=binomial)
+ expect_true(all(rel.diff(sumFun(g3)) < 1e-13))
+ expect_equal(weights(g3, type = "working"), workingWeights(g3),
+              tolerance = 1e-4)  ## FIXME:  why is such a high tolerance required?
+ 
+ 
+ 
+ d.b.2 <- gSim(nperblk = 2, family=binomial())
+ g.b.2 <- glmer(y ~ x + (1|block), data=d.b.2, family=binomial)
+ 
+ expect_true(all(rel.diff(sumFun(g.b.2 )) < 1e-13))
+ 
+ 
+ ## Many blocks of only 2 observations each - (but nicely balanced)
+ ## Want this "as" https://github.com/lme4/lme4/issues/47
+ ## (but it "FAILS" survival already):
+ ##
+ ## n2 = n/2 :
+ n2 <- 2048
+ if(FALSE)
+ n2 <-  100 # for building/testing
+ set.seed(47)
+ dB2 <- gSim(n2, nperblk = 2, x= rep(0:1, each= n2), family=binomial())
+ ##                       --  --     ---  --------
+ gB2 <- glmer(y ~ x + (1|block), data=dB2, family=binomial)
+ expect_true(all(rel.diff(sumFun(gB2)) < 1e-13))
+ 
+ ## NB: Finite sample bias of  \hat\sigma_1 and  \hat\beta_1 ("Intercept")
+ ##     tend to zero only slowly for  n2 -> Inf,  e.g., for
+ ## n2 = 2048,  b1 ~= 4.3 (instead of 4);  s1 ~= 1.3 (instead of 1)
+ 
+ ## FAILS -----
+ ## library(survival)
+ ## (gSurv.B2 <- clogit(y ~ x + strata(block), data=dB2))
+ ## ## --> Error in Surv(rep(1, 200L), y) : Time and status are different lengths
+ ## summary(gSurv.B2)
+ ## (SE.surf <- sqrt(diag(vcov(gSurv.B2))))
+ 
+ 
+ 
+ g3 <-  glmer(y ~ x + (1|block),data=gSim(family=binomial(),nbinom=10),
+              family=binomial)
+ expect_equal(var(sumFun(g3)),0)
+ 
+ ## check dispersion parameter
+ ## (lowered tolerance to pass checks on my machine -- SCW)
+ expect_equal(sigma(g0)^2, 0.4888248, tolerance=1e-4)
+ 
+ } ## skip on windows (for speed)
> 
> proc.time()
   user  system elapsed 
   0.26    0.09    0.29