## ----include = FALSE---------------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4, warning = FALSE, message = FALSE ) ## ----install, eval = FALSE---------------------------------------------------- # install.packages("TestREnlme") ## ----load--------------------------------------------------------------------- library(TestREnlme) ## ----theoph-profiles, fig.cap = "Individual concentration profiles with mean (black and bold)."---- d <- as.data.frame(Theoph) p <- plot_profiles(d, group = "Subject", time = "Time", response = "conc", title = "Theophylline: individual profiles") ## Colour by subject using ggplot2 customisation p + ggplot2::aes(colour = .data$`.group`) + ggplot2::scale_colour_manual(values = rainbow(12), guide = "none") ## ----theoph-setup------------------------------------------------------------- Expr <- conc ~ Dose * exp(ai2 + ai3 - ai1) * (exp(-Time * exp(ai3)) - exp(-Time * exp(ai2))) / (exp(ai2) - exp(ai3)) random <- c("ai1 ~ B1 + bi1", "ai2 ~ B2 + bi2", "ai3 ~ B3 + bi3") start <- c(ai1 = -3.22, ai2 = 0.47, ai3 = -2.45) ## ----theoph-vls--------------------------------------------------------------- ## VLS (default: no subject-specific estimation) DVLS <- Dmethod(d, Expr, group = "Subject", random = random, start = start, verbose = 0) summary(DVLS) ## ----theoph-mm---------------------------------------------------------------- mb <- MM_base(d, Expr, group = "Subject", random = random, start = start, kappa_max = 1e4, verbose = 1) DMM <- Dmethod(d, Expr, group = "Subject", random = random, start = start, method = "MM", MM_base_obj = mb, verbose = 0) DMMF <- Dmethod(d, Expr, group = "Subject", random = random, start = start, method = "MMF", MM_base_obj = mb, verbose = 0) summary(DMM) summary(DMMF) ## Compare variance component estimates across methods plot_variance(list(VLS = DVLS, MM = DMM, MMF = DMMF), component = "diagonal") ## ----theoph-boot, eval = FALSE------------------------------------------------ # BVLS <- bootstrap_se(DVLS, nboot = 200, seed = 42) # print(BVLS) ## ----theoph-tests------------------------------------------------------------- ## Test 1: H0: D = 0 (all RE absent) HVLS <- Dhypothesis_test(d, Expr, group = "Subject", random = random, start = start, Dhatt = DVLS, nperm = 200, seed = 1, verbose = 0) cat("Test 1 - p-value:", HVLS$pvalue, "| Decision:", HVLS$Decision, "\n") ## Test 2: H0: d11 = d13 = d33 = 0 given bi2 retained in the model HVLS13 <- Dhypothesis_test(d, Expr, group = "Subject", random = random, start = start, bi_out = c("bi1", "bi3"), Dhatt = DVLS, nperm = 200, seed = 1, verbose = 0) cat("Test 2 - p-value:", HVLS13$pvalue, "| Decision:", HVLS13$Decision, "\n") ## Test 3: H0: d33 = 0 given bi1 and bi2 retained in the model HVLS3 <- Dhypothesis_test(d, Expr, group = "Subject", random = random, start = start, bi_out = "bi3", Dhatt = DVLS, nperm = 200, seed = 1, verbose = 0) cat("Test 3 - p-value:", HVLS3$pvalue, "| Decision:", HVLS3$Decision, "\n") ## ----theoph-hist, fig.width = 10, fig.height = 4, fig.cap = "Permutation null distributions for the three tests."---- library(ggpubr) ggarrange( plot_perm_hist(HVLS, title = "H0: D = 0"), plot_perm_hist(HVLS13, title = "H0: d11=d13=d33=0"), plot_perm_hist(HVLS3, title = "H0: d33=0"), nrow = 1 ) ## ----theoph-reduced-step1----------------------------------------------------- Expr2 <- conc ~ Dose * exp(ai2 + B3 - ai1) * (exp(-Time * exp(B3)) - exp(-Time * exp(ai2))) / (exp(ai2) - exp(B3)) random2 <- c("ai1 ~ B1 + bi1", "ai2 ~ B2 + bi2") start2 <- c(ai1 = -3.22, ai2 = 0.47, B3 = -2.45) ## ----theoph-reduced-step2----------------------------------------------------- DVLS_C <- Dmethod(d, Expr2, group = "Subject", random = random2, start = start2, verbose = 0) summary(DVLS_C) ## ----theoph-fitted, fig.cap = "Observed vs fitted values for all 12 subjects (reduced model)."---- plot_fitted(DVLS_C, time = "Time", subjects = NULL, overlay = TRUE) ## ----s3-dmethod, eval = FALSE------------------------------------------------- # ## Shows: profiles, fitted, residuals (raw + standardised) # ## For MM/MMF also shows condition-number plot # plot(x=DVLS_C, time = "Time") ## ----theoph-cond, fig.cap = "Condition numbers for the reduced MM model. All subjects are below the exclusion threshold."---- DMM_C <- Dmethod(d, Expr2, group = "Subject", random = random2, start = start2, method = "MM", verbose = 0) plot_condition(DMM_C, kappa_max = 1e4) ## ----theoph-reduced-test------------------------------------------------------ ## Test for the need of all random effects in the reduced model HC <- Dhypothesis_test(d, Expr2, group = "Subject", random = random2, start = start2, Dhatt = DVLS_C, nperm = 200, seed = 1, verbose = 0) cat("Testing of all RE (bi1, bi2) - p-value:", HC$pvalue, "| Decision:", HC$Decision, "\n") # ## Test for the need of bi1 given the presence of bi2 in the reduced model HC1 <- Dhypothesis_test(d, Expr2, group = "Subject", random = random2, start = start2,bi_out = "bi1", Dhatt = DVLS_C, nperm = 200, seed = 1, verbose = 0) cat("Testing bi1 give bi2 in the model- p-value:", HC1$pvalue, "| Decision:", HC1$Decision, "\n") ## ----skill-setup-------------------------------------------------------------- ## Reshape from wide to long format library(reshape2) data("SkillAcq") qrt1 <- melt(SkillAcq, id.vars = c("id", "wm2"), variable.name = "ly", value.name = "Y") qrt1$t <- as.numeric(sub("ly", "", qrt1$ly)) #model formulation Expr_learn <- Y ~ ai1 - (ai1 + ai0) * exp(ai2 * t) random_learn <- c("ai0 ~ B0 + B01 * wm2 + bi0", "ai1 ~ B1 + B11 * wm2 + bi1", "ai2 ~ B2 + B21 * wm2 + bi2") ## ----skill-full--------------------------------------------------------------- ## MM estimator DMM_all <- Dmethod(data=qrt1, Expr=Expr_learn, group = "id", random = random_learn, method = "MM") summary(DMM_all) ## ----skill-full-plot , fig.cap="Skill acquisition data: the per-subject condition numbers from MM for the full three-random-effects model, sorted in decreasing order for converged subjects. The subjects above the dashed line (red dots) were excluded from MM second-stage estimation, and the subjects below (blue dots) were retained and used in the analysis."---- ## Condition number plot -- subjects above threshold excluded set.seed(123) NS <- sort(sample(1:204, 20)) plot_condition(DMM_all, kappa_max = 1e4)