######################################################################
# Setting global variables
options(prompt = "R> ", continue = "+  ", width = 70, useFancyQuotes = FALSE)
par(mar = c(5.1, 4.1, 4.1, 2.1))
par.old <- par(no.readonly = TRUE)

#######################################################
# Section entitled 'Structure of the cglasso Package' #

# install.packages("cglasso")
library("cglasso")

###################################################################
# Section entitled 'Network inference from multiple myeloma data' #

data("MM", package = "cglasso")

MM.datacggm <- datacggm(Y = MM$Y, X = MM$X, up = 40)

Y <- getMatrix(MM.datacggm, name = "Y")
X <- getMatrix(MM.datacggm, name = "X")

dim(MM.datacggm)

dimnames(MM.datacggm)$Y |> lapply(head)

dimnames(MM.datacggm)$X |> lapply(head)

rowNames(MM.datacggm)$Y <- paste0("u", seq_len(nobs(MM.datacggm)))
rowNames(MM.datacggm)$X <- paste0("u", seq_len(nobs(MM.datacggm)))
colNames(MM.datacggm)$Y <- paste0("Y", seq_len(nresp(MM.datacggm)))
colNames(MM.datacggm)$X <- c(paste0("X", seq_len(15L)), "CyAb")

print(MM.datacggm, n = 5L)

summary(MM.datacggm, n = 5L)

mylambda <- c(0.30, 0.20)
myrho <- c(0.50, 0.40)
MM.mdl <- cglasso(data = MM.datacggm, lambda = mylambda, rho = myrho)

MM.mdl

coef(MM.mdl) |> lapply(dim)

MM.mdl.fitted <- fitted(MM.mdl) 
MM.mdl.obsRes <- residuals(MM.mdl)
MM.mdl.imputed <- impute(MM.mdl)

B.pred <- predict(MM.mdl, lambda.new = 0.25, rho.new = 0.45)

MM.mdl.mle <- cggm(MM.mdl)

class(MM.mdl.mle)

mylambda <- seq(from = 0.30, to = 0.20, length = 10L)
myrho <- seq(from = 0.50, to = 0.40, length = 10L)
MM.mdl.path <- cglasso(data = MM.datacggm, lambda = mylambda,
   rho = myrho)

MM.mdl.path.eBIC <- BIC(MM.mdl.path, g = 0.5, mle = TRUE)

plot(MM.mdl.path.eBIC, main = "Multiple Myeloma Dataset")

MM.optmdl <- select.cglasso(MM.mdl.path, GoF = MM.mdl.path.eBIC)

summary(MM.optmdl)

MM.optmdl.graph <- to_graph(MM.optmdl)

names(MM.optmdl.graph)

Gyy <- getGraph(MM.optmdl.graph)
Gxy <- getGraph(MM.optmdl.graph, type = "Gxy")


plot(MM.optmdl.graph, xlim = c(-0.85, 0.85), ylim = c(-0.85, 0.85),
   main = "Multiple Myeloma Dataset", layout = layout_with_kk)

plot(MM.optmdl.graph, xlim = c(-0.85, 0.85), ylim = c(-0.85, 0.85),
   main = "Multiple Myeloma Dataset", type = "Gyy",
   layout = layout_with_kk)

plot(MM.optmdl.graph, xlim = c(-0.85, 0.85), ylim = c(-0.85, 0.85),
   main = "Multiple Myeloma Dataset", type = "Gxy",
   layout = layout_with_kk)
par(par.old)

####################################################################################
# Section entitled 'Network inference from megakaryocyte-erythroid progenitors data'

data("MKMEP", package = "cglasso")

MKMEP.datacggm <- datacggm(Y = MKMEP, up = 40)

rowNames(MKMEP.datacggm)$Y <- paste0("u", seq_len(nobs(MKMEP.datacggm)))
colNames(MKMEP.datacggm)$Y <- paste0("Y", seq_len(nresp(MKMEP.datacggm)))

print(MKMEP.datacggm, n = 5L)

myrho <- seq(from = 6, to = 3, length = 100L)
MKMEP.mdl.path <- cglasso(data = MKMEP.datacggm, rho = myrho)

MKMEP.mdl.path.eBIC <- BIC(MKMEP.mdl.path, g = 0.5, mle = TRUE)

plot(MKMEP.mdl.path.eBIC,
   main = "Megakaryocyte-Erythroid Progenitors Dataset")

MKMEP.optmdl.mle <- cggm(MKMEP.mdl.path, GoF = MKMEP.mdl.path.eBIC)
MKMEP.optmdl.mle

MKMEP.optmdl.graph <- to_graph(MKMEP.optmdl.mle)

plot(MKMEP.optmdl.graph, xlim = c(-0.85, 0.85), ylim = c(-0.85, 0.85),
   main = "Megakaryocyte-Erythroid Progenitors Dataset",
   layout = layout_with_kk)
par(par.old)

################################################################################
# Section entitled 'Simulating `datacggm' objects and plotting functions'

set.seed(1234)
n <- 1000L
p <- 4L
b0 <- rep(10, time = p)
Sigma <- outer(1:p, 1:p, function(i, j) 0.3^abs(i - j))

