jaromilfrossard / permuco

R package for permutation test in linear model with nuisances variables
https://jaromilfrossard.github.io/permuco/
11 stars 2 forks source link

`compute_tfce()` would fail to calculate the correct TFCE values when the original statistics are too small #15

Closed psychelzh closed 2 months ago

psychelzh commented 2 months ago

For example, I generated some random numbers. At first it could give the correct values:

set.seed(42)
stat_map <- as.matrix(rnorm(100))
permuco::compute_tfce(stat_map, ndh = 100)
#> $main
#>      statistic pvalue
#> [1,] 0.5066432   0.17
#> 
#> $distribution
#>   [1] 0.5066432 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.5066432
#>   [8] 0.0000000 1.9848038 0.0000000 0.5066432 1.9848038 0.0000000 0.0000000
#>  [15] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.5066432 0.0000000
#>  [22] 0.0000000 0.0000000 0.5066432 0.5066432 0.0000000 0.0000000 0.0000000
#>  [29] 0.0000000 0.0000000 0.0000000 0.0000000 0.5066432 0.0000000 0.0000000
#>  [36] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
#>  [43] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.5066432 0.0000000
#>  [50] 0.0000000 0.0000000 0.0000000 0.5066432 0.0000000 0.0000000 0.0000000
#>  [57] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
#>  [64] 0.5066432 0.0000000 0.5066432 0.0000000 0.5066432 0.0000000 0.0000000
#>  [71] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
#>  [78] 0.0000000 0.0000000 0.0000000 0.5066432 0.0000000 0.0000000 0.0000000
#>  [85] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.5066432
#>  [92] 0.0000000 0.0000000 0.5066432 0.0000000 0.0000000 0.0000000 0.0000000
#>  [99] 0.0000000 0.0000000
#> 
#> $E
#> [1] 0.5
#> 
#> $H
#> [1] 1
#> 
#> $dhi
#>   [1] 0.03612261 0.06599097 0.09585932 0.12572768 0.15559604 0.18546440
#>   [7] 0.21533276 0.24520112 0.27506947 0.30493783 0.33480619 0.36467455
#>  [13] 0.39454291 0.42441127 0.45427962 0.48414798 0.51401634 0.54388470
#>  [19] 0.57375306 0.60362142 0.63348977 0.66335813 0.69322649 0.72309485
#>  [25] 0.75296321 0.78283157 0.81269992 0.84256828 0.87243664 0.90230500
#>  [31] 0.93217336 0.96204172 0.99191007 1.02177843 1.05164679 1.08151515
#>  [37] 1.11138351 1.14125187 1.17112022 1.20098858 1.23085694 1.26072530
#>  [43] 1.29059366 1.32046202 1.35033037 1.38019873 1.41006709 1.43993545
#>  [49] 1.46980381 1.49967217 1.52954052 1.55940888 1.58927724 1.61914560
#>  [55] 1.64901396 1.67888232 1.70875067 1.73861903 1.76848739 1.79835575
#>  [61] 1.82822411 1.85809247 1.88796082 1.91782918 1.94769754 1.97756590
#>  [67] 2.00743426 2.03730262 2.06717097 2.09703933 2.12690769 2.15677605
#>  [73] 2.18664441 2.21651277 2.24638112 2.27624948 2.30611784 2.33598620
#>  [79] 2.36585456 2.39572292 2.42559127 2.45545963 2.48532799 2.51519635
#>  [85] 2.54506471 2.57493307 2.60480142 2.63466978 2.66453814 2.69440650
#>  [91] 2.72427486 2.75414322 2.78401157 2.81387993 2.84374829 2.87361665
#>  [97] 2.90348501 2.93335337 2.96322172 2.99309008

Then I rescale the numbers to smaller values, and the calculations failed.

