tests/testthat/test-auto-practical.R
In bozenne/repeated: Repeated Measurement Models for Discrete Times
### test-auto-practical.R ---
##----------------------------------------------------------------------
## Author: Brice Ozenne
## Created: Jun 7 2021 (17:03)
## Version:
## Last-Updated: maj 7 2024 (15:18)
## By: Brice Ozenne
## Update #: 131
##----------------------------------------------------------------------
##
### Commentary:
##
### Change Log:
##----------------------------------------------------------------------
##
### Code:
if(FALSE){
library(testthat)
library(numDeriv)
library(lava)
library(multcomp)
library(nlme)
library(qqtest)
library(ggplot2)
library(lme4)
library(lmerTest)
library(LMMstar)
}
context("Check lmm on the example from the course with Julie")
LMMstar.options(optimizer = "FS", method.numDeriv = "simple", precompute.moments = TRUE,
columns.confint = c("estimate","se","df","lower","upper","p.value"))
## * Practical 1
test_that("practical 1 - gastricbypass",{
## ** data
data(gastricbypassL, package = "LMMstar")
## ** summarize
summarize(glucagonAUC~time, data = gastricbypassL, na.rm=TRUE)
summarize(glucagonAUC~time|id, data = gastricbypassL, na.rm=TRUE)
## ** compound symmetry
eCS.gls <- gls(glucagonAUC~visit,
data=gastricbypassL,
correlation=corCompSymm(form=~visit|id),
na.action=na.exclude,
control=glsControl(opt="optim"))
eCS.lmm <- lmm(glucagonAUC~visit,
data=gastricbypassL,
repetition = ~visit|id,
structure = "CS")
expect_equal(as.double(logLik(eCS.gls)), as.double(logLik(eCS.lmm)), tol = 1e-6)
## ** unstructured with missing data
eUN.gls <- gls(glucagonAUC~visit,
data=gastricbypassL,
correlation=corSymm(form=~as.numeric(visit)|id),
weights=varIdent(form=~1|time),
na.action=na.exclude,
control=glsControl(opt="optim"),
method = "REML")
eUN.lmm <- lmm(glucagonAUC~visit,
data=gastricbypassL,
repetition = ~visit|id,
structure = "UN",
method.fit = "REML", trace = 0)
## check moments
expect_equal(as.double(logLik(eUN.gls)), as.double(logLik(eUN.lmm)), tol = 1e-6)
## GS <- numDeriv::jacobian(func = function(p){logLik(eUN.lmm, p = p, transform.sigma = "none", transform.k = "none", transform.rho = "none")}, x = coef(eUN.lmm, transform.sigma = "none", transform.k = "none", transform.rho = "none", effects = "all"))
GS <- rbind(c("(Intercept)" = 0, "visit2" = 0, "visit3" = 0, "visit4" = 0, "sigma" = 2e-08, "k.2" = 7.85e-06, "k.3" = 5.37e-06, "k.4" = -2.109e-05, "rho(1,2)" = -0.00011073, "rho(1,3)" = -7.107e-05, "rho(1,4)" = 9.35e-06, "rho(2,3)" = 7.852e-05, "rho(2,4)" = -2.407e-05, "rho(3,4)" = 4.945e-05))
expect_equal(as.double(GS), as.double(score(eUN.lmm, transform.sigma = "none", transform.k = "none", transform.rho = "none", effects = "all")), tol = 1e-6)
## GS <- -numDeriv::jacobian(func = function(p){score(eUN.lmm, p = p, transform.sigma = "none", transform.k = "none", transform.rho = "none",effects = "all")}, x = coef(eUN.lmm, transform.sigma = "none", transform.k = "none", transform.rho = "none", effects = "all"))
GS <- matrix(c(0.2000435, 0.10068008, -0.00122474, 0.0667917, 0, 0.01960874, -0.00398498, 0.00220684, 0.015223, 0.02694338, -0.00181844, -0.02394494, -0.00163313, 0.00522044, 0.10068008, 0.39828548, 0.05620903, -0.04323161, 0, 0.06528472, -0.03616399, 0.07545552, 0.34097272, 0.31291727, -0.03538287, -0.31640674, 0.0569075, -0.1180618, -0.00122474, 0.05620903, 0.05070213, -0.04276096, 0, 0.00060657, -0.01939405, 0.0573912, 0.24479209, 0.18869963, -0.02474055, -0.19994584, 0.05242174, -0.1138332, 0.0667917, -0.04323161, -0.04276096, 0.08432189, 0, 0.00029111, 0.01617431, -0.04825539, -0.20558724, -0.15785621, 0.02076676, 0.16745253, -0.04422625, 0.09610526, 0, 0, 0, 0, 0.61232751, 2.68368828, 1.30246228, 1.80652473, -7.88842064, -3.41692146, 2.72607489, 2.53245722, -1.48188657, -3.02859195, 0.01960874, 0.06528472, 0.00060657, 0.00029111, 2.68368828, 119.6435905, -1.30957174, 2.67549682, -67.36667611, 2.7100885, -1.4963409, 22.4424789, -13.18941762, -0.05361995, -0.00398498, -0.03616399, -0.01939405, 0.01617431, 1.30246228, -1.30957174, 17.62474246, -5.70108564, 2.50057785, -13.49661757, -1.26046703, 9.8287541, 1.16522916, -13.33546581, 0.00220684, 0.07545552, 0.0573912, -0.04825539, 1.80652473, 2.67549682, -5.70108564, 29.08095721, -4.89725928, -3.62086875, 17.33682661, 3.88337218, -10.663066, -16.12586692, 0.015223, 0.34097272, 0.24479209, -0.20558724, -7.88842064, -67.36667611, 2.50057785, -4.89725928, 494.48230881, 170.86268327, -126.83334603, -160.48986839, 120.86002165, 52.61300133, 0.02694338, 0.31291727, 0.18869963, -0.15785621, -3.41692146, 2.7100885, -13.49661757, -3.62086875, 170.86268327, 195.7773562, -137.85293967, -150.42933562, 105.91037437, 88.19350765, -0.00181844, -0.03538287, -0.02474055, 0.02076676, 2.72607489, -1.4963409, -1.26046703, 17.33682661, -126.83334603, -137.85293967, 182.19196649, 101.08267288, -138.37463781, -84.00686145, -0.02394494, -0.31640674, -0.19994584, 0.16745253, 2.53245722, 22.4424789, 9.8287541, 3.88337218, -160.48986839, -150.42933562, 101.08267288, 155.69993421, -103.48124038, -60.02043103, -0.00163313, 0.0569075, 0.05242174, -0.04422625, -1.48188657, -13.18941762, 1.16522916, -10.663066, 120.86002165, 105.91037437, -138.37463781, -103.48124038, 144.14781848, 59.2350391, 0.00522044, -0.1180618, -0.1138332, 0.09610526, -3.02859195, -0.05361995, -13.33546581, -16.12586692, 52.61300133, 88.19350765, -84.00686145, -60.02043103, 59.2350391, 99.27770758),
nrow = 14,
ncol = 14,
dimnames = list(c("(Intercept)", "visit2", "visit3", "visit4", "sigma", "k.2", "k.3", "k.4", "rho(1,2)", "rho(1,3)", "rho(1,4)", "rho(2,3)", "rho(2,4)", "rho(3,4)"),c("(Intercept)", "visit2", "visit3", "visit4", "sigma", "k.2", "k.3", "k.4", "rho(1,2)", "rho(1,3)", "rho(1,4)", "rho(2,3)", "rho(2,4)", "rho(3,4)"))
)
expect_equal(as.double(GS), as.double(information(eUN.lmm, transform.sigma = "none", transform.k = "none", transform.rho = "none", effects = "all")), tol = 1e-6)
## ** extract information
capture.output(model.tables(eUN.lmm))
capture.output(summary(eUN.lmm, columns = "statistic"))
capture.output(confint(eUN.lmm, effects = "all", backtransform = TRUE))
expect_equal(unname(coef(eUN.lmm, effects = "variance", transform.k = "sd")),
unname(coef(eUN.lmm, effects = "variance")[1]*c(1,coef(eUN.lmm, effects = "variance")[-1])),
tol = 1e-6)
autoplot(eUN.lmm)
apply(residuals(eUN.lmm, type = "normalized", format = "wide")[,-1],2,sd,na.rm=TRUE)
cor(residuals(eUN.lmm, type = "normalized", format = "wide")[,-1], use = "pairwise")
qqtest::qqtest(na.omit(residuals(eUN.lmm, type = "normalized")))
eUN.lmm_anova <- anova(eUN.lmm, effects = c("visit3-visit2=0"), ci = TRUE)
capture.output(eUN.lmm_anova)
capture.output(summary(eUN.lmm_anova))
expect_equal(eUN.lmm_anova$multivariate$df.denom, 19.20982, tol = 1e-1) ## Richardson
})
## * Practical 2
test_that("practical 2 - ncgs",{
data(ncgsL, package = "LMMstar")
ncgsL$visit <- as.numeric(ncgsL$visit)
ncgsL$highdose.time <- ncgsL$time
ncgsL$highdose.time[ncgsL$group=="placebo"] <- "0"
ncgsL$time <- relevel(as.factor(ncgsL$time), ref="0")
ncgsL$highdose.time <- relevel(as.factor(ncgsL$highdose.time), ref="0")
ncgsL$treatment <- factor(ncgsL$group, c("none","placebo","highdose"))
ncgsL$treatment[ncgsL$time=="0"] <- "none"
ncgsL$treatment2 <- factor(ncgsL$group, c("placebo","highdose"))
ncgsL$treatment2[ncgsL$time=="0"] <- "placebo"
## ** unstructured with missing data
e.gls <- gls(cholest~time+highdose.time,
data=ncgsL,
correlation=corSymm(form=~visit|id),
weights=varIdent(form=~1|time),
na.action=na.exclude,
control=glsControl(opt='optim'))
e.lmm <- lmm(cholest~time+highdose.time,
data=ncgsL,
repetition = ~time|id,
structure = "UN",
df = FALSE)
e2.lmm <- suppressMessages(lmm(cholest~treatment*time,
data=ncgsL,
repetition = ~time|id,
structure = "UN"))
e3.lmm <- suppressMessages(lmm(cholest~treatment2*time,
data=ncgsL,
repetition = ~time|id,
structure = "UN"))
expect_equal(as.double(logLik(e.gls)), as.double(logLik(e.lmm)), tol = 1e-6)
expect_equal(as.double(logLik(e.gls)), as.double(logLik(e2.lmm)), tol = 1e-6)
expect_equal(as.double(logLik(e.gls)), as.double(logLik(e3.lmm)), tol = 1e-6)
## ** extract information
confint(e2.lmm, effects = "all", backtransform = TRUE)[,c("estimate","lower","upper")]
e2.lmm_anova <- anova(e2.lmm, effects = c("treatmenthighdose:time6-treatmentplacebo:time6=0","treatmenthighdose:time12-treatmentplacebo:time12=0"), ci = TRUE)
expect_equal(e2.lmm_anova$multivariate$df.denom, 100.0411, tol = 1e-1) ## Richardson
autoplot(e2.lmm)
autoplot(e2.lmm, color = "group")
autoplot(e2.lmm, color = "group", ci.alpha = 0.25)
})
test_that("practical 2 - vitamin",{
data(vitaminL, package = "LMMstar")
vitaminL$visit.num <- as.numeric(vitaminL$visit)
## define treatment
vitaminL$treatment <- factor(vitaminL$group, c("none","control","vitamin"))
vitaminL$treatment[vitaminL$time<=4] <- "none"
## define interaction manually
vitaminL$vita.time <- vitaminL$time
vitaminL$vita.time[vitaminL$group=="control" | vitaminL$time<=4] <- "1"
## ** fit unstructured
e.gls <- gls(weight~visit+vita.time,
data=vitaminL,
correlation=corSymm(form=~visit.num|animal),
weights=varIdent(form=~1|visit),
na.action=na.exclude,
control=glsControl(opt='optim'))
e0.lmm <- suppressWarnings(lmm(weight~visit+vita.time,
data=vitaminL,
repetition = ~visit|animal,
structure = "UN"))
e.lmm <- suppressMessages(suppressWarnings(lmm(weight~treatment*visit,
data=vitaminL,
repetition = ~visit|animal,
structure = "UN")))
expect_equal(as.double(logLik(e0.lmm)), as.double(logLik(e.lmm)), tol = 1e-6)
expect_equal(as.double(logLik(e.lmm)), -234.28331772, tol = 1e-6)
expect_equal(as.double(logLik(e.gls)), -232.08183621, tol = 1e-6)
## ** extract information
autoplot(e.lmm, color = "group",ci =FALSE)
qqtest::qqtest(residuals(e.lmm))
confint(e.lmm)
e.lmm_anova <- anova(e.lmm, effects = "treatmentvitamin:visit6 - treatmentcontrol:visit6 = 0", ci = TRUE)
## expect_equal(e.lmm_anova$all$df.denom, 0.4972829, tol = 1e-1) ## Richardson
})
## * Practical 3
test_that("practical 3 - swabsL",{
data(swabsL, package = "LMMstar")
## ** unstructured
eUN.gls <- gls(swabs ~ crowding + name,
data = swabsL,
correlation = corSymm(form = ~ as.numeric(name) | family),
weights = varIdent(form = ~ 1 | name)
)
eUN.lmm <- lmm(swabs ~ crowding + name, data = swabsL, structure = "UN", repetition = ~name|family)
expect_equal(as.double(logLik(eUN.lmm)),as.double(logLik(eUN.gls)), tol = 1e-6)
## summary(eUN.lmm)
sigma(eUN.lmm)
## ** compound symmetry
eCS.lme <- lme(swabs ~ crowding + name,
data = swabsL,
random =~ 1 | family)
eCS.lmm <- lmm(swabs ~ crowding + name, data = swabsL, structure = "CS", repetition = ~name|family)
expect_equal(as.double(logLik(eCS.lmm)), as.double(logLik(eCS.lme)), tol = 1e-6)
capture.output(summary(eCS.lmm))
sigma(eCS.lmm)
autoplot(eCS.lmm)
## with interaction
eCSI.lmm <- lmm(swabs ~ crowding * name, data = swabsL, structure = "CS", repetition = ~name|family)
eRI.lmm <- lmm(swabs ~ crowding * name + (1|family), data = swabsL)
eCSI.lme <- lme(swabs ~ crowding * name, random =~ 1 | family, data = swabsL)
GS <- as.data.frame(ranef(eCSI.lme, augFrame = TRUE))
expect_equal(as.double(GS[,1]),as.double(ranef(eRI.lmm)), tol = 1e-4)
})
## * Practical 4
test_that("practical 4 - bloodpressureL",{
data(bloodpressureL, package = "LMMstar")
bloodpressureL$period.num <- as.numeric(bloodpressureL$period)
bloodpressureL$treatment.num <- as.numeric(bloodpressureL$treatment)
## ** compound symmetry
eCS.gls <- gls(duration ~ period + treatment,
data = bloodpressureL,
correlation = corCompSymm(form=~ 1 | id))
eCS.lmm <- lmm(duration ~ period + treatment,
data = bloodpressureL,
structure = "CS", repetition = ~ period | id)
expect_equal(as.double(logLik(eCS.lmm)),as.double(logLik(eCS.gls)), tol = 1e-6)
plot(eCS.lmm, color = "sequence", obs.alpha = 0.25, ci.alpha = 0.05)
confint(eCS.gls)
confint(eCS.lmm)
## ** unstructured
eUNP.gls <- gls(duration ~ period + treatment,
data = bloodpressureL,
correlation = corSymm(form=~ period.num | id),
weights = varIdent(form=~ 1|period),
)
eUNP.lmm <- lmm(duration ~ period + treatment,
data = bloodpressureL,
structure = "UN", repetition = ~ period | id)
expect_equal(as.double(logLik(eUNP.lmm)),as.double(logLik(eUNP.gls)), tol = 1e-6)
capture.output(summary(eUNP.lmm))
eUNT.gls <- gls(duration ~ period + treatment,
data = bloodpressureL,
correlation = corSymm(form=~ treatment.num | id),
weights = varIdent(form=~ 1|treatment),
)
eUNT.lmm <- lmm(duration ~ period + treatment,
data = bloodpressureL,
structure = "UN", repetition = ~ treatment | id)
capture.output(summary(eUNT.lmm))
expect_equal(as.double(logLik(eUNT.lmm)),as.double(logLik(eUNT.gls)), tol = 1e-6)
})
## * Practical 6
test_that("practical 6 - vasscoresL",{
data(vasscoresL, package = "LMMstar")
summarize(vas ~ treatment, data = vasscoresL, na.rm = TRUE)
summarize(vas ~ treatment|id, data = vasscoresL, na.rm = TRUE)
summarize(vas ~ treatment+group|id, data = vasscoresL, na.rm = TRUE)
## ** model on all pairs: AB, AC, BC
fit.CS <- lmm(vas~-1+treatment, data=vasscoresL,
repetition=~treatment|id, structure="CS")
capture.output(summary(fit.CS))
GS <- matrix(c(2099.98852585, 1675.79243303, 1675.79243303, 1675.79243303, 2099.98852585, 1675.79243303, 1675.79243303, 1675.79243303, 2099.98852585),
nrow = 3,
ncol = 3,
dimnames = list(c("A", "B", "C"),c("A", "B", "C"))
)
expect_equivalent(sigma(fit.CS), GS, tol = 1e-5)
## GS <- lmer(vas~-1+treatment+(1|id), data=vasscoresL)
GS <- lme(vas~-1+treatment, random =~ 1|id, data=vasscoresL, na.action = na.omit)
expect_equivalent(unname(coef(fit.CS)), unname(fixef(GS)), tol = 1e-5)
## expect_equivalent(unname(model.tables(fit.CS)$se), unname(summary(GS)$coef[,"Std. Error"]), tol = 1e-2)
## expect_equivalent(unname(model.tables(fit.CS)$df), unname(summary(GS)$coef[,"df"]), tol = 1e-2)
expect_equivalent(unname(model.tables(fit.CS)$se), unname(summary(GS)$tTable[,"Std.Error"]), tol = 1e-2)
expect_equivalent(unname(model.tables(fit.CS)$df), c(43.64236203, 43.64236203, 43.64236203), tol = 1e-2)
## autoplot(fit.CS)
suppressWarnings(plot(fit.CS, obs.alpha = 0.1))
fit.UN <- lmm(vas~-1+treatment, data=vasscoresL,
repetition=~treatment|id, structure="UN")
capture.output(summary(fit.UN))
## GS <- lmm(vas~-1+treatment, data=vasscoresL,
## repetition=~treatment|id, structure="UN",
## control = list(optimizer = "FS"))
GS <- matrix(c(2038.74222977, 1892.18702336, 1797.44635389, 1892.18702336, 2111.43949953, 1369.79020341, 1797.44635389, 1369.79020341, 2157.87873809),
nrow = 3,
ncol = 3,
dimnames = list(c("A", "B", "C"),c("A", "B", "C"))
)
expect_equivalent(sigma(fit.UN), GS, tol = 1e-5)
vasscoreL.imputed <- fitted(fit.UN, type = "outcome", keep.data = TRUE)
set.seed(11)
gg <- ggplot(vasscoreL.imputed, aes(x = treatment, y = vas, group = id))
gg + geom_point(aes(color = impute)) + geom_line()
## ** model some pairs: AB, AC
## IMPORTANT: check extraction of the residual variance-covariance matrix when not all pairs of time are analyzed
vasscoreL.red <- vasscoresL[vasscoresL$group %in% c("AB","AC"),]
fit.CS.red <- lmm(vas~-1+treatment, data=vasscoreL.red,
repetition=~treatment|id, structure="CS")
capture.output(summary(fit.CS.red))
GS <- matrix(c(2006.9354323, 1763.57527658, 1763.57527658, 1763.57527658, 2006.9354323, 1763.57527658, 1763.57527658, 1763.57527658, 2006.9354323),
nrow = 3,
ncol = 3,
dimnames = list(c("A", "B", "C"),c("A", "B", "C"))
)
expect_equivalent(sigma(fit.CS.red), GS, tol = 1e-5)
fit.UN.red <- lmm(vas~-1+treatment, data=vasscoreL.red,
repetition=~treatment|id, structure="UN", control = list(optimizer = "FS")) ## FS optimizer mandatory, GLS does something weird as it is able to estimate the correlation (B,C)
capture.output(summary(fit.UN.red))
GS <- matrix(c(1996.68088251, 1926.29216631, 1573.59329163, 1926.29216631, 2231.4711806, NA, 1573.59329163, NA, 1753.88563313),
nrow = 3,
ncol = 3,
dimnames = list(c("A", "B", "C"),c("A", "B", "C"))
)
expect_equivalent(sigma(fit.UN.red), GS, tol = 1e-5)
})
##----------------------------------------------------------------------
### test-auto-practical.R ends here
bozenne/repeated documentation built on July 16, 2025, 11:16 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
### test-auto-practical.R ---
##----------------------------------------------------------------------
## Author: Brice Ozenne
## Created: Jun 7 2021 (17:03)
## Version:
## Last-Updated: maj 7 2024 (15:18)
## By: Brice Ozenne
## Update #: 131
##----------------------------------------------------------------------
##
### Commentary:
##
### Change Log:
##----------------------------------------------------------------------
##
### Code:
if(FALSE){
library(testthat)
library(numDeriv)
library(lava)
library(multcomp)
library(nlme)
library(qqtest)
library(ggplot2)
library(lme4)
library(lmerTest)
library(LMMstar)
}
context("Check lmm on the example from the course with Julie")
LMMstar.options(optimizer = "FS", method.numDeriv = "simple", precompute.moments = TRUE,
columns.confint = c("estimate","se","df","lower","upper","p.value"))
## * Practical 1
test_that("practical 1 - gastricbypass",{
## ** data
data(gastricbypassL, package = "LMMstar")
## ** summarize
summarize(glucagonAUC~time, data = gastricbypassL, na.rm=TRUE)
summarize(glucagonAUC~time|id, data = gastricbypassL, na.rm=TRUE)
## ** compound symmetry
eCS.gls <- gls(glucagonAUC~visit,
data=gastricbypassL,
correlation=corCompSymm(form=~visit|id),
na.action=na.exclude,
control=glsControl(opt="optim"))
eCS.lmm <- lmm(glucagonAUC~visit,
data=gastricbypassL,
repetition = ~visit|id,
structure = "CS")
expect_equal(as.double(logLik(eCS.gls)), as.double(logLik(eCS.lmm)), tol = 1e-6)
## ** unstructured with missing data
eUN.gls <- gls(glucagonAUC~visit,
data=gastricbypassL,
correlation=corSymm(form=~as.numeric(visit)|id),
weights=varIdent(form=~1|time),
na.action=na.exclude,
control=glsControl(opt="optim"),
method = "REML")
eUN.lmm <- lmm(glucagonAUC~visit,
data=gastricbypassL,
repetition = ~visit|id,
structure = "UN",
method.fit = "REML", trace = 0)
## check moments
expect_equal(as.double(logLik(eUN.gls)), as.double(logLik(eUN.lmm)), tol = 1e-6)
## GS <- numDeriv::jacobian(func = function(p){logLik(eUN.lmm, p = p, transform.sigma = "none", transform.k = "none", transform.rho = "none")}, x = coef(eUN.lmm, transform.sigma = "none", transform.k = "none", transform.rho = "none", effects = "all"))
GS <- rbind(c("(Intercept)" = 0, "visit2" = 0, "visit3" = 0, "visit4" = 0, "sigma" = 2e-08, "k.2" = 7.85e-06, "k.3" = 5.37e-06, "k.4" = -2.109e-05, "rho(1,2)" = -0.00011073, "rho(1,3)" = -7.107e-05, "rho(1,4)" = 9.35e-06, "rho(2,3)" = 7.852e-05, "rho(2,4)" = -2.407e-05, "rho(3,4)" = 4.945e-05))
expect_equal(as.double(GS), as.double(score(eUN.lmm, transform.sigma = "none", transform.k = "none", transform.rho = "none", effects = "all")), tol = 1e-6)
## GS <- -numDeriv::jacobian(func = function(p){score(eUN.lmm, p = p, transform.sigma = "none", transform.k = "none", transform.rho = "none",effects = "all")}, x = coef(eUN.lmm, transform.sigma = "none", transform.k = "none", transform.rho = "none", effects = "all"))
GS <- matrix(c(0.2000435, 0.10068008, -0.00122474, 0.0667917, 0, 0.01960874, -0.00398498, 0.00220684, 0.015223, 0.02694338, -0.00181844, -0.02394494, -0.00163313, 0.00522044, 0.10068008, 0.39828548, 0.05620903, -0.04323161, 0, 0.06528472, -0.03616399, 0.07545552, 0.34097272, 0.31291727, -0.03538287, -0.31640674, 0.0569075, -0.1180618, -0.00122474, 0.05620903, 0.05070213, -0.04276096, 0, 0.00060657, -0.01939405, 0.0573912, 0.24479209, 0.18869963, -0.02474055, -0.19994584, 0.05242174, -0.1138332, 0.0667917, -0.04323161, -0.04276096, 0.08432189, 0, 0.00029111, 0.01617431, -0.04825539, -0.20558724, -0.15785621, 0.02076676, 0.16745253, -0.04422625, 0.09610526, 0, 0, 0, 0, 0.61232751, 2.68368828, 1.30246228, 1.80652473, -7.88842064, -3.41692146, 2.72607489, 2.53245722, -1.48188657, -3.02859195, 0.01960874, 0.06528472, 0.00060657, 0.00029111, 2.68368828, 119.6435905, -1.30957174, 2.67549682, -67.36667611, 2.7100885, -1.4963409, 22.4424789, -13.18941762, -0.05361995, -0.00398498, -0.03616399, -0.01939405, 0.01617431, 1.30246228, -1.30957174, 17.62474246, -5.70108564, 2.50057785, -13.49661757, -1.26046703, 9.8287541, 1.16522916, -13.33546581, 0.00220684, 0.07545552, 0.0573912, -0.04825539, 1.80652473, 2.67549682, -5.70108564, 29.08095721, -4.89725928, -3.62086875, 17.33682661, 3.88337218, -10.663066, -16.12586692, 0.015223, 0.34097272, 0.24479209, -0.20558724, -7.88842064, -67.36667611, 2.50057785, -4.89725928, 494.48230881, 170.86268327, -126.83334603, -160.48986839, 120.86002165, 52.61300133, 0.02694338, 0.31291727, 0.18869963, -0.15785621, -3.41692146, 2.7100885, -13.49661757, -3.62086875, 170.86268327, 195.7773562, -137.85293967, -150.42933562, 105.91037437, 88.19350765, -0.00181844, -0.03538287, -0.02474055, 0.02076676, 2.72607489, -1.4963409, -1.26046703, 17.33682661, -126.83334603, -137.85293967, 182.19196649, 101.08267288, -138.37463781, -84.00686145, -0.02394494, -0.31640674, -0.19994584, 0.16745253, 2.53245722, 22.4424789, 9.8287541, 3.88337218, -160.48986839, -150.42933562, 101.08267288, 155.69993421, -103.48124038, -60.02043103, -0.00163313, 0.0569075, 0.05242174, -0.04422625, -1.48188657, -13.18941762, 1.16522916, -10.663066, 120.86002165, 105.91037437, -138.37463781, -103.48124038, 144.14781848, 59.2350391, 0.00522044, -0.1180618, -0.1138332, 0.09610526, -3.02859195, -0.05361995, -13.33546581, -16.12586692, 52.61300133, 88.19350765, -84.00686145, -60.02043103, 59.2350391, 99.27770758),
nrow = 14,
ncol = 14,
dimnames = list(c("(Intercept)", "visit2", "visit3", "visit4", "sigma", "k.2", "k.3", "k.4", "rho(1,2)", "rho(1,3)", "rho(1,4)", "rho(2,3)", "rho(2,4)", "rho(3,4)"),c("(Intercept)", "visit2", "visit3", "visit4", "sigma", "k.2", "k.3", "k.4", "rho(1,2)", "rho(1,3)", "rho(1,4)", "rho(2,3)", "rho(2,4)", "rho(3,4)"))
)
expect_equal(as.double(GS), as.double(information(eUN.lmm, transform.sigma = "none", transform.k = "none", transform.rho = "none", effects = "all")), tol = 1e-6)
## ** extract information
capture.output(model.tables(eUN.lmm))
capture.output(summary(eUN.lmm, columns = "statistic"))
capture.output(confint(eUN.lmm, effects = "all", backtransform = TRUE))
expect_equal(unname(coef(eUN.lmm, effects = "variance", transform.k = "sd")),
unname(coef(eUN.lmm, effects = "variance")[1]*c(1,coef(eUN.lmm, effects = "variance")[-1])),
tol = 1e-6)
autoplot(eUN.lmm)
apply(residuals(eUN.lmm, type = "normalized", format = "wide")[,-1],2,sd,na.rm=TRUE)
cor(residuals(eUN.lmm, type = "normalized", format = "wide")[,-1], use = "pairwise")
qqtest::qqtest(na.omit(residuals(eUN.lmm, type = "normalized")))
eUN.lmm_anova <- anova(eUN.lmm, effects = c("visit3-visit2=0"), ci = TRUE)
capture.output(eUN.lmm_anova)
capture.output(summary(eUN.lmm_anova))
expect_equal(eUN.lmm_anova$multivariate$df.denom, 19.20982, tol = 1e-1) ## Richardson
})
## * Practical 2
test_that("practical 2 - ncgs",{
data(ncgsL, package = "LMMstar")
ncgsL$visit <- as.numeric(ncgsL$visit)
ncgsL$highdose.time <- ncgsL$time
ncgsL$highdose.time[ncgsL$group=="placebo"] <- "0"
ncgsL$time <- relevel(as.factor(ncgsL$time), ref="0")
ncgsL$highdose.time <- relevel(as.factor(ncgsL$highdose.time), ref="0")
ncgsL$treatment <- factor(ncgsL$group, c("none","placebo","highdose"))
ncgsL$treatment[ncgsL$time=="0"] <- "none"
ncgsL$treatment2 <- factor(ncgsL$group, c("placebo","highdose"))
ncgsL$treatment2[ncgsL$time=="0"] <- "placebo"
## ** unstructured with missing data
e.gls <- gls(cholest~time+highdose.time,
data=ncgsL,
correlation=corSymm(form=~visit|id),
weights=varIdent(form=~1|time),
na.action=na.exclude,
control=glsControl(opt='optim'))
e.lmm <- lmm(cholest~time+highdose.time,
data=ncgsL,
repetition = ~time|id,
structure = "UN",
df = FALSE)
e2.lmm <- suppressMessages(lmm(cholest~treatment*time,
data=ncgsL,
repetition = ~time|id,
structure = "UN"))
e3.lmm <- suppressMessages(lmm(cholest~treatment2*time,
data=ncgsL,
repetition = ~time|id,
structure = "UN"))
expect_equal(as.double(logLik(e.gls)), as.double(logLik(e.lmm)), tol = 1e-6)
expect_equal(as.double(logLik(e.gls)), as.double(logLik(e2.lmm)), tol = 1e-6)
expect_equal(as.double(logLik(e.gls)), as.double(logLik(e3.lmm)), tol = 1e-6)
## ** extract information
confint(e2.lmm, effects = "all", backtransform = TRUE)[,c("estimate","lower","upper")]
e2.lmm_anova <- anova(e2.lmm, effects = c("treatmenthighdose:time6-treatmentplacebo:time6=0","treatmenthighdose:time12-treatmentplacebo:time12=0"), ci = TRUE)
expect_equal(e2.lmm_anova$multivariate$df.denom, 100.0411, tol = 1e-1) ## Richardson
autoplot(e2.lmm)
autoplot(e2.lmm, color = "group")
autoplot(e2.lmm, color = "group", ci.alpha = 0.25)
})
test_that("practical 2 - vitamin",{
data(vitaminL, package = "LMMstar")
vitaminL$visit.num <- as.numeric(vitaminL$visit)
## define treatment
vitaminL$treatment <- factor(vitaminL$group, c("none","control","vitamin"))
vitaminL$treatment[vitaminL$time<=4] <- "none"
## define interaction manually
vitaminL$vita.time <- vitaminL$time
vitaminL$vita.time[vitaminL$group=="control" | vitaminL$time<=4] <- "1"
## ** fit unstructured
e.gls <- gls(weight~visit+vita.time,
data=vitaminL,
correlation=corSymm(form=~visit.num|animal),
weights=varIdent(form=~1|visit),
na.action=na.exclude,
control=glsControl(opt='optim'))
e0.lmm <- suppressWarnings(lmm(weight~visit+vita.time,
data=vitaminL,
repetition = ~visit|animal,
structure = "UN"))
e.lmm <- suppressMessages(suppressWarnings(lmm(weight~treatment*visit,
data=vitaminL,
repetition = ~visit|animal,
structure = "UN")))
expect_equal(as.double(logLik(e0.lmm)), as.double(logLik(e.lmm)), tol = 1e-6)
expect_equal(as.double(logLik(e.lmm)), -234.28331772, tol = 1e-6)
expect_equal(as.double(logLik(e.gls)), -232.08183621, tol = 1e-6)
## ** extract information
autoplot(e.lmm, color = "group",ci =FALSE)
qqtest::qqtest(residuals(e.lmm))
confint(e.lmm)
e.lmm_anova <- anova(e.lmm, effects = "treatmentvitamin:visit6 - treatmentcontrol:visit6 = 0", ci = TRUE)
## expect_equal(e.lmm_anova$all$df.denom, 0.4972829, tol = 1e-1) ## Richardson
})
## * Practical 3
test_that("practical 3 - swabsL",{
data(swabsL, package = "LMMstar")
## ** unstructured
eUN.gls <- gls(swabs ~ crowding + name,
data = swabsL,
correlation = corSymm(form = ~ as.numeric(name) | family),
weights = varIdent(form = ~ 1 | name)
)
eUN.lmm <- lmm(swabs ~ crowding + name, data = swabsL, structure = "UN", repetition = ~name|family)
expect_equal(as.double(logLik(eUN.lmm)),as.double(logLik(eUN.gls)), tol = 1e-6)
## summary(eUN.lmm)
sigma(eUN.lmm)
## ** compound symmetry
eCS.lme <- lme(swabs ~ crowding + name,
data = swabsL,
random =~ 1 | family)
eCS.lmm <- lmm(swabs ~ crowding + name, data = swabsL, structure = "CS", repetition = ~name|family)
expect_equal(as.double(logLik(eCS.lmm)), as.double(logLik(eCS.lme)), tol = 1e-6)
capture.output(summary(eCS.lmm))
sigma(eCS.lmm)
autoplot(eCS.lmm)
## with interaction
eCSI.lmm <- lmm(swabs ~ crowding * name, data = swabsL, structure = "CS", repetition = ~name|family)
eRI.lmm <- lmm(swabs ~ crowding * name + (1|family), data = swabsL)
eCSI.lme <- lme(swabs ~ crowding * name, random =~ 1 | family, data = swabsL)
GS <- as.data.frame(ranef(eCSI.lme, augFrame = TRUE))
expect_equal(as.double(GS[,1]),as.double(ranef(eRI.lmm)), tol = 1e-4)
})
## * Practical 4
test_that("practical 4 - bloodpressureL",{
data(bloodpressureL, package = "LMMstar")
bloodpressureL$period.num <- as.numeric(bloodpressureL$period)
bloodpressureL$treatment.num <- as.numeric(bloodpressureL$treatment)
## ** compound symmetry
eCS.gls <- gls(duration ~ period + treatment,
data = bloodpressureL,
correlation = corCompSymm(form=~ 1 | id))
eCS.lmm <- lmm(duration ~ period + treatment,
data = bloodpressureL,
structure = "CS", repetition = ~ period | id)
expect_equal(as.double(logLik(eCS.lmm)),as.double(logLik(eCS.gls)), tol = 1e-6)
plot(eCS.lmm, color = "sequence", obs.alpha = 0.25, ci.alpha = 0.05)
confint(eCS.gls)
confint(eCS.lmm)
## ** unstructured
eUNP.gls <- gls(duration ~ period + treatment,
data = bloodpressureL,
correlation = corSymm(form=~ period.num | id),
weights = varIdent(form=~ 1|period),
)
eUNP.lmm <- lmm(duration ~ period + treatment,
data = bloodpressureL,
structure = "UN", repetition = ~ period | id)
expect_equal(as.double(logLik(eUNP.lmm)),as.double(logLik(eUNP.gls)), tol = 1e-6)
capture.output(summary(eUNP.lmm))
eUNT.gls <- gls(duration ~ period + treatment,
data = bloodpressureL,
correlation = corSymm(form=~ treatment.num | id),
weights = varIdent(form=~ 1|treatment),
)
eUNT.lmm <- lmm(duration ~ period + treatment,
data = bloodpressureL,
structure = "UN", repetition = ~ treatment | id)
capture.output(summary(eUNT.lmm))
expect_equal(as.double(logLik(eUNT.lmm)),as.double(logLik(eUNT.gls)), tol = 1e-6)
})
## * Practical 6
test_that("practical 6 - vasscoresL",{
data(vasscoresL, package = "LMMstar")
summarize(vas ~ treatment, data = vasscoresL, na.rm = TRUE)
summarize(vas ~ treatment|id, data = vasscoresL, na.rm = TRUE)
summarize(vas ~ treatment+group|id, data = vasscoresL, na.rm = TRUE)
## ** model on all pairs: AB, AC, BC
fit.CS <- lmm(vas~-1+treatment, data=vasscoresL,
repetition=~treatment|id, structure="CS")
capture.output(summary(fit.CS))
GS <- matrix(c(2099.98852585, 1675.79243303, 1675.79243303, 1675.79243303, 2099.98852585, 1675.79243303, 1675.79243303, 1675.79243303, 2099.98852585),
nrow = 3,
ncol = 3,
dimnames = list(c("A", "B", "C"),c("A", "B", "C"))
)
expect_equivalent(sigma(fit.CS), GS, tol = 1e-5)
## GS <- lmer(vas~-1+treatment+(1|id), data=vasscoresL)
GS <- lme(vas~-1+treatment, random =~ 1|id, data=vasscoresL, na.action = na.omit)
expect_equivalent(unname(coef(fit.CS)), unname(fixef(GS)), tol = 1e-5)
## expect_equivalent(unname(model.tables(fit.CS)$se), unname(summary(GS)$coef[,"Std. Error"]), tol = 1e-2)
## expect_equivalent(unname(model.tables(fit.CS)$df), unname(summary(GS)$coef[,"df"]), tol = 1e-2)
expect_equivalent(unname(model.tables(fit.CS)$se), unname(summary(GS)$tTable[,"Std.Error"]), tol = 1e-2)
expect_equivalent(unname(model.tables(fit.CS)$df), c(43.64236203, 43.64236203, 43.64236203), tol = 1e-2)
## autoplot(fit.CS)
suppressWarnings(plot(fit.CS, obs.alpha = 0.1))
fit.UN <- lmm(vas~-1+treatment, data=vasscoresL,
repetition=~treatment|id, structure="UN")
capture.output(summary(fit.UN))
## GS <- lmm(vas~-1+treatment, data=vasscoresL,
## repetition=~treatment|id, structure="UN",
## control = list(optimizer = "FS"))
GS <- matrix(c(2038.74222977, 1892.18702336, 1797.44635389, 1892.18702336, 2111.43949953, 1369.79020341, 1797.44635389, 1369.79020341, 2157.87873809),
nrow = 3,
ncol = 3,
dimnames = list(c("A", "B", "C"),c("A", "B", "C"))
)
expect_equivalent(sigma(fit.UN), GS, tol = 1e-5)
vasscoreL.imputed <- fitted(fit.UN, type = "outcome", keep.data = TRUE)
set.seed(11)
gg <- ggplot(vasscoreL.imputed, aes(x = treatment, y = vas, group = id))
gg + geom_point(aes(color = impute)) + geom_line()
## ** model some pairs: AB, AC
## IMPORTANT: check extraction of the residual variance-covariance matrix when not all pairs of time are analyzed
vasscoreL.red <- vasscoresL[vasscoresL$group %in% c("AB","AC"),]
fit.CS.red <- lmm(vas~-1+treatment, data=vasscoreL.red,
repetition=~treatment|id, structure="CS")
capture.output(summary(fit.CS.red))
GS <- matrix(c(2006.9354323, 1763.57527658, 1763.57527658, 1763.57527658, 2006.9354323, 1763.57527658, 1763.57527658, 1763.57527658, 2006.9354323),
nrow = 3,
ncol = 3,
dimnames = list(c("A", "B", "C"),c("A", "B", "C"))
)
expect_equivalent(sigma(fit.CS.red), GS, tol = 1e-5)
fit.UN.red <- lmm(vas~-1+treatment, data=vasscoreL.red,
repetition=~treatment|id, structure="UN", control = list(optimizer = "FS")) ## FS optimizer mandatory, GLS does something weird as it is able to estimate the correlation (B,C)
capture.output(summary(fit.UN.red))
GS <- matrix(c(1996.68088251, 1926.29216631, 1573.59329163, 1926.29216631, 2231.4711806, NA, 1573.59329163, NA, 1753.88563313),
nrow = 3,
ncol = 3,
dimnames = list(c("A", "B", "C"),c("A", "B", "C"))
)
expect_equivalent(sigma(fit.UN.red), GS, tol = 1e-5)
})
##----------------------------------------------------------------------
### test-auto-practical.R ends here
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.