probl <- 0.1
probr <- 0.1
probna <- 0.1

Z.sim <- rcggm(n = n, b0 = b0, Sigma = Sigma, probl = probl,
   probr = probr, probna = probna)

summary(Z.sim)

hist(Z.sim, max.hist = 4L, breaks = 5L)

qqcnorm(Z.sim, max.plot = 4L)

####################################################
# Section entitled 'The objects of class datacggm' #

names(MM.datacggm)

names(MM.datacggm$Info)

MKMEP.datacggm$X; MKMEP.datacggm$Info$q


#########################################
# Section entitled 'Accessor functions' #

Y <- getMatrix(MM.datacggm, name = "Y", ordered = TRUE)

R <- event(MM.datacggm,  ordered = TRUE)
cbind(Y = Y[, 1L], R = R[, 1L]) |> tail()

ColMeans(MM.datacggm) |> lapply(head, n = 5L)

########################################################
# Section entitled 'Fitting conditional glasso models' #

MKMEP.mdl1 <- cglasso(data = MKMEP.datacggm)

MKMEP.mdl1


##########################################################
# Section entitled 'Plotting an object of class cglasso' #

set.seed(1234)
sample(nresp(MM.datacggm), 5L) |> sort() -> Y.id
colNames(MM.datacggm)$Y[Y.id]

mylambda <- c(0.2, 0.1)
MM.submdl <- cglasso(cbind(Y16, Y22, Y28, Y37, Y44) ~ .,
   data = MM.datacggm, nrho = 100L, lambda = mylambda)

plot(MM.submdl, what = "Theta", penalty = "rho", given = 1L,
   matplot.arg1 = list(ylab = "Precision Values",
   main = "Multiple Myeloma Dataset"), GoF = AIC)

plot(MM.submdl, what = Theta ~ rho | 1,
   matplot.arg1 = list(ylab = "Precision Values",
   main = "Multiple Myeloma Dataset"), GoF = AIC)

set.seed(1234)
sample(nresp(MKMEP.datacggm), 5L) |> sort() -> Y.id
colNames(MKMEP.datacggm)$Y[Y.id]

MKMEP.submdl <- cglasso(cbind(Y16, Y22, Y28, Y37, Y58) ~ ., 
   data = MKMEP.datacggm, nrho = 100L)

plot(MKMEP.submdl, matplot.arg1 = list(ylab = "Precision Values",
   main = "Megakaryocyte-Erythroid Progenitors Dataset"),
   GoF = AIC)

###################################################################
# Section entitled 'Further method functions for cglasso objects' #

c(lambda = MM.submdl$lambda[1L], rho = MM.submdl$rho[1L])

Theta_hat <- coef(MM.submdl, type = "Theta", lambda.id = 1L,
   rho.id = 1L)
dim(Theta_hat)

FittedValues <- fitted(MM.submdl, type = "Theta", lambda.id = 1L,
   rho.id = 1L)

WorkingResiduals <- residuals(MM.submdl, type = "working",
   lambda.id = 1L, rho.id = 1L)

####################################################################################
# Section entitled 'Prediction and missing data imputation using a cglasso object' #

lambda.new <- mean(MM.submdl$lambda)
rho.new <- mean(MM.submdl$rho)
mu_pred <- predict(MM.submdl, type = "mu", lambda.new = lambda.new,
   rho.new = rho.new)

#######################################################################
# Section entitled 'Post-hoc maximum likelihood refitting of a model' #

MM.mdl2 <- cglasso(data = MM.datacggm, nlambda = 20L, nrho = 20L, 
   lambda.min.ratio = 0.2, rho.min.ratio = 0.2)

MM.mdl2.mle <- cggm(MM.mdl2, GoF = BIC)
MM.mdl2.mle

class(MM.mdl2.mle)

Theta.MM.mdl2.mle <- coef(MM.mdl2.mle, type = "Theta")

##########################################
# Section entitled 'Measuring model fit' #

mylambda <- c(1, 0.5)
MM.mdl1 <- cglasso(data = MM.datacggm, lambda = mylambda, nrho = 5L)
MM.mdl1.Qvalue <- QFun(MM.mdl1, mle = FALSE)
MM.mdl1.Qvalue

####################################################
# Section entitled 'The goodness-of-fit functions' #

MM.mdl1.AIC <- AIC(MM.mdl1, mle = FALSE)
class(MM.mdl1.AIC)

names(MM.mdl1.AIC)

MM.mdl1.AIC

MKMEP.mdl1.eBIC <- BIC(MKMEP.mdl1, g = 0.5, type = "FD", mle = TRUE)
MKMEP.mdl1.eBIC

MM.mdl1.eBIC <- BIC(MM.mdl1, g = 0.5, type = "CC", mle = TRUE)
MM.mdl1.eBIC

###################################################
# Section entitled 'Summary of the fitted models' #

summary(MM.mdl1, GoF = MM.mdl1.eBIC)

summary(MM.mdl1, GoF = BIC, g = 0.5, mle = TRUE)

sessionInfo()