repeated: Repeated Measurement Models for Discrete Times

### test-auto-pattern-regression.R --- 
##----------------------------------------------------------------------
## Author: Brice Ozenne
## Created: jul 31 2023 (09:54) 
## Version: 
## Last-Updated: May 18 2024 (13:20) 
##           By: Brice Ozenne
##     Update #: 22
##----------------------------------------------------------------------
## 
### Commentary: 
## 
### Change Log:
##----------------------------------------------------------------------
## 
### Code:

if(FALSE){
    library(testthat)
    library(numDeriv)
    library(lava)
    library(multcomp)
    library(nlme)

    library(LMMstar)
}

context("Check lmm on mixed models parametrized with covariance patterns")
LMMstar.options(optimizer = "FS", method.numDeriv = "simple", precompute.moments = TRUE, # "Richardson"
                columns.confint = c("estimate","se","df","lower","upper","p.value"))

## * Simulate data
m <- lvm(c(Y1,Y2,Y3,Y4) ~ 0.05*age + gender)
categorical(m, labels = c("male","female")) <- ~gender
transform(m, id~gender) <- function(x){1:NROW(x)}
distribution(m, ~age) <- gaussian.lvm(mean = 50, sd = 10)

set.seed(10)
dW <- lava::sim(m, 1e2)

## move to the long format
name.varying <- paste0("Y",1:4)
dL <- reshape(dW, direction  = "long",
              idvar = c("id","age","gender"),
              varying = name.varying,
              v.names = "Y",
              timevar = "visit")
rownames(dL) <- NULL
dL$visit <- factor(dL$visit,
                   levels = 1:length(name.varying),
                   labels = name.varying)
 

## * Compound symmetry structure

test_that("Compound symmetry structure (REML)",{

    ## ** fit
    eCS.lmm <- lmm(Y ~ visit + age + gender, repetition = ~visit|id, structure = "CS", data = dL, trace = 0, method = "REML")
    eCS.gls <- gls(Y ~ visit + age + gender, correlation = corCompSymm(form=~1|id), data = dL, method = "REML")
    
    ## ** iteration
    expect_equal(eCS.lmm$opt$n.iter,4)

    ## ** coef
    expect_equal(coef(eCS.lmm, effects = "mean"), coef(eCS.gls), tol = 1e-6)
    coef(eCS.lmm, transform.rho = "cov")
    coef(eCS.lmm, transform.sigma = "square")

    ## ** logLikelihood
    expect_equal(logLik(eCS.lmm), as.double(logLik(eCS.gls)), tol = 1e-6)

    ## ** score
    GS <- jacobian(func = function(p){logLik(eCS.lmm, p = p)}, x = coef(eCS.lmm, effects = "all", transform.sigma = "none", transform.k = "none", transform.rho = "none"))
    GS <- rbind(c(4e-08, 3e-08, 0, -2e-08, -2e-08, 0, 5.27e-06, 1.952e-05))
    test <- score(eCS.lmm, effects = "all")
    expect_true(all(abs(test) < LMMstar.options()$tol.score))
    expect_equal(as.double(GS), as.double(test), tol = 1e-5)

    ## no transformation
    newp <- coef(eCS.lmm, effects = "all", transform.sigma = "none", transform.k = "none", transform.rho = "none")+0.1
    ## GS <- jacobian(func = function(p){logLik(eCS.lmm, p = p)}, x = newp)
    GS <- rbind(c(-1434.54797024, -360.96098759, -360.9609876, -360.96098776, -75945.93430486, -597.3281072, 7235.68958315, 9425.0785116))
    test <- score(eCS.lmm, effects = "all", p = newp, transform.sigma = "none", transform.k = "none", transform.rho = "none")
    expect_equal(as.double(test), as.double(GS), tol = 1e-6)

    ## log transformation
    newp.log <- newp; newp.log["sigma"] <- log(newp["sigma"])
    ## GS <- jacobian(func = function(p){p["sigma"] <- exp(p["sigma"]); logLik(eCS.lmm, p = p)}, x = newp.log)
    GS <- rbind(c(-1434.54797024, -360.96098759, -360.9609876, -360.96098776, -75945.93430486, -597.3281072, 7844.35575598, 9425.0785116))
    test <- score(eCS.lmm, effects = "all", p = newp , transform.sigma = "log", transform.k = "none", transform.rho = "none")
    expect_equal(as.double(test), as.double(GS), tol = 1e-6)

    ## lava transformation
    newp.2 <- newp; newp.2["sigma"] <- newp["sigma"]^2
    ## GS <- jacobian(func = function(p){p["sigma"] <- sqrt(p["sigma"]); logLik(eCS.lmm, p = p)}, x = newp.2)
    GS <- rbind(c(-1434.54797024, -360.96098759, -360.9609876, -360.96098776, -75945.93430486, -597.3281072, 3337.12578619, 9425.0785116))
    test <- score(eCS.lmm, effects = "all", p = newp, transform.sigma = "square", transform.k = "none", transform.rho = "none")
    expect_equal(as.double(test), as.double(GS), tol = 1e-6)

    ## cov transformation
    newp.cov <- newp; newp.cov[c("sigma","rho")] <- c(newp["sigma"]^2,newp["rho"]*newp["sigma"]^2)
    ## GS <- jacobian(func = function(p){p[c("sigma","rho")] <- c(sqrt(p["sigma"]),p["rho"]/p["sigma"]); logLik(eCS.lmm, p = p, transform.sigma = "none", transform.k = "none", transform.rho = "cov")}, x = newp.cov)
    GS <- rbind(c(-1434.54797024, -360.96098759, -360.9609876, -360.96098776, -75945.93430486, -597.3281072, 2657.67539179, 8019.18564457))
    test <- score(eCS.lmm, effects = "all", p = newp, transform.rho = "cov")
    expect_equal(as.double(test), as.double(GS), tol = 1e-6)

    ## ** information
    ## no transformation
    test <- information(eCS.lmm, effects = "all", transform.sigma = "none", transform.k = "none", transform.rho = "none",
                        p = coef(eCS.lmm, effects = "all", transform.sigma = "none", transform.k = "none", transform.rho = "none"))
    ## GS <- -hessian(func = function(p){logLik(eCS.lmm, p = p)},
    ##                x = coef(eCS.lmm, effects = "all", transform.sigma = "none", transform.k = "none", transform.rho = "none"))
    GS <- cbind(c(432.84420535, 108.21105139, 108.21105134, 108.21105134, 21943.19088006, 181.79456625, -1e-08, 5e-08), 
                c(108.21105139, 103.32786316, 1.62772943, 1.6277295, 5485.79771987, 45.44864157, 0, 3.7e-07), 
                c(108.21105134, 1.62772943, 103.32786307, 1.62772943, 5485.79772009, 45.44864156, 0, 1.6e-07), 
                c(108.21105134, 1.6277295, 1.62772943, 103.32786307, 5485.79772003, 45.44864156, 0, 5e-08), 
                c(21943.19088006, 5485.79771987, 5485.79772009, 5485.79772003, 1163984.34440481, 9028.04063345, -6.7e-07, -2.5e-07), 
                c(181.79456625, 45.44864157, 45.44864156, 45.44864156, 9028.04063345, 181.79456625, 0, 0), 
                c(-1e-08, 0, 0, 0, -6.7e-07, 0, 813.63592184, 12.64082829), 
                c(5e-08, 3.7e-07, 1.6e-07, 5e-08, -2.5e-07, 0, 12.64082829, 623.49010272)
                )
    expect_equal(as.double(test), as.double(GS), tol = 1e-6)

    ## log transformation
    newp.log <- newp; newp.log["sigma"] <- log(newp["sigma"])
    ## GS <- -hessian(func = function(p){p["sigma"] <- exp(p["sigma"]); logLik(eCS.lmm, p = p)}, x = newp.log)
    GS <- cbind(c(271.35880659, 67.83969975, 67.83970126, 67.83970148, 13756.6311032, 113.97069813, -2869.09594257, -3431.42859628), 
                c(67.83969975, 86.67977273, -6.28002394, -6.28002374, 3439.15777569, 28.49267455, -721.9219748, -855.3180185), 
                c(67.83970126, -6.28002394, 86.67977266, -6.28002379, 3439.15777583, 28.49267453, -721.92197541, -855.31801867), 
                c(67.83970148, -6.28002374, -6.28002379, 86.67977256, 3439.15777592, 28.49267454, -721.92197526, -855.31801859), 
                c(13756.6311032, 3439.15777569, 3439.15777583, 3439.15777592, 729725.3769578, 5659.861652, -151891.8686083, -181662.13752763), 
                c(113.97069813, 28.49267455, 28.49267453, 28.49267454, 5659.861652, 113.97069816, -1194.6562144, -1428.80460625), 
                c(-2869.09594257, -721.9219748, -721.92197541, -721.92197526, -151891.8686083, -1194.6562144, 16476.71151466, 18757.68659853), 
                c(-3431.42859628, -855.3180185, -855.31801867, -855.31801859, -181662.13752763, -1428.80460625, 18757.68659853, 45449.68536078)
                )
    test <- information(eCS.lmm, p = newp, effects = "all", transform.sigma = "log", transform.k = "none", transform.rho = "none")
    expect_equal(as.double(test), as.double(GS), tol = 1e-6)

    ## lava transformation
    newp.2 <- newp; newp.2["sigma"] <- newp["sigma"]^2
    ## GS <- -hessian(func = function(p){p["sigma"] <- sqrt(p["sigma"]); logLik(eCS.lmm, p = p)}, x = newp.2)
    GS <- cbind(c(271.35880659, 67.83969975, 67.83970126, 67.83970148, 13756.6311032, 113.97069813, -1220.56346615, -3431.42859628), 
                c(67.83969975, 86.67977273, -6.28002394, -6.28002374, 3439.15777569, 28.49267455, -307.11820248, -855.3180185), 
                c(67.83970126, -6.28002394, 86.67977266, -6.28002379, 3439.15777583, 28.49267453, -307.11820248, -855.31801867), 
                c(67.83970148, -6.28002374, -6.28002379, 86.67977256, 3439.15777592, 28.49267454, -307.11820247, -855.31801859), 
                c(13756.6311032, 3439.15777569, 3439.15777583, 3439.15777592, 729725.3769578, 5659.861652, -64617.45070718, -181662.13752763), 
                c(113.97069813, 28.49267455, 28.49267453, 28.49267454, 5659.861652, 113.97069816, -508.22759476, -1428.80460625), 
                c(-1220.56346615, -307.11820248, -307.11820248, -307.11820247, -64617.45070718, -508.22759476, 5821.29834877, 7979.8471127), 
                c(-3431.42859628, -855.3180185, -855.31801867, -855.31801859, -181662.13752763, -1428.80460625, 7979.8471127, 45449.68536078)
                )
    test <- information(eCS.lmm, p = newp, effects = "all", transform.sigma = "square", transform.k = "none", transform.rho = "none")
    expect_equal(as.double(test), as.double(GS), tol = 1e-6)

    ## cov transformation
    newp.cov <- newp; newp.cov[c("sigma","rho")] <- c(newp["sigma"]^2,newp["rho"]*newp["sigma"]^2)
    ## GS <- -hessian(func = function(p){p[c("sigma","rho")] <- c(sqrt(p["sigma"]),p["rho"]/p["sigma"]); logLik(eCS.lmm, p = p)}, x = newp.cov)
    GS <- cbind(c(271.35880659, 67.83969975, 67.83970126, 67.83970148, 13756.6311032, 113.97069813, -973.19305846, -2919.57917442), 
                c(67.83969975, 86.67977273, -6.28002394, -6.28002374, 3439.15777569, 28.49267455, -245.45864557, -727.73441273), 
                c(67.83970126, -6.28002394, 86.67977266, -6.28002379, 3439.15777583, 28.49267453, -245.45864555, -727.73441286), 
                c(67.83970148, -6.28002374, -6.28002379, 86.67977256, 3439.15777592, 28.49267454, -245.45864554, -727.7344128), 
                c(13756.6311032, 3439.15777569, 3439.15777583, 3439.15777592, 729725.3769578, 5659.861652, -51521.49497161, -154564.4849145), 
                c(113.97069813, 28.49267455, 28.49267453, 28.49267454, 5659.861652, 113.97069816, -405.22560363, -1215.67681085), 
                c(-973.19305846, -245.45864557, -245.45864555, -245.45864554, -51521.49497161, -405.22560363, 3750.76786163, 10824.81657287), 
                c(-2919.57917442, -727.73441273, -727.73441286, -727.7344128, -154564.4849145, -1215.67681085, 10824.81657287, 32901.93673129)
                )
    test <- information(eCS.lmm, p = newp, effects = "all", transform.rho = "cov")
    expect_equal(as.double(test), as.double(GS), tol = 1e-6)

    ## no transformation 
    newp <- coef(eCS.lmm, effects = "all", transform.sigma = "none", transform.k = "none", transform.rho = "none")+0.1
    ## GS <- -hessian(func = function(p){logLik(eCS.lmm, p = p)}, x = newp)
    GS <- cbind(c(271.35880659, 67.83969975, 67.83970126, 67.83970148, 13756.6311032, 113.97069813, -2646.47451812, -3431.42859628), 
                c(67.83969975, 86.67977273, -6.28002394, -6.28002374, 3439.15777569, 28.49267455, -665.9059684, -855.3180185), 
                c(67.83970126, -6.28002394, 86.67977266, -6.28002379, 3439.15777583, 28.49267453, -665.90596843, -855.31801867), 
                c(67.83970148, -6.28002374, -6.28002379, 86.67977256, 3439.15777592, 28.49267454, -665.90596843, -855.31801859), 
                c(13756.6311032, 3439.15777569, 3439.15777583, 3439.15777592, 729725.3769578, 5659.861652, -140106.14071148, -181662.13752763), 
                c(113.97069813, 28.49267455, 28.49267453, 28.49267454, 5659.861652, 113.97069816, -1101.95939524, -1428.80460625), 
                c(-2646.47451812, -665.9059684, -665.90596843, -665.90596843, -140106.14071148, -1101.95939524, 20693.21261074, 17302.22362826), 
                c(-3431.42859628, -855.3180185, -855.31801867, -855.31801859, -181662.13752763, -1428.80460625, 17302.22362826, 45449.68536078)
                )
    test <- information(eCS.lmm, p = newp, effects = "all", transform.sigma = "none", transform.k = "none", transform.rho = "none")
    expect_equal(as.double(test), as.double(GS), tol = 1e-6) 

    ## ** degree of freedom
    test <- model.tables(eCS.lmm, effects = "all", transform.sigma = "none")
    expect_equal(test[c("visitY2","visitY3","visitY4"),"df"], rep(297,3), tol = 1e-3)
    expect_equal(test[c("age","genderfemale"),"df"], rep(97,2), tol = 1e-3)

    ## ** anova
    eCS.lmm_anova <- anova(eCS.lmm, effects = "all")
    expect_equal(eCS.lmm_anova$multivariate[eCS.lmm_anova$multivariate$type=="mu","df.denom"], c(297.03106,  97.05084,  97.05084), tol = 1e-1)
    expect_equal(eCS.lmm_anova$multivariate[eCS.lmm_anova$multivariate$type=="rho","df.denom"], c(14.7493), tol = 1e-1)

    ## ** getVarCov
    sigma(eCS.lmm)

    ## ** prediction
    test <- predict(eCS.lmm, newdata = dL, se = TRUE)
    index <- sample.int(NROW(dL))
    test2 <- predict(eCS.lmm, newdata = dL[index,,drop=FALSE], se = TRUE)
    ## GS <- AICcmodavg::predictSE(eCS.gls, newdata = dL)
    GS <- list(fit = c("1" = 2.78132573, "2" = 3.04660727, "3" = 3.47291713, "4" = 3.30670866, "5" = 2.74336483, "6" = 3.19347129) ,
               se.fit = c("1" = 0.11128311, "2" = 0.11700474, "3" = 0.11383802, "4" = 0.11360494, "5" = 0.11036375, "6" = 0.12667589) )
    expect_equivalent(head(test$estimate),GS$fit, tol = 1e-7)
    expect_equivalent(head(test$se),GS$se.fit, tol = 1e-7)
    expect_equivalent(test[index,,drop=FALSE],test2, tol = 1e-7)

})

