0
$\begingroup$

I am looking for a way to analyse the differences in vegetation carbon stocks (named after 'OC_TOT') among different vegetation types (named after 'Facies') and islands (Guadeloupe and Martinique).

The dataset was assembled from several published and unpublished data obtained either from Master Internship, PhD Thesis or different kind of surveys. The sampling design have changed a lot during the years depending on subsequent scientific questions and operators, thus leading to unequal sampling plot areas within the facies and unequal sampling effort (i.e. unbalanced design). The boxplots and a summary of the values is presented below :

enter image description here

FSA::Summarize(OC_TOT~Facies, data = MgHa)
    Facies  n      mean       sd    min       Q1  median       Q3    max
1      DSF  4  97.24250 31.67354  68.58  80.5725  89.130 105.8000 142.13
2   Fringe 11 101.69818 56.02387  25.22  55.5050 112.830 133.0300 218.71
3 High Avi 10  87.50100 30.11720  44.81  67.0375  92.270 100.8925 150.78
4  Low Avi 16  62.10500 27.29377  22.38  45.6025  56.920  77.1950 121.70
5    Mixed 18 123.91833 37.68447  54.17  99.8700 113.695 154.5625 183.96
6      OSF 13 248.32385 45.63600 161.26 221.2600 247.460 289.6400 307.00
7    Scrub 15  55.92733 23.95107  21.43  38.2350  49.790  71.4700  99.96
8      SSF  8 150.14875 12.51794 130.87 142.8100 147.395 161.8250 167.09

Distribution of my values within each facies and their corresponding shapiro.test responses.

I tried several ways that are presented below :

  1. PERMANOVA : PERMANOVA was considered as a way to bypass the heteroskedasticity and the existence of some extraordinary Ctot values considered as outliers for R but that represent true patches of forests. I considered first the function aovp from the lmPerm package (1) that gave me these results :
> ModFac <- aovp(OC_TOT ~ Facies,data=MgHa,perm="Prob")
[2] "Settings:  unique SS "
> summary(ModFac)
Component 1 :
            Df R Sum Sq R Mean Sq Iter  Pr(Prob)    
Facies       7   350247     50035 5000 < 2.2e-16 ***
Residuals   87   111996      1287                   

> ModFac <- glht(ModFac,linfct=mcp(Facies="Tukey"))
> cld(ModFac)
     DSF   Fringe High Avi  Low Avi    Mixed      OSF    Scrub      SSF 
    "ad"    "bcd"     "ac"     "ab"     "cd"      "e"      "a"      "d" 
Warning messages:
1: In RET$pfunction("adjusted", ...) : Completion with error > abseps
2: In RET$pfunction("adjusted", ...) : Completion with error > abseps
3: In RET$pfunction("adjusted", ...) : Completion with error > abseps
4: In RET$pfunction("adjusted", ...) : Completion with error > abseps

If I understand right, this method takes into account ~75% of the variability of the OC_TOT values which is quite good. However, I don't understand the warning message and how DSF can be "ad" and not being "b" in the same time... Quite disturbing...

  1. Basic linear regression : I have created the following model test <- lm(sqrt(OC_TOT) ~ Facies*Ile, data = MgHa) and its analysis gives me the following results :
> shapiro.test(test$residuals)

    Shapiro-Wilk normality test

data:  test$residuals
W = 0.99051, p-value = 0.7372

> leveneTest(sqrt(OC_TOT)~Facies, data = MgHa)
Levene's Test for Homogeneity of Variance (center = median)
      Df F value  Pr(>F)  
group  7  2.1583 0.04573 *
      87                  

> leveneTest(sqrt(OC_TOT)~Ile, data = MgHa)
Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  1  1.8466 0.1775
      93               

> leveneTest(sqrt(OC_TOT)~Facies*Ile, data = MgHa)
Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group 13  0.9603 0.4969
      81     

enter image description here

> Comp <- aov(test))
> summary(Comp)
            Df Sum Sq Mean Sq F value Pr(>F)    
Facies       7  674.4   96.34  31.634 <2e-16 ***
Ile          1    4.0    4.00   1.315  0.255    
Facies:Ile   5   19.4    3.88   1.275  0.283    
Residuals   81  246.7    3.05                   

TukeyHSD(Comp, which = "Facies")
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = test)

