## ----include = FALSE---------------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ## ----setup-------------------------------------------------------------------- library(SmoothPLS) library(ggplot2) library(pls) ## ----------------------------------------------------------------------------- N_states = 3 nind = 50 # number of individuals (train set) start = 0 # First time end = 100 # end time curve_type = 'cat' TTRatio = 0.2 # Train Test Ratio means we have floor(nind*TTRatio/(1-TTRatio)) NotS_ratio = 0.2 # noise variance over total variance for Y beta_0_real=65.4321 # Intercept value for the link between X(t) and Y nbasis = 7 # number of basis functions norder = 4 # 4 for cubic splines basis regul_time = seq(start, end, 1) # regularisation time sequence regul_time_0 = seq(start, end, 1) all_curves_in_pls = TRUE #Keep all or (all - 1) states in PLS method if FALSE. int_mode = 1 # in case of integration errors. ## ----------------------------------------------------------------------------- df_cfd_os = generate_X_df(nind = nind, start = start, end = end, curve_type = 'cat') ## ----------------------------------------------------------------------------- head(df_cfd_os) ## ----------------------------------------------------------------------------- beta_func_cfd_os = beta_5_real_func Y_df_cfd_os = generate_Y_df(df_cfd_os, curve_type = 'cat', beta_real_func_or_list = beta_func_cfd_os, beta_0_real = beta_0_real, NotS_ratio, id_col = 'id', time_col = 'time') Y_cfd_os = Y_df_cfd_os$Y_noised ## ----------------------------------------------------------------------------- plot(regul_time, beta_5_real_func(regul_time)) ## ----------------------------------------------------------------------------- # Initialized the lambdas values lambdas = lambda_determination(N_states) lambdas ## ----------------------------------------------------------------------------- # Initialized the transition matrix transition_df = transfer_probabilities(N_states) transition_df ## ----------------------------------------------------------------------------- df_cfd = generate_X_df_multistates(nind = nind, N_states, start, end, lambdas, transition_df) head(df_cfd) ## ----------------------------------------------------------------------------- plot_CFD_individuals(df_cfd) ## ----------------------------------------------------------------------------- plot_CFD_individuals(df_cfd, by_cfda = TRUE) ## ----------------------------------------------------------------------------- ##### beta_real ##### ###### beta_0_real ###### beta_0_real ## ----------------------------------------------------------------------------- beta_func_list = beta_list_generation(N_states = N_states) ## ----------------------------------------------------------------------------- for(i in 1:N_states){ plot(0:end, beta_func_list[[i]](0:end, end_time = end), ylab=paste0("Beta(t) n°=", i), type = 'l') title(paste0("Beta(t) n°=", i)) } ## ----------------------------------------------------------------------------- df_processed = cat_data_to_indicator(data = df_cfd) ## ----------------------------------------------------------------------------- Y_df_cfd = generate_Y_df(df = df_processed, curve_type = curve_type, beta_real_func_or_list = beta_func_list, beta_0_real = beta_0_real, NotS_ratio = NotS_ratio) Y_cfd = Y_df_cfd$Y_noised names(Y_df_cfd) ## ----------------------------------------------------------------------------- head(Y_df_cfd) ## ----------------------------------------------------------------------------- df_sfd = generate_X_df(nind = nind, start = start, end = end, curve_type = 'num', noise_sd = 0.15, seed = 123) head(df_sfd) ## ----------------------------------------------------------------------------- # Visualisation ggplot(df_sfd, aes(x = time, y = value, group = id, color = factor(id))) + geom_line(alpha = 0.8) + labs(title = "Noised cosinus curves", x = "Time", y = "Value", color = "Individual") + theme_minimal() + theme(legend.position = "none") ## ----------------------------------------------------------------------------- df_sfd_regul = regularize_time_series(df_sfd, time_seq = regul_time, curve_type = 'num') df_sfd_regul ## ----------------------------------------------------------------------------- df_sfd_long = convert_to_wide_format(df_sfd_regul) df_sfd_long ## ----------------------------------------------------------------------------- beta_func_sfd = beta_1_real_func Y_df_sfd = generate_Y_df(df = df_sfd_regul, curve_type = 'num', beta_real_func_or_list = beta_func_sfd, beta_0_real = beta_0_real, NotS_ratio = NotS_ratio, seed = 123) Y_sfd = Y_df_sfd$Y_noised head(Y_df_sfd) ## ----------------------------------------------------------------------------- beta_func_list_tot = append(beta_list_generation(N_states = N_states), append(beta_func_sfd, beta_func_cfd_os)) ## ----hist Y noised------------------------------------------------------------ oldpar <- par(mfrow=c(1, 3)) hist(Y_df_sfd$Y_real) hist(Y_df_sfd$Y_noised) plot(regul_time, beta_4_real_func(regul_time)) par(oldpar) ## ----------------------------------------------------------------------------- Y = Y_cfd + Y_sfd + Y_cfd_os ## ----------------------------------------------------------------------------- nind_test = floor(nind*TTRatio/(1-TTRatio)) df_cfd_os_test = generate_X_df(nind = nind_test, start = start, end = end, curve_type = 'cat') ## ----------------------------------------------------------------------------- Y_df_cfd_os_test = generate_Y_df(df_cfd_os_test, curve_type = 'cat', beta_real_func_or_list = beta_5_real_func, beta_0_real = beta_0_real, NotS_ratio = NotS_ratio) Y_cfd_os_test = Y_df_cfd_os_test$Y_noised ## ----------------------------------------------------------------------------- nind_test = floor(nind*TTRatio/(1-TTRatio)) df_cfd_test = generate_X_df_multistates(nind = nind_test, N_states, start, end, lambdas, transition_df) ## ----------------------------------------------------------------------------- df_cfd_test_processed = cat_data_to_indicator(df_cfd_test) Y_df_cfd_test = generate_Y_df(df_cfd_test_processed, curve_type = 'cat', beta_real_func_or_list = beta_func_list, beta_0_real = beta_0_real, NotS_ratio= NotS_ratio) Y_cfd_test = Y_df_cfd_test$Y_noised ## ----------------------------------------------------------------------------- df_sfd_test = generate_X_df(nind = nind_test, start = start, end = end, curve_type = 'num', noise_sd = 0.15, seed = 123) ## ----------------------------------------------------------------------------- df_sfd_regul_test = regularize_time_series(df_sfd_test, time_seq = regul_time, curve_type = 'num') Y_df_sfd_test = generate_Y_df(df = df_sfd_regul_test, curve_type = 'num', beta_real_func_or_list = beta_func_sfd, beta_0_real = beta_0_real, NotS_ratio = NotS_ratio, seed = 123) Y_sfd_test = Y_df_sfd_test$Y_noised ## ----------------------------------------------------------------------------- basis = create_bspline_basis(start, end, nbasis, norder) #basis = fda::create.fourier.basis(c(start,end), nbasis = nbasis) # All the states will share the same basis. basis_list = obj_list_creation(N_rep = N_states, obj = basis) plot(basis, main=paste0(nbasis, " ", basis$type," functions basis")) ## ----------------------------------------------------------------------------- spls_obj = smoothPLS(df_list = list(df_cfd, df_sfd, df_cfd_os), Y = Y, basis_obj = basis, regul_time_obj = regul_time, curve_type_obj = list('cat', 'num', 'cat'), orth_obj = c(TRUE, TRUE, TRUE), id_col_obj = 'id', time_col_obj = 'time', print_steps = TRUE, plot_rmsep = TRUE, print_nbComp = TRUE, plot_reg_curves = FALSE) for(k in 1:(length(spls_obj$reg_obj)-1)){ y_lim = eval_max_min_y(f_list = list(beta_func_list_tot[[k]], spls_obj$reg_obj[[k+1]]), regul_time = regul_time_0) plot(regul_time_0, beta_func_list_tot[[k]](regul_time_0), ylim = y_lim) plot(spls_obj$reg_obj[[k+1]], add = TRUE, col = 'blue') title(paste0("Delta - ", names(spls_obj$reg_obj)[k+1])) } ## ----------------------------------------------------------------------------- basis_2 = create_bspline_basis(start, end, (2*nbasis), norder) spls_obj_0 = smoothPLS(df_list = list(df_cfd, df_sfd, df_cfd_os), Y = Y_cfd + Y_sfd, basis_obj = list(basis_2, basis, basis), regul_time_obj = start:end, curve_type_obj = list('cat', 'num', 'cat'), orth_obj = c(TRUE, FALSE, TRUE), id_col_obj = 'id', time_col_obj = 'time', plot_reg_curves = FALSE) for(k in 1:(length(spls_obj_0$reg_obj)-1)){ y_lim = eval_max_min_y(f_list = list(beta_func_list_tot[[k]], spls_obj_0$reg_obj[[k+1]]), regul_time = regul_time_0) plot(regul_time_0, beta_func_list_tot[[k]](regul_time_0), ylim = y_lim) plot(spls_obj_0$reg_obj[[k+1]], add = TRUE, col = 'blue') title(paste0("Delta - ", names(spls_obj_0$reg_obj)[k+1])) } ## ----------------------------------------------------------------------------- Y_hat = smoothPLS_predict(df_predict_list = list(df_cfd, df_sfd, df_cfd_os), delta_list = spls_obj_0$reg_obj, curve_type_obj = list('cat', 'num', 'cat'), id_col_obj = 'id', time_col_obj = 'time', regul_time_obj = 0:100, int_mode = 1, nb_pt = 10, subdivision = 100) hist(Y_hat) ## ----------------------------------------------------------------------------- spls_obj_2 = smoothPLS(df_list = df_sfd, Y = Y_df_sfd$Y_noised, basis_obj = basis, regul_time_obj = start:end, curve_type_obj = 'num', orth_obj = c(TRUE), id_col_obj = 'id', time_col_obj = 'time') for(k in 1:(length(spls_obj_2$reg_obj)-1)){ y_lim = eval_max_min_y(f_list = list(beta_func_sfd, spls_obj_2$reg_obj[[k+1]]), regul_time = regul_time_0) plot(regul_time_0, beta_func_sfd(regul_time_0), ylim = y_lim) plot(spls_obj_2$reg_obj[[k+1]], add = TRUE, col = 'blue') title(paste0("delta - ", names(spls_obj_2$reg_obj)[k+1])) } ## ----------------------------------------------------------------------------- Y_hat_2 = smoothPLS_predict(df_predict_list = df_sfd, delta_list = spls_obj_2$reg_obj, curve_type_obj = 'num', id_col_obj = 'id', time_col_obj = 'time', regul_time_obj = regul_time, int_mode = 1, nb_pt = 10, subdivision = 100) hist(Y_hat_2) ## ----------------------------------------------------------------------------- spls_obj_3 = smoothPLS(df_list = df_cfd, Y = Y_df_cfd$Y_noised, basis_obj = basis, regul_time_obj = regul_time, curve_type_obj = 'cat', orth_obj = c(TRUE), id_col_obj = 'id', time_col_obj = 'time') for(k in 1:(length(spls_obj_3$reg_obj)-1)){ y_lim = eval_max_min_y(f_list = list(beta_func_list_tot[[k]], spls_obj_3$reg_obj[[k+1]]), regul_time = regul_time_0) plot(regul_time_0, beta_func_list_tot[[k]](regul_time_0), ylim = y_lim) plot(spls_obj_3$reg_obj[[k+1]], add = TRUE, col = 'blue') title(paste0("delta - ", names(spls_obj_3$reg_obj)[k+1])) } ## ----------------------------------------------------------------------------- Y_hat_3 = smoothPLS_predict(df_predict_list = df_sfd, delta_list = spls_obj_3$reg_obj, curve_type_obj = 'num', id_col_obj = 'id', time_col_obj = 'time', regul_time_obj = 0:100, int_mode = 1, nb_pt = 10, subdivision = 100) hist(Y_hat_3) ## ----------------------------------------------------------------------------- spls_obj_4 = smoothPLS(df_list = df_cfd_os, Y = Y_df_cfd_os$Y_noised, basis_obj = basis, regul_time_obj = regul_time, curve_type_obj = 'cat', orth_obj = c(TRUE), id_col_obj = 'id', time_col_obj = 'time') for(k in 1:(length(spls_obj_4$reg_obj)-1)){ y_lim = eval_max_min_y(f_list = list(beta_func_cfd_os, spls_obj_4$reg_obj[[k+1]]), regul_time = regul_time_0) plot(regul_time_0, beta_func_cfd_os(regul_time_0), ylim = y_lim) plot(spls_obj_4$reg_obj[[k+1]], add = TRUE, col = 'blue') title(paste0("delta - ", names(spls_obj_4$reg_obj)[k+1])) } ## ----------------------------------------------------------------------------- Y_hat_4 = smoothPLS_predict(df_predict_list = df_sfd, delta_list = spls_obj_4$reg_obj, curve_type_obj = 'num', id_col_obj = 'id', time_col_obj = 'time', regul_time_obj = 0:100, int_mode = 1, nb_pt = 10, subdivision = 100) hist(Y_hat_4) ## ----------------------------------------------------------------------------- naive_pls_obj = naivePLS(df_list = list(df_cfd, df_sfd, df_cfd_os), Y = Y, regul_time_obj = regul_time, curve_type_obj = list('cat', 'num', 'cat'), id_col_obj = 'id', time_col_obj = 'time', print_steps = TRUE, plot_rmsep = TRUE, print_nbComp = TRUE, plot_reg_curves = TRUE) ## ----------------------------------------------------------------------------- fpls_obj = funcPLS(df_list = list(df_cfd, df_sfd, df_cfd_os), Y = Y, basis_obj = basis, curve_type_obj = list('cat', 'num', 'cat'), regul_time_obj = regul_time, id_col_obj = 'id', time_col_obj = 'time', print_steps = TRUE, plot_rmsep = TRUE, print_nbComp = TRUE, plot_reg_curves = TRUE) ## ----------------------------------------------------------------------------- spls_obj = smoothPLS(df_list = list(df_cfd, df_sfd, df_cfd_os), Y = Y, basis_obj = basis, int_mode = 1, regul_time_obj = regul_time, curve_type_obj = list('cat', 'num', 'cat'), orth_obj = c(TRUE, TRUE, TRUE), id_col_obj = 'id', time_col_obj = 'time', print_steps = TRUE, plot_rmsep = TRUE, print_nbComp = TRUE, plot_reg_curves = TRUE) ## ----------------------------------------------------------------------------- # Warning ms_spls_obj$delta_ms_list[[1]] is the intercept! cat("curve_1 : smooth PLS regression curve.\n") cat("curve_2 : functional PLS regression curve.\n") cat("curve_3 : naive PLS regression coefficients\n") N_curves = N_states + 2 # CFD multistate + sfd + cfd_one_state for(i in 1:N_curves){ start = 0 print(paste0("Curve_", i, " : ", names(spls_obj$reg_obj)[i+1])) evaluate_curves_distances(real_f = beta_func_list_tot[[i]], regul_time = regul_time, fun_fd_list = list(spls_obj$reg_obj[[i+1]], fpls_obj$reg_obj[[i+1]], approxfun( x = regul_time, y = naive_pls_obj$opti_reg_coef[ start:(start+length(regul_time) )]) ) ) start = start + length(regul_time) } ## ----------------------------------------------------------------------------- for(i in 1:N_curves){ start = 0 y_lim = eval_max_min_y(f_list = list(spls_obj$reg_ob[[i+1]], fpls_obj$reg_ob[[i+1]], approxfun( x = regul_time, y = naive_pls_obj$opti_reg_coef[ start:(start+length(regul_time))]), beta_func_list_tot[[i]] ), regul_time = regul_time_0) plot(regul_time_0, beta_func_list_tot[[i]](regul_time_0), col = 'black', ylim = y_lim, xlab = 'Time', ylab = 'Value', type = 'l') lines(regul_time_0, approxfun(x = regul_time, y = naive_pls_obj$opti_reg_coef[ start:(start+ length(regul_time))])(regul_time_0), col = 'green') title(paste0(names(spls_obj$reg_obj)[i+1], " regression curves")) plot(spls_obj$reg_obj[[i+1]], col = 'blue', add = TRUE) plot(fpls_obj$reg_obj[[i+1]], col = 'red', add = TRUE) legend("topleft", legend = c("Real curve", "NaivePLS coef", "SmoothPLS reg curve", "FunctionalPLS reg curve"), col = c("black", "green", "blue", "red"), lty = 1, lwd = 1) start = start + length(regul_time) } ## ----------------------------------------------------------------------------- train_results = data.frame(matrix(ncol = 5, nrow = 3)) colnames(train_results) = c("PRESS", "RMSE", "MAE", "R2", "var_Y") rownames(train_results) = c("NaivePLS", "FPLS", "SmoothPLS") test_results = train_results ## ----------------------------------------------------------------------------- print(paste0("There is ", 100*NotS_ratio, "% of noised in Y")) ## ----------------------------------------------------------------------------- Y_train = Y # Naive Y_hat = predict(naive_pls_obj$plsr_model, ncomp = naive_pls_obj$nbCP_opti, newdata = naive_pls_obj$plsr_model$model$`as.matrix(df_mod_wide)`) train_results["NaivePLS", ] = evaluate_results(Y_train, Y_hat) # FPLS Y_hat_fpls = (predict(fpls_obj$plsr_model, ncomp = fpls_obj$nbCP_opti, newdata = fpls_obj$trans_alphas) + fpls_obj$reg_obj$Intercept + mean(Y)) Y_hat_fpls = smoothPLS_predict(df_predict_list = list(df_cfd, df_sfd, df_cfd_os), delta_list = fpls_obj$reg_obj, curve_type_obj = list('cat', 'num', 'cat'), int_mode = int_mode, regul_time_obj = regul_time) train_results["FPLS", ] = evaluate_results(Y_train, Y_hat_fpls) # Smooth PLS Y_hat_spls = smoothPLS_predict(df_predict_list = list(df_cfd, df_sfd, df_cfd_os), delta_list = spls_obj$reg_obj, curve_type_obj = list('cat', 'num', 'cat'), int_mode = int_mode, regul_time_obj = regul_time) train_results["SmoothPLS", ] = evaluate_results(Y_train, Y_hat_spls) train_results["NaivePLS", "nb_cp"] = naive_pls_obj$nbCP_opti train_results["FPLS", "nb_cp"] = fpls_obj$nbCP_opti train_results["SmoothPLS", "nb_cp"] = spls_obj$nbCP_opti ## ----------------------------------------------------------------------------- train_results ## ----------------------------------------------------------------------------- Y_test = Y_cfd_test + Y_sfd_test + Y_cfd_os_test # Naive df_test_wide = naivePLS_formatting(df_list = list(df_cfd_test, df_sfd_test, df_cfd_os_test), regul_time_obj = regul_time, curve_type_obj = list('cat', 'num', 'cat'), id_col_obj = 'id', time_col_obj = 'time') Y_hat = predict(naive_pls_obj$plsr_model, ncomp = naive_pls_obj$nbCP_opti, newdata = as.matrix(df_test_wide)) test_results["NaivePLS", ] = evaluate_results(Y_test, Y_hat) # FPLS Y_hat_fpls = smoothPLS_predict(df_predict_list = list(df_cfd_test, df_sfd_test, df_cfd_os_test), delta_list = fpls_obj$reg_obj, curve_type_obj = list('cat', 'num', 'cat'), int_mode = int_mode, regul_time_obj = regul_time) test_results["FPLS", ] = evaluate_results(Y_test, Y_hat_fpls) # Smooth PLS Y_hat_spls = smoothPLS_predict(df_predict_list = list(df_cfd_test, df_sfd_test, df_cfd_os_test), delta_list = spls_obj$reg_obj, curve_type_obj = list('cat', 'num', 'cat'), int_mode = int_mode, regul_time_obj = regul_time) test_results["SmoothPLS", ] = evaluate_results(Y_test, Y_hat_spls) test_results["NaivePLS", "nb_cp"] = naive_pls_obj$nbCP_opti test_results["FPLS", "nb_cp"] = fpls_obj$nbCP_opti test_results["SmoothPLS", "nb_cp"] = spls_obj$nbCP_opti ## ----------------------------------------------------------------------------- test_results ## ----------------------------------------------------------------------------- train_results test_results ## ----------------------------------------------------------------------------- plot_model_metrics_base(train_results, test_results) ## ----------------------------------------------------------------------------- plot_model_metrics_base(train_results, test_results, models_to_plot = c('FPLS', 'SmoothPLS'))