## * Unstructed covariance matrix
test_that("Unstructured covariance matrix (REML)",{

    ## ** fit
    eUNexp.lmm <- lmm(Y ~ visit + age + gender, repetition = ~visit|id, structure = "UN", data = dL, trace = 0, method = "REML", type.information = "expected")
    eUN.lmm <- lmm(Y ~ visit + age + gender, repetition = ~visit|id, structure = "UN", data = dL, trace = 0, method = "REML", type.information = "observed")
    eUN.gls <- gls(Y ~ visit + age + gender, correlation = corSymm(form=~1|id), weights = varIdent(form=~1|visit), data = dL, method = "REML")

    ## ** iteration
    expect_equal(eUNexp.lmm$opt$n.iter,4)
    expect_equal(eUN.lmm$opt$n.iter,4)

    ## ** coef
    expect_equal(coef(eUN.lmm, effects = "mean"), coef(eUN.gls), tol = 1e-6)
    expect_equal(coef(eUNexp.lmm, effects = "mean"), coef(eUN.gls), tol = 1e-6)

    ## ** logLikelihood
    expect_equal(logLik(eUN.lmm), as.double(logLik(eUN.gls)), tol = 1e-6)
    expect_equal(logLik(eUNexp.lmm), as.double(logLik(eUN.gls)), tol = 1e-6)

    ## ** score
    test <- score(eUN.lmm, effects = "all", transform.sigma = "none", transform.k = "none", transform.rho = "none")
    ## GS <- jacobian(func = function(p){logLik(eUN.lmm, p = p)},
    ##                x = coef(eUN.lmm, effects = "all", transform.sigma = "none", transform.k = "none", transform.rho = "none"))
    GS <- rbind(c(-4e-08, -1.6e-07, 0, 0, 0, 0, 4.95e-06, 1.533e-05, -1.655e-05, 7.73e-06, 8.99e-06, -1.1e-05, 6.85e-06, -1e-08, 1.745e-05, 0))
    ## expect_true(all(abs(test) < 1e-3))
    expect_equal(as.double(GS), as.double(test), tol = 1e-5)

    ## no transformation
    newp <- coef(eUN.lmm, effects = "all", transform.sigma = "none", transform.k = "none", transform.rho = "none")+0.1
    ## GS <- jacobian(func = function(p){logLik(eUN.lmm, p = p)}, x = newp)
    GS <- rbind(c(-1274.98904333, -342.88438277, -334.63254795, -215.90335215, -67511.19682982, -530.78211036, 6631.4375721, 1641.16355148, 1614.71521605, 962.06391598, 1678.30491401, 1640.1332732, 1126.74501451, 1657.87479712, 1141.3752835, 1114.1769838))
    test <- score(eUN.lmm, effects = "all", p = newp, transform.sigma = "none", transform.k = "none", transform.rho = "none")
    expect_equal(as.double(test), as.double(GS), tol = 1e-6)

    ## log transformation
    newp.log <- newp; newp.log["sigma"] <- log(newp["sigma"])
    ## GS <- jacobian(func = function(p){p["sigma"] <- exp(p["sigma"]); logLik(eUN.lmm, p = p)}, x = newp.log)
    GS <- rbind(c(-1274.98904333, -342.88438277, -334.63254795, -215.90335215, -67511.19682982, -530.78211036, 6916.83142601, 1641.16355148, 1614.71521605, 962.06391598, 1678.30491401, 1640.1332732, 1126.74501451, 1657.87479712, 1141.3752835, 1114.1769838))
    test <- score(eUN.lmm, effects = "all", p = newp, transform.sigma = "log", transform.k = "none", transform.rho = "none")
    expect_equal(as.double(test), as.double(GS), tol = 1e-6)

    ## lava transformation
    newp.2 <- newp; newp.2["sigma"] <- newp["sigma"]^2
    ## GS <- jacobian(func = function(p){p["sigma"] <- sqrt(p["sigma"]); logLik(eUN.lmm, p = p)}, x = newp.2)
    GS <- rbind(c(-1274.98904333, -342.88438277, -334.63254795, -215.90335215, -67511.19682982, -530.78211036, 3178.90964642, 1641.16355148, 1614.71521605, 962.06391598, 1678.30491401, 1640.1332732, 1126.74501451, 1657.87479712, 1141.3752835, 1114.1769838))
    test <- score(eUN.lmm, effects = "all", p = newp, transform.sigma = "square", transform.k = "none", transform.rho = "none")
    expect_equal(as.double(test), as.double(GS), tol = 1e-6)

    ## cov transformation
    ## name.sigma <- eUN.lmm$design$param[eUN.lmm$design$param$type=="rho","name"]
    ## name.rho <- eUN.lmm$design$param[eUN.lmm$design$param$type=="rho","name"]
    ## name.k1 <- eUN.lmm$design$param[eUN.lmm$design$param$type=="rho","k.x"]
    ## name.k2 <- eUN.lmm$design$param[eUN.lmm$design$param$type=="rho","k.y"]
    ## newp.cov <- newp
    ## newp.cov[name.rho] <- newp[name.rho]*ifelse(is.na(newp[name.k1]),1,newp[name.k1])*newp[name.k2]*newp["sigma"]^2
    ## newp.cov[name.rho] <- newp[name.rho]*ifelse(is.na(newp[name.k1]),1,newp[name.k1])*newp[name.k2]*newp["sigma"]^2
    ## GS <- jacobian(func = function(p){p[name.rho] <- c(p[name.rho]*); logLik(eUN.lmm, p = p)}, x = newp.cov)
    ## test <- score(eUN.lmm, p = newp, transform.rho = "cov")
    ## expect_equal(as.double(test), as.double(GS), tol = 1e-6)

    ## ** information
    expect_equal(vcov(eUNexp.lmm, effects = "mean"), vcov(eUN.gls), tol = 1e-6)

    ## no transformation
    newp <- coef(eUN.lmm, effects = "all", transform.sigma = "none", transform.k = "none", transform.rho = "none")+0.1
    ## GS <- -hessian(func = function(p){logLik(eUN.lmm, p = p)}, x = newp)
    GS <- cbind(c(241.4276261, 64.50068304, 62.84798106, 40.45548526, 12239.25904978, 101.39960062, -2444.76400233, -605.09410326, -591.34394345, -356.6684903, -614.45166963, -597.89210276, -411.70248495, -600.92192108, -413.7827081, -402.6264939), 
                c(64.50068304, 73.21443141, 3.49099427, -9.33761434, 3269.88504359, 27.09028677, -657.47341151, -645.9497367, -16.4036038, 41.06398553, -332.18747695, -4.58242056, 54.95083336, -358.00042572, -187.64104164, 35.21654172), 
                c(62.84798106, 3.49099427, 73.4157247, -5.92613477, 3186.10076536, 26.39615222, -641.65069748, -16.36801713, -641.00001661, 26.06133703, 17.95060609, -311.03354514, 52.78431089, -366.97750658, 19.02772528, -205.98042191), 
                c(40.45548526, -9.33761434, -5.92613477, 65.9452637, 2050.9052171, 16.99130389, -413.99002328, 43.78071542, 27.84592494, -468.76477714, 87.12601096, 69.89098536, -301.62532094, 71.84055748, -306.57808259, -308.77047073), 
                c(12239.25904978, 3269.88504359, 3186.10076536, 2050.9052171, 649235.8381143, 5035.57247599, -129451.26452782, -31952.35290987, -31341.07385733, -18837.71167947, -32546.69377367, -31728.40145323, -21802.76288715, -31788.23696115, -21844.0631514, -21294.2253707), 
                c(101.39960062, 27.09028677, 26.39615222, 16.99130389, 5035.57247599, 101.39960056, -1017.76325408, -257.40463535, -238.49688193, -150.98532049, -258.3743199, -244.46171114, -173.04102013, -248.67924617, -175.90487759, -166.54859374), 
                c(-2444.76400233, -657.47341151, -641.65069748, -413.99002328, -129451.26452782, -1017.76325408, 19797.77242215, 3313.55987297, 3263.23133819, 2001.37516506, 3210.47594884, 3125.03071949, 2134.38017671, 3187.27455335, 2161.55979863, 2119.10105462), 
                c(-605.09410326, -645.9497367, -16.36801713, 43.78071542, -31952.35290987, -257.40463535, 3313.55987297, 4698.45923152, 79.52612717, -200.35755175, 1618.75890508, 22.29548051, -267.99319578, 1744.82398992, 899.12913933, -171.78219472), 
                c(-591.34394345, -16.4036038, -641.00001661, 27.84592494, -31341.07385733, -238.49688193, 3263.23133819, 79.52612717, 4680.85329314, -127.68986192, -87.29462845, 1519.58952368, -258.26729124, 1792.94480904, -93.0270774, 996.52858573), 
                c(-356.6684903, 41.06398553, 26.06133703, -468.76477714, -18837.71167947, -150.98532049, 2001.37516506, -200.35755175, -127.68986192, 3035.28503039, -397.77172353, -319.5137104, 1373.54307929, -327.56863472, 1385.87089749, 1403.891274), 
                c(-614.45166963, -332.18747695, 17.95060609, 87.12601096, -32546.69377367, -258.3743199, 3210.47594884, 1618.75890508, -87.29462845, -397.77172353, 3397.00621427, 1570.20497911, 669.40246524, 1521.79966982, 640.13813427, -509.51606881), 
                c(-597.89210276, -4.58242056, -311.03354514, 69.89098536, -31728.40145323, -244.46171114, 3125.03071949, 22.29548051, 1519.58952368, -319.5137104, 1570.20497911, 3105.35116462, 663.7000477, 1566.17475066, -355.5590459, 666.96803064), 
                c(-411.70248495, 54.95083336, 52.78431089, -301.62532094, -21802.76288715, -173.04102013, 2134.38017671, -267.99319578, -258.26729124, 1373.54307929, 669.40246524, 663.7000477, 2294.20243311, -527.6926834, 1418.45732112, 1384.39801595), 
                c(-600.92192108, -358.00042572, -366.97750658, 71.84055748, -31788.23696115, -248.67924617, 3187.27455335, 1744.82398992, 1792.94480904, -327.56863472, 1521.79966982, 1566.17475066, -527.6926834, 3602.07056013, 807.0376631, 836.8177412), 
                c(-413.7827081, -187.64104164, 19.02772528, -306.57808259, -21844.0631514, -175.90487759, 2161.55979863, 899.12913933, -93.0270774, 1385.87089749, 640.13813427, -355.5590459, 1418.45732112, 807.0376631, 2296.85216508, 1521.87542022), 
                c(-402.6264939, 35.21654172, -205.98042191, -308.77047073, -21294.2253707, -166.54859374, 2119.10105462, -171.78219472, 996.52858573, 1403.891274, -509.51606881, 666.96803064, 1384.39801595, 836.8177412, 1521.87542022, 2334.12612433)
                )
    test <- information(eUN.lmm, effects = "all", p = newp, transform.sigma = "none", transform.k = "none", transform.rho = "none")
    expect_equal(as.double(test), as.double(GS), tol = 1e-6)

    ## log transformation
    newp.log <- newp; newp.log["sigma"] <- log(newp["sigma"])
    ## GS <- -hessian(func = function(p){p["sigma"] <- exp(p["sigma"]); logLik(eUN.lmm, p = p)}, x = newp.log)
    GS <- cbind(c(241.4276261, 64.50068304, 62.84798106, 40.45548526, 12239.25904978, 101.39960062, -2549.97808671, -605.09410326, -591.34394345, -356.6684903, -614.45166963, -597.89210276, -411.70248495, -600.92192108, -413.7827081, -402.6264939), 
                c(64.50068304, 73.21443141, 3.49099427, -9.33761434, 3269.88504359, 27.09028677, -685.76876562, -645.9497367, -16.4036038, 41.06398553, -332.18747695, -4.58242056, 54.95083336, -358.00042572, -187.64104164, 35.21654172), 
                c(62.84798106, 3.49099427, 73.4157247, -5.92613477, 3186.10076536, 26.39615222, -669.26509654, -16.36801713, -641.00001661, 26.06133703, 17.95060609, -311.03354514, 52.78431089, -366.97750658, 19.02772528, -205.98042191), 
                c(40.45548526, -9.33761434, -5.92613477, 65.9452637, 2050.9052171, 16.99130389, -431.80670431, 43.78071542, 27.84592494, -468.76477714, 87.12601096, 69.89098536, -301.62532094, 71.84055748, -306.57808259, -308.77047073), 
                c(12239.25904978, 3269.88504359, 3186.10076536, 2050.9052171, 649235.8381143, 5035.57247599, -135022.39366066, -31952.35290987, -31341.07385733, -18837.71167947, -32546.69377367, -31728.40145323, -21802.76288715, -31788.23696115, -21844.0631514, -21294.2253707), 
                c(101.39960062, 27.09028677, 26.39615222, 16.99130389, 5035.57247599, 101.39960056, -1061.56422075, -257.40463535, -238.49688193, -150.98532049, -258.3743199, -244.46171114, -173.04102013, -248.67924617, -175.90487759, -166.54859374), 
                c(-2549.97808671, -685.76876562, -669.26509654, -431.80670431, -135022.39366066, -1061.56422075, 14621.66285348, 3456.16388747, 3403.66938875, 2087.50734443, 3348.64358899, 3259.52109889, 2226.23642464, 3324.44368914, 2254.58576331, 2210.29974341), 
                c(-605.09410326, -645.9497367, -16.36801713, 43.78071542, -31952.35290987, -257.40463535, 3456.16388747, 4698.45923152, 79.52612717, -200.35755175, 1618.75890508, 22.29548051, -267.99319578, 1744.82398992, 899.12913933, -171.78219472), 
                c(-591.34394345, -16.4036038, -641.00001661, 27.84592494, -31341.07385733, -238.49688193, 3403.66938875, 79.52612717, 4680.85329314, -127.68986192, -87.29462845, 1519.58952368, -258.26729124, 1792.94480904, -93.0270774, 996.52858573), 
                c(-356.6684903, 41.06398553, 26.06133703, -468.76477714, -18837.71167947, -150.98532049, 2087.50734443, -200.35755175, -127.68986192, 3035.28503039, -397.77172353, -319.5137104, 1373.54307929, -327.56863472, 1385.87089749, 1403.891274), 
                c(-614.45166963, -332.18747695, 17.95060609, 87.12601096, -32546.69377367, -258.3743199, 3348.64358899, 1618.75890508, -87.29462845, -397.77172353, 3397.00621427, 1570.20497911, 669.40246524, 1521.79966982, 640.13813427, -509.51606881), 
                c(-597.89210276, -4.58242056, -311.03354514, 69.89098536, -31728.40145323, -244.46171114, 3259.52109889, 22.29548051, 1519.58952368, -319.5137104, 1570.20497911, 3105.35116462, 663.7000477, 1566.17475066, -355.5590459, 666.96803064), 
                c(-411.70248495, 54.95083336, 52.78431089, -301.62532094, -21802.76288715, -173.04102013, 2226.23642464, -267.99319578, -258.26729124, 1373.54307929, 669.40246524, 663.7000477, 2294.20243311, -527.6926834, 1418.45732112, 1384.39801595), 
                c(-600.92192108, -358.00042572, -366.97750658, 71.84055748, -31788.23696115, -248.67924617, 3324.44368914, 1744.82398992, 1792.94480904, -327.56863472, 1521.79966982, 1566.17475066, -527.6926834, 3602.07056013, 807.0376631, 836.8177412), 
                c(-413.7827081, -187.64104164, 19.02772528, -306.57808259, -21844.0631514, -175.90487759, 2254.58576331, 899.12913933, -93.0270774, 1385.87089749, 640.13813427, -355.5590459, 1418.45732112, 807.0376631, 2296.85216508, 1521.87542022), 
                c(-402.6264939, 35.21654172, -205.98042191, -308.77047073, -21294.2253707, -166.54859374, 2210.29974341, -171.78219472, 996.52858573, 1403.891274, -509.51606881, 666.96803064, 1384.39801595, 836.8177412, 1521.87542022, 2334.12612433)
                )
    test <- information(eUN.lmm, effects = "all", p = newp, transform.sigma = "log", transform.k = "none", transform.rho = "none")
    expect_equal(as.double(test), as.double(GS), tol = 1e-6)

    ## lava transformation
    newp.2 <- newp; newp.2["sigma"] <- newp["sigma"]^2
    ## GS <- -hessian(func = function(p){p["sigma"] <- sqrt(p["sigma"]); logLik(eUN.lmm, p = p)}, x = newp.2)
    GS <- cbind(c(241.4276261, 64.50068304, 62.84798106, 40.45548526, 12239.25904978, 101.39960062, -1171.9455677, -605.09410326, -591.34394345, -356.6684903, -614.45166963, -597.89210276, -411.70248495, -600.92192108, -413.7827081, -402.6264939), 
                c(64.50068304, 73.21443141, 3.49099427, -9.33761434, 3269.88504359, 27.09028677, -315.17277323, -645.9497367, -16.4036038, 41.06398553, -332.18747695, -4.58242056, 54.95083336, -358.00042572, -187.64104164, 35.21654172), 
                c(62.84798106, 3.49099427, 73.4157247, -5.92613477, 3186.10076536, 26.39615222, -307.58784494, -16.36801713, -641.00001661, 26.06133703, 17.95060609, -311.03354514, 52.78431089, -366.97750658, 19.02772528, -205.98042191), 
                c(40.45548526, -9.33761434, -5.92613477, 65.9452637, 2050.9052171, 16.99130389, -198.45423622, 43.78071542, 27.84592494, -468.76477714, 87.12601096, 69.89098536, -301.62532094, 71.84055748, -306.57808259, -308.77047073), 
                c(12239.25904978, 3269.88504359, 3186.10076536, 2050.9052171, 649235.8381143, 5035.57247599, -62055.00226189, -31952.35290987, -31341.07385733, -18837.71167947, -32546.69377367, -31728.40145323, -21802.76288715, -31788.23696115, -21844.0631514, -21294.2253707), 
                c(101.39960062, 27.09028677, 26.39615222, 16.99130389, 5035.57247599, 101.39960056, -487.88477474, -257.40463535, -238.49688193, -150.98532049, -258.3743199, -244.46171114, -173.04102013, -248.67924617, -175.90487759, -166.54859374), 
                c(-1171.9455677, -315.17277323, -307.58784494, -198.45423622, -62055.00226189, -487.88477474, 6010.43004435, 1588.41990599, 1564.29393588, 959.39843355, 1539.00460522, 1498.04475893, 1023.15699407, 1527.88256107, 1036.18607886, 1015.83264726), 
                c(-605.09410326, -645.9497367, -16.36801713, 43.78071542, -31952.35290987, -257.40463535, 1588.41990599, 4698.45923152, 79.52612717, -200.35755175, 1618.75890508, 22.29548051, -267.99319578, 1744.82398992, 899.12913933, -171.78219472), 
                c(-591.34394345, -16.4036038, -641.00001661, 27.84592494, -31341.07385733, -238.49688193, 1564.29393588, 79.52612717, 4680.85329314, -127.68986192, -87.29462845, 1519.58952368, -258.26729124, 1792.94480904, -93.0270774, 996.52858573), 
                c(-356.6684903, 41.06398553, 26.06133703, -468.76477714, -18837.71167947, -150.98532049, 959.39843355, -200.35755175, -127.68986192, 3035.28503039, -397.77172353, -319.5137104, 1373.54307929, -327.56863472, 1385.87089749, 1403.891274), 
                c(-614.45166963, -332.18747695, 17.95060609, 87.12601096, -32546.69377367, -258.3743199, 1539.00460522, 1618.75890508, -87.29462845, -397.77172353, 3397.00621427, 1570.20497911, 669.40246524, 1521.79966982, 640.13813427, -509.51606881), 
                c(-597.89210276, -4.58242056, -311.03354514, 69.89098536, -31728.40145323, -244.46171114, 1498.04475893, 22.29548051, 1519.58952368, -319.5137104, 1570.20497911, 3105.35116462, 663.7000477, 1566.17475066, -355.5590459, 666.96803064), 
                c(-411.70248495, 54.95083336, 52.78431089, -301.62532094, -21802.76288715, -173.04102013, 1023.15699407, -267.99319578, -258.26729124, 1373.54307929, 669.40246524, 663.7000477, 2294.20243311, -527.6926834, 1418.45732112, 1384.39801595), 
                c(-600.92192108, -358.00042572, -366.97750658, 71.84055748, -31788.23696115, -248.67924617, 1527.88256107, 1744.82398992, 1792.94480904, -327.56863472, 1521.79966982, 1566.17475066, -527.6926834, 3602.07056013, 807.0376631, 836.8177412), 
                c(-413.7827081, -187.64104164, 19.02772528, -306.57808259, -21844.0631514, -175.90487759, 1036.18607886, 899.12913933, -93.0270774, 1385.87089749, 640.13813427, -355.5590459, 1418.45732112, 807.0376631, 2296.85216508, 1521.87542022), 
                c(-402.6264939, 35.21654172, -205.98042191, -308.77047073, -21294.2253707, -166.54859374, 1015.83264726, -171.78219472, 996.52858573, 1403.891274, -509.51606881, 666.96803064, 1384.39801595, 836.8177412, 1521.87542022, 2334.12612433)
                )
    test <- information(eUN.lmm, effects = "all", p = newp, transform.sigma = "square", transform.k = "none", transform.rho = "none")
    expect_equal(as.double(test), as.double(GS), tol = 1e-6)

    ## cov transformation

    ## no transformation 
    newp <- coef(eUN.lmm, effects = "all", transform.sigma = "none", transform.k = "none", transform.rho = "none")+0.1
    ## GS <- -hessian(func = function(p){logLik(eUN.lmm, p = p)}, x = newp)
    GS <- cbind(c(241.4276261, 64.50068304, 62.84798106, 40.45548526, 12239.25904978, 101.39960062, -2444.76400233, -605.09410326, -591.34394345, -356.6684903, -614.45166963, -597.89210276, -411.70248495, -600.92192108, -413.7827081, -402.6264939), 
                c(64.50068304, 73.21443141, 3.49099427, -9.33761434, 3269.88504359, 27.09028677, -657.47341151, -645.9497367, -16.4036038, 41.06398553, -332.18747695, -4.58242056, 54.95083336, -358.00042572, -187.64104164, 35.21654172), 
                c(62.84798106, 3.49099427, 73.4157247, -5.92613477, 3186.10076536, 26.39615222, -641.65069748, -16.36801713, -641.00001661, 26.06133703, 17.95060609, -311.03354514, 52.78431089, -366.97750658, 19.02772528, -205.98042191), 
                c(40.45548526, -9.33761434, -5.92613477, 65.9452637, 2050.9052171, 16.99130389, -413.99002328, 43.78071542, 27.84592494, -468.76477714, 87.12601096, 69.89098536, -301.62532094, 71.84055748, -306.57808259, -308.77047073), 
                c(12239.25904978, 3269.88504359, 3186.10076536, 2050.9052171, 649235.8381143, 5035.57247599, -129451.26452782, -31952.35290987, -31341.07385733, -18837.71167947, -32546.69377367, -31728.40145323, -21802.76288715, -31788.23696115, -21844.0631514, -21294.2253707), 
                c(101.39960062, 27.09028677, 26.39615222, 16.99130389, 5035.57247599, 101.39960056, -1017.76325408, -257.40463535, -238.49688193, -150.98532049, -258.3743199, -244.46171114, -173.04102013, -248.67924617, -175.90487759, -166.54859374), 
                c(-2444.76400233, -657.47341151, -641.65069748, -413.99002328, -129451.26452782, -1017.76325408, 19797.77242215, 3313.55987297, 3263.23133819, 2001.37516506, 3210.47594884, 3125.03071949, 2134.38017671, 3187.27455335, 2161.55979863, 2119.10105462), 
                c(-605.09410326, -645.9497367, -16.36801713, 43.78071542, -31952.35290987, -257.40463535, 3313.55987297, 4698.45923152, 79.52612717, -200.35755175, 1618.75890508, 22.29548051, -267.99319578, 1744.82398992, 899.12913933, -171.78219472), 
                c(-591.34394345, -16.4036038, -641.00001661, 27.84592494, -31341.07385733, -238.49688193, 3263.23133819, 79.52612717, 4680.85329314, -127.68986192, -87.29462845, 1519.58952368, -258.26729124, 1792.94480904, -93.0270774, 996.52858573), 
                c(-356.6684903, 41.06398553, 26.06133703, -468.76477714, -18837.71167947, -150.98532049, 2001.37516506, -200.35755175, -127.68986192, 3035.28503039, -397.77172353, -319.5137104, 1373.54307929, -327.56863472, 1385.87089749, 1403.891274), 
                c(-614.45166963, -332.18747695, 17.95060609, 87.12601096, -32546.69377367, -258.3743199, 3210.47594884, 1618.75890508, -87.29462845, -397.77172353, 3397.00621427, 1570.20497911, 669.40246524, 1521.79966982, 640.13813427, -509.51606881), 
                c(-597.89210276, -4.58242056, -311.03354514, 69.89098536, -31728.40145323, -244.46171114, 3125.03071949, 22.29548051, 1519.58952368, -319.5137104, 1570.20497911, 3105.35116462, 663.7000477, 1566.17475066, -355.5590459, 666.96803064), 
                c(-411.70248495, 54.95083336, 52.78431089, -301.62532094, -21802.76288715, -173.04102013, 2134.38017671, -267.99319578, -258.26729124, 1373.54307929, 669.40246524, 663.7000477, 2294.20243311, -527.6926834, 1418.45732112, 1384.39801595), 
                c(-600.92192108, -358.00042572, -366.97750658, 71.84055748, -31788.23696115, -248.67924617, 3187.27455335, 1744.82398992, 1792.94480904, -327.56863472, 1521.79966982, 1566.17475066, -527.6926834, 3602.07056013, 807.0376631, 836.8177412), 
                c(-413.7827081, -187.64104164, 19.02772528, -306.57808259, -21844.0631514, -175.90487759, 2161.55979863, 899.12913933, -93.0270774, 1385.87089749, 640.13813427, -355.5590459, 1418.45732112, 807.0376631, 2296.85216508, 1521.87542022), 
                c(-402.6264939, 35.21654172, -205.98042191, -308.77047073, -21294.2253707, -166.54859374, 2119.10105462, -171.78219472, 996.52858573, 1403.891274, -509.51606881, 666.96803064, 1384.39801595, 836.8177412, 1521.87542022, 2334.12612433)
                )
    test <- information(eUN.lmm, effects = "all", p = newp, transform.sigma = "none", transform.k = "none", transform.rho = "none")
    expect_equal(as.double(test), as.double(GS), tol = 1e-6) 

    ## ** degree of freedom
    test <- model.tables(eUN.lmm, effects = "all")
    expect_equal(test$df,
                 c(118.415, 99.001, 99, 98.999, 96.184, 94.178, 99.452, 198.732, 199.05, 198.625, 49.522, 49.276, 49.593, 50.304, 49.359, 49.47),
                 tol = 1e-2)

    ## ** anova
    eUN.lmm_anova <- anova(eUN.lmm, effects = "all", ci = TRUE)$multivariate
    expect_equal(eUN.lmm_anova[eUN.lmm_anova$type=="mu","df.denom"], c(99.00144, 96.18348, 94.17755), tol = 1e-1)
    expect_equal(eUN.lmm_anova[eUN.lmm_anova$type=="k","df.denom"], c(189.236), tol = 1e-1)
    expect_equal(eUN.lmm_anova[eUN.lmm_anova$type=="rho","df.denom"], c(20.66968), tol = 1e-1)

    ## ** getVarCov
    Omega.GS <- matrix(c(0.88931784, -0.04684082, 0.0074039, 0.04333597, -0.04684082, 0.94870934, -0.12406129, 0.037128, 0.0074039, -0.12406129, 0.94419983, -0.00623928, 0.04333597, 0.037128, -0.00623928, 1.09138882), 
                       nrow = 4, 
                       ncol = 4, 
                       dimnames = list(c("Y1", "Y2", "Y3", "Y4"),c("Y1", "Y2", "Y3", "Y4")) 
                       ) 
    expect_equivalent(sigma(eUN.lmm), Omega.GS, tol = 1e-5)
})

