################################################### ### 4.1 Hematologic parameters from Sprague-Dawley rats ################################################### library("CLME") data("rat.blood", package = "CLME") const <- list(order = "simple", node = 1, decreasing = FALSE) hct1 <- clme(hct ~ time + temp + sex + (1|id), data = rat.blood, constraints = const, levels = list(2, levels(rat.blood$time))) hct1b <- summary(hct1, seed = 42) hct1b plot(hct1b, ci = TRUE, legendx = "bottomright", inset = 0.08) hct2 <- clme(hct ~ time + temp + sex + (1|id), data = rat.blood, levels = list(2, levels(rat.blood$time))) summary(hct2, seed = 42) hct3 <- clme(hct ~ time + temp + sex + (1|id), data = rat.blood, gfix = rat.blood$time, constraints = const, levels = list(2, levels(rat.blood$time))) summary(hct3, seed = 42) const <- list(order = "simple.tree", node = 1, decreasing = FALSE) wbc <- summary(clme(wbc ~ time + temp + sex + (1|id), data = rat.blood, constraints = const, levels = list(2, levels(rat.blood$time)), tsf = w.stat), seed = 42) wbc plot(wbc, legend = "topleft", inset = 0.08) ################################################### ### 4.2. Fibroid growth rates ################################################### data("fibroid", package = "CLME") race.age <- factor(paste0(fibroid$race, ".", fibroid$age), levels = c("Black.Yng", "Black.Mid", "Black.Old", "White.Yng", "White.Mid", "White.Old") ) fibroid$race.age <- race.age fibroid$initVol <- cut(fibroid$vol, c(0, 14000, 65000, Inf), c("small", "medium", "large"), right = FALSE) const <- list(A = cbind(2:6, 1:5)[-3, ], B = rbind(c(3, 1), c(6, 4))) const w.blk.wht <- function (theta, cov.theta, B, A, ...) { stats <- vector("numeric", length = nrow(B)) ctd <- diag(cov.theta) stats <- apply(B, 1, FUN = function(a, theta, cov, ctd) { std <- sqrt(ctd[a[1]] + ctd[a[2]] - 2 * cov.theta[a[1], a[2]]) (theta[a[2]] - theta[a[1]]) / std }, theta = theta, cov = cov.theta, ctd = ctd) names(stats) <- c("Black.Yng - Black.Old", "White.Yng - White.Old" ) return(stats) } fib <- summary(clme(lfgr ~ race.age + initVol + (1|id), data = fibroid, constraints = const, tsf = w.blk.wht, levels = list(10, levels(race.age))), seed = 42) fib plot(x = 1, y = 0, col = 0, ylim = c(-6, 22), xlim = c(0.9, 3.1), xlab = "", ylab = "Estimated Coefficient", xaxt = "n") axis(side = 1, at = 1:3, labels = c("Young (<35)", "Middle aged (35-44)", "Older (>45)") ) for (y1 in seq(-15, 25, 5)) { lines( x = c(0, 7), y = c(y1, y1), col = "gray", lty = 2 ) } lty1 <- 1 + (fib$p.value.ind[1] < 0.05) lty2 <- 1 + (fib$p.value.ind[2] < 0.05) lty3 <- 1 + (fib$p.value.ind[3] < 0.05) lty4 <- 1 + (fib$p.value.ind[4] < 0.05) points(c(1, 2), fib$theta[1:2], col = 1, type = "l", lwd = 2 , lty = lty1) points(c(2, 3), fib$theta[2:3], col = 1, type = "l", lwd = 2 , lty = lty2) points(c(1, 2), fib$theta[4:5], col = 3, type = "l", lwd = 2 , lty = lty3) points(c(2, 3), fib$theta[5:6], col = 3, type = "l", lwd = 2 , lty = lty4) points(1:3, fib$theta[1:3], col = 1, cex = 1.5, pch = 21, bg = "white") points(1:3, fib$theta[4:6], col = 3, cex = 1.5, pch = 24, bg = "white") legend("bottom", lty = c(1, 1), pch = c(21, 24), col = c(1, 3), pt.bg = 0, pt.cex = 1.1, lwd = 2, inset = 0.03, cex = 0.9, legend = c("Blacks ", "Whites")) library("lme4") cons <- list(order = "simple", decreasing = FALSE, node = 1) clme1 <- clme(hct ~ time + temp + sex + (1|id), data = rat.blood, constraints = cons) lmer1 <- lmer(hct ~ time + temp + sex + (1|id) - 1, data = rat.blood) lmer2 <- lmer(hct ~ time + temp + sex + (1|id) - 1, data = rat.blood, REML = FALSE) rbind(fixef(clme1), fixef(lmer1), fixef(lmer2)) VarCorr(clme1) VarCorr(lmer1) VarCorr(lmer2) c(AIC(clme1), AIC(lmer1), AIC(lmer2)) library("rbenchmark") set.seed(42) BMclme <- function(x = 0){ clme(hct ~ time + temp + sex + (1|id), data = rat.blood, constraints = cons) } BMlmer <- function(x = 0){ lmer(hct ~ time + temp + sex + (1|id), data = rat.blood ) } benchmark(BMlmer(), BMclme(), replications = 100 )[, 1:4] set.seed(42) beta_sum <- function(.) { beta = fixef(.) } BMclme <- function(x = 0) { summary(clme1, nsim = 100, seed = 42) } BMlmer1 <- function(x = 0) { bootMer(lmer1, FUN = beta_sum, nsim = 100, seed = 42, type = "parametric") } BMlmer2 <- function(x = 0) { bootMer(lmer1, FUN = beta_sum, nsim = 100, seed = 42, use.u = TRUE, type = "semiparametric") } benchmark(BMclme(), BMlmer1(), BMlmer2(), replications = 10)[, 1:4]