Scaling of Canonical Coefficients in prc() using pyrifos data set

[Hansam]

Hello,

I am looking for some help with an issue I encountered while reproducing the PRC from Szoecs 2015 (as by the tutorial in the Supplementary Material: https://www.researchgate.net/publication/272094723_Supplemental_Material_for_the_paper_Analysing_chemical-induced_changes_in_macroinvertebrate_communities_in_aquatic_mesocosm_experiments_A_comparison_of_methods).

I am using the pyrifos data set that is supplied by the vegan package and also in Canoco 5 to validate my code for the use in GLP studies (similar to this closed issue [https://github.com/vegandevs/vegan/issues/17].

I can reproduce the results with the exception of the canonical coefficients which are approximately 10x lower. As a result, when plotting the data the shape of the curves is the same, but not the values on the y-axis.

All the distances and results are proportional and the statistical tests on the PRC also match. However, as this exercise is to validate the code, I will need to convince the GLP auditors that the code and its result are reliable. As I cannot explain the 10 fold difference in the coefficients this will be difficult.

Any help/explanation is very much appreciated!

Thank you, Hanna

require(vegan)
data(pyrifos)
head(pyrifos[, c(1:11)])

week <- gl(11, 12, labels = c(-4, -1, 0.1, 1, 2, 4, 8, 12, 15, 19, 24))
dose <- factor(rep(c(0.1, 0, 0, 0.9, 0, 44, 6, 0.1, 44, 0.9, 0, 6), 11))
ditch <- gl(12, 1, length=132)

pyrifos_prc <- prc(response = pyrifos, treatment = dose, time = week)
pyrifos_prc_sum <- summary(pyrifos_prc, scaling = 1)

plot(pyrifos_prc, select = abs(pyrifos_prc_sum$sp) > 1, scaling=1) 

pyrifos_prc_sum

Call:
prc(response = pyrifos, treatment = dose, time = week) 
Species scores:
    Simve     Daplo     Cerpu     Alogu     Aloco     Alore     Aloaf     Copsp     Ostsp     Slyla     Acrha     Aloex     Chysp 
-2.688099 -1.464566 -0.542739 -0.280040 -0.177019 -0.315038 -0.426524 -1.169368 -2.312186  0.556899 -0.105535 -0.228092 -0.095042 
    Alona     Plead     Oxyte     Grate     Copdi     NauLa     CilHa     Strvi     amosp     Ascmo     Synsp     Squro     Squmu 
-0.063689 -0.138397 -0.025401 -0.096840 -1.428854 -4.847070 -0.895241 -3.069709  1.357663 -0.069736  0.026494 -0.264390  0.452667 
    Polar     Kerqu     Anufi     Mytve     Mytvi     Mytmu     Lepsp     Leppa     Colob     Colbi     Colun     Lecsp     Lecqu 
-0.461989 -0.495348 -0.432767 -0.074372 -0.090928 -0.105891 -0.998286 -0.084809  0.723051  0.139569  0.828338 -0.472866  0.088860 
    Lecco     Leclu     Lecfl     Tripo     Cepsp     Monlo     Monae     Scalo     Trilo   Tripo.1     Tricy     Trisp     Tepat 
 0.290531 -0.049788  0.408035 -0.215234  0.809597  0.527913  0.089948  0.077192  0.039086 -0.246435 -0.335400 -0.078423 -0.007719 
    Rotne     Notla     Filsp     Lopox  hydrspec  bothrosp  olchaeta  erpoocto  glsicomp  alglhete  hebdstag   sphidae  ansuvote 
-0.143730  0.114301  0.168356  0.030990  0.048698 -0.398665  1.165154  0.901030  0.144389  0.073049 -0.902223 -1.463655 -0.140685 
 armicris  bathcont  binitent  gyraalbu  hippcomp  lymnstag  lymnaes7  physfont  plbacorn  popyanti  radiovat  radipere  valvcris 
-1.680010 -0.073282  1.950500 -0.033051 -0.404473  0.263679 -0.135150  0.026383 -0.084761 -1.272223  0.019815  0.625468 -0.010579 
 valvpisc  hycarina  gammpule  aselaqua  proameri  collembo  caenhora  caenluct  caenrobu  cloedipt  cloesimi  aeshniae  libellae 
 0.267577 -1.044034 -1.526450 -1.578743 -0.116577 -0.029906 -5.767844 -2.376188 -0.126181 -4.734035 -1.242207 -0.212699  0.081867 
 conagrae  corident  coripanz  coripunc  cymabons  hesplinn  hespsahl  notoglau  notomacu  notoobli  notoviri  pacoconc  pleaminu 
-1.630574 -0.013761 -0.120439  0.176746 -0.046360  0.069465 -0.033240 -0.555201 -0.050060 -0.082356 -0.215863 -0.016709 -0.071168 
 sigadist  sigafall  sigastri  sigarasp  colyfusc  donacis6  gyrimari  haliconf  haliflav  haligruf  haliobli  herubrev  hya_herm 
-0.076479  0.018364 -0.060206  0.277312 -0.036493  0.078608 -0.010579 -0.446911 -0.044538 -0.450146 -0.131472  0.128486  0.322379 
 hyglpusi  hyhyovat  hypoplan  hyporusp  hytuinae  hytuvers  laphminu  noteclav  rhantusp  sialluta  ablalong  ablaphmo  cltanerv 
 0.011775 -0.024196 -0.014976 -0.232042 -2.316179 -1.772351 -0.632897 -0.007912 -0.067618 -1.109341 -0.014976 -2.992695 -0.075631 
 malopisp  mopetenu  prdiussp  pstavari  chironsp  crchirsp  crclglat  ditendsp  mitegchl  pachgarc  pachgvit  popegnub  popedisp 
-0.047658 -0.008716 -0.554911 -0.082842 -1.889916 -0.016709 -0.028953 -0.083481 -0.230629  0.012187  0.029907 -0.224272  0.069649 
 acriluce  chclpige  conescut  cricotsp  liesspec  psclbarb  psclgsli  psclobvi  psclplat  psclpsil  pscladsp  cladotsp  laa_spec 
 0.007950  0.008744 -0.821036 -0.121530 -0.107387  0.028639 -0.601568  0.362378  0.052054 -0.007339 -0.005674 -0.539894 -0.034105 
 patanysp  tatarssp  zaa_spec  anopmacu  cepogoae  chaoobsc  cucidae4  tabanusp  agdasphr  athrater  cyrncren  holodubi  holopici 
-0.146807 -0.669430 -0.049949  0.163731 -2.555403 -2.442310 -0.033240  0.011601 -0.271815 -0.067618 -0.071168 -0.094754 -0.611618 
 leceriae  lilurhom  monaangu  mystazur  mystloni  oecefurv  oecelacu  triabico  paponysp 
-0.298633 -0.009063 -0.644295 -0.033240 -2.998460 -0.536628 -0.259064 -0.088915 -0.097788 

Coefficients for dose + week:dose interaction
which are contrasts to dose 0 
rows are dose, columns are week
          -4      -1     0.1        1       2       4        8      12      15      19       24
0.1 0.007262 0.01383 0.01026 0.004093 0.02114 0.01373 0.008126 0.01545 0.01130 0.02176 0.007883
0.9 0.008155 0.01947 0.01948 0.047992 0.05008 0.04333 0.013615 0.03569 0.02423 0.01767 0.015537
6   0.016718 0.01240 0.04567 0.116348 0.10902 0.11582 0.056456 0.04727 0.03097 0.03313 0.018339
44  0.014064 0.01970 0.07353 0.126863 0.13057 0.14716 0.129933 0.10143 0.07867 0.05804 0.031213

pyrifos_prc

Call: prc(response = pyrifos, treatment = dose, time = week)

               Inertia Proportion Rank
Total         288.9920     1.0000     
Conditional    63.3493     0.2192   10
Constrained    96.6837     0.3346   44
Unconstrained 128.9589     0.4462   77
Inertia is variance 

Eigenvalues for constrained axes:
  RDA1   RDA2   RDA3   RDA4   RDA5   RDA6   RDA7   RDA8   RDA9  RDA10  RDA11  RDA12  RDA13  RDA14  RDA15  RDA16  RDA17  RDA18 
25.282  8.297  6.044  4.766  4.148  3.857  3.587  3.334  3.087  2.551  2.466  2.209  2.129  1.941  1.799  1.622  1.579  1.440 
 RDA19  RDA20  RDA21  RDA22  RDA23  RDA24  RDA25  RDA26  RDA27  RDA28  RDA29  RDA30  RDA31  RDA32  RDA33  RDA34  RDA35  RDA36 
 1.398  1.284  1.211  1.133  1.001  0.923  0.862  0.788  0.750  0.712  0.685  0.611  0.584  0.537  0.516  0.442  0.417  0.404 
 RDA37  RDA38  RDA39  RDA40  RDA41  RDA42  RDA43  RDA44 
 0.368  0.340  0.339  0.306  0.279  0.271  0.205  0.179 

Eigenvalues for unconstrained axes:
   PC1    PC2    PC3    PC4    PC5    PC6    PC7    PC8 
17.156  9.189  7.585  6.064  5.730  4.843  4.518  4.105 
(Showing 8 of 77 unconstrained eigenvalues)

[jarioksa]

Which version of vegan you are using? The version in CRAN and in github have different kind of scores.

[Hansam]

Hello :),

I'm using the CRAN version. Could you explain why the two versions have different scores? And should the scores from the github package match the Szoecs paper?

Thanks a lot, Hanna

[jarioksa]

Returning back after long delays: I had now time to have a look at the prc() function and compare it against other alternatives.

I asked about version, because CRAN versions since vegan 2.5-1 use wrong kind of scores. The github version was fixed in November 2022 and since then it is similar as before 2.5-1 release of vegan. So older prc (before 2.5-1) and github version give similar results, but in between the result diverged a bit.

About scaling: the numeric scaling of vegan results is not fixed, but it depends on the options you choose. Arguments scaling and const will change numeric values. However, the results will differ only by scaling, that is, a constant multiplier. The relative differences will be unchanged.