SugiharaLab / rEDM

Applications of Empirical Dynamic Modeling from Time Series
Other
117 stars 43 forks source link

The Multiview function does not select the best views for the average prediction #43

Closed duriah closed 3 years ago

duriah commented 3 years ago

There following issue is shown here: The Multiview function does not select the k best views for the average prediction

The issue is shown in the example code below:

library(rEDM)
data(block_3sp)

L.3views = Multiview(dataFrame = block_3sp, lib = "1 99", pred = "105 190", E=3,
                D = 3, columns = "x_t y_t z_t", target = "x_t", multiview = 3)

L.10views = Multiview(dataFrame = block_3sp, lib = "1 99", pred = "105 190", E=3,
                             D = 3, columns = "x_t y_t z_t", target = "x_t", multiview = 10)

L.3views$View[,1:7]

  col_1 col_2 col_3   name_1   name_2   name_3    rho
1     1     2     7 x_t(t-0) x_t(t-1) z_t(t-0) 0.9208
2     1     2     6 x_t(t-0) x_t(t-1) y_t(t-2) 0.8677
3     1     2     3 x_t(t-0) x_t(t-1) x_t(t-2) 0.9319

L.10views$View[,1:7]

   col_1 col_2 col_3   name_1   name_2   name_3    rho
1      1     2     7 x_t(t-0) x_t(t-1) z_t(t-0) 0.9208
2      1     2     6 x_t(t-0) x_t(t-1) y_t(t-2) 0.8677
3      1     2     3 x_t(t-0) x_t(t-1) x_t(t-2) 0.9319
4      1     2     8 x_t(t-0) x_t(t-1) z_t(t-1) 0.9183
5      1     7     9 x_t(t-0) z_t(t-0) z_t(t-2) 0.8858
6      1     4     9 x_t(t-0) y_t(t-0) z_t(t-2) 0.7774
7      1     3     7 x_t(t-0) x_t(t-2) z_t(t-0) 0.8724
8      1     2     5 x_t(t-0) x_t(t-1) y_t(t-1) 0.9017
9      1     2     4 x_t(t-0) x_t(t-1) y_t(t-0) 0.8805
10     1     3     6 x_t(t-0) x_t(t-2) y_t(t-2) 0.8463

If we compare the two View dataframes above, we see that in the model in which only 3 views are used not the best 3 views are selected (based on rho as described here https://science.sciencemag.org/content/353/6302/922.abstract), because in the model with 10 views we see that there would have been better views to select.

For instance, row 4 should have been selected over row 2.

It follows that the multiview function does no actually select the k best views to perform the average forecasting. It is unclear how it selects them.

Thank you and best regards, Uriah

sessionInfo(package = "rEDM")

R version 4.0.3 (2020-10-10)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 18363)

Matrix products: default

locale:
[1] LC_COLLATE=English_Switzerland.1252  LC_CTYPE=English_Switzerland.1252    LC_MONETARY=English_Switzerland.1252
[4] LC_NUMERIC=C                         LC_TIME=English_Switzerland.1252    

attached base packages:
character(0)

other attached packages:
[1] rEDM_1.7.3

loaded via a namespace (and not attached):
 [1] compiler_4.0.3  graphics_4.0.3  htmltools_0.5.0 tools_4.0.3     utils_4.0.3     yaml_2.2.1      grDevices_4.0.3
 [8] Rcpp_1.0.5      stats_4.0.3     datasets_4.0.3  rmarkdown_2.6   knitr_1.30      methods_4.0.3   xfun_0.20      
[15] digest_0.6.27   rlang_0.4.10    base_4.0.3      evaluate_0.14  
SoftwareLiteracy commented 3 years ago

Dear Uriah,

Thank you for your use and analysis of rEDM. I have reproduced your examples, and I think I see what's happening. I see it as a disconnect between the documentation and function behavior. Thank you for pointing this out!

I believe the function is properly reporting the top multiview combinations, the disconnect comes from the value of rho that is reported for each view. More precisely, the value of rho that is used to select the top views may not be the same as that used in the final evaluation of the view, which is reported in the View output table.

This is because the default behavior is to select the top views based on in-sample library and prediction forecast. That is, the default sets pred = lib when evaluating/selecting the top views. I believe this is what is described in the Ye paper. To be pedantic, if the user specifies lib = "1 10" pred = "11 20", then the function evaluates the views with lib = "1 10" pred = "1 10", in-sample prediction, but the final output rho is forecast using the specified pred = "11 20".

You rightly say: why the disconnect? I don't have a good answer. We need to think about this and how to address it. I suspect it is a good thing to have this flexibility, that the documentation/example should be clarified.

