tests/testthat/test-confint.R

In spaMM: Mixed-Effect Models, with or without Spatial Random Effects

cat(crayon::yellow("\ntest confint() for HLfit, corrHLfit, and fitme (with sparse_precision=",
                   capture.output(spaMM.getOption("sparse_precision")),"):\n"))
data("wafers")
wfit <- HLfit(y ~X1+(1|batch),family=Gamma(log),data=wafers,HLmethod="ML")
ci <- confint(wfit,"X1",verbose=FALSE)
testthat::expect_equal(ci$interval[[1]],0.01828659,tolerance=1e-4)
testthat::expect_equal(ci$interval[[2]],0.17271333,tolerance=1e-4)

# HGLM... DHGLM!
wfit <- HLfit(y ~ X1+X2+X1*X3+X2*X3+I(X2^2)+(1|batch), family=Gamma(log),HLmethod="ML", 
              rand.family=inverse.Gamma(log),
              resid.model = ~ X3+I(X3^2) , data=wafers)
ci <- confint(wfit,"X1",verbose=FALSE)
# values originally set for "ML","exp" fit but sufficiently closer to the "obs" ones, as many others in this script. 
testthat::expect_equal(ci$interval[[1]],0.0361157,tolerance=1e-4)
testthat::expect_equal(ci$interval[[2]],0.1313484,tolerance=1e-4)


#### Checks of consistency of procedures for profiling out one or two parameters (with fixed phi and lambda only for a faster test).
## The CI's of the more constrained model should be within the other (even if both fits coincide at the ML, the additional constraint may matter at the bounds)

## with corrHLfit
data("blackcap")
fitobject <- corrHLfit(migStatus ~ means + Matern(1|longitude+latitude),data=blackcap,HLmethod="ML",
                       init.corrHLfit=list(nu=0.535929513,rho=0.007485725),ranFix=list(phi=0.05,lambda=2))
ci <- confint(fitobject,"means",verbose=FALSE)
testthat::expect_equal(ci$interval[[1]],0.1606797,tolerance=1e-4) 
testthat::expect_equal(ci$interval[[2]],0.9558826,tolerance=1e-4)
refitobject <- corrHLfit(migStatus ~ means + Matern(1|longitude+latitude),data=blackcap,HLmethod="ML",
                         init.corrHLfit=list(rho=0.007485725),ranFix=list(nu=0.535929513,phi=0.05,lambda=2))
ci <- confint(refitobject,"means",verbose=FALSE)
testthat::expect_equal(ci$interval[[1]],0.1586536,tolerance=1e-4)
testthat::expect_equal(ci$interval[[2]],0.9570798,tolerance=1e-4)

## with fitme (derived from gentle intro)
rSample <- function(nb,rho,sigma2_u,resid,intercept,slope,pairs=TRUE) {
  ## sample pairs of adjacent locations
  if (pairs) {
    x <- rnorm(nb/2); x <- c(x,x+0.001)
    y <- rnorm(nb/2); y <- c(y,y+0.001)
  } else {x <- rnorm(nb);y <- rnorm(nb)}
  dist <- dist(cbind(x,y)) ## distance matrix between locations
  m <- exp(-rho*as.matrix(dist)) ## correlation matrix
  b <- MASS::mvrnorm(1,rep(0,nb),m*sigma2_u) ## correlated random effects
  pred <- sample(nb) ##  some predictor variable
  obs <- intercept+slope*pred + b +rnorm(nb,0,sqrt(resid)) ## response
  data.frame(obs=obs,x,y,pred=pred)
}

rngcheck <- ("sample.kind" %in% names (formals(RNGkind)))
if (rngcheck) suppressWarnings(RNGkind("Mersenne-Twister", "Inversion", "Rounding"  )) 
set.seed(123)
d1 <- rSample(nb=40,rho=3,sigma2_u=0.5,resid=0.5,intercept=-1,slope=0.1)
if (rngcheck) RNGkind("Mersenne-Twister", "Inversion", "Rejection"  )
# Rnews: The included LAPACK sources have been updated to 3.10.1. ... 
# ... Using 3.10.x may give some different signs from earlier versions in SVDs or eigendecompositions... 
# => MASS::mvrnorm affected.
d1$obs <- c(1.56404438213562,0.838761745570047,1.18539152814119,2.77609094294087,-0.362117933355623,1.05566922884628,1.04891131933423,-1.36129005102859,
            3.14126410723178,-0.267510545755898,1.68242492955138,-0.294886401760225,1.70572380683716,2.25124059611453,-0.850815921177352,0.439615106586318,
            2.0392252364533,2.42453461515297,-0.53309120933375,-0.554803498441966,1.71344222448971,1.43407925526473,1.27070104772916,-0.557869677646131,
            1.29371981088591,1.98292535259867,-0.440298987269969,0.493786122749392,0.0552051443241781,0.943700535986965,-0.0148351370346089,0.260693893323206,
            2.45333962421848,0.276596082599056,-1.53290736864249,1.88815293467526,-1.90295848044005,2.84485580265504,-0.514747120202024,-2.44239024049091)

