tests/testthat/test-auto-predict.R

In bozenne/repeated: Repeated Measurement Models for Discrete Times

### test-auto-predict.R --- 
##----------------------------------------------------------------------
## Author: Brice Ozenne
## Created: feb 22 2024 (10:15) 
## Version: 
## Last-Updated: jul 15 2024 (10:21) 
##           By: Brice Ozenne
##     Update #: 34
##----------------------------------------------------------------------
## 
### Commentary: 
## 
### Change Log:
##----------------------------------------------------------------------
## 
### Code:

if(FALSE){
    ## library(pracma) ## trapz
    library(nlmeU)
    library(testthat)

    library(LMMstar)
}

context("Check predict and fitted functions")
LMMstar.options(optimizer = "FS", precompute.moments = TRUE,
                columns.confint = c("estimate","se","df","lower","upper","p.value"))

## * load and reshape data
data("armd.wide", package = "nlmeU")
armd.long <- reshape(armd.wide, direction = "long", idvar = "subject",
                     timevar = "week.num", times = c(0,4,12,24,52), varying = paste0("visual",c(0,4,12,24,52)), 
                     v.names = "visual")
armd.long$week <- as.factor(armd.long$week.num)
armd.long$lesion[is.na(armd.long$lesion)] <- 1
armd.long$week2 <- armd.long$week

armd.longNNA <- armd.long[!is.na(armd.long$lesion) & !is.na(armd.long$visual),]
## sum(is.na(armd.longNNA))
## 0

nd <- armd.long[armd.long$subject %in% c(2,10),]
nd.X <- model.matrix( ~ lesion + week + week:treat.f, data = nd)
nd2 <- armd.long[armd.long$subject %in% 1:3,]
nd2[nd2$subject == 1, "visual"] <- c(59,55,45,40,38)
nd2[nd2$subject == 2, "visual"] <- NA
nd2[nd2$subject == 3, "visual"] <- c(40,40,37,NA,NA)

ndNNA <- armd.longNNA[armd.longNNA$subject %in% c(1,2,10),]
ndNNA.X <- model.matrix( ~ lesion + week + week:treat.f, data = ndNNA)

## * predict - lm
test_that("predict (lmm)", {
    e.lm <- lm(visual ~ lesion + week + week:treat.f,
               data = armd.long)
    e.lm2 <- suppressWarnings(lmm(visual ~ lesion + week + week:treat.f,
                                  data = armd.long))
    
    expect_equivalent(predict(e.lm, newdata = armd.long, type = "terms"),
                      predict(e.lm2, newdata = armd.long, type = "terms"),
                      tol = 1e-6)

    GS <- predict(e.lm, newdata = armd.long, type = "response", se.fit = TRUE)
    test <- predict(e.lm2, newdata = armd.long, se = TRUE)
    expect_equivalent(GS$fit, test$estimate, tol = 1e-6)
    expect_equivalent(GS$se.fit, test$se, tol = 1e-6)
    expect_equivalent(GS$df, unique(round(test$df)))
})


## * predict/fitted - lmm
e.lmm <- lmm(visual ~ lesion + week + week:treat.f,
             repetition = ~ week | subject, structure = "UN",
             data = armd.long)
e.lmmNNA <- lmm(visual ~ lesion + week + week:treat.f,
             repetition = ~ week | subject, structure = "UN",
             data = armd.longNNA)