$Facies
                        diff         lwr        upr     p adj
Fringe-DSF       -0.06113884  -3.2316523  3.1093746 1.0000000
High Avi-DSF     -0.53866288  -3.7511717  2.6738459 0.9995136
Low Avi-DSF      -2.06552155  -5.1010570  0.9700139 0.4129741
Mixed-DSF         1.23624359  -1.7653742  4.2378614 0.9030854
OSF-DSF           5.92365832   2.8188622  9.0284545 0.0000019
Scrub-DSF        -2.45141982  -5.5071252  0.6042856 0.2119837
SSF-DSF           2.47449595  -0.8507666  5.7997585 0.2981994
High Avi-Fringe  -0.47752404  -2.8501190  1.8950710 0.9984077
Low Avi-Fringe   -2.00438271  -4.1312278  0.1224624 0.0793799
Mixed-Fringe      1.29738243  -0.7807667  3.3755316 0.5264751
OSF-Fringe        5.98479717   3.7602178  8.2093765 0.0000000
Scrub-Fringe     -2.39028098  -4.5458156 -0.2347463 0.0192162
SSF-Fringe        2.53563479   0.0124691  5.0588005 0.0479851
Low Avi-High Avi -1.52685867  -3.7158144  0.6620971 0.3803973
Mixed-High Avi    1.77490647  -0.3667661  3.9165790 0.1790112
OSF-High Avi      6.46232121   4.1782873  8.7463551 0.0000000
Scrub-High Avi   -1.91275694  -4.1295986  0.3040847 0.1423476
SSF-High Avi      3.01315883   0.4374215  5.5888961 0.0107882
Mixed-Low Avi     3.30176514   1.4360148  5.1675154 0.0000113
OSF-Low Avi       7.98917987   5.9616016 10.0167582 0.0000000
Scrub-Low Avi    -0.38589827  -2.3374746  1.5656780 0.9985784
SSF-Low Avi       4.54001750   2.1887018  6.8913332 0.0000014
OSF-Mixed         4.68741474   2.7109766  6.6638529 0.0000000
Scrub-Mixed      -3.68766341  -5.5860532 -1.7892736 0.0000012
SSF-Mixed         1.23825236  -1.0691096  3.5456143 0.7062913
Scrub-OSF        -8.37507815 -10.4327306 -6.3174257 0.0000000
SSF-OSF          -3.44916238  -5.8892380 -1.0090868 0.0008311
SSF-Scrub         4.92591577   2.5486179  7.3032136 0.0000002

> rcompanion::cldList(p.adj ~ Comp, data = PT, threshold = 0.05)

    Group Letter MonoLetter
1  Fringe     ab      ab   
2 HighAvi    abc      abc  
3  LowAvi     ac      a c  
4   Mixed     bd       b d 
5     OSF      e          e
6   Scrub      c        c  
7     SSF      d         d 
8     DSF   abcd      abcd

Here again I can't figure out why the DSF facies (97.24 +/- 31.67) isn't different from the SSF facies while the Fringe facies (101.70 +/- 56.02) is...

  1. Linear Mixed Model with nested random factors : I made the following model with Site nested in Facies and plot nested in Site
> test2 <- lmer(sqrt(OC_TOT) ~ Facies*Ile + (1|Facies:Site) + (1|Site:Parcelle), data = MgHa)
fixed-effect model matrix is rank deficient so dropping 2 columns / coefficients
boundary (singular) fit: see ?isSingular 
> summary(test2)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: sqrt(OC_TOT) ~ Facies * Ile + (1 | Facies:Site) + (1 | Site:Parcelle)
   Data: MgHa

REML criterion at convergence: 337.3

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-0.70503 -0.00812 -0.00052  0.00790  0.70913 

Random effects:
 Groups        Name        Variance  Std.Dev. 
 Site:Parcelle (Intercept) 3.080e+00 1.755e+00
 Facies:Site   (Intercept) 2.918e-12 1.708e-06
 Residual                  5.188e-04 2.278e-02
Number of obs: 95, groups:  Site:Parcelle, 94; Facies:Site, 16

Fixed effects:
                             Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)                    9.7697     0.8775 80.0005  11.133  < 2e-16 ***
