#======================================================================= # Loading data and packages #======================================================================= if (!require("pacman")) { install.packages("pacman") } pacman::p_load(lcc, cccrm, dplyr, reshape2, tidyverse, ICC) #Datasets: data(simulated_hue_block) data("simulated_hue", package = "lcc") data(bfat, package = "cccrm") data(hue) data(bdaw, package = "cccrm") #======================================================================= # Specifying models in the lcc function #======================================================================= attach(simulated_hue_block) m1 <- lcc(data = simulated_hue_block, subject = "Fruit", resp = "Hue", method = "Method",time = "Time", qf = 2, qr = 1) m2 <- update(m1, covar ="Block") m3 <- update(m2, var.class = varIdent, weights.form = "method", lme.control = list(opt="optim")) summary(m3) lccPlot(m3,scale="free") AIC(m2, m3); BIC(m2, m3); anova(m2, m3) detach(simulated_hue_block) #======================================================================= # Example 1: Percentage body fat dataset #======================================================================= bfat <- bfat %>% mutate(VISITNO = replace(VISITNO, VISITNO == 2, 12.5)) %>% mutate(VISITNO = replace(VISITNO, VISITNO == 3, 13)) %>% mutate(VISITNO = replace(VISITNO, VISITNO == 4, 13.5)) %>% mutate(SUBJECT = factor(SUBJECT)) %>% mutate(MET = factor(MET, labels = c("skinfold", "DEXA"))) bfat$TIME <- 12 * (bfat$VISITNO - 12) #----------------------------------------------------------------------- # Scatter plot of body data #----------------------------------------------------------------------- bfat_wide <- dcast(bfat, SUBJECT + VISITNO ~ MET, value.var="BF") bfat_wide colnames(bfat_wide) bfat_wide$VISITNO <- as.factor(bfat_wide$VISITNO) levels(bfat_wide$VISITNO) <- c("12.5 (Visit 2)", "13 (Visit 3)", "13.5 (Visit 4)") p1 <- ggplot(bfat_wide, aes(x = DEXA, y = skinfold)) + facet_wrap(~VISITNO) + geom_point(alpha = 0.7, colour = "black", fill = "gray", size = 1) + geom_abline(intercept = 0, slope = 1) + geom_smooth(se = FALSE, method = "lm") p1 +labs(y = "Skinfold caliper", x = "DEXA")+ theme(legend.position = "none", axis.line.x = element_line(color="black", size = 0.5), axis.line.y = element_line(color="black", size = 0.5), aspect.ratio = 1)+ theme_bw() #----------------------------------------------------------------------- # lcc fit #----------------------------------------------------------------------- set.seed(134) m.bfat.1 <- lcc(data = bfat, subject = "SUBJECT", resp = "BF", method = "MET", time = "TIME", qf = 1, qr = 1, components = TRUE, ci = T, nboot = 10000) set.seed(134) m.bfat.2 <- update(m.bfat.1, lme.control = list(opt = "optim")) summary(m.bfat.2, type = "lcc") lccPlot(m.bfat.2, type = "lcc") + ylim(0,1) + geom_hline(yintercept = 1, linetype = "dashed") + theme_bw() lccPlot(m.bfat.2, type = "lpc")+ ylim(0,1) + geom_hline(yintercept = 1, linetype = "dashed") + theme_bw() lccPlot(m.bfat.2, type = "la")+ ylim(0,1) + geom_hline(yintercept = 1, linetype = "dashed") + theme_bw() #----------------------------------------------------------------------- # Diagnostic plots #----------------------------------------------------------------------- plot(m.bfat.2) #======================================================================= # Example 2: hue data #======================================================================= #----------------------------------------------------------------------- # Scatter Plots #----------------------------------------------------------------------- hue_wide <- dcast(hue, Fruit + Time ~ Method, value.var="H_mean") #----------------------------------------------------------------------- # Scanner versus Colorimeter #----------------------------------------------------------------------- p1 <- ggplot(hue_wide, aes(x = Colorimeter, y = Scanner)) + geom_point(alpha = 0.7, colour = "black", fill = "gray", size = 1) + geom_abline(intercept = 0, slope = 1) + geom_smooth(se = FALSE, method = "lm") p1 +labs(y = "Scanner", x = "Colorimeter")+ theme_bw()+ theme(legend.position = "none", aspect.ratio = 1, axis.line.x = element_line(color="black", size = 0.5), axis.line.y = element_line(color="black", size = 0.5), axis.title.x = element_text(size=14), axis.title.y = element_text(size=14), axis.text.x = element_text(size = 14, face = "plain"), axis.text.y = element_text(size = 14, face = "plain")) #----------------------------------------------------------------------- # Individual Profiles #----------------------------------------------------------------------- p2 <- ggplot(hue, aes(y = H_mean, x = Time, group = Fruit)) + facet_wrap(~ Method) + geom_line(aes(color = Fruit)) p2 +labs(y = "Mean Hue", x = "Time (Day)")+ theme_bw()+ theme(legend.position = "none", aspect.ratio = 1, axis.line.x = element_line(color="black", size = 0.5), axis.line.y = element_line(color="black", size = 0.5), axis.title.x = element_text(size=14), axis.title.y = element_text(size=14), axis.text.x = element_text(size = 14, face = "plain"), axis.text.y = element_text(size = 14, face = "plain"), strip.text.x = element_text(size = 14)) #----------------------------------------------------------------------- # Individual confidence intervals #----------------------------------------------------------------------- # Model m.hue.1 <- lmList(H_mean ~ poly(Time, 2, raw = TRUE)| Fruit, hue) summary(m.hue.1) # box-plots of residuals by Subject plot(m.hue.1, Fruit ~ resid(.)) # observed versus fitted values by Subject plot(m.hue.1, H_mean ~ fitted(.) | Fruit, abline = c(0,1)) # Confidence interval plot(intervals(m.hue.1)) #----------------------------------------------------------------------- # lcc model 1 #----------------------------------------------------------------------- set.seed(6836) m.hue.2 <- lcc(data = hue, subject = "Fruit", resp = "H_mean", method = "Method", time = "Time", qf = 2, qr = 2, ci = TRUE, nboot = 10000, components = TRUE) lccPlot(m.hue.2)+ ylim(0,1) + geom_hline(yintercept = 1, linetype = "dashed") + scale_x_continuous(breaks = seq(1,max(hue$Time),2))+ theme_bw() + theme(legend.position = "none", aspect.ratio = 1, axis.line.x = element_line(color="black", size = 0.5), axis.line.y = element_line(color="black", size = 0.5), axis.title.x = element_text(size=14), axis.title.y = element_text(size=14), axis.text.x = element_text(size = 14, face = "plain"), axis.text.y = element_text(size = 14, face = "plain")) lccPlot(m.hue.2, type = "lpc")+ ylim(0,1) + geom_hline(yintercept = 1, linetype = "dashed") + scale_x_continuous(breaks = seq(1,max(hue$Time),2))+ theme_bw() + theme(legend.position = "none", aspect.ratio = 1, axis.line.x = element_line(color="black", size = 0.5), axis.line.y = element_line(color="black", size = 0.5), axis.title.x = element_text(size=14), axis.title.y = element_text(size=14), axis.text.x = element_text(size = 14, face = "plain"), axis.text.y = element_text(size = 14, face = "plain")) lccPlot(m.hue.2, type = "la")+ ylim(0,1) + geom_hline(yintercept = 1, linetype = "dashed") + scale_x_continuous(breaks = seq(1,max(hue$Time),2))+ theme_bw() + theme(legend.position = "none", aspect.ratio = 1, axis.line.x = element_line(color="black", size = 0.5), axis.line.y = element_line(color="black", size = 0.5), axis.title.x = element_text(size=14), axis.title.y = element_text(size=14), axis.text.x = element_text(size = 14, face = "plain"), axis.text.y = element_text(size = 14, face = "plain")) #----------------------------------------------------------------------- # Residuals #----------------------------------------------------------------------- plot(m.hue.2, which = c(5, 6)) plot(m.hue.2$model, resid(., type = "p") ~ time | method, abline = 0) hue$res <- residuals(m.hue.2) res <- with(hue, tapply(res,list(Time,Method),var)) plot(res[, 1] ~ seq(0, 14, 1), type = "b", ylab = "Residual variance", xlab = "Time") points(res[, 2] ~ seq(0, 14, 1), col = "red", pch = 2, type = "b", lty = 2) legend(y = 2.5, x = 1.5, legend = c("Colorimeter", "Scanner"), pch = c(1, 2), lty = c(1, 2), col = c("black", "red")) #----------------------------------------------------------------------- # lcc model 2 #----------------------------------------------------------------------- lcc_time <- with(hue, list(time = Time, from =min(Time), to=max(Time), n=50)) set.seed(6836) m.hue.3 <- update(m.hue.2, var.class = varIdent, weights.form = "method", time_lcc = lcc_time, lme.control = lmeControl(opt = "optim")) #----------------------------------------------------------------------- # Likelihood ratio test #----------------------------------------------------------------------- anova(m.hue.2, m.hue.3) summary(m.hue.3) summary(m.hue.3, type = "lcc")$gof #----------------------------------------------------------------------- # Plots #----------------------------------------------------------------------- lccPlot(m.hue.3)+ ylim(0,1) + geom_hline(yintercept = 1, linetype = "dashed") + scale_x_continuous(breaks = seq(1,max(hue$Time),2))+ theme_bw() + theme(legend.position = "none", aspect.ratio = 1, axis.line.x = element_line(color="black", size = 0.5), axis.line.y = element_line(color="black", size = 0.5), axis.title.x = element_text(size=14), axis.title.y = element_text(size=14), axis.text.x = element_text(size = 14, face = "plain"), axis.text.y = element_text(size = 14, face = "plain")) lccPlot(m.hue.3, type = "lpc")+ ylim(0,1) + geom_hline(yintercept = 1, linetype = "dashed") + scale_x_continuous(breaks = seq(1,max(hue$Time),2))+ theme_bw() + theme(legend.position = "none", aspect.ratio = 1, axis.line.x = element_line(color="black", size = 0.5), axis.line.y = element_line(color="black", size = 0.5), axis.title.x = element_text(size=14), axis.title.y = element_text(size=14), axis.text.x = element_text(size = 14, face = "plain"), axis.text.y = element_text(size = 14, face = "plain")) lccPlot(m.hue.3, type = "la")+ ylim(0,1) + geom_hline(yintercept = 1, linetype = "dashed") + scale_x_continuous(breaks = seq(1,max(hue$Time),2))+ theme_bw() + theme(legend.position = "none", aspect.ratio = 1, axis.line.x = element_line(color="black", size = 0.5), axis.line.y = element_line(color="black", size = 0.5), axis.title.x = element_text(size=14), axis.title.y = element_text(size=14), axis.text.x = element_text(size = 14, face = "plain"), axis.text.y = element_text(size = 14, face = "plain")) #----------------------------------------------------------------------- # Residuals #----------------------------------------------------------------------- plot(m.hue.3, which = 5, ask = F) plot(m.hue.3$model, resid(., type = "p") ~ time | method, abline = 0) plot(Variogram(m.hue.3$model, form = ~ hue$Time, maxDist = 42), ylim = c(0, 1.4)) #======================================================================= # Example 3: the blood draw dataset #======================================================================= #----------------------------------------------------------------------- # Selecting the individual profiles #----------------------------------------------------------------------- bdaw$SUBJ <- as.factor(bdaw$SUBJ) bdaw$MET <- as.factor(bdaw$MET) levels(bdaw$MET) <- c("1 hour", "2 hours") length(unique(bdaw$SUBJ)) # Model fit_list <- lmList(AUC ~ poly(VNUM, 4) | SUBJ, data = bdaw) int <- intervals(fit_list) zero_included <- function(x) { flag <- min(x) < 0 & max(x) > 0 return(flag) } selected_subj<- names( which(apply(int[,,4], 1, zero_included) & apply(int[,,5], 1, zero_included))) bdaw_subset <- subset(bdaw, SUBJ %in% selected_subj) #----------------------------------------------------------------------- # Plots #----------------------------------------------------------------------- bdaw1 <- bdaw_subset %>% filter(MET == "1 hour") bdaw2 <- bdaw_subset %>% filter(MET != "1 hour") p1 <- ggplot(bdaw_subset, aes(y = AUC, x = VNUM, group = SUBJ, color = SUBJ)) + geom_line(size = 0.6) + facet_wrap(~MET) p1 +labs(y = "AUC", x = "Visits")+ theme_bw()+ theme(legend.position = "none", aspect.ratio = 1, axis.line.x = element_line(color="black", size = 0.5), axis.line.y = element_line(color="black", size = 0.5), axis.title.x = element_text(size=14), axis.title.y = element_text(size=14), axis.text.x = element_text(size = 14, face = "plain"), axis.text.y = element_text(size = 14, face = "plain"), strip.text.x = element_text(size = 14)) p2 <- ggplot(bdaw1, aes(y = bdaw2$AUC, x = bdaw1$AUC)) + geom_point() + geom_smooth(method = "lm", se = FALSE) + geom_abline(intercept = 0, slope = 1) p2 +labs(y = "AUC every 2 hours", x = "AUC every 1 hour")+ theme_bw()+ theme(legend.position = "none", aspect.ratio = 1, axis.line.x = element_line(color="black", size = 0.5), axis.line.y = element_line(color="black", size = 0.5), axis.title.x = element_text(size=14), axis.title.y = element_text(size=14), axis.text.x = element_text(size = 14, face = "plain"), axis.text.y = element_text(size = 14, face = "plain"), strip.text.x = element_text(size = 14)) #----------------------------------------------------------------------- # Individual confidence intervals #----------------------------------------------------------------------- fit_list2 <- lmList(AUC ~ poly(VNUM, 4) | SUBJ, data = bdaw_subset) plot(intervals(fit_list2)) length(unique(bdaw_subset$SUBJ)) #----------------------------------------------------------------------- # lcc model #----------------------------------------------------------------------- m.bw.1 <- lcc(data = bdaw_subset, subject = "SUBJ", resp = "AUC", method = "MET", time = "VNUM", qf = 1, qr = 1) summary(m.bw.1, type = "lcc")$gof #----------------------------------------------------------------------- # Plot #----------------------------------------------------------------------- lccPlot(m.bw.1) #----------------------------------------------------------------------- # Residuals #----------------------------------------------------------------------- plot(m.bw.1, which = c(1, 2, 4, 5, 6)) #----------------------------------------------------------------------- # lcc model 2 #----------------------------------------------------------------------- set.seed(542361) m.bw.2 <- update(m.bw.1, qf = 2, qr = 2, components = TRUE, time_lcc = list(from = 3, to = 7, n = 50), ci = TRUE, nboot = 10000, show.warnings = TRUE, numCore = 4, lme.control = lmeControl(msMaxIter = 200, msMaxEval = 600, maxIter = 200)) summary(m.bw.2) #----------------------------------------------------------------------- # Likelihood ratio test for qr=2 #----------------------------------------------------------------------- m.bw.3 <- update(m.bw.1, qf = 2) anova(m.bw.3, m.bw.2) #----------------------------------------------------------------------- summary(m.bw.1, type = "lcc")$gof summary(m.bw.2, type = "lcc")$gof summary(m.bw.3, type = "lcc")$gof #----------------------------------------------------------------------- # Plots #----------------------------------------------------------------------- lccPlot(m.bw.2)+ ylim(0.85, 1)+ geom_hline(yintercept = 1, linetype = "dashed") + scale_x_continuous(breaks = seq(1,max(hue$Time),2))+ theme_bw() + theme(legend.position = "none", aspect.ratio = 1, axis.line.x = element_line(color="black", size = 0.5), axis.line.y = element_line(color="black", size = 0.5), axis.title.x = element_text(size=14), axis.title.y = element_text(size=14), axis.text.x = element_text(size = 14, face = "plain"), axis.text.y = element_text(size = 14, face = "plain")) lccPlot(m.bw.2, type = "lpc")+ ylim(0.85, 1)+ geom_hline(yintercept = 1, linetype = "dashed") + scale_x_continuous(breaks = seq(1,max(hue$Time),2))+ theme_bw() + theme(legend.position = "none", aspect.ratio = 1, axis.line.x = element_line(color="black", size = 0.5), axis.line.y = element_line(color="black", size = 0.5), axis.title.x = element_text(size=14), axis.title.y = element_text(size=14), axis.text.x = element_text(size = 14, face = "plain"), axis.text.y = element_text(size = 14, face = "plain")) lccPlot(m.bw.2, type = "la")+ ylim(0.85, 1) + geom_hline(yintercept = 1, linetype = "dashed") + scale_x_continuous(breaks = seq(1,max(hue$Time),2))+ theme_bw() + theme(legend.position = "none", aspect.ratio = 1, axis.line.x = element_line(color="black", size = 0.5), axis.line.y = element_line(color="black", size = 0.5), axis.title.x = element_text(size=14), axis.title.y = element_text(size=14), axis.text.x = element_text(size = 14, face = "plain"), axis.text.y = element_text(size = 14, face = "plain")) #----------------------------------------------------------------------- # Residuals #----------------------------------------------------------------------- plot(m.bw.2) #----------------------------------------------------------------------- # Testing Interaction #----------------------------------------------------------------------- m.bw.4 <- lcc(data = bdaw_subset, subject = "SUBJ", resp = "AUC", method = "MET", time = "VNUM", qf = 2, qr = 2, REML = FALSE, interaction = FALSE) m.bw.5 <- update(m.bw.4, interaction = TRUE) anova(m.bw.4, m.bw.5) #-----------------------------------------------------------------------