test_that("predict/fitted (lmm)", {

    ##-- static
    testthat::expect_visible(predict(e.lmm, newdata = nd)) ## check it outputs something
    expect_equal(predict(e.lmm, newdata = nd, type = "static", format = "long", se = FALSE),
                 as.double(nd.X %*% coef(e.lmm)), tol = 1e-6)
    expect_equivalent(predict(e.lmm, newdata = nd, type = "static", format = "long", se = TRUE)[,c("estimate","se")],
                      data.frame(estimate = nd.X %*% coef(e.lmm),
                                 se = sqrt(diag(nd.X %*% vcov(e.lmm) %*% t(nd.X)))),
                      tol = 1e-6)
    expect_visible(predict(e.lmm, newdata = nd, type = "static", format = "long", se = TRUE, keep.data = TRUE))
    expect_visible(predict(e.lmm, newdata = nd, type = "static", format = "wide"))
    GS <- data.frame("subject" = c("2", "10"), 
                     "lesion" = c(1, 1), 
                     "treat.f" = as.factor(c("Active", "Placebo")), 
                     "estimate_0" = c(57.47046, 58.19500), 
                     "estimate_4" = c(53.99483, 56.91638), 
                     "estimate_12" = c(51.62498, 55.84074), 
                     "estimate_24" = c(48.38935, 52.16911), 
                     "estimate_52" = c(41.27514, 46.88251))
    expect_equivalent(predict(e.lmm, newdata = nd, type = "static", format = "wide", keep.data = TRUE), GS, tol = 1e-6)

    ##-- time ordering
    test1 <- predict(e.lmm, newdata = nd[NROW(nd):1,], se = c(TRUE,TRUE))
    test2 <- predict(e.lmm, newdata = nd, se = c(TRUE,TRUE))
    expect_equivalent(test1, test2[NROW(nd):1,], tol = 1e-6)
    ## sqrt(diag(attr(predict(e.lmm, newdata = nd, se = "total2", simplify = FALSE),"vcov")))

    ##-- static NA
    expect_equivalent(predict(e.lmmNNA, newdata = nd, type = "static", format = "wide", keep.data = TRUE), GS, tol = 1e-6)
    GSNNA <- GS
    GSNNA[2,"estimate_52"] <- NA
    expect_equivalent(predict(e.lmmNNA, newdata = ndNNA[ndNNA$subject %in% GSNNA$subject,,drop=FALSE], type = "static", format = "wide", keep.data = TRUE),
                      GSNNA, tol = 1e-6)

    ##-- unique static
    Ufit <- fitted(e.lmm, newdata = "unique")
    
    ##-- dynamic
    predL <- predict(e.lmm, newdata = nd2, type = "dynamic", keep.data = TRUE)
    expect_equal(is.na(nd2$visual),!is.na(predL$estimate))
    expect_equal(c(NA, 1.5957344, NA, NA, 1.66930219, NA, NA, 1.78849944, NA, NA, 1.92606744, 1.07271763, NA, 2.00302196, 1.50857931),predL$se, tol = 1e-6)
    expect_equal(c(NA, 243.9075, NA, NA, 257.7478, NA, NA, 272.432, NA, NA, 275.4818, 174.074, NA, 263.3198, 198.7918), predL$df, tol = 1e-1)
    ## similar to not using a transformation
    ## predL2 <- predict(e.lmm, newdata = nd2, type = "dynamic", transform.sigma = "none", transform.k = "none", transform.rho = "none",
    ##                   keep.data = TRUE)
    ## expect_equal(predL2$estimate,predL$estimate, tol = 1e-6)
    ## expect_equal(predL2$se,predL$se, tol = 1e-3)
    ## expect_equal(predL2$df,predL$df, tol = 1e-1)

    predLse2 <- predict(e.lmm, newdata = nd2, type = "dynamic", se = c(TRUE,TRUE))
    predW <- predict(e.lmm, newdata = nd2, type = "dynamic", format = "wide", keep.data = TRUE)
    expect_equal(colnames(predW), c("subject", "lesion", "treat.f", "estimate_0", "estimate_4", "estimate_12", "estimate_24", "estimate_52"))
    fitL <- fitted(e.lmm, newdata = nd2, type = "outcome", format = "long", keep.data = TRUE, export.vcov = TRUE)
    ## GS <- estimate(e.lmm, function(p){ ## p <- NULL
    ##     predict(e.lmm, p = p, newdata = nd2, type = "dynamic")
    ## })
    expect_equal(fitL$se, c(NA, 1.5957344, NA, NA, 1.66930219, NA, NA, 1.78849944, NA, NA, 1.92606744, 1.07271669, NA, 2.00302196, 1.50857614), tol = 1e-6)
    expect_equal(fitL$df, c(NA, 243.784, NA, NA, 257.264, NA, NA, 271.656, NA, NA, 275.089, 137.352, NA, 262.783, 155.409), tol = 1e-1)

    fitW <- fitted(e.lmm, newdata = nd2, type = "outcome", format = "wide")
    
    ##-- impute
    imputeW <- fitted(e.lmm, newdata = nd2, type = "impute", format = "wide")
    imputeL <- fitted(e.lmm, newdata = nd2, type = "impute", format = "long")
   
    ##-- change
    changeL <- fitted(e.lmm, newdata = nd2, type = "change", format = "long", keep.data = TRUE)
    expect_equal(fitL[fitL$subject==3,"se"], changeL[changeL$subject==3,"se"])
    expect_equal(fitL[fitL$subject==3,"df"], changeL[changeL$subject==3,"df"])

    M.change <- matrix(c(-1,rep(0,4)), byrow = TRUE, nrow = 5, ncol = 5) + diag(1, 5, 5)
    GS.est <- M.change %*% fitL[fitL$subject==2,"visual"]
    GS.se <- sqrt(diag(M.change %*% attr(fitL,"vcov")[fitL$subject==2,fitL$subject==2] %*% t(M.change)))

    expect_equal(changeL[changeL$subject==2,"visual"], GS.est[,1], tol = 1e-6)
    expect_equal(changeL[changeL$subject==2,"se"], c(NA,GS.se[-1]), tol = 1e-6)
   
    changeW <- fitted(e.lmm, newdata = nd2, type = "change", format = "wide")
    expect_equal(changeW$visual_0, rep(0,3), tol = 1e-6)
    expect_equal(changeW$visual_4, fitW$visual_4-fitW$visual_0, tol = 1e-6)
    expect_equal(changeW$visual_12, fitW$visual_12-fitW$visual_0, tol = 1e-6)
   
    ## auc
    aucL <- fitted(e.lmm, newdata = nd2, type = "auc", format = "long", keep.data = TRUE)
    aucbL <- fitted(e.lmm, newdata = nd2, type = "auc-b", format = "long", keep.data = TRUE)
})



##----------------------------------------------------------------------
### test-auto-predict.R ends here