#=======================================================================
# 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)
#-----------------------------------------------------------------------