In fact, this disconnect (imo) is one reason the trainLib argument was added to Multiview. When True (default), in-sample prediction sets are used in the selection of the top views, as described above. When trainLib = False, then the multiview evaluations are performed according to the specified lib and pred arguments. This allows one to have multiview select views based on out-of-sample predictions, rather than purely in-sample predictions.

To address/affirm this we can use the default trainLib = True and set lib = pred, ala:

> Multiview( dataFrame = block_3sp, lib = "1 99", pred = "1 99", E = 3, D = 3, columns = "x_t y_t z_t", target = "x_t", multiview = 3, trainLib = TRUE ) $ View
  col_1 col_2 col_3   name_1   name_2   name_3    rho    MAE   RMSE
1     1     2     7 x_t(t-0) x_t(t-1) z_t(t-0) 0.8875 0.3303 0.4302
2     1     2     6 x_t(t-0) x_t(t-1) y_t(t-2) 0.8834 0.3343 0.4385
3     1     2     3 x_t(t-0) x_t(t-1) x_t(t-2) 0.8692 0.3329 0.4635

which yields the same result as trainLib = FALSE iflib = pred`

> Multiview( dataFrame = block_3sp, lib = "1 99", pred = "1 99", E = 3, D = 3, columns = "x_t y_t z_t", target = "x_t", multiview = 3, trainLib = FALSE ) $ View
  col_1 col_2 col_3   name_1   name_2   name_3    rho    MAE   RMSE
1     1     2     7 x_t(t-0) x_t(t-1) z_t(t-0) 0.8875 0.3303 0.4302
2     1     2     6 x_t(t-0) x_t(t-1) y_t(t-2) 0.8834 0.3343 0.4385
3     1     2     3 x_t(t-0) x_t(t-1) x_t(t-2) 0.8692 0.3329 0.4635

Both of these results select the same top views as the original example, which used default trainLib = TRUE, thereby evaluated the views on lib = pred = 1:99, but reported final rho using different lib = "1 99" and pred = "105 190":

> Multiview( dataFrame = block_3sp, lib = "1 99", pred = "105 190", E = 3, D = 3, columns = "x_t y_t z_t", target = "x_t", multiview = 3, trainLib = TRUE ) $ View
  col_1 col_2 col_3   name_1   name_2   name_3    rho    MAE   RMSE
1     1     2     7 x_t(t-0) x_t(t-1) z_t(t-0) 0.9208 0.2485 0.3164
2     1     2     6 x_t(t-0) x_t(t-1) y_t(t-2) 0.8677 0.3294 0.4113
3     1     2     3 x_t(t-0) x_t(t-1) x_t(t-2) 0.9319 0.2277 0.2934

As long as lib == pred, the reported values of rho in the output table will reflect the values used in view selection:

 Multiview( dataFrame = block_3sp, lib = "20 100", pred = "20 100", E = 3, D = 3, columns = "x_t y_t z_t", target = "x_t", multiview = 3, trainLib = TRUE ) $ View
  col_1 col_2 col_3   name_1   name_2   name_3    rho    MAE   RMSE
1     1     2     3 x_t(t-0) x_t(t-1) x_t(t-2) 0.8820 0.3104 0.4226
2     1     2     7 x_t(t-0) x_t(t-1) z_t(t-0) 0.8749 0.3309 0.4306
3     1     3     7 x_t(t-0) x_t(t-2) z_t(t-0) 0.8682 0.3363 0.4411
> 
> 
> Multiview( dataFrame = block_3sp, lib = "20 100", pred = "20 100", E = 3, D = 3, columns = "x_t y_t z_t", target = "x_t", multiview = 7, trainLib = TRUE ) $ View
  col_1 col_2 col_3   name_1   name_2   name_3    rho    MAE   RMSE
1     1     2     3 x_t(t-0) x_t(t-1) x_t(t-2) 0.8820 0.3104 0.4226
2     1     2     7 x_t(t-0) x_t(t-1) z_t(t-0) 0.8749 0.3309 0.4306
3     1     3     7 x_t(t-0) x_t(t-2) z_t(t-0) 0.8682 0.3363 0.4411
4     1     2     9 x_t(t-0) x_t(t-1) z_t(t-2) 0.8525 0.3514 0.4645
5     1     2     8 x_t(t-0) x_t(t-1) z_t(t-1) 0.8518 0.3362 0.4657
6     1     2     4 x_t(t-0) x_t(t-1) y_t(t-0) 0.8394 0.3485 0.4896
7     1     4     7 x_t(t-0) y_t(t-0) z_t(t-0) 0.8365 0.3749 0.4900
> 

I agree that this situation is not stated in the documents.

duriah commented 3 years ago

Thank you for your detailed answer!

The matter is now pretty clear. I think the disconnect makes the function a bit less transparent, so it might be a good idea to store the in-sample rho-values of the views somewhere accessible or to mention this behaviour in the documentation. But it is by no means anything urgent.

Thanks! Uriah