## * Stratified compound symmetry structure
test_that("Stratified compound symmetry structure (REML)",{

    ## ** fit
    eSCS0.lmm <- lmm(Y ~ gender, repetition = gender~visit|id, structure = "CS", data = dL, trace = 0, method = "REML") ## just to check that it runs
    eSCS.lmm <- lmm(Y ~ (visit + age)*gender, repetition = gender~visit|id, structure = "CS", data = dL, trace = 0, method = "REML")
    eSCS.gls <- list(male=gls(Y ~ visit + age, correlation = corCompSymm(form=~1|id), data = dL[dL$gender=="male",], method = "REML"),
                     female=gls(Y ~ visit + age, correlation = corCompSymm(form=~1|id), data = dL[dL$gender=="female",], method = "REML"))

    ## ** iteration
    expect_equal(eSCS0.lmm$opt$n.iter,3)
    expect_equal(eSCS.lmm$opt$n.iter,4)

    ## ** coef
    coef(eSCS.lmm, transform.rho = "cov")
    coef(eSCS.lmm, transform.sigma = "square")

    ## ** logLikelihood
    expect_equal(logLik(eSCS.lmm), as.double(logLik(eSCS.gls$male)+logLik(eSCS.gls$female)), tol = 1e-6)

    ## ** score
    test <- score(eSCS.lmm, effects = "all", transform.sigma = "none", transform.k = "none", transform.rho = "none")
    ## GS <- jacobian(func = function(p){logLik(eSCS.lmm, p = p)},
    ##                x = coef(eSCS.lmm, effects = "all", transform.sigma = "none", transform.k = "none", transform.rho = "none"))
    GS <- rbind(c(0, -1e-08, 0, 0, 0, 0, 0, 0, 0, -7.2e-07, 4.57e-06, 1.351e-05, 2.924e-05, 4.737e-05))
    expect_true(all(abs(test) < 1e-3))
    expect_equal(as.double(GS), as.double(test), tol = 1e-5)

    ## no transformation
    newp <- coef(eSCS.lmm, effects = "all", transform.sigma = "none", transform.k = "none", transform.rho = "none")+0.1
    ## GS <- jacobian(func = function(p){logLik(eSCS.lmm, p = p)}, x = newp)
    GS <- rbind(c(-1955.49435178, -492.45269369, -492.45269438, -492.45269347, -103153.13334405, -971.75283595, -245.32412161, -245.32412097, -245.32411962, -50018.28657107, 4921.78267544, 9821.6312136, 8008.6677653, 9633.44513632))
    test <- score(eSCS.lmm, effects = "all", p = newp, transform.sigma = "none", transform.k = "none", transform.rho = "none")
    expect_equal(as.double(test), as.double(GS), tol = 1e-6)

    ## log transformation
    newp.log <- newp; newp.log[c("sigma:male","sigma:female")] <- log(newp[c("sigma:male","sigma:female")])
    ## GS <- jacobian(func = function(p){p[c("sigma:male","sigma:female")] <- exp(p[c("sigma:male","sigma:female")]); logLik(eSCS.lmm, p = p)}, x = newp.log)
    GS <- rbind(c(-1955.49435178, -492.45269369, -492.45269438, -492.45269347, -103153.13334405, -971.75283595, -245.32412161, -245.32412097, -245.32411962, -50018.28657107, 5447.52101479, 10322.19190469, 8008.6677653, 9633.44513632)
                )
    test <- score(eSCS.lmm, effects = "all", p = newp, transform.sigma = "log", transform.k = "none", transform.rho = "none")
    expect_equal(as.double(test), as.double(GS), tol = 1e-6)

    ## lava transformation
    newp.2 <- newp; newp.2[c("sigma:male","sigma:female")] <- newp[c("sigma:male","sigma:female")]^2
    ## GS <- jacobian(func = function(p){p[c("sigma:male","sigma:female")] <- sqrt(p[c("sigma:male","sigma:female")]); logLik(eSCS.lmm, p = p)}, x = newp.2)
    GS <- rbind(c(-1955.49435178, -492.45269369, -492.45269438, -492.45269347, -103153.13334405, -971.75283595, -245.32412161, -245.32412097, -245.32411962, -50018.28657107, 2223.39157933, 4672.67226777, 8008.6677653, 9633.44513632)
                )
    test <- score(eSCS.lmm, effects = "all", p = newp, transform.sigma = "square", transform.k = "none", transform.rho = "none")
    expect_equal(as.double(test), as.double(GS), tol = 1e-6)

    ## cov transformation
    newp.cov <- newp.2; newp.cov[c("rho:male","rho:female")] <- c(newp["rho:male"]*newp["sigma:male"]^2, newp["rho:female"]*newp["sigma:female"]^2)
    ## GS <- jacobian(func = function(p){
    ##     p[c("sigma:male","rho:male","sigma:female","rho:female")] <- c(sqrt(p["sigma:male"]),p["rho:male"]/p["sigma:male"],sqrt(p["sigma:female"]),p["rho:female"]/p["sigma:female"])
    ##     logLik(eSCS.lmm, p = p)
    ## }, x = newp.cov)
    GS <- rbind(c(-1955.49435178, -492.45269369, -492.45269438, -492.45269347, -103153.13334405, -971.75283595, -245.32412161, -245.32412097, -245.32411962, -50018.28657107, 2170.98071956, 2901.02998954, 6537.43397205, 8721.77776846)
                )
    test <- score(eSCS.lmm, p = newp, effects = "all", transform.rho = "cov")
    expect_equal(as.double(test), as.double(GS), tol = 1e-6)

    ## ** information
    test <- information(eSCS.lmm, effects = "all", transform.sigma = "none", transform.k = "none", transform.rho = "none")
    ## GS <- -hessian(func = function(p){logLik(eSCS.lmm, p = p)},
    ##                x = coef(eSCS.lmm, effects = "all", transform.sigma = "none", transform.k = "none", transform.rho = "none"))
    GS <- cbind(c(457.97114389, 114.49278603, 114.49278679, 114.49278621, 23307.04647392, 141.87716384, 35.46929092, 35.46929103, 35.46929067, 7045.71553977, 0, 2e-08, 1e-08, 9e-08), 
                c(114.49278603, 106.75872674, 2.57801984, 2.57801987, 5826.76161842, 35.46929096, 47.70481528, -4.07850809, -4.07850814, 1761.42888422, 0, 0, -2e-08, -1e-08), 
                c(114.49278679, 2.57801984, 106.75872668, 2.57801985, 5826.76161867, 35.46929096, -4.0785081, 47.70481536, -4.07850813, 1761.42888504, 0, 0, -6e-08, 2e-08), 
                c(114.49278621, 2.57801987, 2.57801985, 106.75872673, 5826.76161857, 35.46929095, -4.07850812, -4.07850811, 47.70481527, 1761.42888419, 0, 0, 0, 2e-08), 
                c(23307.04647392, 5826.76161842, 5826.76161867, 5826.76161857, 1242852.77049857, 7045.71553626, 1761.42888404, 1761.42888415, 1761.42888407, 363108.48111926, -4.4e-07, -1.9e-07, -5e-08, 1.1e-07), 
                c(141.87716384, 35.46929096, 35.46929096, 35.46929095, 7045.71553626, 141.8771638, 35.46929095, 35.46929096, 35.46929095, 7045.71553625, 0, 0, 0, 1e-08), 
                c(35.46929092, 47.70481528, -4.0785081, -4.07850812, 1761.42888404, 35.46929095, 47.70481529, -4.07850811, -4.07850812, 1761.42888388, 0, 0, 0, 0), 
                c(35.46929103, -4.07850809, 47.70481536, -4.07850811, 1761.42888415, 35.46929096, -4.07850811, 47.70481526, -4.07850809, 1761.42888438, 0, 0, 0, 3e-08), 
                c(35.46929067, -4.07850814, -4.07850813, 47.70481527, 1761.42888407, 35.46929095, -4.07850812, -4.07850809, 47.70481531, 1761.42888374, 0, 0, 0, 3e-08), 
                c(7045.71553977, 1761.42888422, 1761.42888504, 1761.42888419, 363108.48111926, 7045.71553625, 1761.42888388, 1761.42888438, 1761.42888374, 363108.48111231, 0, 1e-08, 0, 5.4e-07), 
                c(0, 0, 0, 0, -4.4e-07, 0, 0, 0, 0, 0, 447.87140481, 0, 74.92116221, 1e-08), 
                c(2e-08, 0, 0, 0, -1.9e-07, 0, 0, 0, 0, 1e-08, 0, 360.48605438, 0, -47.84326993), 
                c(1e-08, -2e-08, -6e-08, 0, -5e-08, 0, 0, 0, 0, 0, 74.92116221, 0, 552.38940001, 1e-08), 
                c(9e-08, -1e-08, 2e-08, 2e-08, 1.1e-07, 1e-08, 0, 3e-08, 3e-08, 5.4e-07, 1e-08, -47.84326993, 1e-08, 181.4441314)
                )
    expect_equal(as.double(GS), as.double(test), tol = 1e-5)

    ## no transformation
    newp <- coef(eSCS.lmm, effects = "all", transform.sigma = "none", transform.k = "none", transform.rho = "none")+0.1
    ## GS <- -hessian(func = function(p){logLik(eSCS.lmm, p = p)}, x = newp)
    GS <- cbind(c(279.44142308, 69.86035586, 69.86035606, 69.86035566, 14207.15684851, 94.50886027, 23.62721448, 23.62721507, 23.6272162, 4693.37367605, -1777.60193561, -1849.25800177, -2881.91137511, -1811.41069509), 
                c(69.86035586, 89.04957008, -6.39640479, -6.3964049, 3551.78921234, 23.62721507, 41.69548432, -6.02275625, -6.02275504, 1173.34341889, -446.55656124, -466.85492136, -719.27500286, -449.85857459), 
                c(69.86035606, -6.39640479, 89.04957038, -6.39640505, 3551.7892123, 23.62721507, -6.02275701, 41.69548455, -6.02275647, 1173.34341887, -446.55656131, -466.85492138, -719.27501045, -449.8585746), 
                c(69.86035566, -6.3964049, -6.39640505, 89.04957033, 3551.78921194, 23.62721504, -6.02275664, -6.02275649, 41.69548428, 1173.34341854, -446.5565613, -466.85492139, -719.27500744, -449.85857472), 
                c(14207.15684851, 3551.78921234, 3551.7892123, 3551.78921194, 756577.3760879, 4693.37367455, 1173.34341811, 1173.34341859, 1173.34341892, 241878.028924, -96013.64275974, -95185.4352719, -155660.7262142, -93237.34996608), 
                c(94.50886027, 23.62721507, 23.62721507, 23.62721504, 4693.37367455, 94.50886034, 23.62721498, 23.62721509, 23.62721535, 4693.37367477, 1e-08, -1849.25800171, -1.8e-06, -1811.41069515), 
                c(23.62721448, 41.69548432, -6.02275701, -6.02275664, 1173.34341811, 23.62721498, 41.69548512, -6.02275649, -6.02275528, 1173.34341866, 0, -466.8549214, -9.15e-06, -449.85857465), 
                c(23.62721507, -6.02275625, 41.69548455, -6.02275649, 1173.34341859, 23.62721509, -6.02275649, 41.69548439, -6.02275627, 1173.3434187, 0, -466.85492137, 0, -449.85857465), 
                c(23.6272162, -6.02275504, -6.02275647, 41.69548428, 1173.34341892, 23.62721535, -6.02275528, -6.02275627, 41.69547097, 1173.34342243, 4.7e-07, -466.8549213, -4e-08, -449.85857172), 
                c(4693.37367605, 1173.34341889, 1173.34341887, 1173.34341854, 241878.028924, 4693.37367477, 1173.34341866, 1173.3434187, 1173.34342243, 241878.02892262, 6e-08, -95185.43527173, -7.8e-07, -93237.34996538), 
                c(-1777.60193561, -446.55656124, -446.55656131, -446.5565613, -96013.64275974, 1e-08, 0, 0, 4.7e-07, 6e-08, 13710.94732216, 0, 14463.98413491, 7e-08), 
                c(-1849.25800177, -466.85492136, -466.85492138, -466.85492139, -95185.4352719, -1849.25800171, -466.8549214, -466.85492137, -466.8549213, -95185.43527173, 0, 28331.18238373, 1e-06, 18256.64690681), 
                c(-2881.91137511, -719.27500286, -719.27501045, -719.27500744, -155660.7262142, -1.8e-06, -9.15e-06, 0, -4e-08, -7.8e-07, 14463.98413491, 1e-06, 47191.56357782, -1e-08), 
                c(-1811.41069509, -449.85857459, -449.8585746, -449.85857472, -93237.34996608, -1811.41069515, -449.85857465, -449.85857465, -449.85857172, -93237.34996538, 7e-08, 18256.64690681, -1e-08, 36048.89557121)
                )
    test <- information(eSCS.lmm, effects = "all", p = newp, transform.sigma = "none", transform.k = "none", transform.rho = "none")
    expect_equal(as.double(test), as.double(GS), tol = 1e-6)

    ## log transformation
    newp.log <- newp; newp.log[c("sigma:male","sigma:female")] <- log(newp[c("sigma:male","sigma:female")])
    ## GS <- -hessian(func = function(p){p[c("sigma:male","sigma:female")] <- exp(p[c("sigma:male","sigma:female")]); logLik(eSCS.lmm, p = p)}, x = newp.log)
    GS <- cbind(c(279.44142308, 69.86035586, 69.86035606, 69.86035566, 14207.15684851, 94.50886027, 23.62721448, 23.62721507, 23.6272162, 4693.37367605, -1967.48303217, -1943.50567329, -2881.91137511, -1811.41069509), 
                c(69.86035586, 89.04957008, -6.39640479, -6.3964049, 3551.78921234, 23.62721507, 41.69548432, -6.02275625, -6.02275504, 1173.34341889, -494.2571446, -490.64824157, -719.27500286, -449.85857459), 
                c(69.86035606, -6.39640479, 89.04957038, -6.39640505, 3551.7892123, 23.62721507, -6.02275701, 41.69548455, -6.02275647, 1173.34341887, -494.25714472, -490.64824354, -719.27501045, -449.8585746), 
                c(69.86035566, -6.3964049, -6.39640505, 89.04957033, 3551.78921194, 23.62721504, -6.02275664, -6.02275649, 41.69548428, 1173.34341854, -494.25714488, -490.64824263, -719.27500744, -449.85857472), 
                c(14207.15684851, 3551.78921234, 3551.7892123, 3551.78921194, 756577.3760879, 4693.37367455, 1173.34341811, 1173.34341859, 1173.34341892, 241878.028924, -106269.69354432, -100036.57314084, -155660.7262142, -93237.34996608), 
                c(94.50886027, 23.62721507, 23.62721507, 23.62721504, 4693.37367455, 94.50886034, 23.62721498, 23.62721509, 23.62721535, 4693.37367477, -1e-08, -1943.50567227, -1.8e-06, -1811.41069515), 
                c(23.62721448, 41.69548432, -6.02275701, -6.02275664, 1173.34341811, 23.62721498, 41.69548512, -6.02275649, -6.02275528, 1173.34341866, -7.2e-07, -490.64824247, -9.15e-06, -449.85857465), 
                c(23.62721507, -6.02275625, 41.69548455, -6.02275649, 1173.34341859, 23.62721509, -6.02275649, 41.69548439, -6.02275627, 1173.3434187, 2e-08, -490.64824221, 0, -449.85857465), 
                c(23.6272162, -6.02275504, -6.02275647, 41.69548428, 1173.34341892, 23.62721535, -6.02275528, -6.02275627, 41.69547097, 1173.34342243, 0, -490.64824042, -4e-08, -449.85857172), 
                c(4693.37367605, 1173.34341889, 1173.34341887, 1173.34341854, 241878.028924, 4693.37367477, 1173.34341866, 1173.3434187, 1173.34342243, 241878.02892262, -6e-08, -100036.57314116, -7.8e-07, -93237.34996538), 
                c(-1967.48303217, -494.2571446, -494.25714472, -494.25714488, -106269.69354432, -1e-08, -7.2e-07, 2e-08, 0, -6e-08, 11349.04203083, -3.34e-06, 16009.00784061, -7e-08), 
                c(-1943.50567329, -490.64824157, -490.64824354, -490.64824263, -100036.57314084, -1943.50567227, -490.64824247, -490.64824221, -490.64824042, -100036.57314116, -3.34e-06, 20970.38381485, -2.11e-05, 19187.09924985), 
                c(-2881.91137511, -719.27500286, -719.27501045, -719.27500744, -155660.7262142, -1.8e-06, -9.15e-06, 0, -4e-08, -7.8e-07, 16009.00784061, -2.11e-05, 47191.56357782, -1e-08), 
                c(-1811.41069509, -449.85857459, -449.8585746, -449.85857472, -93237.34996608, -1811.41069515, -449.85857465, -449.85857465, -449.85857172, -93237.34996538, -7e-08, 19187.09924985, -1e-08, 36048.89557121)
                )
    test <- information(eSCS.lmm, effects = "all", p = newp, transform.sigma = "log", transform.k = "none", transform.rho = "none")
    expect_equal(as.double(test), as.double(GS), tol = 1e-6)

    ## lava transformation
    newp.2 <- newp; newp.2[c("sigma:male","sigma:female")] <- newp[c("sigma:male","sigma:female")]^2
    ## GS <- -hessian(func = function(p){p[c("sigma:male","sigma:female")] <- sqrt(p[c("sigma:male","sigma:female")]); logLik(eSCS.lmm, p = p)}, x = newp.2)
    GS <- cbind(c(279.44142308, 69.86035586, 69.86035606, 69.86035566, 14207.15684851, 94.50886027, 23.62721448, 23.62721507, 23.6272162, 4693.37367605, -803.02309881, -879.7903722, -2881.91137511, -1811.41069509), 
                c(69.86035586, 89.04957008, -6.39640479, -6.3964049, 3551.78921234, 23.62721507, 41.69548432, -6.02275625, -6.02275504, 1173.34341889, -201.72977221, -222.10771275, -719.27500286, -449.85857459), 
                c(69.86035606, -6.39640479, 89.04957038, -6.39640505, 3551.7892123, 23.62721507, -6.02275701, 41.69548455, -6.02275647, 1173.34341887, -201.72977227, -222.10771285, -719.27501045, -449.8585746), 
                c(69.86035566, -6.3964049, -6.39640505, 89.04957033, 3551.78921194, 23.62721504, -6.02275664, -6.02275649, 41.69548428, 1173.34341854, -201.72977226, -222.10771274, -719.27500744, -449.85857472), 
                c(14207.15684851, 3551.78921234, 3551.7892123, 3551.78921194, 756577.3760879, 4693.37367455, 1173.34341811, 1173.34341859, 1173.34341892, 241878.028924, -43373.69991689, -45284.77337606, -155660.7262142, -93237.34996608), 
                c(94.50886027, 23.62721507, 23.62721507, 23.62721504, 4693.37367455, 94.50886034, 23.62721498, 23.62721509, 23.62721535, 4693.37367477, 0, -879.79037216, -1.8e-06, -1811.41069515), 
                c(23.62721448, 41.69548432, -6.02275701, -6.02275664, 1173.34341811, 23.62721498, 41.69548512, -6.02275649, -6.02275528, 1173.34341866, -6e-08, -222.10771277, -9.15e-06, -449.85857465), 
                c(23.62721507, -6.02275625, 41.69548455, -6.02275649, 1173.34341859, 23.62721509, -6.02275649, 41.69548439, -6.02275627, 1173.3434187, 0, -222.10771274, 0, -449.85857465), 
                c(23.6272162, -6.02275504, -6.02275647, 41.69548428, 1173.34341892, 23.62721535, -6.02275528, -6.02275627, 41.69547097, 1173.34342243, 0, -222.10771293, -4e-08, -449.85857172), 
                c(4693.37367605, 1173.34341889, 1173.34341887, 1173.34341854, 241878.028924, 4693.37367477, 1173.34341866, 1173.3434187, 1173.34342243, 241878.02892262, 6e-08, -45284.77337601, -7.8e-07, -93237.34996538), 
                c(-803.02309881, -201.72977221, -201.72977227, -201.72977226, -43373.69991689, 0, -6e-08, 0, 0, 6e-08, 3705.51531106, -1e-08, 6534.03505258, -3e-08), 
                c(-879.7903722, -222.10771275, -222.10771285, -222.10771274, -45284.77337606, -879.79037216, -222.10771277, -222.10771274, -222.10771293, -45284.77337601, -1e-08, 8527.74605114, -8e-08, 8685.65779444), 
                c(-2881.91137511, -719.27500286, -719.27501045, -719.27500744, -155660.7262142, -1.8e-06, -9.15e-06, 0, -4e-08, -7.8e-07, 6534.03505258, -8e-08, 47191.56357782, -1e-08), 
                c(-1811.41069509, -449.85857459, -449.8585746, -449.85857472, -93237.34996608, -1811.41069515, -449.85857465, -449.85857465, -449.85857172, -93237.34996538, -3e-08, 8685.65779444, -1e-08, 36048.89557121)
                )
    test <- information(eSCS.lmm, effects = "all", p = newp, transform.sigma = "square", transform.k = "none", transform.rho = "none")
    expect_equal(as.double(test), as.double(GS), tol = 1e-6)

    ## cov transformation
    newp.cov <- newp.2; newp.cov[c("rho:male","rho:female")] <- c(newp["rho:male"]*newp["sigma:male"]^2, newp["rho:female"]*newp["sigma:female"]^2)
    ## GS <- -hessian(func = function(p){
    ##     p[c("sigma:male","rho:male","sigma:female","rho:female")] <- c(sqrt(p["sigma:male"]),p["rho:male"]/p["sigma:male"],sqrt(p["sigma:female"]),p["rho:female"]/p["sigma:female"])
    ##     logLik(eSCS.lmm, p = p)
    ## }, x = newp.cov)
    GS <- cbind(c(279.44142308, 69.86035586, 69.86035606, 69.86035566, 14207.15684851, 94.50886027, 23.62721448, 23.62721507, 23.6272162, 4693.37367605, -784.1631014, -546.66222064, -2352.48931285, -1639.98666169), 
                c(69.86035586, 89.04957008, -6.39640479, -6.3964049, 3551.78921234, 23.62721507, 41.69548432, -6.02275625, -6.02275504, 1173.34341889, -197.02264455, -139.3763057, -587.14045541, -407.28591478), 
                c(69.86035606, -6.39640479, 89.04957038, -6.39640505, 3551.7892123, 23.62721507, -6.02275701, 41.69548455, -6.02275647, 1173.34341887, -197.02264462, -139.37630569, -587.14046161, -407.28591479), 
                c(69.86035566, -6.3964049, -6.39640505, 89.04957033, 3551.78921194, 23.62721504, -6.02275664, -6.02275649, 41.69548428, 1173.34341854, -197.0226446, -139.37630576, -587.14045915, -407.2859149), 
                c(14207.15684851, 3551.78921234, 3551.7892123, 3551.78921194, 756577.3760879, 4693.37367455, 1173.34341811, 1173.34341859, 1173.34341892, 241878.028924, -42355.01457061, -28137.92416658, -127065.04371116, -84413.77249962), 
                c(94.50886027, 23.62721507, 23.62721507, 23.62721504, 4693.37367455, 94.50886034, 23.62721498, 23.62721509, 23.62721535, 4693.37367477, 1e-08, -546.66222061, -1.47e-06, -1639.98666178), 
                c(23.62721448, 41.69548432, -6.02275701, -6.02275664, 1173.34341811, 23.62721498, 41.69548512, -6.02275649, -6.02275528, 1173.34341866, 0, -139.37630563, -7.47e-06, -407.28591483), 
                c(23.62721507, -6.02275625, 41.69548455, -6.02275649, 1173.34341859, 23.62721509, -6.02275649, 41.69548439, -6.02275627, 1173.3434187, 0, -139.3763057, 0, -407.28591483), 
                c(23.6272162, -6.02275504, -6.02275647, 41.69548428, 1173.34341892, 23.62721535, -6.02275528, -6.02275627, 41.69547097, 1173.34342243, 1.9e-07, -139.37630543, -3e-08, -407.28591218), 
                c(4693.37367605, 1173.34341889, 1173.34341887, 1173.34341854, 241878.028924, 4693.37367477, 1173.34341866, 1173.3434187, 1173.34342243, 241878.02892262, 6e-08, -28137.92416645, -6.3e-07, -84413.77249899), 
                c(-784.1631014, -197.02264455, -197.02264462, -197.0226446, -42355.01457061, 1e-08, 0, 0, 1.9e-07, 6e-08, 3536.45003887, 0, 10418.07264568, 3e-08), 
                c(-546.66222064, -139.3763057, -139.37630569, -139.37630576, -28137.92416658, -546.66222061, -139.37630563, -139.3763057, -139.37630543, -28137.92416645, 0, 3344.32069652, 7.1e-07, 9757.88125453), 
                c(-2352.48931285, -587.14045541, -587.14046161, -587.14045915, -127065.04371116, -1.47e-06, -7.47e-06, 0, -3e-08, -6.3e-07, 10418.07264568, 7.1e-07, 31445.49551933, 0), 
                c(-1639.98666169, -407.28591478, -407.28591479, -407.2859149, -84413.77249962, -1639.98666178, -407.28591483, -407.28591483, -407.28591218, -84413.77249899, 3e-08, 9757.88125453, 0, 29548.72459835)
                )
    test <- information(eSCS.lmm, p = newp, effects = "all", transform.rho = "cov")
    expect_equal(as.double(test), as.double(GS), tol = 1e-6)

    ## ** degree of freedom
    test <- model.tables(eSCS.lmm, effects = "all")
    GS <- c(82.36345237, 171.0218177, 171.0218177, 171.0218177, 56.03889211, 81.14993151, 285.63027224, 285.63027224, 285.63027224, 63.97167501, 178.18701642, 88.36291507, 6.93663225, 8.08291417)
    expect_equal(test$df,GS, tol = 1e-3)
    
    ## ** anova
    eSCS.lmm_anova <- anova(eSCS.lmm, effects = "all")$multivariate
    expect_equal(eSCS.lmm_anova[eSCS.lmm_anova$type=="mu","df.denom"], c(171.0218177, 56.03889211, 81.14993152, 285.63027224, 63.97167502), tol = 1e-1)
    expect_equal(eSCS.lmm_anova[eSCS.lmm_anova$type=="rho","df.denom"], c(7.450136), tol = 1e-1)

    ## ** getVarCov
    Omega.GS <- list("female" = matrix(c(0.90433467, 0.09326275, 0.09326275, 0.09326275, 0.09326275, 0.90433467, 0.09326275, 0.09326275, 0.09326275, 0.09326275, 0.90433467, 0.09326275, 0.09326275, 0.09326275, 0.09326275, 0.90433467), 
                                  nrow = 4, 
                                  ncol = 4, 
                                  dimnames = list(c("Y1", "Y2", "Y3", "Y4"),c("Y1", "Y2", "Y3", "Y4")) 
                                  ),
                     "male" = matrix(c(1.01368386, -0.09324164, -0.09324164, -0.09324164, -0.09324164, 1.01368386, -0.09324164, -0.09324164, -0.09324164, -0.09324164, 1.01368386, -0.09324164, -0.09324164, -0.09324164, -0.09324164, 1.01368386), 
                                  nrow = 4, 
                                  ncol = 4, 
                                  dimnames = list(c("Y1", "Y2", "Y3", "Y4"),c("Y1", "Y2", "Y3", "Y4")) 
                                  ) 
                     )
    expect_equivalent(sigma(eSCS.lmm)[names(Omega.GS)], Omega.GS, tol = 1e-5)
})

