inst/doc/bsfadgp_regular_data_example.R

In DGP4LCF: Dependent Gaussian Processes for Longitudinal Correlated Factors

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