permuco::compute_tfce(stat_map * 1e-3, ndh = 100)
#> $main
#>      statistic pvalue
#> [1,]         0      1
#> 
#> $distribution
#>   [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
#>  [38] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
#>  [75] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
#> 
#> $E
#> [1] 0.5
#> 
#> $H
#> [1] 1
#> 
#> $dhi
#>   [1] 3.612261e-05 6.599097e-05 9.585932e-05 1.257277e-04 1.555960e-04
#>   [6] 1.854644e-04 2.153328e-04 2.452011e-04 2.750695e-04 3.049378e-04
#>  [11] 3.348062e-04 3.646745e-04 3.945429e-04 4.244113e-04 4.542796e-04
#>  [16] 4.841480e-04 5.140163e-04 5.438847e-04 5.737531e-04 6.036214e-04
#>  [21] 6.334898e-04 6.633581e-04 6.932265e-04 7.230948e-04 7.529632e-04
#>  [26] 7.828316e-04 8.126999e-04 8.425683e-04 8.724366e-04 9.023050e-04
#>  [31] 9.321734e-04 9.620417e-04 9.919101e-04 1.021778e-03 1.051647e-03
#>  [36] 1.081515e-03 1.111384e-03 1.141252e-03 1.171120e-03 1.200989e-03
#>  [41] 1.230857e-03 1.260725e-03 1.290594e-03 1.320462e-03 1.350330e-03
#>  [46] 1.380199e-03 1.410067e-03 1.439935e-03 1.469804e-03 1.499672e-03
#>  [51] 1.529541e-03 1.559409e-03 1.589277e-03 1.619146e-03 1.649014e-03
#>  [56] 1.678882e-03 1.708751e-03 1.738619e-03 1.768487e-03 1.798356e-03
#>  [61] 1.828224e-03 1.858092e-03 1.887961e-03 1.917829e-03 1.947698e-03
#>  [66] 1.977566e-03 2.007434e-03 2.037303e-03 2.067171e-03 2.097039e-03
#>  [71] 2.126908e-03 2.156776e-03 2.186644e-03 2.216513e-03 2.246381e-03
#>  [76] 2.276249e-03 2.306118e-03 2.335986e-03 2.365855e-03 2.395723e-03
#>  [81] 2.425591e-03 2.455460e-03 2.485328e-03 2.515196e-03 2.545065e-03
#>  [86] 2.574933e-03 2.604801e-03 2.634670e-03 2.664538e-03 2.694406e-03
#>  [91] 2.724275e-03 2.754143e-03 2.784012e-03 2.813880e-03 2.843748e-03
#>  [96] 2.873617e-03 2.903485e-03 2.933353e-03 2.963222e-03 2.993090e-03

Created on 2024-08-25 with reprex v2.1.0

jaromilfrossard commented 2 months ago

Hello,

Thanks you for reporting the bug. I think there is indeed something wrong with the approximation of the integral of the TFCE, when the values of the statistics are below 1. The issue might come from the cpp code.

#zero
permuco:::tfce_distribution(distribution = matrix(.9),E=0.5,H=1,dh=.1,dhi=c(0.6,.7,.8,.9)) 

#not zero
permuco:::tfce_distribution(distribution = matrix(1.1),E=0.5,H=1,dh=.1,dhi=c(0.6,.7,.8,.9,1,1.1))
permuco:::tfce_distribution(distribution = matrix(1.1),E=0.5,H=1,dh=.1,dhi=c(.8,.9,1,1.1))

#also zero !?
permuco:::tfce_distribution(distribution = matrix(1.1),E=0.5,H=1,dh=.1,dhi=c(1,1.1))

I double check the cpp code, and I think:

https://github.com/jaromilfrossard/permuco/blob/19ec0f916f5a090a532cf423adb584d6a15369b3/src/tfce.cpp#L61C1-L68C4

jaromilfrossard commented 2 months ago

In line

jaromilfrossard commented 2 months ago

Fixed in 1.1.3 This bug should not have affected the results of TFCE with the usually statistics as they usually have greater than 1.