## * Stratified unstructed covariance matrix
test_that("Stratified unstructured (REML)",{

    ## ** fit
    eSUN.lmm <- lmm(Y ~ (visit + age)*gender, repetition = gender~visit|id, structure = "UN", data = dL, trace = 0, method = "REML")
    eSUN.gls <- list(male=gls(Y ~ visit + age, correlation = corSymm(form=~1|id), weights = varIdent(form=~1|visit), data = dL[dL$gender=="male",], method = "REML"),
                     female=gls(Y ~ visit + age, correlation = corSymm(form=~1|id), weights = varIdent(form=~1|visit), data = dL[dL$gender=="female",], method = "REML"))
    coef(eSUN.lmm, transform.sigma = "none", transform.k = "sd", effects = "variance")
    coef(eSUN.lmm, transform.sigma = "none", effects = "variance")
    capture.output(summary(eSUN.lmm))
    ## anova(eSUN.lmm, effects = "all")

    ## ** iteration
    expect_equal(eSUN.lmm$opt$n.iter,5)

    ## ** coef
    coef(eSUN.lmm, transform.rho = "cov")
    coef(eSUN.lmm, transform.k = "var")
    coef(eSUN.lmm, transform.k = "log")

    ## lava trans
    fct.trans <- function(p, inverse = FALSE){
        newp <- p
        if(inverse){
            newp[c("sigma:male","k.Y2:male","k.Y3:male","k.Y4:male")] <- c(sqrt(newp["sigma:male"]), sqrt(newp[c("k.Y2:male","k.Y3:male","k.Y4:male")]/newp["sigma:male"]));
            newp[c("sigma:female","k.Y2:female","k.Y3:female","k.Y4:female")] <- c(sqrt(newp["sigma:female"]), sqrt(newp[c("k.Y2:female","k.Y3:female","k.Y4:female")]/newp["sigma:female"]));
        }else{
            newp[c("sigma:male","k.Y2:male","k.Y3:male","k.Y4:male")] <- p["sigma:male"]^2 * c(1,p[c("k.Y2:male","k.Y3:male","k.Y4:male")]^2);
            newp[c("sigma:female","k.Y2:female","k.Y3:female","k.Y4:female")] <- p["sigma:female"]^2 * c(1,p[c("k.Y2:female","k.Y3:female","k.Y4:female")]^2);
        }
        return(newp)
    }
    expect_equal(fct.trans(coef(eSUN.lmm, effects = "all", transform.sigma = "none", transform.k = "none", transform.rho = "none")),
                 coef(eSUN.lmm, effects = "all", transform.k = "var", transform.rho = "none", transform.names=FALSE))

    ## ** logLikelihood
    expect_equal(logLik(eSUN.lmm), as.double(logLik(eSUN.gls$male)+logLik(eSUN.gls$female)), tol = 1e-6)

    ## ** score
    test <- score(eSUN.lmm, effects = "all", transform.sigma = "none", transform.k = "none", transform.rho = "none")
    ## GS <- jacobian(func = function(p){logLik(eSUN.lmm, p = p)}, x = coef(eSUN.lmm, effects = "all", transform.sigma = "none", transform.k = "none", transform.rho = "none"))
    GS <- rbind(c(0, 0, 0, 0, 1.2e-07, 0, -3e-08, 0, 0, 0, 2.1e-07, 1.08e-06, -9e-08, -1e-08, 4.5e-07, -5.123e-05, 9.4e-06, 2.556e-05, -2.096e-05, 1.581e-05, 2.913e-05, -4.119e-05, 2.76e-06, 3.397e-05, -2.4e-07, -2.2e-07, 8.8e-07, -4e-08, 8.5e-07, 5.8e-07))
    expect_true(all(abs(test) < 1e-3))
    expect_equal(as.double(GS), as.double(test), tol = 1e-5)

    ## no transformation
    newp <- coef(eSUN.lmm, effects = "all", transform.sigma = "none", transform.k = "none", transform.rho = "none")+0.1
    ## GS <- jacobian(func = function(p){logLik(eSUN.lmm, p = p)}, x = newp)
    GS <- c("(Intercept)" = -1783.17328079, "visitY2" = -494.19593683, "visitY3" = -489.74409497, "visitY4" = -280.85528465, "age" = -94052.55029456, "genderfemale" = -895.1100431, "visitY2:genderfemale" = -228.57059312, "visitY3:genderfemale" = -305.42867252, "visitY4:genderfemale" = -112.07212363, "age:genderfemale" = -46076.161193, "sigma:male" = 4635.77219007, "sigma:female" = 9146.88022718, "k.Y2:male" = 1300.93756056, "k.Y3:male" = 876.88257525, "k.Y4:male" = 748.51536421, "k.Y2:female" = 2186.33054223, "k.Y3:female" = 3054.55980002, "k.Y4:female" = 1010.82864051, "rho(Y1,Y2):male" = 1625.93441559, "rho(Y1,Y3):male" = 1155.44499441, "rho(Y1,Y4):male" = 1113.83275036, "rho(Y2,Y3):male" = 1285.80100653, "rho(Y2,Y4):male" = 1238.78723381, "rho(Y3,Y4):male" = 881.61137318, "rho(Y1,Y2):female" = 1648.82180272, "rho(Y1,Y3):female" = 2154.97282953, "rho(Y1,Y4):female" = 892.39161387, "rho(Y2,Y3):female" = 2147.6597046, "rho(Y2,Y4):female" = 892.35942774, "rho(Y3,Y4):female" = 1162.00385284)
    test <- score(eSUN.lmm, effects = "all", p = newp, transform.sigma = "none", transform.k = "none", transform.rho = "none")
    expect_equal(test, GS, tol = 1e-6)

    ## log transformation
    newp.log <- newp; newp.log[eSUN.lmm$design$param$type %in% "sigma"] <- log(newp[eSUN.lmm$design$param$type %in% "sigma"])
    ## GS <- jacobian(func = function(p){p[eSUN.lmm$design$param$type %in% "sigma"] <- exp(p[eSUN.lmm$design$param$type %in% "sigma"]); logLik(eSUN.lmm, p = p)}, x = newp.log)
    GS <- c("(Intercept)" = -1783.17328079, "visitY2" = -494.19593683, "visitY3" = -489.74409497, "visitY4" = -280.85528465, "age" = -94052.55029456, "genderfemale" = -895.1100431, "visitY2:genderfemale" = -228.57059312, "visitY3:genderfemale" = -305.42867252, "visitY4:genderfemale" = -112.07212363, "age:genderfemale" = -46076.161193, "log(sigma):male" = 4885.74713566, "log(sigma):female" = 9483.46984917, "k.Y2:male" = 1300.93756056, "k.Y3:male" = 876.88257525, "k.Y4:male" = 748.51536421, "k.Y2:female" = 2186.33054223, "k.Y3:female" = 3054.55980002, "k.Y4:female" = 1010.82864051, "rho(Y1,Y2):male" = 1625.93441559, "rho(Y1,Y3):male" = 1155.44499441, "rho(Y1,Y4):male" = 1113.83275036, "rho(Y2,Y3):male" = 1285.80100653, "rho(Y2,Y4):male" = 1238.78723381, "rho(Y3,Y4):male" = 881.61137318, "rho(Y1,Y2):female" = 1648.82180272, "rho(Y1,Y3):female" = 2154.97282953, "rho(Y1,Y4):female" = 892.39161387, "rho(Y2,Y3):female" = 2147.6597046, "rho(Y2,Y4):female" = 892.35942774, "rho(Y3,Y4):female" = 1162.00385284)
    test <- score(eSUN.lmm, effects = "all", p = newp, transform.sigma = "log", transform.k = "none", transform.rho = "none")
    expect_equal(as.double(test), as.double(GS), tol = 1e-6)


    ## lava transformation
    newp.2 <- fct.trans(newp); 
    ## GS <- jacobian(func = function(p){logLik(eSUN.lmm, p = fct.trans(p,inverse = TRUE))}, x = newp.2)
    GS <- c("(Intercept)" = -1783.17328079, "visitY2" = -494.19593683, "visitY3" = -489.74409497, "visitY4" = -280.85528465, "age" = -94052.55029456, "genderfemale" = -895.1100431, "visitY2:genderfemale" = -228.57059312, "visitY3:genderfemale" = -305.42867252, "visitY4:genderfemale" = -112.07212363, "age:genderfemale" = -46076.161193, "sigma:male" = 665.79608398, "sigma:female" = 1218.87880933, "k.Y2:male" = 517.53206249, "k.Y3:male" = 340.71006511, "k.Y4:male" = 274.52072123, "k.Y2:female" = 924.93858688, "k.Y3:female" = 1332.7320395, "k.Y4:female" = 395.13143761, "rho(Y1,Y2):male" = 1625.93441559, "rho(Y1,Y3):male" = 1155.44499441, "rho(Y1,Y4):male" = 1113.83275036, "rho(Y2,Y3):male" = 1285.80100653, "rho(Y2,Y4):male" = 1238.78723381, "rho(Y3,Y4):male" = 881.61137318, "rho(Y1,Y2):female" = 1648.82180272, "rho(Y1,Y3):female" = 2154.97282953, "rho(Y1,Y4):female" = 892.39161387, "rho(Y2,Y3):female" = 2147.6597046, "rho(Y2,Y4):female" = 892.35942774, "rho(Y3,Y4):female" = 1162.00385284)
    test <- score(eSUN.lmm, effects = "all", p = newp, transform.k = "var", transform.rho = "none", transform.names = FALSE)
    expect_equal(as.double(test), as.double(GS), tol = 1e-6)

    ## cov transformation
    ## newp.cov <- newp; newp.cov[c("sigma","rho")] <- c(newp["sigma"]^2,newp["rho"]*newp["sigma"]^2)
    ## GS <- jacobian(func = function(p){p[c("sigma","rho")] <- c(sqrt(p["sigma"]),p["rho"]/p["sigma"]); logLik(eSUN.lmm, p = p)}, x = newp.cov)
    ## test <- score(eSUN.lmm, p = newp, transform.rho = "cov")
    ## expect_equal(as.double(test), as.double(GS), tol = 1e-6)

    ## ** information
    name.all <- names(coef(eSUN.lmm, effects = "all"))
    name.mu <- names(coef(eSUN.lmm, effects = "mean"))
    name.vcov <- names(coef(eSUN.lmm, effects = c("variance","covariance")))
    name.vcovM <- grep(":male",name.vcov,value=TRUE)
    name.vcovF <- grep(":female",name.vcov,value=TRUE)

    ## expect_equal(as.double(vcov(eSUN.lmm, effects = "mean")), as.double(Matrix::bdiag(lapply(eSUN.gls,vcov))), tol = 1e-6)
    test <- information(eSUN.lmm, effects = "all",
                        p = coef(eSUN.lmm, effects = "all", transform.sigma = "none", transform.k = "none", transform.rho = "none"),
                        transform.sigma = "none", transform.k = "none", transform.rho = "none")
    ## GS <- -hessian(func = function(p){logLik(eSUN.lmm, p = p)},
    ##                x = coef(eSUN.lmm, effects = "all", transform.sigma = "none", transform.k = "none", transform.rho = "none"))
    
    
    GS.mu <- matrix(c(481.58902581, 141.41457601, 126.49375658, 90.87132755, 24505.0293415, 151.42333094, 39.89451563, 53.34143654, 22.82307691, 7519.78462744, 141.41457601, 123.53334861, 16.52169518, -6.89363657, 7203.84822284, 39.89451563, 57.11914725, 7.0825439, -18.90835881, 1981.18852306, 126.49375658, 16.52169518, 110.82914455, -0.53727869, 6412.26409582, 53.34143654, 7.0825439, 52.49669252, -1.50781629, 2648.9716743, 90.87132755, -6.89363657, -0.53727869, 99.56459117, 4634.12494769, 22.82307691, -18.90835881, -1.50781629, 47.23106187, 1133.40937505, 24505.0293415, 7203.84822284, 6412.26409582, 4634.12494769, 1306448.44293862, 7519.78462744, 1981.18852306, 2648.9716743, 1133.40937505, 387540.13845115, 151.42333094, 39.89451563, 53.34143654, 22.82307691, 7519.78462744, 151.42333094, 39.89451563, 53.34143654, 22.82307691, 7519.78462744, 39.89451563, 57.11914725, 7.0825439, -18.90835881, 1981.18852306, 39.89451563, 57.11914725, 7.0825439, -18.90835881, 1981.18852306, 53.34143654, 7.0825439, 52.49669252, -1.50781629, 2648.9716743, 53.34143654, 7.0825439, 52.49669252, -1.50781629, 2648.9716743, 22.82307691, -18.90835881, -1.50781629, 47.23106187, 1133.40937505, 22.82307691, -18.90835881, -1.50781629, 47.23106187, 1133.40937505, 7519.78462744, 1981.18852306, 2648.9716743, 1133.40937505, 387540.13845115, 7519.78462744, 1981.18852306, 2648.9716743, 1133.40937505, 387540.13845115), 
                    nrow = 10, ncol = 10, dimnames = list(c("(Intercept)", "visitY2", "visitY3", "visitY4", "age", "genderfemale", "visitY2:genderfemale", "visitY3:genderfemale", "visitY4:genderfemale", "age:genderfemale"),c("(Intercept)", "visitY2", "visitY3", "visitY4", "age", "genderfemale", "visitY2:genderfemale", "visitY3:genderfemale", "visitY4:genderfemale", "age:genderfemale")) ) 
    expect_equal(as.double(test[name.mu,name.mu]),as.double(GS.mu), tol = 1e-5)

    GS.muvcov <- matrix(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -55.65356149, 0, 0, 0, 0, 0, 0, 0, 0, 0, -16.11419454, 0, 0, 0, 0, 0, 0, 0, 0, 0, 140.80528561, 0, 0, 0, 0, 0, 0, 0, 0, 0, 199.48144377, 0, 0, 0, 0, 199.48144377, 0, 0, 0, 0, -38.99126336, 0, 0, 0, 0, -38.99126336, 0, 0, 0, 0, -94.48017644, 0, 0, 0, 0, -94.48017644), 
                        nrow = 10, ncol = 8, dimnames = list(c("(Intercept)", "visitY2", "visitY3", "visitY4", "age", "genderfemale", "visitY2:genderfemale", "visitY3:genderfemale", "visitY4:genderfemale", "age:genderfemale"),c("sigma:male", "sigma:female", "k.Y2:male", "k.Y3:male", "k.Y4:male", "k.Y2:female", "k.Y3:female", "k.Y4:female"))) 
    expect_equal(as.double(test[name.mu,name.vcov]),as.double(GS.muvcov), tol = 1e-5)

    GS.vcovF <- matrix(c(371.47143055, 87.01549659, 89.82765268, 80.01553877, 87.01549659, 89.83790539, -0.89130901, -6.1318067, 89.82765268, -0.89130901, 88.63803396,
                         -0.09068914, 80.01553877, -6.1318067, -0.09068914, 74.93849286), 
                       nrow = 4, ncol = 4, dimnames = list(c("sigma:female", "k.Y2:female", "k.Y3:female", "k.Y4:female"),c("sigma:female", "k.Y2:female", "k.Y3:female", "k.Y4:female")) ) 
    expect_equal(as.double(test[name.vcovF,name.vcovF]),as.double(GS.vcovF), tol = 1e-5)

    GS.vcovM <- matrix(c(498.91800085, 115.22711759, 112.45930461, 105.6207779, 115.22711759, 112.7824544, -1.31024125, -2.24418408, 112.45930461, -1.31024125, 102.94836181,
                         -0.06439274, 105.6207779, -2.24418408, -0.06439274, 91.58135279),
                       nrow = 4, ncol = 4, dimnames = list(c("sigma:male", "k.Y2:male", "k.Y3:male", "k.Y4:male"),c("sigma:male", "k.Y2:male", "k.Y3:male", "k.Y4:male")))
    expect_equal(as.double(test[name.vcovM,name.vcovM]),as.double(GS.vcovM), tol = 1e-5)

    expect_true(all(abs(as.double(test[name.vcovF,name.vcovM])<1e-10)))

    test_logSigma <- information(eSUN.lmm, effects = "all",
                        p = coef(eSUN.lmm, effects = "all", transform.sigma = "none", transform.k = "none", transform.rho = "none"),
                        transform.sigma = "log", transform.k = "none", transform.rho = "none")
    ## GS <- -hessian(func = function(p){p[eSUN.lmm$design$param$type %in% "sigma"]<-exp(p[eSUN.lmm$design$param$type %in% "sigma"]);logLik(eSUN.lmm, p = p)},
    ##                x = coef(eSUN.lmm, effects = "all", transform.sigma = "log", transform.k = "none", transform.rho = "none", transform.names = FALSE))

    GS.vcovF <- matrix(c(326.00000202, 81.51596908, 84.15039212, 74.95842051, 81.51596908, 89.83790539, -0.89130901, -6.1318067, 84.15039212, -0.89130901, 88.63803396, -0.09068914, 74.95842051, -6.1318067, -0.09068914, 74.93849286), 
                       nrow = 4, ncol = 4, dimnames = list(c("log(sigma):female", "k.Y2:female", "k.Y3:female", "k.Y4:female"),c("log(sigma):female", "k.Y2:female", "k.Y3:female", "k.Y4:female")) )
    expect_equal(as.double(test_logSigma[name.all %in% name.vcovF,name.all %in% name.vcovF]),as.double(GS.vcovF), tol = 1e-5)

    GS.vcovM <- matrix(c(454.00000041, 109.91780271, 107.27752213, 100.75409391, 109.91780271, 112.7824544, -1.31024125, -2.24418408, 107.27752213, -1.31024125, 102.94836181, -0.06439274, 100.75409391, -2.24418408, -0.06439274, 91.58135279), nrow = 4, ncol = 4, dimnames = list(c("log(sigma):male", "k.Y2:male", "k.Y3:male", "k.Y4:male"),c("log(sigma):male", "k.Y2:male", "k.Y3:male", "k.Y4:male")) )
    expect_equal(as.double(test_logSigma[name.all %in% name.vcovM,name.all %in% name.vcovM]),as.double(GS.vcovM), tol = 1e-5)

    expect_true(all(abs(as.double(test_logSigma[name.all %in% name.vcovF,name.all %in% name.vcovM])<1e-10)))

    test_lava <- information(eSUN.lmm, effects = "all",
                             p = coef(eSUN.lmm, effects = "all", transform.sigma = "none", transform.k = "none", transform.rho = "none"),
                             transform.k = "var", transform.rho = "none", transform.names = FALSE) 
    ## GS <- -hessian(func = function(p){logLik(eSUN.lmm, p = fct.trans(p,inverse = TRUE))},
    ##                x = coef(eSUN.lmm, effects = "all", transform.k = "var", transform.rho = "none", transform.names = FALSE))
 
    GS.vcovF <- matrix(c(27.03458543, -0.2372027, -0.16490875, -0.14958512, -0.2372027, 29.19263486, -0.29964208, -1.8271702, -0.16490875, -0.29964208, 30.82870221, -0.02795802, -0.14958512, -1.8271702, -0.02795802, 20.4772797),nrow = 4,ncol = 4,dimnames = list(c("sigma:female", "k.Y2:female", "k.Y3:female", "k.Y4:female"),c("sigma:female", "k.Y2:female", "k.Y3:female", "k.Y4:female")) ) 
    expect_equal(as.double(test_lava[name.all %in% name.vcovF,name.all %in% name.vcovF]),as.double(GS.vcovF), tol = 1e-5)

    GS.vcovM <- matrix(c(35.16315581, -0.73325294, -0.07775833, -0.02930618, -0.73325294, 32.00020684, -0.36228096, -0.58262412, -0.07775833, -0.36228096, 27.73937621, -0.01629108, -0.02930618, -0.58262412, -0.01629108, 21.75483409),nrow = 4,ncol = 4,dimnames = list(c("sigma:male", "k.Y2:male", "k.Y3:male", "k.Y4:male"),c("sigma:male", "k.Y2:male", "k.Y3:male", "k.Y4:male")) ) 
    expect_equal(as.double(test_lava[name.all %in% name.vcovM,name.all %in% name.vcovM]),as.double(GS.vcovM), tol = 1e-5)

    expect_true(all(abs(as.double(test_lava[name.all %in% name.vcovF,name.all %in% name.vcovM])<1e-10)))

    ## test <- information(eSUN.lmm, p = coef(eSUN.lmm, transform.sigma = "none"), transform.rho = "cov") 
    ## GS <- -hessian(func = function(p){p[c("sigma","rho")] <- c(sqrt(p["sigma"]),p["rho"]/p["sigma"]); logLik(eSUN.lmm, p = p)}, x = coef(eSUN.lmm, transform.rho = "cov", transform.names = FALSE))
    ## expect_equal(as.double(test), as.double(GS), tol = 1e-6)

    ## no transformation 
    newp <- coef(eSUN.lmm, effects = "all", transform.sigma = "none", transform.k = "none", transform.rho = "none")+0.1
    ## GS <- -hessian(func = function(p){logLik(eSUN.lmm, p = p)}, x = newp)

    test2 <- information(eSUN.lmm, effects = "all", p = newp, transform.sigma = "none", transform.k = "none", transform.rho = "none")

    GS.mu <- matrix(c(254.24318203, 71.76445765, 63.87468527, 42.26944907, 12924.03081837, 87.11643724, 21.99098334, 29.42732688, 10.75450474, 4326.26096328, 71.76445765, 84.43958159, 4.6800939, -14.44164788, 3652.66391742, 21.99098334, 42.93868075, 3.64466597, -17.80841204, 1092.08704759, 63.87468527, 4.6800939, 75.46467934, -7.23359469, 3233.51172331, 29.42732688, 3.64466597, 35.92910321, -3.23861115, 1461.38087768, 42.26944907, -14.44164788, -7.23359469, 72.24811721, 2155.3499046, 10.75450474, -17.80841204, -3.23861115, 36.57718605, 534.07595056, 12924.03081837, 3652.66391742, 3233.51172331, 2155.3499046, 688101.15548784, 4326.26096328, 1092.08704759, 1461.38087768, 534.07595056, 222958.48295566, 87.11643724, 21.99098334, 29.42732688, 10.75450474, 4326.26096328, 87.11643724, 21.99098334, 29.42732688, 10.75450474, 4326.26096328, 21.99098334, 42.93868075, 3.64466597, -17.80841204, 1092.08704759, 21.99098334, 42.93868075, 3.64466597, -17.80841204, 1092.08704759, 29.42732688, 3.64466597, 35.92910321, -3.23861115, 1461.38087768, 29.42732688, 3.64466597, 35.92910321, -3.23861115, 1461.38087768, 10.75450474, -17.80841204, -3.23861115, 36.57718605, 534.07595056, 10.75450474, -17.80841204, -3.23861115, 36.57718605, 534.07595056, 4326.26096328, 1092.08704759, 1461.38087768, 534.07595056, 222958.48295566, 4326.26096328, 1092.08704759, 1461.38087768, 534.07595056, 222958.48295566), 
                    nrow = 10, ncol = 10, dimnames = list(c("(Intercept)", "visitY2", "visitY3", "visitY4", "age", "genderfemale", "visitY2:genderfemale", "visitY3:genderfemale", "visitY4:genderfemale", "age:genderfemale"),c("(Intercept)", "visitY2", "visitY3", "visitY4", "age", "genderfemale", "visitY2:genderfemale", "visitY3:genderfemale", "visitY4:genderfemale", "age:genderfemale")) ) 
    expect_equal(as.double(test2[name.mu,name.mu]),as.double(GS.mu), tol = 1e-5)

    GS.muvcov <- matrix(c(-1685.25252986, -504.06971427, -349.77017266, -320.29503872, -91043.43888514, 0, 0, 0, 0, 0, -1726.68115881, -440.9162199, -589.17664735, -216.18886504, -88881.62971225, -1726.68115881, -440.9162199, -589.17664735, -216.18886504, -88881.62971225, -469.83394988, -430.76127601, -4.89049, -15.90175999, -25401.77112182, 0, 0, 0, 0, 0, -318.00367211, -4.77655838, -341.4762561, 18.42935809, -17182.13567777, 0, 0, 0, 0, 0, -274.74336896, -14.66016298, 17.39566748, -292.8401056, -14759.89276707, 0, 0, 0, 0, 0, -414.54823988, -611.40121182, -34.25022396, 167.35198939, -21203.33350299, -414.54823988, -611.40121182, -34.25022396, 167.35198939, -21203.33350299, -571.70133223, -35.32324114, -634.71501176, 31.38785378, -29437.85032793, -571.70133223, -35.32324114, -634.71501176, 31.38785378, -29437.85032793, -187.56664228, 154.63145744, 28.1210453, -411.78645691, -9715.60136534, -187.56664228, 154.63145744, 28.1210453, -411.78645691, -9715.60136534), 
                        nrow = 10, ncol = 8, dimnames = list(c("(Intercept)", "visitY2", "visitY3", "visitY4", "age", "genderfemale", "visitY2:genderfemale", "visitY3:genderfemale", "visitY4:genderfemale", "age:genderfemale"),c("sigma:male", "sigma:female", "k.Y2:male", "k.Y3:male", "k.Y4:male", "k.Y2:female", "k.Y3:female", "k.Y4:female"))) 
    expect_equal(as.double(test2[name.mu,name.vcov]),as.double(GS.muvcov), tol = 1e-5)

    GS.vcovF <- matrix(c(26769.98049729, 4288.95680591, 5965.86978904, 2016.17040192, 4288.95680591, 7939.27669307, 342.27188376, -1506.40817923, 5965.86978904, 342.27188376, 9291.01643238, -282.90000023, 2016.17040192, -1506.40817923, -282.90000023, 4618.83374848), 
                       nrow = 4, ncol = 4, dimnames = list(c("sigma:female", "k.Y2:female", "k.Y3:female", "k.Y4:female"),c("sigma:female", "k.Y2:female", "k.Y3:female", "k.Y4:female")) ) 
    expect_equal(as.double(test2[name.vcovF,name.vcovF]),as.double(GS.vcovF), tol = 1e-5)

    GS.vcovM <- matrix(c(13604.49242326, 2563.84567473, 1757.06336823, 1508.27378162, 2563.84567473, 3347.35271581, 23.47824893, 71.74218771, 1757.06336823, 23.47824893, 2467.93754138, -83.83308856, 1508.27378162, 71.74218771, -83.83308856, 1986.67629708),
                       nrow = 4, ncol = 4, dimnames = list(c("sigma:male", "k.Y2:male", "k.Y3:male", "k.Y4:male"),c("sigma:male", "k.Y2:male", "k.Y3:male", "k.Y4:male")))
    expect_equal(as.double(test2[name.vcovM,name.vcovM]),as.double(GS.vcovM), tol = 1e-5)

    expect_true(all(abs(as.double(test2[name.vcovF,name.vcovM])<1e-10)))
   
    ## ** degree of freedom
    test <- confint(eSUN.lmm, effects = "all")
    GS <- c(75.02445251, 57.00262333, 57.00100521, 57.00082937, 55.02679116, 67.96985133, 96.26931669, 94.16508864, 93.47141533, 54.92384753, 56.85151867, 41.11400945, 109.52022023, 113.40272184, 114.03144185, 81.54478741, 81.63717886, 79.9545627, 29.54320962, 28.55308578, 28.49261999, 29.06453625, 29.60876751, 28.48658921, 20.70521888, 20.51343872, 20.7126253, 20.32470404, 23.65317661, 20.31145426)
    expect_equal(test$df, GS, tol = 1e-3)
    
    ## ** anova
    eSUN.lmm_anova <- anova(eSUN.lmm, effects = "all", ci = TRUE)$multivariate
    expect_equal(eSUN.lmm_anova[eSUN.lmm_anova$type=="mu","df.denom"], c(57.00206139, 55.02679117, 67.96985137, 93.51660476, 54.92384757), tol = 1e-1)
    expect_equal(eSUN.lmm_anova[eSUN.lmm_anova$type=="k","df.denom"], c(86.07826), tol = 1e-1)
    expect_equal(eSUN.lmm_anova[eSUN.lmm_anova$type=="rho","df.denom"], c(9.100919), tol = 1e-1)


    ## ** getVarCov
    Omega.GS <- list("female" = matrix(c(0.87759105, 0.1151575, 0.06705109, 0.12241344, 0.1151575, 0.87666546, -0.09762837, 0.35757795, 0.06705109, -0.09762837, 0.81905439, -0.00726966, 0.12241344, 0.35757795, -0.00726966, 1.04251101), 
                                    nrow = 4, 
                                    ncol = 4, 
                                    dimnames = list(c("Y1", "Y2", "Y3", "Y4"),c("Y1", "Y2", "Y3", "Y4"))),
                     "male" = matrix(c(0.90996917, -0.18047956, -0.03950488, -0.00529246, -0.18047956, 0.9682825, -0.1395401, -0.21029739, -0.03950488, -0.1395401, 1.01961468, 0.01518698, -0.00529246, -0.21029739, 0.01518698, 1.15655052), 
                                    nrow = 4, 
                                    ncol = 4, 
                                    dimnames = list(c("Y1", "Y2", "Y3", "Y4"),c("Y1", "Y2", "Y3", "Y4"))) 
                     )

    expect_equivalent(sigma(eSUN.lmm)[names(Omega.GS)], Omega.GS, tol = 1e-5)
})