FaciesFringe                  -1.2636     1.1329 80.0003  -1.115   0.2680    
FaciesHigh Avi                -0.2170     1.1329 80.0003  -0.192   0.8486    
FaciesLow Avi                 -2.1890     1.1329 80.0003  -1.932   0.0569 .  
FaciesMixed                    1.1365     1.0133 80.0004   1.122   0.2654    
FaciesOSF                      5.8066     1.0133 80.0004   5.731  1.7e-07 ***
FaciesScrub                   -2.4036     1.0133 80.0002  -2.372   0.0201 *  
FaciesSSF                      2.4052     1.1000 79.9996   2.186   0.0317 *  
IleMartinique                 -0.2390     1.1329 79.9997  -0.211   0.8335    
FaciesFringe:IleMartinique     2.8845     1.5533 79.9999   1.857   0.0670 .  
FaciesHigh Avi:IleMartinique  -0.5651     1.6021 79.9999  -0.353   0.7252    
FaciesLow Avi:IleMartinique    0.4366     1.4508 79.9998   0.301   0.7642    
FaciesMixed:IleMartinique      0.5382     1.4330 79.9998   0.376   0.7082    
FaciesOSF:IleMartinique        1.7603     2.1495 79.9999   0.819   0.4152    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation matrix not shown by default, as p = 14 > 12.
Use print(x, correlation=TRUE)  or
    vcov(x)        if you need it

fit warnings:
fixed-effect model matrix is rank deficient so dropping 2 columns / coefficients
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see ?isSingular

> summary(glht(test2, linfct = mcp(Facies = "Tukey")), test = adjusted("holm"))

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: lmer(formula = sqrt(OC_TOT) ~ Facies * Ile + (1 | Facies:Site) + 
    (1 | Site:Parcelle), data = MgHa)

Linear Hypotheses:
                        Estimate Std. Error z value Pr(>|z|)    
Fringe - DSF == 0        -1.2636     1.1329  -1.115 1.000000    
High Avi - DSF == 0      -0.2170     1.1329  -0.192 1.000000    
Low Avi - DSF == 0       -2.1890     1.1329  -1.932 0.567957    
Mixed - DSF == 0          1.1365     1.0133   1.122 1.000000    
OSF - DSF == 0            5.8066     1.0133   5.731 2.20e-07 ***
Scrub - DSF == 0         -2.4036     1.0133  -2.372 0.229924    
SSF - DSF == 0            2.4052     1.1000   2.186 0.345349    
High Avi - Fringe == 0    1.0466     1.0133   1.033 1.000000    
Low Avi - Fringe == 0    -0.9254     1.0133  -0.913 1.000000    
Mixed - Fringe == 0       2.4002     0.8775   2.735 0.099769 .  
OSF - Fringe == 0         7.0703     0.8775   8.057 2.31e-14 ***
Scrub - Fringe == 0      -1.1400     0.8775  -1.299 1.000000    
SSF - Fringe == 0         3.6689     0.9764   3.757 0.002918 ** 
Low Avi - High Avi == 0  -1.9720     1.0133  -1.946 0.567957    
Mixed - High Avi == 0     1.3535     0.8775   1.542 1.000000    
OSF - High Avi == 0       6.0236     0.8775   6.864 1.67e-10 ***
Scrub - High Avi == 0    -2.1866     0.8775  -2.492 0.177943    
SSF - High Avi == 0       2.6222     0.9764   2.686 0.108610    
Mixed - Low Avi == 0      3.3256     0.8775   3.790 0.002715 ** 
OSF - Low Avi == 0        7.9957     0.8775   9.112  < 2e-16 ***
Scrub - Low Avi == 0     -0.2146     0.8775  -0.245 1.000000    
SSF - Low Avi == 0        4.5943     0.9764   4.705 5.07e-05 ***
OSF - Mixed == 0          4.6701     0.7165   6.518 1.71e-09 ***
Scrub - Mixed == 0       -3.5401     0.7165  -4.941 1.63e-05 ***
SSF - Mixed == 0          1.2687     0.8347   1.520 1.000000    
Scrub - OSF == 0         -8.2103     0.7165 -11.459  < 2e-16 ***
SSF - OSF == 0           -3.4014     0.8347  -4.075 0.000874 ***
SSF - Scrub == 0          4.8088     0.8347   5.761 1.92e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- holm method)

Warning message:
In mcp2matrix(model, linfct = linfct) :
  covariate interactions found -- default contrast might be inappropriate
