## ----setup--------------------------------------------------------------------
library(DGP4LCF)
## ----eval = FALSE-------------------------------------------------------------
#
# # for reproducibility purpose
# set.seed(456)
#
# # create objects required in the function input
# q<- 8
# n<- 17
#
# obs_time_num<- rep(q, times = n)
# obs_time_index<- list()
# a_person<- list()
# col_person_index<- list()
# observed_x_train_regular_8<- list()
#
# for (person_index in 1:n){
# obs_time_index[[person_index]]<- (1:q)
# a_person[[person_index]]<- sim_fcs_truth$a_full[(1:q)]
# col_person_index[[person_index]]<- ((person_index-1)*q + 1):(person_index*q)
# observed_x_train_regular_8[[person_index]]<- sim_fcs_truth$observed_x_train[,,person_index]
# }
#
# mcem_parameter_setup_result<-
# mcem_parameter_setup(p = 100, k = 4, n = 17, q = 8,
# obs_time_num = obs_time_num,
# obs_time_index = obs_time_index,
# a_person = a_person,
# col_person_index = col_person_index,
# y_init = sim_fcs_init$y_init,
# a_init = sim_fcs_init$a_init_2,
# z_init = sim_fcs_init$z_init_2,
# phi_init = sim_fcs_init$phi_init,
# a_full = sim_fcs_truth$a_full,
# x = observed_x_train_regular_8,
# train_index = (1:8))
#
## ----eval = FALSE-------------------------------------------------------------
#
# mcem_algorithm_result<-
# mcem_algorithm(ind_x = TRUE,
# x = observed_x_train_regular_8,
# mcem_parameter_setup_result = mcem_parameter_setup_result)
#
## ----fig1, fig.height = 4, fig.width=7----------------------------------------
old_par <- par(no.readonly = TRUE)
par(mfrow = c(2,4))
for (em_index in 1:sim_fcs_results_regular_8$mcem_algorithm_result$index_used){
mcem_cov_plot(sim_fcs_results_regular_8$mcem_algorithm_result$sigmay_record[em_index,,], k = 4, q = 8, title = paste0("MCEM Iteration ", em_index))
}
mcem_cov_plot(sim_fcs_truth$gp_sigmay_truth, k = 4, q = 10, title = "Truth: Correlated Factors")
par(old_par)
## ----eval = FALSE-------------------------------------------------------------
#
# gibbs_after_mcem_diff_initials_result<-
# gibbs_after_mcem_diff_initials(ind_x = TRUE,
# tot_chain = 5,
# mcem_parameter_setup_result = mcem_parameter_setup_result,
# mcem_algorithm_result = mcem_algorithm_result)
#
## ----eval = FALSE-------------------------------------------------------------
#
# tot_chain<- 5
#
# for (chain_index in 1:tot_chain){
#
# gibbs_after_mcem_algorithm(chain_index = chain_index,
# mc_num = 10000,
# burnin = 3000,
# thin_step = 10 ,
# pathname = "path",
# pred_indicator = TRUE,
# pred_time_index = (9:10),
# x = observed_x_train_regular_8,
# gibbs_after_mcem_diff_initials_result = gibbs_after_mcem_diff_initials_result,
# mcem_algorithm_result = mcem_algorithm_result,
# mcem_parameter_setup_result = mcem_parameter_setup_result)
# }
#
## ----eval = FALSE-------------------------------------------------------------
#
# constant_list<- list(num_time_test = 2,
# mc_num = 10000,
# thin_step = 10,
# burnin = 3000,
# pathname = "path",
# p = 100,
# k = 4,
# n = 17,
# q = 8,
# ind_x = TRUE,
# pred_indicator = TRUE)
#
# for (chain_index in 1:tot_chain){
#
# gibbs_after_mcem_load_chains_result<- gibbs_after_mcem_load_chains(chain_index = chain_index,
# gibbs_after_mcem_algorithm_result = constant_list)
#
# save(gibbs_after_mcem_load_chains_result,
# file = paste0("path/chain_", chain_index,"_result.RData"))
# }
#
# gibbs_after_mcem_combine_chains_result<- gibbs_after_mcem_combine_chains(tot_chain = 5,
# gibbs_after_mcem_algorithm_result = constant_list)
#
## ----eval = FALSE-------------------------------------------------------------
# numerics_summary_do_not_need_alignment_result<-
# numerics_summary_do_not_need_alignment(pred_x_truth = sim_fcs_truth$observed_x_pred_reformat,
# pred_x_truth_indicator = TRUE,
# gibbs_after_mcem_combine_chains_result = gibbs_after_mcem_combine_chains_result)
#
## -----------------------------------------------------------------------------
sim_fcs_results_regular_8$numerics_summary_do_not_need_alignment_result$convergence_summary
## -----------------------------------------------------------------------------
pred_result_overview<- matrix(c(sim_fcs_results_regular_8$numerics_summary_do_not_need_alignment_result$pred_x_result$mae_using_median_est,
sim_fcs_results_regular_8$numerics_summary_do_not_need_alignment_result$pred_x_result$mean_width_interval,
sim_fcs_results_regular_8$numerics_summary_do_not_need_alignment_result$pred_x_result$proportion_of_within_interval_biomarkers),
nrow = 3, ncol = 1)
rownames(pred_result_overview)<- c("Mean Absolute Error", "Mean Width of Intervals", "Proportion of Coverage")
colnames(pred_result_overview)<- "Value"
pred_result_overview
## -----------------------------------------------------------------------------
time_index<- 1
person_index<- 2
gene_index<- 3
pred_result_specific_gene<-
matrix(c(sim_fcs_results_regular_8$numerics_summary_do_not_need_alignment_result$pred_x_result$pred_lower_quantile[person_index, gene_index, time_index],
sim_fcs_results_regular_8$numerics_summary_do_not_need_alignment_result$pred_x_result$pred_median_quantile[person_index, gene_index, time_index],
sim_fcs_results_regular_8$numerics_summary_do_not_need_alignment_result$pred_x_result$pred_upper_quantile[person_index, gene_index, time_index],
sim_fcs_truth$observed_x_pred_reformat[person_index, gene_index, time_index]),
nrow = 4, ncol = 1)
rownames(pred_result_specific_gene)<- c("2.5% Quantile", "Point Estimate", "97.5% Quantile", "Truth")
colnames(pred_result_specific_gene)<- "Value"
pred_result_specific_gene
## ----eval = FALSE-------------------------------------------------------------
# numerics_summary_need_alignment_result<-
# numerics_summary_need_alignment(gibbs_after_mcem_combine_chains_result = gibbs_after_mcem_combine_chains_result)
## -----------------------------------------------------------------------------
sim_fcs_results_regular_8$numerics_summary_need_alignment_result$convergence_summary
## ----fig2, fig.height = 6, fig.width=7----------------------------------------
factor_loading_heatmap(sim_fcs_results_regular_8$numerics_summary_need_alignment_result$reordered_summary$big_l,
heatmap_title = "Estimated Factor Loadings")
## ----fig3, fig.height = 4, fig.width=7----------------------------------------
# only display trajectories at training time points
factor_score_trajectory(sim_fcs_results_regular_8$numerics_summary_need_alignment_result$reordered_summary$latent_y[(1:8),,],
factor_index = 1,
person_index = 1,
trajectory_title = paste0("Estimated Trajectory of Factor 1 for Person 1"))
## -----------------------------------------------------------------------------
q<- 8
k<- 4
n<- 17
a_train<- sim_fcs_truth$a_full[(1:q)]
h3n2_data<- list()
list_temp <- vector("list", k)
for (list_index in 1:k){
list_temp[[list_index]]<- a_train
}
h3n2_data$input<- list_temp
fcs_sigma_y_init_truth_for_train_data<- GPFDA::mgpCovMat(Data=h3n2_data, hp=sim_fcs_truth$gp_hp_truth)
# rescale truth
d_matrix<- diag(sqrt(diag(fcs_sigma_y_init_truth_for_train_data)))
d_matrix_inv<- solve(d_matrix)
fcs_real_y_rescaled<- array(0, dim = c(q,k,n))
for (person_index in 1:n){
fcs_real_y_rescaled[,,person_index]<- matrix(d_matrix_inv%*%as.numeric(sim_fcs_truth$real_y[1:q,,person_index]),
nrow = q,
ncol = k)
}
# compare
mae_regular_8<- mean(abs(sim_fcs_results_regular_8$numerics_summary_need_alignment_result$reordered_summary$latent_y[(1:q),,] - fcs_real_y_rescaled))
mae_regular_8