## * Missing data
test_that("missing values",{

    ## ** full cluster missing
    dL$Ymiss <- dL$Y
    dL$Ymiss[dL$id==1] <- NA

    eCS.lmm <- suppressWarnings(lmm(Ymiss ~ visit + age + gender, repetition = ~visit|id, structure = "CS", data = dL, trace = 0, method.fit = "ML"))
    eCS2.lmm <- lmm(Ymiss ~ visit + age + gender, repetition = ~visit|id, structure = "CS", data = dL[dL$id!=1,], trace = 0, method.fit = "ML")
    expect_equal(logLik(eCS2.lmm), logLik(eCS.lmm))
    
    ## logLik(eCS.lmm, indiv = TRUE)
    ## score(eCS.lmm, indiv = TRUE)
    ## information(eCS.lmm, indiv = TRUE)
    
    ## ** only part of the cluster is missing
    set.seed(11)
    dL$Ymiss <- dL$Y
    dL$Ymiss[which(rbinom(NROW(dL), size = 1, prob = 0.1)==1)] <- NA

    eCS.lmm <- lmm(Ymiss ~ visit + age + gender, repetition = ~visit|id, structure = "CS", data = dL, trace = 0, method.fit = "ML")
    eCS.gls <- gls(Ymiss ~ visit + age + gender, correlation = corCompSymm(form=~1|id), data = dL, method = "ML", na.action = na.omit)
    expect_equal(as.double(logLik(eCS.lmm)), as.double(logLik(eCS.gls)))
})