> multComp::cld(summary(glht(test2, linfct = mcp(Facies = "Tukey")), test = adjusted("holm")))
     DSF   Fringe High Avi  Low Avi    Mixed      OSF    Scrub      SSF 
    "ab"     "ac"     "ab"      "a"     "bc"      "d"      "a"      "b" 

Here again, I can't figure out why Fringe isn't "abc" at least and, if SSF isn't "c", what does "c" represents... enter image description here

  1. Non-parametric approach : Kruskall-Wallis followed by a Dunn Test. I guess that's all I can do with this kind of experimental design and dataset.
> FSA::dunnTest(MgHa$OC_TOT, MgHa$Facies, method = "bonferroni", two.sided=TRUE)
Dunn (1964) Kruskal-Wallis multiple comparison
  p-values adjusted with the Bonferroni method.

           Comparison           Z      P.unadj        P.adj
1        DSF - Fringe  0.02965118 9.763452e-01 1.000000e+00
2      DSF - High Avi  0.25138801 8.015141e-01 1.000000e+00
3   Fringe - High Avi  0.30075791 7.635991e-01 1.000000e+00
4       DSF - Low Avi  1.28966481 1.971671e-01 1.000000e+00
5    Fringe - Low Avi  1.79647023 7.241976e-02 1.000000e+00
6  High Avi - Low Avi  1.41950658 1.557514e-01 1.000000e+00
7         DSF - Mixed -0.91506269 3.601587e-01 1.000000e+00
8      Fringe - Mixed -1.36692691 1.716482e-01 1.000000e+00
9    High Avi - Mixed -1.65956962 9.700106e-02 1.000000e+00
10    Low Avi - Mixed -3.57040908 3.564242e-04 9.979876e-03
11          DSF - OSF -2.75847564 5.807163e-03 1.626006e-01
12       Fringe - OSF -3.89220280 9.933815e-05 2.781468e-03
13     High Avi - OSF -4.10330634 4.072874e-05 1.140405e-03
14      Low Avi - OSF -6.15479455 7.517489e-10 2.104897e-08
15        Mixed - OSF -2.94359638 3.244228e-03 9.083838e-02
16        DSF - Scrub  1.53953118 1.236747e-01 1.000000e+00
17     Fringe - Scrub  2.13884024 3.244861e-02 9.085611e-01
18   High Avi - Scrub  1.75780145 7.878130e-02 1.000000e+00
19    Low Avi - Scrub  0.40456036 6.858007e-01 1.000000e+00
20      Mixed - Scrub  3.92491689 8.675969e-05 2.429271e-03
21        OSF - Scrub  6.44854214 1.129312e-10 3.162073e-09
22          DSF - SSF -1.59194318 1.113975e-01 1.000000e+00
23       Fringe - SSF -2.13526941 3.273901e-02 9.166923e-01
24     High Avi - SSF -2.36872574 1.784948e-02 4.997856e-01
25      Low Avi - SSF -3.91629775 8.991918e-05 2.517737e-03
26        Mixed - SSF -1.10384090 2.696622e-01 1.000000e+00
27          OSF - SSF  1.34048125 1.800889e-01 1.000000e+00
28        Scrub - SSF -4.20560792 2.603813e-05 7.290676e-04

> cldList(P.adj ~ Comparison, data = PT, threshold = 0.05)
    Group Letter MonoLetter
1     DSF    abc        abc
2  Fringe     ab        ab 
3 HighAvi     ab        ab 
4  LowAvi      a        a  
5   Mixed     bc         bc
6     OSF      c          c
7   Scrub      a        a  
8     SSF     bc         bc  

I still cannot figure out how DSF can be grouped with OSF and but not Fringe.

What about pairwise.t.test if I have multinormality accros all my facies ?


(1) : I have learned later on that the package had been deprecated in 2018. So I considered the vegan package and the adonis function with an euclidean distance. I also started to read Anderson 2001 and Anderson 2014 and found this quote : « PERMANOVA makes no explicit assumptions regarding either the distributions of original variables in Y or the distributions of dissimilarities in D. For a given test, PERMANOVA assumes only exchangeability of permutable units under a true null hypothesis. ». This is further discussed in Anderson 2001 : Observational tests: validity through exchangeability but I guess I cannot use this type of approach.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.