HLMf <- fitme(obs~pred+Matern(1|x+y),init=list(rho=59.11287),fixed=list(nu=48.96201,phi=0.447761,lambda=0.3697),data=d1,method="ML")
ci <- confint(HLMf,"pred",verbose=FALSE) 
testthat::expect_equal(ci$interval[[1]],0.06483438,tolerance=1e-4)  
testthat::expect_equal(ci$interval[[2]],0.10567544,tolerance=1e-4)  
HLM <- fitme(obs~pred+Matern(1|x+y),init=list(nu=48.96201,rho=59.11287),fixed=list(phi=0.447761,lambda=0.3697),data=d1,method="ML")
ci <- confint(HLM,"pred",verbose=FALSE) ## practically identical in the two fits (+ nu drifts to higher values whatever the initial one)
testthat::expect_equal(ci$interval[[1]],0.06483437,tolerance=1e-4)  
testthat::expect_equal(ci$interval[[2]],0.10567543,tolerance=1e-4)

# compar to lme4, and tests with ranCoefs
if(requireNamespace("lme4", quietly = TRUE)) {
  data("sleepstudy",package = "lme4")
  (mlfit <- fitme(Reaction ~ Days + (1|Subject), data = sleepstudy))
  (fitci <- confint(mlfit, parm = "Days",verbose=FALSE)$interval)
  rl_mer <- lme4::lmer(Reaction ~ Days + (1|Subject), data = sleepstudy)
  ci_mer <- suppressMessages(suppressWarnings(confint(rl_mer)))[4, ] ## suppress bobyqa conv warning in lme4:::optwrap, and message "Computing profile confidence intervals ..."
  testthat::expect_true(diff(range(ci_mer-fitci))<1e-5) 
  
  if (spaMM.getOption("example_maxtime")>15) { # indeed a bit slow
    (mlfit <- fitme(Reaction ~ Days + (Days|Subject), data = sleepstudy, method="ML"))
    (fitci <- (zut <- confint(mlfit, parm = "Days",verbose=FALSE))$interval) # Now OK, (failed before v4.1.41)
    testthat::expect_true(diff(range(c(7.358652, 13.575920)-fitci))<1e-5) 
    (mlfit <- fitme(Reaction ~ Days + (Days|Subject), data = sleepstudy, method="ML", fixed=list(ranCoefs=list("1"=c(600,-0.1,140)))))
    (fitci <- (zut <- confint(mlfit, parm = "Days",verbose=FALSE))$interval) # seems correct
    testthat::expect_true(diff(range(c(4.849028, 16.085544)-fitci))<1e-5) 
    if (TRUE) {
      ## not sure what to think about these checks for partially fixed ranCoefs but there is no obvious error: 
      # The checked values are those from v4.1.41 and 4.1.50. In between, a stopping bug was introduced.
      (mlfit <- fitme(Reaction ~ Days + (Days|Subject), data = sleepstudy, method="ML", fixed=list(ranCoefs=list("1"=c(NA,-0.1,140)))))
      (fitci <- (zut <- confint(mlfit, parm = "Days",verbose=FALSE))$interval) 
      testthat::expect_true(diff(range(c(4.848926, 16.085646)-fitci))<1e-5) 
      (mlfit <- fitme(Reaction ~ Days + (Days|Subject), data = sleepstudy, method="ML", fixed=list(ranCoefs=list("1"=c(600,NA,140)))))
      (fitci <- (zut <- confint(mlfit, parm = "Days",verbose=FALSE))$interval) 
      testthat::expect_true(diff(range(c(4.848961, 16.085611)-fitci))<1e-5) 
      (mlfit <- fitme(Reaction ~ Days + (Days|Subject), data = sleepstudy, method="ML", fixed=list(ranCoefs=list("1"=c(200,NA,140)))))
      (fitci <- (zut <- confint(mlfit, parm = "Days",verbose=FALSE))$interval) 
      testthat::expect_true(diff(range(c(4.840307, 16.094265)-fitci))<1e-5) 
    }
    (mlfit <- HLfit(Reaction ~ Days + (Days|Subject), data = sleepstudy, method="ML", family=Gamma(log)))
    cat(crayon::yellow("\n[try()-]error EXPECTED here for confint(<*HLfit() return value with ranCoefs*>):"))
    try(zut <- confint(mlfit, parm = "Days",verbose=FALSE)) # catch inappropriate request.
    (mlfit <- fitme(Reaction ~ Days + (Days|Subject), data = sleepstudy, method="ML", family=Gamma(log)))
    (fitci <- (zut <- confint(mlfit, parm = "Days",verbose=FALSE))$interval) 
    testthat::expect_true(diff(range(c(0.02427526, 0.04343179)-fitci))<1e-5) 
  }
  
}

# GLM with a 'method'...
library(spaMM)
set.seed(123)
iris$bidon <- rbinom(nrow(iris), 1, 0.5)
fML <- fitme(bidon ~ Species, data = iris, family = binomial(link = "logit"))
fPQL <- fitme(bidon ~ Species, data = iris, family = binomial(link = "logit"), method = "PQL/L")
testthat::expect_true(identical(confint(fML, "Speciesversicolor", verbose = FALSE)$interval, 
                                confint(fPQL, "Speciesversicolor", verbose = FALSE)$interval))