## * Baseline constraint
test_that("Baseline constraint",{

    dL$group <- as.factor(dL$id %% 2)
    dL$treat <- (dL$group==1)*(dL$visit!="Y1")
    table(dL$treat, baselineAdjustment(dL, variable = "group", repetition = ~visit|id, constrain = "Y1", new.level = "0"))
    dL$treat.visit <- baselineAdjustment(dL, variable = "group", repetition = ~visit|id, constrain = "Y1", collapse.time = ".")

    ## eUN.lmm <- lmm(Y ~ group*visit, repetition = ~group*visit|id, structure = "UN", data = dL, trace = 0, method = "REML")
    ## logLik(eUN.lmm)
    eCUN.lmm <- suppressMessages(lmm(Y ~ treat*visit, repetition = ~treat*visit|id, structure = "UN", data = dL, trace = 0, method = "REML", df = FALSE, control = list(optimizer = "FS")))
    eCUN2.lmm <- lmm(Y ~ treat.visit, repetition = ~treat.visit|id, structure = "UN", data = dL, trace = 0, method = "REML", df = FALSE, control = list(optimizer = "FS"))
    
    expect_equal(logLik(eCUN2.lmm), logLik(eCUN.lmm), tol = 1e-5)
    expect_equal(logLik(eCUN2.lmm), -618.14359397, tol = 1e-5)

    capture.output(summary(eCUN2.lmm))
    capture.output(summary(anova(eCUN2.lmm), method = "none"))
    ## plot(eCUN2.lmm, color = "group", time.var = "visit")

    ## baseline constrain for order 3 interaction
    eCUN.I2.lmm <- suppressMessages(lmm(Y ~ gender*treat*visit, repetition = ~treat*visit|id, structure = "UN", data = dL, trace = 0, method = "REML", df = FALSE, control = list(optimizer = "FS")))
    eCUN2.I2.lmm <- suppressMessages(lmm(Y ~ gender:treat.visit, repetition = ~treat.visit|id, structure = "UN", data = dL, trace = 0, method = "REML", df = FALSE, control = list(optimizer = "FS")))
    expect_equal(logLik(eCUN.I2.lmm), logLik(eCUN2.I2.lmm), tol = 1e-5)
    expect_equal(logLik(eCUN.I2.lmm), -598.96051963, tol = 1e-5)

})


##----------------------------------------------------------------------
### test-auto-pattern-regression.R ends here
bozenne/repeated documentation built on July 16, 2025, 11:16 p.m.
rdrr.io home R language documentation Run R code online
CRAN packages Bioconductor packages R-Forge packages GitHub packages
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
bozenne/repeated
Repeated Measurement Models for Discrete Times

tests/testthat/test-auto-pattern-regression.R
In bozenne/repeated: Repeated Measurement Models for Discrete Times

R Package Documentation

Browse R Packages

We want your feedback!

bozenne/repeated Repeated Measurement Models for Discrete Times

tests/testthat/test-auto-pattern-regression.R In bozenne/repeated: Repeated Measurement Models for Discrete Times

R Package Documentation

Browse R Packages

We want your feedback!

bozenne/repeated
Repeated Measurement Models for Discrete Times

tests/testthat/test-auto-pattern-regression.R
In bozenne/repeated: Repeated Measurement Models for Discrete Times
