This notebook present a way of performing the SmoothPLS algorithm for
a Scalar Functional Data.
It shows that Smooth PLS over perform FPLS with a time vector
regul_time without a lot of points, and even with a lot of points.
Generate_X_df
For SFD for X_df two new arguments are important: the noise added to
the signal and the seed for repeatability.
df = generate_X_df(nind = nind, start = start,end = end, curve_type = 'num',
noise_sd = 0.15, seed = 123)
head(df)
#> id time value
#> 1 1 0 1.2591356
#> 2 1 1 0.7193661
#> 3 1 2 1.0470795
#> 4 1 3 0.7297674
#> 5 1 4 0.6105385
#> 6 1 5 0.5351208
# Visualisation
ggplot(df, 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")
generate_Y_df
plot(regul_time_0, beta_real_func(regul_time_0))
WARNING! Here we have to regularize X(t) before evaluating Y!
Y_df = generate_Y_df(df = df_regul, curve_type = 'num',
beta_real_func_or_list = beta_real_func,
beta_0_real = beta_0_real, NotS_ratio = NotS_ratio,
seed = 123)
Y = Y_df$Y_noised
head(Y_df)
#> id Y_real Y_noised
#> 1 1 -201.7150 -203.2424
#> 2 2 -197.9957 -198.6230
#> 3 3 -198.3660 -194.1184
#> 4 4 -199.5417 -199.3496
#> 5 5 -197.7755 -197.4232
#> 6 6 -198.4119 -193.7382
oldpar <- par(mfrow=c(1, 2))
hist(Y_df$Y_real)
hist(Y_df$Y_noised)
Naive PLS
naive_pls_obj = naivePLS(df_list = df, Y = Y_df$Y_noised,
regul_time_obj = regul_time,
curve_type_obj = 'num',
id_col_obj = 'id', time_col_obj = 'time',
print_steps = TRUE, plot_rmsep = TRUE,
print_nbComp = TRUE, plot_reg_curves = TRUE)
#> ### Naive PLS ###
#> => Input format assertions.
#> => Input format assertions OK.
#> => Data formatting.
#> => PLS model.
#> [1] "Optimal number of PLS components : 4"
Regression coefficients
for(i in 2:length(naive_pls_obj$reg_obj)){
plot(naive_pls_obj$reg_obj[[i]])
title(naive_pls_obj$curves_names[i-1])
}
plot(naive_pls_obj$reg_obj$Num_1,
xlab='regul_time', ylab='beta coefficient value',
main="naive pls regression coefficients")
lines(x = regul_time, y = naive_pls_obj$reg_obj$Num_1[, 2], col='blue')
lines(x = regul_time_0, y=beta_real_func(start:end, end), col='red')
FPLS
fpls_obj = funcPLS(df_list = df, Y = Y_df$Y_noised,
basis_obj = basis,
curve_type_obj = 'num',
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)
#> ### Functional PLS ###
#> => Input format assertions.
#> => Input format assertions OK.
#> => Building alpha matrix.
#> => Building curve names.
#> ==> Evaluating alpha for : NumFD_1.
#> => Evaluate metrix and root_metric.
#> => plsr(Y ~ trans_alphas).
#> Optimal number of PLS components : 2 .
#> => Build Intercept and regression curves for optimal number of components.
#> ==> Build regression curve for : NumFD_1
#Smooth PLS
spls_obj = smoothPLS(df_list = df, Y = Y_df$Y_noised, basis_obj = basis,
orth_obj = TRUE, curve_type_obj = curve_type,
regul_time_obj = regul_time,
int_mode = int_mode, print_steps = TRUE, plot_rmsep = TRUE,
print_nbComp = TRUE, plot_reg_curves = TRUE)
#> ### Smooth PLS ###
#> ## Use parallelization in case of heavy computational load. ##
#> ## Threshold can be manualy adjusted : (default 2500) ##
#> ## >options(SmoothPLS.parallel_threshold = 500) ##
#> => Input format assertions.
#> => Input format assertions OK.
#> => Orthonormalize basis.
#> => Data objects formatting.
#> => Evaluate Lambda matrix.
#> ==> Lambda for : NumFD_1_value.
#> => PLSR model.
#> => Optimal number of PLS components : 9
#> => Evaluate SmoothPLS functions and <w_i, p_j> coef.
#> => Build regression functions and intercept.
#> ==> Build regression curve for : NumFD_1_value
Delta
delta_0 = spls_obj$reg_obj$Intercept
print(delta_0)
#> [1] -13.72161
delta = spls_obj$reg_obj$NumFD_1
y_lim = eval_max_min_y(f_list = list(beta_real_func, delta),
regul_time = regul_time)
plot(regul_time, beta_real_func(regul_time), type='l', xlab="Beta_t",
ylim = y_lim)
plot(delta, add=TRUE, col='blue')
#> [1] "done"
Curve comparison
cat("curve_1 : smooth PLS regression curve.\n")
#> curve_1 : smooth PLS regression curve.
cat("curve_2 : functional PLS regression curve.\n")
#> curve_2 : functional PLS regression curve.
cat("curve_3 : naive PLS regression coefficients\n")
#> curve_3 : naive PLS regression coefficients
evaluate_curves_distances(real_f = beta_real_func,
regul_time = regul_time,
fun_fd_list = list(
delta,
fpls_obj$reg_obj$NumFD_1,
approxfun(x = regul_time,
y = naive_pls_obj$opti_reg_coef)
)
)
#> [1] "real_f -> curve_1 / INPROD : 5.42764621642823 / DIST : 3.50987893386841"
#> [1] "real_f -> curve_2 / INPROD : 6.66383089607853 / DIST : 3.40650741798858"
#> [1] "real_f -> curve_3 / INPROD : 247.075818435331 / DIST : 68.8619172938924"
y_lim = eval_max_min_y(f_list = list(spls_obj$reg_obj$NumFD_1,
fpls_obj$reg_obj$NumFD_1,
approxfun(x = regul_time,
y = naive_pls_obj$opti_reg_coef),
beta_real_func ),
regul_time = regul_time_0)
plot(regul_time_0, beta_real_func(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)(regul_time_0),
col = 'green')
title(paste0(names(spls_obj$reg_obj)[2], " regression curves"))
plot(spls_obj$reg_obj$NumFD_1, col = 'blue', add = TRUE)
#> [1] "done"
plot(fpls_obj$reg_obj$NumFD_1, col = 'red', add = TRUE)
#> [1] "done"
legend("bottomleft",
legend = c("Real curve", "NaivePLS coef",
"SmoothPLS reg curve", "FunctionalPLS reg curve"),
col = c("black", "green", "blue", "red"),
lty = 1,
lwd = 1)
y_lim = eval_max_min_y(f_list = list(beta_real_func,
spls_obj$reg_obj$NumFD_1,
fpls_obj$reg_obj$NumFD_1),
regul_time = regul_time)
plot(regul_time_0, beta_real_func(regul_time_0), type='l', xlab="Beta_t",
ylim = y_lim, lwd = 2)
plot(spls_obj$reg_obj$NumFD_1, add=TRUE, col='blue', lwd = 1)
#> [1] "done"
plot(fpls_obj$reg_obj$NumFD_1, add=TRUE, col='red', lwd = 1)
#> [1] "done"
legend("bottomleft",
legend = c("Real curve",
"SmoothPLS reg curve", "FunctionalPLS reg curve"),
col = c("black", "blue", "red"),
lty = 1,
lwd = 1)
evaluate_curves_distances(real_f = beta_real_func,
regul_time = regul_time,
fun_fd_list = list(spls_obj$reg_obj$NumFD_1,
fpls_obj$reg_obj$NumFD_1))
#> [1] "real_f -> curve_1 / INPROD : 5.42764621642823 / DIST : 3.50987893386841"
#> [1] "real_f -> curve_2 / INPROD : 6.66383089607853 / DIST : 3.40650741798858"
Plot the second derivatives for the smoothness :
plot((fda::deriv.fd(expr = spls_obj$reg_obj$NumFD_1, Lfdobj = 2)),
col='blue')
#> [1] "done"
plot((fda::deriv.fd(expr = fpls_obj$reg_obj$NumFD_1, Lfdobj = 2)),
add=TRUE, col='red')
#> [1] "done"
title("Second derivatives")
legend("bottomleft",
legend = c("delta_SmoothPLS''", "delta_FPLS''"),
col = c("blue", "red"),
lty = 1,
lwd = 1)
Results
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"))
#> [1] "There is 20% of noised in Y"
Train Set
Y_train = Y_df$Y_noised
# Naive
Y_hat = predict(naive_pls_obj$plsr_model,
ncomp = naive_pls_obj$nbCP_opti,
newdata = as.matrix(df_wide[, -c(1)]))
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 = df,
delta_list = fpls_obj$reg_obj,
curve_type_obj = curve_type,
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 = df,
delta_list = spls_obj$reg_obj,
curve_type_obj = curve_type,
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
#> PRESS RMSE MAE R2 var_Y nb_cp
#> NaivePLS 447.7783 2.116077 1.658857 0.8713332 87.13332 4
#> FPLS 2958.4590 5.439172 4.598839 0.1499018 14.96341 2
#> SmoothPLS 15533.5802 12.463378 11.223560 -3.4634956 13.79539 9
Test Set
Test set creation
nind_test = number_of_test_id(TTRatio = TTRatio, nind = nind)
df_test = generate_X_df(nind = nind_test, start = start, end = end,
curve_type = curve_type, lambda_0 = lambda_0, lambda_1 = lambda_1,
prob_start = prob_start)
Y_df_test = generate_Y_df(df_test, curve_type = curve_type,
beta_real_func_or_list = beta_real_func,
beta_0_real, NotS_ratio,
id_col = 'id', time_col = 'time')
df_test_regul = regularize_time_series(df_test, time_seq = regul_time,
curve_type = curve_type)
df_test_wide = convert_to_wide_format(df_test_regul)
Evaluation
Y_test = Y_df_test$Y_noised
# Naive
Y_hat = predict(naive_pls_obj$plsr_model,
ncomp = naive_pls_obj$nbCP_opti,
newdata = as.matrix(df_test_wide[, -c(1)]))
test_results["NaivePLS", ] = evaluate_results(Y_test, Y_hat)
# FPLS
Y_hat_fpls = smoothPLS_predict(df_predict_list = df_test,
delta_list = fpls_obj$reg_obj,
curve_type_obj = curve_type,
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 = df_test,
delta_list = list(delta_0, delta),
curve_type_obj = curve_type,
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
#> PRESS RMSE MAE R2 var_Y nb_cp
#> NaivePLS 222.84261 2.985583 2.3244464 -2.813555 232.68890 4
#> FPLS 31.38556 1.120456 0.9331762 0.462892 81.55013 2
#> SmoothPLS 3881.27654 12.459978 12.4399294 -65.421150 76.32215 9
