### test-auto-pattern-regression.R --- 
##----------------------------------------------------------------------
## Author: Brice Ozenne
## Created: jul 31 2023 (09:54) 
## Version: 
## Last-Updated: jul 31 2023 (10:29) 
##           By: Brice Ozenne
##     Update #: 10
##----------------------------------------------------------------------
## 
### Commentary: 
## 
### Change Log:
##----------------------------------------------------------------------
## 
### Code:

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

    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.sigma = "none", transform.k = "none", 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.sigma = "none", transform.k = "none", 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)
    index <- sample.int(NROW(dL))
    test2 <- predict(eCS.lmm, newdata = dL[index,,drop=FALSE])
    if(require(AICcmodavg)){
        GS <- AICcmodavg::predictSE(eCS.gls, newdata = dL)
        expect_equivalent(test$estimate,GS$fit, tol = 1e-7)
        expect_equivalent(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("2:2" = 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")) 
                                  ),
                     "1:1" = 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), 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 <- rbind(c(-1783.17328003, -494.19593699, -489.74409476, -280.85528479, -94052.55029524, -895.11004272, -228.57059199, -305.42867251, -112.07212518, -46076.16119356, 4635.77219006, 9146.88022716, 1300.93756054, 876.88257529, 748.51536418, 2186.33054224, 3054.55980008, 1010.82864048, 1648.82180282, 2154.97282945, 892.39161399, 2147.65971004, 892.35942784, 1162.00385316, 1625.93441573, 1155.44499505, 1113.83275099, 1285.80100671, 1238.78723302, 881.61137276))
    test <- score(eSUN.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[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 <- rbind(c(-1783.17328003, -494.19593699, -489.74409476, -280.85528479, -94052.55029524, -895.11004272, -228.57059199, -305.42867251, -112.07212518, -46076.16119356, 4885.74713667, 9483.4698487, 1300.93756054, 876.88257529, 748.51536418, 2186.33054224, 3054.55980008, 1010.82864048, 1648.82180282, 2154.97282945, 892.39161399, 2147.65971004, 892.35942784, 1162.00385316, 1625.93441573, 1155.44499505, 1113.83275099, 1285.80100671, 1238.78723302, 881.61137276))    
    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 <- rbind(c(-1783.17328003, -494.19593699, -489.74409476, -280.85528479, -94052.55029524, -895.11004272, -228.57059199, -305.42867251, -112.07212518, -46076.16119356, 665.79608401, 1218.87880924, 517.53206247, 340.71006511, 274.52072123, 924.9385869, 1332.73203954, 395.13143762, 1648.82180282, 2154.97282945, 892.39161399, 2147.65971004, 892.35942784, 1162.00385316, 1625.93441573, 1155.44499505, 1113.83275099, 1285.80100671, 1238.78723302, 881.61137276)
                )
    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
    ## 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 <- cbind(c(481.58902557, 141.4145761, 126.4937567, 90.87132761, 24505.0293416, 151.42333095, 39.89451566, 53.34143661, 22.82307698, 7519.7846284, 0, 0, -1e-08, 1e-08, 0, 1e-08, 1e-08, 2e-08, 4e-08, 9e-08, 9e-08, 1.7e-07, 3e-08, 4e-07, 2e-08, 2.6e-07, 8.2e-07, 9e-08, 5e-08, 1.6e-07), 
                c(141.4145761, 123.53334857, 16.52169525, -6.89363657, 7203.84822301, 39.89451564, 57.11914728, 7.08254392, -18.90835878, 1981.18852358, 0, 0, 0, 0, 1e-08, 0, 0, 0, -1e-08, 7e-08, -3e-08, 1e-08, 1e-08, -1.8e-07, 2e-08, 1.9e-07, 2.3e-07, 3e-08, 0, 4.2e-07), 
                c(126.4937567, 16.52169525, 110.82914443, -0.53727867, 6412.26409583, 53.34143653, 7.08254388, 52.49669251, -1.50781625, 2648.97167493, 1e-08, -1e-08, -2e-08, 1e-08, 1e-08, 0, 0, 0, -6e-08, 8e-08, 3e-08, -2e-08, 0, -4e-08, -5e-08, 2.8e-07, 0, -3e-08, -5e-08, 1.2e-07), 
                c(90.87132761, -6.89363657, -0.53727867, 99.56459117, 4634.12494773, 22.82307691, -18.90835881, -1.50781628, 47.23106187, 1133.40937526, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2e-08, 0, 1e-08, 0, 9e-08, 1e-08, 3e-08, -7e-08, 1e-08, 0, 8e-08), 
                c(24505.0293416, 7203.84822301, 6412.26409583, 4634.12494773, 1306448.44293833, 7519.78462743, 1981.18852307, 2648.97167428, 1133.40937504, 387540.13845235, -4.9e-07, -2.2e-07, -55.6535616, -16.1141946, 140.80528556, 199.48144369, -38.99126342, -94.48017649, 62.7205861, -45.39488044, -82.81165625, 126.92432422, -7.39505442, -98.5050688, -86.80408649, -60.47376097, 56.0232937, -39.89494287, 99.31782595, 76.6938026), 
                c(151.42333095, 39.89451564, 53.34143653, 22.82307691, 7519.78462743, 151.42333094, 39.89451563, 53.34143654, 22.8230769, 7519.7846276, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1e-08, 0, -1e-08, 0, -1e-07, 0, 0, -2e-08, 0, 0, -1e-08), 
                c(39.89451566, 57.11914728, 7.08254388, -18.90835881, 1981.18852307, 39.89451563, 57.11914725, 7.08254391, -18.90835882, 1981.18852349, 0, 0, 0, 0, 0, 0, 0, 0, -1e-08, 2e-08, -1e-08, 1e-08, 0, -1.2e-07, 0, 0, 0, 0, 0, -1e-08), 
                c(53.34143661, 7.08254392, 52.49669251, -1.50781628, 2648.97167428, 53.34143654, 7.08254391, 52.49669251, -1.50781629, 2648.97167472, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2e-08, 0, 0, 0, -1.4e-07, 0, 5e-08, 0, 0, 0, 0), 
                c(22.82307698, -18.90835878, -1.50781625, 47.23106187, 1133.40937504, 22.8230769, -18.90835882, -1.50781629, 47.23106187, 1133.40937519, 0, 0, 0, 0, 0, 0, 0, 0, -2e-08, 2e-08, -4e-08, -2e-08, -1e-08, -8e-08, 0, 0, 0, 0, 0, 0), 
                c(7519.7846284, 1981.18852358, 2648.97167493, 1133.40937526, 387540.13845235, 7519.7846276, 1981.18852349, 2648.97167472, 1133.40937519, 387540.13844065, 5e-08, -3e-08, 0, 8e-08, 8e-08, 199.48144381, -38.99126326, -94.48017637, 62.72058643, -45.39487981, -82.811656, 126.92432487, -7.39505416, -98.50506137, 4.5e-07, 2.12e-06, 1.681e-05, 0, 2.2e-07, 6.2e-06), 
                c(0, 0, 1e-08, 0, -4.9e-07, 0, 0, 0, 0, 5e-08, 498.91800084, 0, 115.22711757, 112.45930459, 105.62077788, 0, 0, 0, -1e-08, -1e-08, 0, -1e-08, 0, -1.3e-07, 25.49189912, 8.07827461, 5.10112289, 18.51692731, 25.39714368, 1.57974633), 
                c(0, 0, -1e-08, 0, -2.2e-07, 0, 0, 0, 0, -3e-08, 0, 371.47143055, 0, 0, 0, 87.01549657, 89.82765267, 80.01553875, -10.28461852, -8.89442284, -8.21656242, 11.9023847, -37.96514389, -3.28134592, -1e-08, -1e-08, -2.1e-07, -1e-08, -1e-08, -8e-08), 
                c(-1e-08, 0, -2e-08, 0, -55.6535616, 0, 0, 0, 0, 0, 115.22711757, 0, 112.7824544, -1.31024125, -2.24418408, 0, 0, 0, 0, 0, -1e-08, -1e-08, 0, -1.2e-07, 11.94927774, -0.19182148, -0.19004235, 8.70562069, 11.88553603, -0.16828589), 
                c(1e-08, 0, 1e-08, 0, -16.1141946, 0, 0, 0, 0, 8e-08, 112.45930459, 0, -1.31024125, 102.94836181, -0.06439274, 0, 0, 0, -1e-08, 0, 0, -1e-08, 0, -1.7e-07, -0.18216424, 3.80923856, -0.13344723, 8.54620445, -0.1598135, 0.86031718), 
                c(0, 1e-08, 1e-08, 0, 140.80528556, 0, 0, 0, 0, 8e-08, 105.62077788, 0, -2.24418408, -0.06439274, 91.58135278, 0, 0, 0, 0, 0, 0, 0, 0, 1e-08, -0.15910624, -0.11764696, 2.32434774, -0.14089147, 10.94386634, 0.81552409), 
                c(1e-08, 0, 0, 0, 199.48144369, 0, 0, 0, 0, 199.48144381, 0, 87.01549657, 0, 0, 0, 89.83790539, -0.89130901, -6.1318067, -4.72773497, -0.1325759, -0.06399443, 5.71221553, -17.72741051, -0.09324222, 0, -1e-08, 0, -1e-08, 0, -4e-08), 
                c(1e-08, 0, 0, 0, -38.99126342, 0, 0, 0, 0, -38.99126326, 0, 89.82765267, 0, 0, 0, -0.89130901, 88.63803396, -0.09068914, -0.14189633, -4.23553146, -0.08853111, 5.85757284, -0.09980804, -1.53685315, 0, 0, -8e-08, 0, 0, 1e-08), 
                c(2e-08, 0, 0, 0, -94.48017649, 0, 0, 0, 0, -94.48017637, 0, 80.01553875, 0, 0, 0, -6.1318067, -0.09068914, 74.93849286, -0.05382399, -0.06955407, -3.45331919, -0.07842752, -16.22805835, -1.29681765, 0, 0, -9e-08, 0, 0, 0), 
                c(4e-08, -1e-08, -6e-08, 0, 62.7205861, 0, -1e-08, 0, -2e-08, 62.72058643, -1e-08, -10.28461852, 0, -1e-08, 0, -4.72773497, -0.14189633, -0.05382399, 50.9890118, 6.40124366, -17.87273437, -5.43814218, -2.50747803, 1.08339806, -4e-08, 0, -1.49e-06, -3e-08, -4e-08, 0), 
                c(9e-08, 7e-08, 8e-08, 2e-08, -45.39488044, 1e-08, 2e-08, 2e-08, 2e-08, -45.39487981, -1e-08, -8.89442284, 0, 0, 0, -0.1325759, -4.23553146, -0.06955407, 6.40124366, 43.64783296, -1.11149819, -5.43451324, -0.44808313, -3.77487948, 0, 1.8e-07, 0, -4e-08, 0, 4.4e-07), 
                c(9e-08, -3e-08, 3e-08, 0, -82.81165625, 0, -1e-08, 0, -4e-08, -82.811656, 0, -8.21656242, -1e-08, 0, 0, -0.06399443, -0.08853111, -3.45331919, -17.87273437, -1.11149819, 50.00966136, 1.76934478, -3.8485412, -4.50868395, -2e-08, -9e-08, -8e-07, -6e-08, -2e-08, -2.9e-07), 
                c(1.7e-07, 1e-08, -2e-08, 1e-08, 126.92432422, -1e-08, 1e-08, 0, -2e-08, 126.92432487, -1e-08, 11.9023847, -1e-08, -1e-08, 0, 5.71221553, 5.85757284, -0.07842752, -5.43814218, -5.43451324, 1.76934478, 50.71957223, -4.22066211, -18.35978405, -2e-08, -1e-07, -8.9e-07, -6e-08, -2e-08, -3.5e-07), 
                c(3e-08, 1e-08, 0, 0, -7.39505442, 0, 0, 0, -1e-08, -7.39505416, 0, -37.96514389, 0, 0, 0, -17.72741051, -0.09980804, -16.22805835, -2.50747803, -0.44808313, -3.8485412, -4.22066211, 64.86434918, 7.40531477, -1e-08, -3e-08, -5e-07, -2e-08, -1e-08, -9e-08), 
                c(4e-07, -1.8e-07, -4e-08, 9e-08, -98.5050688, -1e-07, -1.2e-07, -1.4e-07, -8e-08, -98.50506137, -1.3e-07, -3.28134592, -1.2e-07, -1.7e-07, 1e-08, -0.09324222, -1.53685315, -1.29681765, 1.08339806, -3.77487948, -4.50868395, -18.35978405, 7.40531477, 49.19089192, -6.4e-07, -2.99e-06, -2.275e-05, -8.4e-07, -2.3e-07, -1.355e-05), 
                c(2e-08, 2e-08, -5e-08, 1e-08, -86.80408649, 0, 0, 0, 0, 4.5e-07, 25.49189912, -1e-08, 11.94927774, -0.18216424, -0.15910624, 0, 0, 0, -4e-08, 0, -2e-08, -2e-08, -1e-08, -6.4e-07, 68.49310859, 10.2811001, 13.39724145, 6.3299363, 5.47958227, 1.11448435), 
                c(2.6e-07, 1.9e-07, 2.8e-07, 3e-08, -60.47376097, 0, 0, 5e-08, 0, 2.12e-06, 8.07827461, -1e-08, -0.19182148, 3.80923856, -0.11764696, -1e-08, 0, 0, 0, 1.8e-07, -9e-08, -1e-07, -3e-08, -2.99e-06, 10.2811001, 61.22307964, 1.0795372, 13.34154613, 0.42815576, 2.71608178), 
                c(8.2e-07, 2.3e-07, 0, -7e-08, 56.0232937, -2e-08, 0, 0, 0, 1.681e-05, 5.10112289, -2.1e-07, -0.19004235, -0.13344723, 2.32434774, 0, -8e-08, -9e-08, -1.49e-06, 0, -8e-07, -8.9e-07, -5e-07, -2.275e-05, 13.39724145, 1.0795372, 62.13284815, 0.89834529, 13.50221379, 4.22666712), 
                c(9e-08, 3e-08, -3e-08, 1e-08, -39.89494287, 0, 0, 0, 0, 0, 18.51692731, -1e-08, 8.70562069, 8.54620445, -0.14089147, -1e-08, 0, 0, -3e-08, -4e-08, -6e-08, -6e-08, -2e-08, -8.4e-07, 6.3299363, 13.34154613, 0.89834529, 66.00144213, 2.87473442, 12.76327821), 
                c(5e-08, 0, -5e-08, 0, 99.31782595, 0, 0, 0, 0, 2.2e-07, 25.39714368, -1e-08, 11.88553603, -0.1598135, 10.94386634, 0, 0, 0, -4e-08, 0, -2e-08, -2e-08, -1e-08, -2.3e-07, 5.47958227, 0.42815576, 13.50221379, 2.87473442, 68.39943276, 9.63020519), 
                c(1.6e-07, 4.2e-07, 1.2e-07, 8e-08, 76.6938026, -1e-08, -1e-08, 0, 0, 6.2e-06, 1.57974633, -8e-08, -0.16828589, 0.86031718, 0.81552409, -4e-08, 1e-08, 0, 0, 4.4e-07, -2.9e-07, -3.5e-07, -9e-08, -1.355e-05, 1.11448435, 2.71608178, 4.22666712, 12.76327821, 9.63020519, 60.85664305)
                )
    expect_equal(as.double(test),as.double(GS), tol = 1e-5)

    test <- 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 <- cbind(c(481.58902557, 141.4145761, 126.4937567, 90.87132761, 24505.0293416, 151.42333095, 39.89451566, 53.34143661, 22.82307698, 7519.7846284, 2e-07, 1e-07, -1e-08, 1e-08, 0, 1e-08, 1e-08, 2e-08, 4e-08, 9e-08, 9e-08, 1.7e-07, 3e-08, 4e-07, 2e-08, 2.6e-07, 8.2e-07, 9e-08, 5e-08, 1.6e-07), 
                c(141.4145761, 123.53334857, 16.52169525, -6.89363657, 7203.84822301, 39.89451564, 57.11914728, 7.08254392, -18.90835878, 1981.18852358, 1e-07, 6e-08, 0, 0, 1e-08, 0, 0, 0, -1e-08, 7e-08, -3e-08, 1e-08, 1e-08, -1.8e-07, 2e-08, 1.9e-07, 2.3e-07, 3e-08, 0, 4.2e-07), 
                c(126.4937567, 16.52169525, 110.82914443, -0.53727867, 6412.26409583, 53.34143653, 7.08254388, 52.49669251, -1.50781625, 2648.97167493, 1.2e-07, 9e-08, -2e-08, 1e-08, 1e-08, 0, 0, 0, -6e-08, 8e-08, 3e-08, -2e-08, 0, -4e-08, -5e-08, 2.8e-07, 0, -3e-08, -5e-08, 1.2e-07), 
                c(90.87132761, -6.89363657, -0.53727867, 99.56459117, 4634.12494773, 22.82307691, -18.90835881, -1.50781628, 47.23106187, 1133.40937526, 4e-08, 2e-08, 0, 0, 0, 0, 0, 0, 0, 2e-08, 0, 1e-08, 0, 9e-08, 1e-08, 3e-08, -7e-08, 1e-08, 0, 8e-08), 
                c(24505.0293416, 7203.84822301, 6412.26409583, 4634.12494773, 1306448.44293833, 7519.78462743, 1981.18852307, 2648.97167428, 1133.40937504, 387540.13845235, 3e-08, -2e-07, -55.6535616, -16.1141946, 140.80528556, 199.48144369, -38.99126342, -94.48017649, 62.7205861, -45.39488044, -82.81165625, 126.92432422, -7.39505442, -98.5050688, -86.80408649, -60.47376097, 56.0232937, -39.89494287, 99.31782595, 76.6938026), 
                c(151.42333095, 39.89451564, 53.34143653, 22.82307691, 7519.78462743, 151.42333094, 39.89451563, 53.34143654, 22.8230769, 7519.7846276, 0, 1e-08, 0, 0, 0, 0, 0, 0, 0, 1e-08, 0, -1e-08, 0, -1e-07, 0, 0, -2e-08, 0, 0, -1e-08), 
                c(39.89451566, 57.11914728, 7.08254388, -18.90835881, 1981.18852307, 39.89451563, 57.11914725, 7.08254391, -18.90835882, 1981.18852349, 0, -3e-08, 0, 0, 0, 0, 0, 0, -1e-08, 2e-08, -1e-08, 1e-08, 0, -1.2e-07, 0, 0, 0, 0, 0, -1e-08), 
                c(53.34143661, 7.08254392, 52.49669251, -1.50781628, 2648.97167428, 53.34143654, 7.08254391, 52.49669251, -1.50781629, 2648.97167472, 4e-08, 1e-08, 0, 0, 0, 0, 0, 0, 0, 2e-08, 0, 0, 0, -1.4e-07, 0, 5e-08, 0, 0, 0, 0), 
                c(22.82307698, -18.90835878, -1.50781625, 47.23106187, 1133.40937504, 22.8230769, -18.90835882, -1.50781629, 47.23106187, 1133.40937519, 1e-08, -5e-08, 0, 0, 0, 0, 0, 0, -2e-08, 2e-08, -4e-08, -2e-08, -1e-08, -8e-08, 0, 0, 0, 0, 0, 0), 
                c(7519.7846284, 1981.18852358, 2648.97167493, 1133.40937526, 387540.13845235, 7519.7846276, 1981.18852349, 2648.97167472, 1133.40937519, 387540.13844065, 1.84e-06, 1.13e-06, 0, 8e-08, 8e-08, 199.48144381, -38.99126326, -94.48017637, 62.72058643, -45.39487981, -82.811656, 126.92432487, -7.39505416, -98.50506137, 4.5e-07, 2.12e-06, 1.681e-05, 0, 2.2e-07, 6.2e-06), 
                c(2e-07, 1e-07, 1.2e-07, 4e-08, 3e-08, 0, 0, 4e-08, 1e-08, 1.84e-06, 454.00000034, 1e-08, 109.9178027, 107.27752212, 100.75409393, -1e-08, 0, 0, 1e-08, 1.5e-07, -8e-08, 1e-08, -3e-08, -1.07e-06, 24.31730997, 7.70605243, 4.866079, 17.66372364, 24.22692056, 1.50695768), 
                c(1e-07, 6e-08, 9e-08, 2e-08, -2e-07, 1e-08, -3e-08, 1e-08, -5e-08, 1.13e-06, 1e-08, 326.00000229, -1e-08, -1e-08, 1e-08, 81.51596907, 84.15039213, 74.95842051, -9.63461324, -8.33228005, -7.6972618, 11.1501336, -35.56568222, -3.07396025, -4e-08, 3e-08, -1.43e-06, -5e-08, -3e-08, 5e-08), 
                c(-1e-08, 0, -2e-08, 0, -55.6535616, 0, 0, 0, 0, 0, 109.9178027, -1e-08, 112.7824544, -1.31024125, -2.24418408, 0, 0, 0, 0, 0, -1e-08, -1e-08, 0, -1.2e-07, 11.94927774, -0.19182148, -0.19004235, 8.70562069, 11.88553603, -0.16828589), 
                c(1e-08, 0, 1e-08, 0, -16.1141946, 0, 0, 0, 0, 8e-08, 107.27752212, -1e-08, -1.31024125, 102.94836181, -0.06439274, 0, 0, 0, -1e-08, 0, 0, -1e-08, 0, -1.7e-07, -0.18216424, 3.80923856, -0.13344723, 8.54620445, -0.1598135, 0.86031718), 
                c(0, 1e-08, 1e-08, 0, 140.80528556, 0, 0, 0, 0, 8e-08, 100.75409393, 1e-08, -2.24418408, -0.06439274, 91.58135278, 0, 0, 0, 0, 0, 0, 0, 0, 1e-08, -0.15910624, -0.11764696, 2.32434774, -0.14089147, 10.94386634, 0.81552409), 
                c(1e-08, 0, 0, 0, 199.48144369, 0, 0, 0, 0, 199.48144381, -1e-08, 81.51596907, 0, 0, 0, 89.83790539, -0.89130901, -6.1318067, -4.72773497, -0.1325759, -0.06399443, 5.71221553, -17.72741051, -0.09324222, 0, -1e-08, 0, -1e-08, 0, -4e-08), 
                c(1e-08, 0, 0, 0, -38.99126342, 0, 0, 0, 0, -38.99126326, 0, 84.15039213, 0, 0, 0, -0.89130901, 88.63803396, -0.09068914, -0.14189633, -4.23553146, -0.08853111, 5.85757284, -0.09980804, -1.53685315, 0, 0, -8e-08, 0, 0, 1e-08), 
                c(2e-08, 0, 0, 0, -94.48017649, 0, 0, 0, 0, -94.48017637, 0, 74.95842051, 0, 0, 0, -6.1318067, -0.09068914, 74.93849286, -0.05382399, -0.06955407, -3.45331919, -0.07842752, -16.22805835, -1.29681765, 0, 0, -9e-08, 0, 0, 0), 
                c(4e-08, -1e-08, -6e-08, 0, 62.7205861, 0, -1e-08, 0, -2e-08, 62.72058643, 1e-08, -9.63461324, 0, -1e-08, 0, -4.72773497, -0.14189633, -0.05382399, 50.9890118, 6.40124366, -17.87273437, -5.43814218, -2.50747803, 1.08339806, -4e-08, 0, -1.49e-06, -3e-08, -4e-08, 0), 
                c(9e-08, 7e-08, 8e-08, 2e-08, -45.39488044, 1e-08, 2e-08, 2e-08, 2e-08, -45.39487981, 1.5e-07, -8.33228005, 0, 0, 0, -0.1325759, -4.23553146, -0.06955407, 6.40124366, 43.64783296, -1.11149819, -5.43451324, -0.44808313, -3.77487948, 0, 1.8e-07, 0, -4e-08, 0, 4.4e-07), 
                c(9e-08, -3e-08, 3e-08, 0, -82.81165625, 0, -1e-08, 0, -4e-08, -82.811656, -8e-08, -7.6972618, -1e-08, 0, 0, -0.06399443, -0.08853111, -3.45331919, -17.87273437, -1.11149819, 50.00966136, 1.76934478, -3.8485412, -4.50868395, -2e-08, -9e-08, -8e-07, -6e-08, -2e-08, -2.9e-07), 
                c(1.7e-07, 1e-08, -2e-08, 1e-08, 126.92432422, -1e-08, 1e-08, 0, -2e-08, 126.92432487, 1e-08, 11.1501336, -1e-08, -1e-08, 0, 5.71221553, 5.85757284, -0.07842752, -5.43814218, -5.43451324, 1.76934478, 50.71957223, -4.22066211, -18.35978405, -2e-08, -1e-07, -8.9e-07, -6e-08, -2e-08, -3.5e-07), 
                c(3e-08, 1e-08, 0, 0, -7.39505442, 0, 0, 0, -1e-08, -7.39505416, -3e-08, -35.56568222, 0, 0, 0, -17.72741051, -0.09980804, -16.22805835, -2.50747803, -0.44808313, -3.8485412, -4.22066211, 64.86434918, 7.40531477, -1e-08, -3e-08, -5e-07, -2e-08, -1e-08, -9e-08), 
                c(4e-07, -1.8e-07, -4e-08, 9e-08, -98.5050688, -1e-07, -1.2e-07, -1.4e-07, -8e-08, -98.50506137, -1.07e-06, -3.07396025, -1.2e-07, -1.7e-07, 1e-08, -0.09324222, -1.53685315, -1.29681765, 1.08339806, -3.77487948, -4.50868395, -18.35978405, 7.40531477, 49.19089192, -6.4e-07, -2.99e-06, -2.275e-05, -8.4e-07, -2.3e-07, -1.355e-05), 
                c(2e-08, 2e-08, -5e-08, 1e-08, -86.80408649, 0, 0, 0, 0, 4.5e-07, 24.31730997, -4e-08, 11.94927774, -0.18216424, -0.15910624, 0, 0, 0, -4e-08, 0, -2e-08, -2e-08, -1e-08, -6.4e-07, 68.49310859, 10.2811001, 13.39724145, 6.3299363, 5.47958227, 1.11448435), 
                c(2.6e-07, 1.9e-07, 2.8e-07, 3e-08, -60.47376097, 0, 0, 5e-08, 0, 2.12e-06, 7.70605243, 3e-08, -0.19182148, 3.80923856, -0.11764696, -1e-08, 0, 0, 0, 1.8e-07, -9e-08, -1e-07, -3e-08, -2.99e-06, 10.2811001, 61.22307964, 1.0795372, 13.34154613, 0.42815576, 2.71608178), 
                c(8.2e-07, 2.3e-07, 0, -7e-08, 56.0232937, -2e-08, 0, 0, 0, 1.681e-05, 4.866079, -1.43e-06, -0.19004235, -0.13344723, 2.32434774, 0, -8e-08, -9e-08, -1.49e-06, 0, -8e-07, -8.9e-07, -5e-07, -2.275e-05, 13.39724145, 1.0795372, 62.13284815, 0.89834529, 13.50221379, 4.22666712), 
                c(9e-08, 3e-08, -3e-08, 1e-08, -39.89494287, 0, 0, 0, 0, 0, 17.66372364, -5e-08, 8.70562069, 8.54620445, -0.14089147, -1e-08, 0, 0, -3e-08, -4e-08, -6e-08, -6e-08, -2e-08, -8.4e-07, 6.3299363, 13.34154613, 0.89834529, 66.00144213, 2.87473442, 12.76327821), 
                c(5e-08, 0, -5e-08, 0, 99.31782595, 0, 0, 0, 0, 2.2e-07, 24.22692056, -3e-08, 11.88553603, -0.1598135, 10.94386634, 0, 0, 0, -4e-08, 0, -2e-08, -2e-08, -1e-08, -2.3e-07, 5.47958227, 0.42815576, 13.50221379, 2.87473442, 68.39943276, 9.63020519), 
                c(1.6e-07, 4.2e-07, 1.2e-07, 8e-08, 76.6938026, -1e-08, -1e-08, 0, 0, 6.2e-06, 1.50695768, 5e-08, -0.16828589, 0.86031718, 0.81552409, -4e-08, 1e-08, 0, 0, 4.4e-07, -2.9e-07, -3.5e-07, -9e-08, -1.355e-05, 1.11448435, 2.71608178, 4.22666712, 12.76327821, 9.63020519, 60.85664305)
                )
    expect_equal(as.double(test), as.double(GS), 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.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 <- cbind(c(481.58902557, 141.4145761, 126.4937567, 90.87132761, 24505.0293416, 151.42333095, 39.89451566, 53.34143661, 22.82307698, 7519.7846284, 1e-08, 1e-08, 0, 2e-08, 1e-08, 1e-08, 1e-08, 0, 4e-08, 9e-08, 9e-08, 1.7e-07, 3e-08, 4e-07, 2e-08, 2.6e-07, 8.2e-07, 9e-08, 5e-08, 1.6e-07), 
      c(141.4145761, 123.53334857, 16.52169525, -6.89363657, 7203.84822301, 39.89451564, 57.11914728, 7.08254392, -18.90835878, 1981.18852358, 0, 1e-08, 0, 0, 0, 1e-08, 0, 0, -1e-08, 7e-08, -3e-08, 1e-08, 1e-08, -1.8e-07, 2e-08, 1.9e-07, 2.3e-07, 3e-08, 0, 4.2e-07), 
      c(126.4937567, 16.52169525, 110.82914443, -0.53727867, 6412.26409583, 53.34143653, 7.08254388, 52.49669251, -1.50781625, 2648.97167493, 0, 0, 1e-08, 0, 0, 1e-08, 0, 1e-08, -6e-08, 8e-08, 3e-08, -2e-08, 0, -4e-08, -5e-08, 2.8e-07, 0, -3e-08, -5e-08, 1.2e-07), 
      c(90.87132761, -6.89363657, -0.53727867, 99.56459117, 4634.12494773, 22.82307691, -18.90835881, -1.50781628, 47.23106187, 1133.40937526, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2e-08, 0, 1e-08, 0, 9e-08, 1e-08, 3e-08, -7e-08, 1e-08, 0, 8e-08), 
      c(24505.0293416, 7203.84822301, 6412.26409583, 4634.12494773, 1306448.44293833, 7519.78462743, 1981.18852307, 2648.97167428, 1133.40937504, 387540.13845235, -46.30599547, -33.46207953, -29.64480404, -8.36463948, 68.62673186, 113.71282932, -22.99507592, -49.38832322, 62.7205861, -45.39488044, -82.81165625, 126.92432422, -7.39505442, -98.5050688, -86.80408649, -60.47376097, 56.0232937, -39.89494287, 99.31782595, 76.6938026), 
      c(151.42333095, 39.89451564, 53.34143653, 22.82307691, 7519.78462743, 151.42333094, 39.89451563, 53.34143654, 22.8230769, 7519.7846276, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1e-08, 0, -1e-08, 0, -1e-07, 0, 0, -2e-08, 0, 0, -1e-08), 
      c(39.89451566, 57.11914728, 7.08254388, -18.90835881, 1981.18852307, 39.89451563, 57.11914725, 7.08254391, -18.90835882, 1981.18852349, 0, 0, 0, 0, 0, 0, 0, 0, -1e-08, 2e-08, -1e-08, 1e-08, 0, -1.2e-07, 0, 0, 0, 0, 0, -1e-08), 
      c(53.34143661, 7.08254392, 52.49669251, -1.50781628, 2648.97167428, 53.34143654, 7.08254391, 52.49669251, -1.50781629, 2648.97167472, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2e-08, 0, 0, 0, -1.4e-07, 0, 5e-08, 0, 0, 0, 0), 
      c(22.82307698, -18.90835878, -1.50781625, 47.23106187, 1133.40937504, 22.8230769, -18.90835882, -1.50781629, 47.23106187, 1133.40937519, 0, 0, 0, 0, 0, 0, 0, 0, -2e-08, 2e-08, -4e-08, -2e-08, -1e-08, -8e-08, 0, 0, 0, 0, 0, 0), 
      c(7519.7846284, 1981.18852358, 2648.97167493, 1133.40937526, 387540.13845235, 7519.7846276, 1981.18852349, 2648.97167472, 1133.40937519, 387540.13844065, 1e-07, -33.46207949, 9e-08, 4e-08, 8e-08, 113.71282936, -22.99507584, -49.38832318, 62.72058643, -45.39487981, -82.811656, 126.92432487, -7.39505416, -98.50506137, 4.5e-07, 2.12e-06, 1.681e-05, 0, 2.2e-07, 6.2e-06), 
      c(1e-08, 0, 0, 0, -46.30599547, 0, 0, 0, 0, 1e-07, 35.16315564, 0, -0.73325294, -0.07775833, -0.02930618, 0, 0, 0, 0, 1e-08, 0, 0, 0, -1.3e-07, 6.79325508, 2.2002666, 1.41925623, -0.1121543, -0.11111399, -0.08216028), 
      c(1e-08, 1e-08, 0, 0, -33.46207953, 0, 0, 0, 0, -33.46207949, 0, 27.03459311, 0, 0, 0, -0.2372027, -0.16490876, -0.14958512, -2.68554804, -2.29727168, -2.15586471, -0.0754545, -0.03642425, -0.04707602, -1e-08, 2e-08, -2.2e-07, 0, 0, 0), 
      c(0, 0, 1e-08, 0, -29.64480404, 0, 0, 0, 0, 9e-08, -0.73325294, 0, 32.00020684, -0.36228097, -0.58262412, 0, 0, 0, -1e-08, 0, 0, -1e-08, 0, -5e-08, 6.36498344, -0.10217692, -0.10122932, 4.63719504, 6.33103034, -0.08964022), 
      c(2e-08, 0, 0, 0, -8.36463948, 0, 0, 0, 0, 4e-08, -0.07775833, 0, -0.36228097, 27.73937621, -0.01629108, 0, 0, 0, 0, 0, 0, 0, 0, -5e-08, -0.09455876, 1.97731928, -0.06927055, 4.43620803, -0.08295682, 0.44657793), 
      c(1e-08, 0, 0, 0, 68.62673186, 0, 0, 0, 0, 8e-08, -0.02930618, 0, -0.58262412, -0.01629108, 21.75483409, 0, 0, 0, -1e-08, 0, -1e-08, -1e-08, 0, -1.6e-07, -0.07754639, -0.05733966, 1.13285783, -0.06866875, 5.33390332, 0.3974762), 
      c(1e-08, 1e-08, 1e-08, 0, 113.71282932, 0, 0, 0, 0, 113.71282936, 0, -0.2372027, 0, 0, 0, 29.19263486, -0.29964208, -1.8271702, -2.69500817, -0.07557383, -0.03647952, 3.25620358, -10.10537105, -0.05315193, 0, 2e-08, 1e-08, 0, 0, 0), 
      c(1e-08, 0, 0, 0, -22.99507592, 0, 0, 0, 0, -22.99507584, 0, -0.16490876, 0, 0, 0, -0.29964208, 30.82870221, -0.02795802, -0.08368329, -2.49790233, -0.05221117, 3.45450031, -0.05886174, -0.90635828, 0, 2e-08, -1e-07, 0, 0, -4e-08), 
      c(0, 0, 1e-08, 0, -49.38832322, 0, 0, 0, 0, -49.38832318, 0, -0.14958512, 0, 0, 0, -1.8271702, -0.02795802, 20.47727971, -0.02813582, -0.03635852, -1.80517916, -0.040997, -8.48301328, -0.6778952, 0, -1e-08, -1e-07, -1e-08, 0, -8e-08), 
      c(4e-08, -1e-08, -6e-08, 0, 62.7205861, 0, -1e-08, 0, -2e-08, 62.72058643, 0, -2.68554804, -1e-08, 0, -1e-08, -2.69500817, -0.08368329, -0.02813582, 50.9890118, 6.40124366, -17.87273437, -5.43814218, -2.50747803, 1.08339806, -4e-08, 0, -1.49e-06, -3e-08, -4e-08, 0), 
      c(9e-08, 7e-08, 8e-08, 2e-08, -45.39488044, 1e-08, 2e-08, 2e-08, 2e-08, -45.39487981, 1e-08, -2.29727168, 0, 0, 0, -0.07557383, -2.49790233, -0.03635852, 6.40124366, 43.64783296, -1.11149819, -5.43451324, -0.44808313, -3.77487948, 0, 1.8e-07, 0, -4e-08, 0, 4.4e-07), 
      c(9e-08, -3e-08, 3e-08, 0, -82.81165625, 0, -1e-08, 0, -4e-08, -82.811656, 0, -2.15586471, 0, 0, -1e-08, -0.03647952, -0.05221117, -1.80517916, -17.87273437, -1.11149819, 50.00966136, 1.76934478, -3.8485412, -4.50868395, -2e-08, -9e-08, -8e-07, -6e-08, -2e-08, -2.9e-07), 
      c(1.7e-07, 1e-08, -2e-08, 1e-08, 126.92432422, -1e-08, 1e-08, 0, -2e-08, 126.92432487, 0, -0.0754545, -1e-08, 0, -1e-08, 3.25620358, 3.45450031, -0.040997, -5.43814218, -5.43451324, 1.76934478, 50.71957223, -4.22066211, -18.35978405, -2e-08, -1e-07, -8.9e-07, -6e-08, -2e-08, -3.5e-07), 
      c(3e-08, 1e-08, 0, 0, -7.39505442, 0, 0, 0, -1e-08, -7.39505416, 0, -0.03642425, 0, 0, 0, -10.10537105, -0.05886174, -8.48301328, -2.50747803, -0.44808313, -3.8485412, -4.22066211, 64.86434918, 7.40531477, -1e-08, -3e-08, -5e-07, -2e-08, -1e-08, -9e-08), 
      c(4e-07, -1.8e-07, -4e-08, 9e-08, -98.5050688, -1e-07, -1.2e-07, -1.4e-07, -8e-08, -98.50506137, -1.3e-07, -0.04707602, -5e-08, -5e-08, -1.6e-07, -0.05315193, -0.90635828, -0.6778952, 1.08339806, -3.77487948, -4.50868395, -18.35978405, 7.40531477, 49.19089192, -6.4e-07, -2.99e-06, -2.275e-05, -8.4e-07, -2.3e-07, -1.355e-05), 
      c(2e-08, 2e-08, -5e-08, 1e-08, -86.80408649, 0, 0, 0, 0, 4.5e-07, 6.79325508, -1e-08, 6.36498344, -0.09455876, -0.07754639, 0, 0, 0, -4e-08, 0, -2e-08, -2e-08, -1e-08, -6.4e-07, 68.49310859, 10.2811001, 13.39724145, 6.3299363, 5.47958227, 1.11448435), 
      c(2.6e-07, 1.9e-07, 2.8e-07, 3e-08, -60.47376097, 0, 0, 5e-08, 0, 2.12e-06, 2.2002666, 2e-08, -0.10217692, 1.97731928, -0.05733966, 2e-08, 2e-08, -1e-08, 0, 1.8e-07, -9e-08, -1e-07, -3e-08, -2.99e-06, 10.2811001, 61.22307964, 1.0795372, 13.34154613, 0.42815576, 2.71608178), 
      c(8.2e-07, 2.3e-07, 0, -7e-08, 56.0232937, -2e-08, 0, 0, 0, 1.681e-05, 1.41925623, -2.2e-07, -0.10122932, -0.06927055, 1.13285783, 1e-08, -1e-07, -1e-07, -1.49e-06, 0, -8e-07, -8.9e-07, -5e-07, -2.275e-05, 13.39724145, 1.0795372, 62.13284815, 0.89834529, 13.50221379, 4.22666712), 
      c(9e-08, 3e-08, -3e-08, 1e-08, -39.89494287, 0, 0, 0, 0, 0, -0.1121543, 0, 4.63719504, 4.43620803, -0.06866875, 0, 0, -1e-08, -3e-08, -4e-08, -6e-08, -6e-08, -2e-08, -8.4e-07, 6.3299363, 13.34154613, 0.89834529, 66.00144213, 2.87473442, 12.76327821), 
      c(5e-08, 0, -5e-08, 0, 99.31782595, 0, 0, 0, 0, 2.2e-07, -0.11111399, 0, 6.33103034, -0.08295682, 5.33390332, 0, 0, 0, -4e-08, 0, -2e-08, -2e-08, -1e-08, -2.3e-07, 5.47958227, 0.42815576, 13.50221379, 2.87473442, 68.39943276, 9.63020519), 
      c(1.6e-07, 4.2e-07, 1.2e-07, 8e-08, 76.6938026, -1e-08, -1e-08, 0, 0, 6.2e-06, -0.08216028, 0, -0.08964022, 0.44657793, 0.3974762, 0, -4e-08, -8e-08, 0, 4.4e-07, -2.9e-07, -3.5e-07, -9e-08, -1.355e-05, 1.11448435, 2.71608178, 4.22666712, 12.76327821, 9.63020519, 60.85664305)
)
    expect_equal(as.double(test), as.double(GS), tol = 1e-6)

    ## 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)
    GS <- cbind(c(254.24318231, 71.76445703, 63.87468511, 42.26944905, 12924.03081807, 87.11643718, 21.99098267, 29.42732681, 10.75450504, 4326.26096348, -1685.25252981, -1726.68115882, -469.83394986, -318.00367208, -274.74336896, -414.54823988, -571.70133222, -187.56664227, -314.54812398, -407.83236118, -166.70298638, -403.27731988, -164.83442591, -213.7219395, -586.32171168, -416.0100563, -403.40300041, -460.03702007, -446.09385133, -316.51037674), 
                c(71.76445703, 84.43958198, 4.68009386, -14.44164787, 3652.66391671, 21.99098328, 42.93868024, 3.64466591, -17.80841298, 1092.08704692, -504.06971427, -440.9162199, -430.761276, -4.77655839, -14.66016296, -611.4012118, -35.32324115, 154.63145745, -257.27754453, 31.76949129, 158.22138849, -418.42304079, -24.83696373, 163.33350651, -264.50314921, -22.01509943, -36.66976293, -198.94377374, -210.09147055, -21.65770202), 
                c(63.87468511, 4.68009386, 75.46467938, -7.23359472, 3233.51172295, 29.42732682, 3.64466597, 35.92910323, -3.23861152, 1461.38087738, -349.77017263, -589.17664735, -4.89048996, -341.47625608, 17.39566751, -34.25022399, -634.71501177, 28.12104531, 18.88935919, -186.57280892, 48.14045054, -279.75897635, 11.10961124, -98.61846211, 6.20946658, -227.55291796, 33.73702541, -268.4410888, 23.57608675, -161.66434405), 
                c(42.26944905, -14.44164787, -7.23359472, 72.24811722, 2155.34990444, 10.75450471, -17.80841206, -3.23861116, 36.57718587, 534.07595046, -320.29503872, -216.18886504, -15.90175999, 18.42935809, -292.8401056, 167.35198938, 31.38785378, -411.7864569, 155.51717831, 61.80531573, -261.11741824, 185.37910164, -213.11551281, -350.92953675, 0.65579184, 38.29981353, -211.83431463, 11.06206915, -267.12705195, -160.67894449), 
                c(12924.03081807, 3652.66391671, 3233.51172295, 2155.34990444, 688101.15548783, 4326.26096324, 1092.08704716, 1461.3808776, 534.07595024, 222958.4829556, -91043.43888476, -88881.62971152, -25401.77112173, -17182.13567772, -14759.89276701, -21203.33350271, -29437.85032769, -9715.60136513, -16145.25690429, -21028.24848252, -8643.28684641, -20664.3922041, -8492.85373763, -11060.40761815, -31730.32736475, -22510.76473119, -21758.53782695, -24865.12288714, -24034.7390626, -17051.07701002), 
                c(87.11643718, 21.99098328, 29.42732682, 10.75450471, 4326.26096324, 87.1164373, 21.99098329, 29.42732687, 10.75450498, 4326.26096323, 0, -1726.6811588, 0, 1e-08, 1e-08, -414.54823987, -571.70133223, -187.56664226, -314.54812399, -407.83236156, -166.70298636, -403.27732221, -164.83442597, -213.72193982, 0, 0, 1.6e-07, 3.7e-07, 0, -7e-08), 
                c(21.99098267, 42.93868024, 3.64466597, -17.80841206, 1092.08704716, 21.99098329, 42.93868229, 3.64466582, -17.80841384, 1092.08704707, 0, -440.91621995, 0, 0, 0, -611.40121185, -35.32324119, 154.6314574, -257.27754486, 31.76949117, 158.22138828, -418.42304136, -24.83696382, 163.33350655, 3e-08, 0, 0, -8e-08, -3e-08, 3e-08), 
                c(29.42732681, 3.64466591, 35.92910323, -3.23861116, 1461.3808776, 29.42732687, 3.64466582, 35.92910324, -3.23861168, 1461.38087749, 0, -589.17664737, 0, 0, 1e-08, -34.25022397, -634.71501178, 28.12104529, 18.88935921, -186.57280909, 48.14045059, -279.75897944, 11.10961119, -98.61846248, 1e-08, 2e-07, 1.3e-07, 2.7e-07, -1e-08, 1e-08), 
                c(10.75450504, -17.80841298, -3.23861152, 36.57718587, 534.07595024, 10.75450498, -17.80841384, -3.23861168, 36.57718122, 534.07595186, 0, -216.18886506, 0, 1.1e-07, 1.1e-07, 167.35198937, 31.38785372, -411.78645698, 155.51717801, 61.80531579, -261.11741745, 185.37910791, -213.11551266, -350.92953694, 0, -2.19e-06, 1.36e-06, 3.19e-06, 0, -1.14e-06), 
                c(4326.26096348, 1092.08704692, 1461.38087738, 534.07595046, 222958.4829556, 4326.26096323, 1092.08704707, 1461.38087749, 534.07595186, 222958.48295492, 1e-08, -88881.62971144, 0, 0, 0, -21203.33350265, -29437.85032765, -9715.60136507, -16145.25690448, -21028.2484823, -8643.28684641, -20664.39219991, -8492.85373762, -11060.40761838, 0, 0, 3e-08, -8e-08, 0, 0), 
                c(-1685.25252981, -504.06971427, -349.77017263, -320.29503872, -91043.43888476, 0, 0, 0, 0, 1e-08, 13604.49242314, -1e-08, 2563.84567435, 1757.06336791, 1508.27378135, 0, -1e-08, 0, 2e-08, 3e-08, 2e-08, -2.2e-07, -1e-08, -8e-08, 3093.93260142, 2187.10076688, 2104.28729566, 2442.35130733, 2360.02389702, 1660.9057482), 
                c(-1726.68115882, -440.9162199, -589.17664735, -216.18886504, -88881.62971152, -1726.6811588, -440.91621995, -589.17664737, -216.18886506, -88881.62971144, -1e-08, 26769.98049706, 0, -1e-08, -1e-08, 4288.9568048, 5965.86978793, 2016.17040113, 3165.16281546, 4141.62562364, 1709.75340866, 4151.08646439, 1674.03852124, 2232.97978755, -8e-08, -2.1e-07, -4e-08, -2e-07, -1.2e-07, -3e-08), 
                c(-469.83394986, -430.761276, -4.89048996, -15.90175999, -25401.77112173, 0, 0, 0, 0, 0, 2563.84567435, 0, 3347.35271579, 23.47824893, 71.7421877, -1e-08, 0, 0, -1e-08, -2e-08, -1e-08, -2e-07, -1e-08, -4e-08, 1313.73774478, 108.80106317, 180.97680782, 984.49320892, 1036.33854862, 106.84414197), 
                c(-318.00367208, -4.77655839, -341.47625608, 18.42935809, -17182.13567772, 1e-08, 0, 0, 1.1e-07, 0, 1757.06336791, -1e-08, 23.47824893, 2467.93754137, -83.83308856, 0, -1e-08, 0, 0, 3e-08, 0, -2.2e-07, -1e-08, -4e-08, -30.01433247, 1098.70246188, -162.84627518, 1296.12715248, -114.00793336, 769.68404901), 
                c(-274.74336896, -14.66016296, 17.39566751, -292.8401056, -14759.89276701, 1e-08, 0, 1e-08, 1.1e-07, 0, 1508.27378135, -1e-08, 71.7421877, -83.83308856, 1986.67629707, 0, -1e-08, 0, 0, 0, -1e-08, 2e-08, -1e-08, 0, -2.90197157, -173.7924841, 959.48051191, -50.28802291, 1213.2657455, 724.0476357), 
                c(-414.54823988, -611.4012118, -34.25022399, 167.35198938, -21203.33350271, -414.54823987, -611.40121185, -34.25022397, 167.35198937, -21203.33350265, 0, 4288.9568048, -1e-08, 0, 0, 7939.27669298, 342.27188373, -1506.40817921, 2500.26023072, -307.07604093, -1539.21108571, 4054.43352205, 249.06857644, -1581.20762488, -7e-08, -6e-08, -3e-08, -9e-08, -7e-08, 0), 
                c(-571.70133222, -35.32324115, -634.71501177, 31.38785378, -29437.85032769, -571.70133223, -35.32324119, -634.71501178, 31.38785372, -29437.85032765, -1e-08, 5965.86978793, 0, -1e-08, -1e-08, 342.27188373, 9291.0164323, -282.90000023, -184.99215326, 1875.11176139, -486.32414398, 2789.19059695, -108.0596607, 999.59776839, -3e-08, -1.3e-07, -3e-08, -1.9e-07, -7e-08, 0), 
                c(-187.56664227, 154.63145745, 28.12104531, -411.7864569, -9715.60136513, -187.56664226, 154.6314574, 28.12104529, -411.78645698, -9715.60136507, 0, 2016.17040113, 0, 0, 0, -1506.40817921, -282.90000023, 4618.83374843, -1400.66458758, -556.75547986, 2352.83689877, -1665.76049115, 1871.72040806, 3152.80473172, -7e-08, -1.1e-07, 4e-08, 0, -7e-08, -3e-08), 
                c(-314.54812398, -257.27754453, 18.88935919, 155.51717831, -16145.25690429, -314.54812399, -257.27754486, 18.88935921, 155.51717801, -16145.25690448, 2e-08, 3165.16281546, -1e-08, 0, 0, 2500.26023072, -184.99215326, -1400.66458758, 3437.47328404, 1850.78664073, -428.49771503, 2314.31401937, -208.68372341, -1664.08152245, 3e-08, -2.8e-07, 0, -4e-08, -3.7e-07, 0), 
                c(-407.83236118, 31.76949129, -186.57280892, 61.80531573, -21028.24848252, -407.83236156, 31.76949117, -186.57280909, 61.80531579, -21028.2484823, 3e-08, 4141.62562364, -2e-08, 3e-08, 0, -307.07604093, 1875.11176139, -556.75547986, 1850.78664073, 4067.38830859, 632.26310835, 1011.11985831, -609.85676953, 74.89510741, -2.3e-07, -8e-07, 4e-08, 0, -2.2e-07, -2.1e-07), 
                c(-166.70298638, 158.22138849, 48.14045054, -261.11741824, -8643.28684641, -166.70298636, 158.22138828, 48.14045059, -261.11741745, -8643.28684641, 2e-08, 1709.75340866, -1e-08, 0, -1e-08, -1539.21108571, -486.32414398, 2352.83689877, -428.49771503, 632.26310835, 2509.81200925, -1843.01608217, 1427.19223845, 2495.519007, -3.7e-07, -6e-07, 0, 0, -5.5e-07, -1.5e-07), 
                c(-403.27731988, -418.42304079, -279.75897635, 185.37910164, -20664.3922041, -403.27732221, -418.42304136, -279.75897944, 185.37910791, -20664.39219991, -2.2e-07, 4151.08646439, -2e-07, -2.2e-07, 2e-08, 4054.43352205, 2789.19059695, -1665.76049115, 2314.31401937, 1011.11985831, -1843.01608217, 5947.23203622, 180.56394375, -859.54195269, 4.9e-07, -8.6e-06, 0, 1e-08, -2.57e-06, 2e-07), 
                c(-164.83442591, -24.83696373, 11.10961124, -213.11551281, -8492.85373763, -164.83442597, -24.83696382, 11.10961119, -213.11551266, -8492.85373762, -1e-08, 1674.03852124, -1e-08, -1e-08, -1e-08, 249.06857644, -108.0596607, 1871.72040806, -208.68372341, -609.85676953, 1427.19223845, 180.56394375, 1809.36562987, 2119.81992769, -1.8e-07, -1.3e-07, 2e-08, -2e-07, -1.6e-07, -1e-08), 
                c(-213.7219395, 163.33350651, -98.61846211, -350.92953675, -11060.40761815, -213.72193982, 163.33350655, -98.61846248, -350.92953694, -11060.40761838, -8e-08, 2232.97978755, -4e-08, -4e-08, 0, -1581.20762488, 999.59776839, 3152.80473172, -1664.08152245, 74.89510741, 2495.519007, -859.54195269, 2119.81992769, 4009.54679281, -9.5e-07, -6.5e-07, 5.2e-07, 0, -8.9e-07, 0), 
                c(-586.32171168, -264.50314921, 6.20946658, 0.65579184, -31730.32736475, 0, 3e-08, 1e-08, 0, 0, 3093.93260142, -8e-08, 1313.73774478, -30.01433247, -2.90197157, -7e-08, -3e-08, -7e-08, 3e-08, -2.3e-07, -3.7e-07, 4.9e-07, -1.8e-07, -9.5e-07, 3642.37366676, 1368.52211531, 1355.13022652, 1238.2975598, 1229.46053944, -33.18725375), 
                c(-416.0100563, -22.01509943, -227.55291796, 38.29981353, -22510.76473119, 0, 0, 2e-07, -2.19e-06, 0, 2187.10076688, -2.1e-07, 108.80106317, 1098.70246188, -173.7924841, -6e-08, -1.3e-07, -1.1e-07, -2.8e-07, -8e-07, -6e-07, -8.6e-06, -1.3e-07, -6.5e-07, 1368.52211531, 2344.58234317, 584.88507596, 1701.73311588, -181.37104293, 885.41956388), 
                c(-403.40300041, -36.66976293, 33.73702541, -211.83431463, -21758.53782695, 1.6e-07, 0, 1.3e-07, 1.36e-06, 3e-08, 2104.28729566, -4e-08, 180.97680782, -162.84627518, 959.48051191, -3e-08, -3e-08, 4e-08, 0, 4e-08, 0, 0, 2e-08, 5.2e-07, 1355.13022652, 584.88507596, 2223.19872755, -58.86208849, 1744.68177124, 953.02849857), 
                c(-460.03702007, -198.94377374, -268.4410888, 11.06206915, -24865.12288714, 3.7e-07, -8e-08, 2.7e-07, 3.19e-06, -8e-08, 2442.35130733, -2e-07, 984.49320892, 1296.12715248, -50.28802291, -9e-08, -1.9e-07, 0, -4e-08, 0, 0, 1e-08, -2e-07, 0, 1238.2975598, 1701.73311588, -58.86208849, 2866.35444736, 842.76055259, 1227.3960449), 
                c(-446.09385133, -210.09147055, 23.57608675, -267.12705195, -24034.7390626, 0, -3e-08, -1e-08, 0, 0, 2360.02389702, -1.2e-07, 1036.33854862, -114.00793336, 1213.2657455, -7e-08, -7e-08, -7e-08, -3.7e-07, -2.2e-07, -5.5e-07, -2.57e-06, -1.6e-07, -8.9e-07, 1229.46053944, -181.37104293, 1744.68177124, 842.76055259, 2980.15548721, 1289.98604601), 
                c(-316.51037674, -21.65770202, -161.66434405, -160.67894449, -17051.07701002, -7e-08, 3e-08, 1e-08, -1.14e-06, 0, 1660.9057482, -3e-08, 106.84414197, 769.68404901, 724.0476357, 0, 0, -3e-08, 0, -2.1e-07, -1.5e-07, 2e-07, -1e-08, 0, -33.18725375, 885.41956388, 953.02849857, 1227.3960449, 1289.98604601, 1642.25537941)
                )
    test <- information(eSUN.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 <- confint(eSUN.lmm, effects = "all")
    GS <- c(75.02445251, 57.00262333, 57.00100521, 57.00082937, 55.02679117, 67.96985135, 96.26931669, 94.16508864, 93.47141533, 54.92384754, 56.85151868, 41.11400945, 109.52022024, 113.40272184, 114.03144185, 81.5447874, 81.63717886, 79.9545627, 20.70521888, 20.51343872, 20.7126253, 20.32470404, 23.65317661, 20.31145426, 29.54320962, 28.55308578, 28.49261999, 29.06453625, 29.60876751, 28.48658921)
    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("2:2" = 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"))),
                     "1:1" = 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), 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