R Under development (unstable) (2023-11-03 r85470 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.

> ## Tests of a variety of GLMM families and links
> ## coding: family {g=Gamma, P=Poisson, G=Gaussian, B=binomial}
> ##         link   {l=log, i=inverse, c=cloglog, i=identity}
> ##         model  {1 = intercept-only, 2 = with continuous predictor}
> 
> library("lme4")
Loading required package: Matrix
> 
> source(system.file("testdata/lme-tst-funs.R", package="lme4", mustWork=TRUE))
> ##-> gSim(), a general simulation function ...
> 
> str(gSim)
function (nblk = 26, nperblk = 100, sigma_B = 1, beta = c(4, 3), x = runif(n), 
    shape = 2, nbinom = 10, sd = 1, dInitial = NULL, family = Gamma())  
> ## function (nblk = 26, nperblk = 100, sigma = 1, beta = c(4, 3),
> ##           x = runif(n), shape = 2, nbinom = 10, family = Gamma())
> 
> if (.Platform$OS.type != "windows") withAutoprint({
+ set.seed(101)
+ ## Gamma, inverse link (= default) :
+ d <- gSim()
+ ## Gamma, log link   eta = log(mu) :
+ dgl <- gSim(dInitial = d, family = Gamma(link = log))
+ ## Poisson, log link
+ dP <- gSim(dInitial = d, family = poisson())
+ ##  Gaussian, log link --- need to use a non-identity link, otherwise glmer calls lmer
+ dG <- gSim(dInitial = d, family = gaussian(link = log), sd = 2)
+ ## Gaussian with inverse link :     (sd small enough to avoid negative values) :
+ dGi <- gSim(dInitial = d, family = gaussian(link = inverse), sd = 0.01)
+ ## binomial with cloglog link
+ dBc <- d
+ dBc$eta <- d$eta - 5 # <==> beta intercept 5 less: otherwise  y will be constant
+ dBc <- gSim(dInitial = dBc, ## beta = c(-1, 3),
+             nbinom = 1, family = binomial(link="cloglog"))
+ 
+ ## binomial with identity link
+ dBi <- d
+ dBc$eta <- d$eta / 10 # <==> beta slope / 10 : scale so range goes from 0.2-0.8
+ dBi <- gSim(dInitial = dBc, ## beta = c(4, 3/10),
+             nbinom = 1, family = binomial(link="identity"))
+ 
+ 
+ ############
+ ## Gamma/inverse
+ 
+ ## GLMs
+ gm0 <- glm(y ~ 1,       data=d, family=Gamma)
+ gm1 <- glm(y ~ block-1, data=d, family=Gamma)
+ stopifnot(all.equal(sd(coef(gm1)),1.00753942148611))
+ 
+ gm2  <- glmer(y ~ 1 + (1|block), d, Gamma, nAGQ=0)
+ gm3  <- glmer(y ~ x + (1|block), d, Gamma, nAGQ=0)
+ gm2B <- glmer(y ~ 1 + (1|block), d, Gamma)
+ gm3B <- glmer(y ~ x + (1|block), d, Gamma)
+ 
+ ##   y ~ x + (1|block),  Gamma   is TRUE model
+ summary(gm3)
+ summary(gm3B)# should be better
+ ## Both have "correct" beta ~= (4, 3)  -- but *too* small  (sigma_B, sigma) !!
+ stopifnot(exprs = {
+     all.equal(fixef(gm3 ), c(`(Intercept)` = 4.07253, x = 3.080585), tol = 1e-5) # 1.21e-7
+     all.equal(fixef(gm3B), c(`(Intercept)` = 4.159398, x = 3.058521),tol = 1e-5) # 1.13e-7
+ })
+ VarCorr(gm3)  # both variances / std.dev. should be ~ 1  but are too small
+ 
+ 
+ 
+ ##
+ ## library(hglm)
+ ## h1 <- hglm2(y~x+(1|block), data=d, family=Gamma())
+ ## lme4.0 fails on all of these ...
+ 
+ ## Gamma/log
+ ggl1 <- glmer(y ~ 1 + (1|block), data=dgl, family=Gamma(link="log"))
+ ggl2 <- glmer(y ~ x + (1|block), data=dgl, family=Gamma(link="log"))# true model
+ (h.1.2 <- anova(ggl1, ggl2))
+ stopifnot(
+     all.equal(unlist(h.1.2[2,]),
+               c(npar = 4, AIC = 34216.014, BIC = 34239.467, logLik = -17104.007,
+                 deviance = 34208.014, Chisq = 2458.5792, Df = 1, `Pr(>Chisq)` = 0))
+ )
+ ## "true" model :
+ summary(ggl2)
+ VarCorr(ggl2)
+ 
+ ##
+ ## library(lme4.0)
+ ## ggl1 <- glmer(y ~ 1 + (1|block), data=dgl, family=Gamma(link="log"), verbose= 2)
+ ## fails
+ 
+ ## Poisson/log
+ gP1 <- glmer(y ~ 1 + (1|block), data=dP, family=poisson)
+ gP2 <- glmer(y ~ x + (1|block), data=dP, family=poisson)
+ 
+ ## Gaussian/log
+ gG1 <- glmer(y ~ 1 + (1|block), data=dG, family=gaussian(link="log"))
+ gG2 <- glmer(y ~ x + (1|block), data=dG, family=gaussian(link="log"))
+ 
+ ## works with lme4.0 but AIC/BIC/logLik are crazy, and scale
+ ## parameter is not reported
+ ## glmmML etc. doesn't allow models with scale parameters
+ ## gG1B <- glmmadmb(y ~ 1 + (1|block), data=dG,
+ ##                  family="gaussian",link="log",verbose=TRUE)
+ ## what is the best guess at the estimate of the scale parameter?
+ ## is it the same as sigma?
+ ## gG1B$alpha
+ 
+ ## if(Sys.info()["user"] != "maechler") { # <- seg.faults (MM)
+ 
+ ## Gaussian/inverse
+ gGi1 <- glmer(y ~ 1 + (1|block), data=dGi, family=gaussian(link="inverse"))
+ gGi2 <- glmer(y ~ x + (1|block), data=dGi, family=gaussian(link="inverse"))
+ 
+ 
+ ## Binomial/cloglog
+ gBc1 <- glmer(y ~ 1 + (1|block), data=dBc, family=binomial(link="cloglog"))
+ gBc2 <- glmer(y ~ x + (1|block), data=dBc, family=binomial(link="cloglog"))
+ ## library("glmmADMB")
+ ## glmmadmbfit <- glmmadmb(y ~ x + (1|block), data=dBc,
+ ## family="binomial",link="cloglog")
+ glmmadmbfit <- list(fixef = c("(Intercept)" = -0.717146132730349, x =2.83642900561633),
+                     VarCorr = structure(list(
+                         block = structure(0.79992, .Dim = c(1L, 1L),
+                                           .Dimnames = list("(Intercept)", "(Intercept)"))),
+                         class = "VarCorr"))
+ stopifnot(all.equal(fixef(gBc2), glmmadmbfit$fixef, tolerance=5e-3))
+ ## pretty loose tolerance ...
+ stopifnot(all.equal(unname(unlist(VarCorr(gBc2))),
+                     c(glmmadmbfit$VarCorr$block), tolerance=2e-2))
+ 
+ gBi1 <- glmer(y ~ 1 + (1|block), data=dBi, family=binomial(link="identity"))
+ gBi2 <- glmer(y ~ x + (1|block), data=dBi, family=binomial(link="identity"))
+ 
+ ## FIXME: should test more of the *results* of these efforts, not
+ ##  just that they run without crashing ...
+ }) ## skip on windows (for speed)
> 
> proc.time()
   user  system elapsed 
   1.10    0.14    1.25