I have come across a bug in the updated SIMPER code. While running SIMPER on newer versions of R and RStudio, the resulting p-values for taxon with 0 observations is unreasonable.
For example, we have used the Dune data set, grouping by management. The focus will be on the contrast HF_NM.
As you can see in the results below, using an older version of R (3.6.0) Chenalbu and Cirsarve have zero observations, and therefore have a p value of 1.0.
Contrast: HF_NM
Taxon contr sd contr.sd ava avb p
Achimill 0.016555 0.014900 1.1111 1.2 0.3333 0.669323
Agrostol 0.031915 0.028887 1.1048 1.4 2.1667 0.948207
Airaprae 0.012605 0.018243 0.6910 0.0 0.8333 0.286853
Alopgeni 0.024459 0.032399 0.7549 1.6 0.0000 0.940239
Anthodor 0.028717 0.024799 1.1580 1.8 1.3333 0.286853
Bellpere 0.008801 0.013732 0.6409 0.4 0.3333 0.916335
Bromhord 0.012094 0.015169 0.7973 0.8 0.0000 0.836653
Chenalbu 0.000000 0.000000 NaN 0.0 0.0000 1.000000
Cirsarve 0.000000 0.000000 NaN 0.0 0.0000 1.000000
Comapalu 0.010105 0.014556 0.6942 0.0 0.6667 0.215139
Eleopalu 0.032043 0.032315 0.9916 0.8 2.1667 0.494024
Elymrepe 0.029811 0.038676 0.7708 2.0 0.0000 0.458167
Empenigr 0.004536 0.010325 0.4393 0.0 0.3333 0.581673
Hyporadi 0.017141 0.026548 0.6457 0.0 1.1667 0.442231
Juncarti 0.027414 0.028537 0.9607 1.6 1.1667 0.175299
Juncbufo 0.018181 0.024648 0.7376 1.2 0.0000 0.386454
Lolipere 0.054533 0.029625 1.8408 4.0 0.3333 0.151394
Planlanc 0.041633 0.029560 1.4084 3.0 0.8333 0.007968
Poaprat 0.041750 0.018852 2.2146 3.4 0.6667 0.015936
Poatriv 0.071553 0.013681 5.2302 4.8 0.0000 0.019920
Ranuflam 0.019285 0.019939 0.9672 0.4 1.3333 0.334661
Rumeacet 0.046546 0.030806 1.5109 3.2 0.0000 0.003984
Sagiproc 0.015282 0.016535 0.9243 0.8 0.5000 0.892430
Salirepe 0.025340 0.027291 0.9285 0.0 1.8333 0.135458
Scorautu 0.020703 0.014125 1.4658 2.8 3.1667 0.780876
Trifprat 0.025843 0.025972 0.9951 1.8 0.0000 0.003984
Trifrepe 0.030372 0.022871 1.3280 2.8 1.8333 0.617530
Vicilath 0.002399 0.005461 0.4393 0.0 0.1667 0.836653
Bracruta 0.035340 0.020104 1.7579 2.8 2.8333 0.290837
Callcusp 0.016833 0.024901 0.6760 0.0 1.1667 0.438247
Whereas, in the results below, using a newer version of R (4.3.2) Chenalbu and Cirsarve have a p value of 0.4980 and 0.5976, respectively. Because there are 0 observations for both of these taxon, it is unreasonable to have a p value less than 0.9999.
Contrast: HF_NM
Taxon average sd ratio ava avb p
Achimill 0.01656 0.01490 1.11100 1.20000 0.333 0.6534
Agrostol 0.03192 0.02889 1.10500 1.40000 2.167 0.9442
Airaprae 0.01261 0.01824 0.69100 0.00000 0.833 0.2829
Alopgeni 0.02446 0.03240 0.75500 1.60000 0.000 0.9163
Anthodor 0.02872 0.02480 1.15800 1.80000 1.333 0.2829
Bellpere 0.00880 0.01373 0.64100 0.40000 0.333 0.9203
Bromhord 0.01209 0.01517 0.79700 0.80000 0.000 0.8606
Chenalbu 0.00000 0.00000 NaN 0.00000 0.000 0.4980
Cirsarve 0.00000 0.00000 NaN 0.00000 0.000 0.5976
Comapalu 0.01011 0.01456 0.69400 0.00000 0.667 0.2749
Eleopalu 0.03204 0.03231 0.99200 0.80000 2.167 0.5538
Elymrepe 0.02981 0.03868 0.77100 2.00000 0.000 0.4900
Empenigr 0.00454 0.01033 0.43900 0.00000 0.333 0.5578
Hyporadi 0.01714 0.02655 0.64600 0.00000 1.167 0.4263
Juncarti 0.02741 0.02854 0.96100 1.60000 1.167 0.1912
Juncbufo 0.01818 0.02465 0.73800 1.20000 0.000 0.3307
Lolipere 0.05453 0.02962 1.84100 4.00000 0.333 0.1633
Planlanc 0.04163 0.02956 1.40800 3.00000 0.833 0.0040
Poaprat 0.04175 0.01885 2.21500 3.40000 0.667 0.0199
Poatriv 0.07155 0.01368 5.23000 4.80000 0.000 0.0040
Ranuflam 0.01928 0.01994 0.96700 0.40000 1.333 0.3586
Rumeacet 0.04655 0.03081 1.51100 3.20000 0.000 0.0040
Sagiproc 0.01528 0.01653 0.92400 0.80000 0.500 0.8566
Salirepe 0.02534 0.02729 0.92900 0.00000 1.833 0.1155
Scorautu 0.02070 0.01412 1.46600 2.80000 3.167 0.8526
Trifprat 0.02584 0.02597 0.99500 1.80000 0.000 0.0040
Trifrepe 0.03037 0.02287 1.32800 2.80000 1.833 0.6733
Vicilath 0.00240 0.00546 0.43900 0.00000 0.167 0.8685
Bracruta 0.03534 0.02010 1.75800 2.80000 2.833 0.3386
Callcusp 0.01683 0.02490 0.67600 0.00000 1.167 0.4980
Does anyone know where the error may be arising and how to resolve it? Because the taxon with 0 observations are incorrect, can I have confidence in the other resulting p values?
The SIMPER code produces the correct results when using an older ASUS running windows 10 (with R 3.6.0, Rstudio 1.3.1093, vegan 2.5-7, lattice 0.20-38, permute 0.9-5).
The SIMPER analysis produces unreasonable values when a use a Dell computer running Windows 10 pro (with R 4.3.2, vegan 2.6-4, lattice 0.21-9, permute 0.9-7). I have tried downloading older versions of R and Rstudio, as well as lattice, and it does not resolve the issue.
Below is the code I have used. Note: the suggested SIMPER code does not produce p values (only cumulative sum), thus, I have included an altered code which produces p values.
options(max.print=1000000)
options(scipen = 999)
library(vegan)
data(dune)
data(dune.env)
##As suggested Online -- Rdir.io https://rdrr.io/rforge/vegan/man/simper.html
(sim <- with(dune.env, simper(dune, Management)))
summary(sim)
##altered to output p value
dune2 <- with(dune.env, simper(dune, Management, permutations = 250))
summary(dune2, ordered = F)
I have tried un-installing and re-installing R and RStudio, using various versions of R and Rstudio, running different versions of Vegan, lattice, and permute, and using different computers. Only the Asus using an older version of R and RStudio produces reasonable results. I have tried using the same version of R and R studio on my Dell computer and cannot yield the same results.