I initially used a GLM (because LMA data are neither normal nor heteroskedastic) and on a reduced dataset because of a lot of lines with NAs for LMA. In order to have a larger dataset, we supplemented LMA measurements with average values for some species (variable "Species") having n>5 LMA measurement. The variance of LMA is now highly heterogeneous and reduced for some of the species.

I initially used a GLM (because LMA data are neither normal nor homoskedastic) and on a reduced dataset because of a lot of lines with NAs for LMA. In order to have a larger dataset, we supplemented LMA measurements with average values for some species (variable "Species") having n>5 LMA measurement. The variance of LMA is now highly heterogeneous and reduced for some of the species.

but since we took species-mean LMA, it is quite logical ... if we consider the relation to species is biased because we used species mean values, then the best model would be xm2b: LMA ~ Genus + PI + Locality + specimen_TCT + Pheno

but since we took species-mean LMA, it is quite logical ... if we consider the relation to species is biased because we used species mean values, then the best model would be xm2b: LMA ~ Genus + PI + Locality + specimen_TCT + Pheno

Thus, I switched to GLS to take into account this particularity with "weights = varIdent(form=~1|Species)" and I made several models with method = REML to test the effect of the different variables (as some categorical variables were nested, I did not test them at the same time).

Bests,

Agathe

x <- c("species_LMA","Locality","Family","Genus","Species", "specimen_TCT","taxon_TCT","combined_TCT","cat_TCT","PI","Pheno","Growth_form")
noNA4 <- subset(data2, complete.cases(data2[,x])) # n = 561

Here is are the different models I made:

xm1a <- gls(species_LMA ~ -1 + Family + PI + Locality + cat_TCT + Pheno , weights = varIdent(form=~1|Species),data=noNA4, na.action=na.omit)
xm1b <- gls(species_LMA ~ -1 + Family + PI + Locality + specimen_TCT + Pheno , weights = varIdent(form=~1|Species),data=noNA4, na.action=na.omit)
xm1c <- gls(species_LMA ~ -1 + Family + PI + Locality + combined_TCT + Pheno , weights = varIdent(form=~1|Species),data=noNA4, na.action=na.omit)

xm2a <- gls(species_LMA ~ -1 + Genus + PI + Locality + cat_TCT + Pheno , weights = varIdent(form=~1|Species),data=noNA4, na.action=na.omit)
xm2b <- gls(species_LMA ~ -1 + Genus + PI + Locality + specimen_TCT + Pheno , weights = varIdent(form=~1|Species),data=noNA4, na.action=na.omit)
xm2c <- gls(species_LMA ~ -1 + Genus + PI + Locality + combined_TCT + Pheno , weights = varIdent(form=~1|Species),data=noNA4, na.action=na.omit)

xm3a <- gls(species_LMA ~ -1 + Species + PI + Locality + cat_TCT + Pheno , weights = varIdent(form=~1|Species),data=noNA4, na.action=na.omit)
#xm3b <- gls(species_LMA ~ -1 + Species + PI + Locality + specimen_TCT + Pheno , weights = varIdent(form=~1|Species),data=noNA3, na.action=na.omit) # do not work for some reason... correlation?
xm3c <- gls(species_LMA ~ -1 + Species + PI + Locality + combined_TCT + Pheno , weights = varIdent(form=~1|Species),data=noNA4, na.action=na.omit)

model.sel(xm1a,xm1b,xm1c,xm2a,xm2b,xm2c,xm3a,xm3c)

Model selection table 
     cat_TCT Fml Lcl Phn     PI spc_TCT cmb_TCT Gns Spc df    logLik   AICc  delta weight
xm3c               +   +  8.254               +       + 60 -2051.166 4237.2   0.00      1
xm3a       +       +   +  8.604                       + 56 -2064.354 4253.6  16.37      0
xm2b               +   + 10.380       +           +     52 -2085.641 4286.3  49.10      0
xm2c               +   +  9.439               +   +     49 -2092.155 4292.1  54.84      0
xm2a       +       +   + 10.510                   +     45 -2105.152 4308.5  71.26      0
xm1b           +   +   +  2.933       +                 47 -2111.535 4326.0  88.79      0
xm1c           +   +   +  5.425               +         44 -2115.838 4327.5  90.26      0
xm1a       +   +   +   + 11.680                         40 -2138.920 4364.2 127.04      0
Models ranked by AICc(x) 

# best option are *3c* then 3a, 2b, 2c, 2a, 1c, 1b, 1a but AIC are highly similar
# xm3c : LMA ~ Species + PI + Locality + combined_TCT + Pheno
# but since we took species-mean LMA, it is quite logical... if we consider the relation to species is biased because we used species mean values, then the best model would be xm2b: LMA ~ Genus + PI + Locality + specimen_TCT + Pheno

anova(update(xm3c, . ~ ., method = "ML"), update(xm2b, . ~ ., method = "ML"))

Error in eigen(val, only.values = TRUE) : 
  infinite or missing values in 'x'

> anova(xm1b) 
Denom. DF: 531 
             numDF   F-value p-value
Family           9 2137401.7  <.0001
PI               1       0.1  0.7331
Locality         1     118.4  <.0001
specimen_TCT    10       5.2  <.0001
Pheno            1       0.1  0.7080

> anova(xm2b) # nothing but Species is significant... marginal corr of PI... 
Denom. DF: 526 
             numDF   F-value p-value
Genus           14 247220.46  <.0001
PI               1      2.76  0.0970
Locality         1      0.10  0.7500
specimen_TCT    10      0.56  0.8456
Pheno            1      0.41  0.5226

> anova(xm3c) # nothing but Species is significant
Denom. DF: 518 
             numDF   F-value p-value
Species         25 195602.70  <.0001
PI               1      2.12  0.1463
Locality         1      0.18  0.6691
combined_TCT     7      0.75  0.6304
Pheno            1      0.17  0.6784

``

Thus, I switched to GLS to take into account this particularity with weights = varIdent(form=\~1|Species) and I made several models with method = REML to test the effect of the different variables (as some categorical variables were nested, I did not test them at the same time).

x <- c("species_LMA", "Locality", "Family", 
    "Genus", "Species", "specimen_TCT", 
    "taxon_TCT", "combined_TCT", "cat_TCT", 
    "PI", "Pheno", "Growth_form")  
noNA4 <- subset(data2, complete.cases(data2[, 
                x])) # n = 561

Here are are the different models I made:

xm1a <- gls(species_LMA ~ -1 + Family + PI + 
    Locality + cat_TCT + Pheno , 
    weights = varIdent(form=~1|Species), 
    data=noNA4, na.action=na.omit)  
xm1b <- gls(species_LMA ~ -1 + Family + PI + 
    Locality + specimen_TCT + Pheno , 
    weights = varIdent(form=~1|Species), 
    data=noNA4, na.action=na.omit)  
xm1c <- gls(species_LMA ~ -1 + Family + PI + 
    Locality + combined_TCT + Pheno , 
    weights = varIdent(form=~1|Species), 
    data=noNA4, na.action=na.omit)

xm2a <- gls(species_LMA ~ -1 + Genus + PI + 
    Locality + cat_TCT + Pheno , 
    weights = varIdent(form=~1|Species), 
    data=noNA4, na.action=na.omit)  
xm2b <- gls(species_LMA ~ -1 + Genus + PI + 
    Locality + specimen_TCT + Pheno , 
    weights = varIdent(form=~1|Species), 
    data=noNA4, na.action=na.omit)     
xm2c <- gls(species_LMA ~ -1 + Genus + PI +  
    Locality + combined_TCT + Pheno , 
    weights = varIdent(form=~1|Species), 
    data=noNA4, na.action=na.omit)

xm3a <- gls(species_LMA ~ -1 + Species + PI + 
    Locality + cat_TCT + Pheno , 
    weights = varIdent(form=~1|Species), 
    data=noNA4, na.action=na.omit)    
#xm3b <- gls(species_LMA ~ -1 + Species + PI + 
    Locality + specimen_TCT + Pheno , 
    weights = varIdent(form=~1|Species), 
    data=noNA3, na.action=na.omit) 
# do not work for some reason ... correlation?`  
xm3c <- gls(species_LMA ~ -1 + Species + PI + 
    Locality + combined_TCT + Pheno , 
    weights = varIdent(form=~1|Species), 
    data=noNA4, na.action=na.omit)

model.sel(xm1a, xm1b, xm1c, xm2a, xm2b, xm2c, 
          xm3a, xm3c) 

    Model selection table 
         cat_TCT Fml Lcl Phn     PI spc_TCT cmb_TCT Gns Spc df    logLik   AICc  delta weight
    xm3c               +   +  8.254               +       + 60 -2051.166 4237.2   0.00      1
    xm3a       +       +   +  8.604                       + 56 -2064.354 4253.6  16.37      0
    xm2b               +   + 10.380       +           +     52 -2085.641 4286.3  49.10      0
    xm2c               +   +  9.439               +   +     49 -2092.155 4292.1  54.84      0
    xm2a       +       +   + 10.510                   +     45 -2105.152 4308.5  71.26      0
    xm1b           +   +   +  2.933       +                 47 -2111.535 4326.0  88.79      0
    xm1c           +   +   +  5.425               +         44 -2115.838 4327.5  90.26      0
    xm1a       +   +   +   + 11.680                         40 -2138.920 4364.2 127.04      0
    Models ranked by AICc(x) 

# best option are *3c* then 3a, 2b, 2c, 2a, 1c, 1b, 1a but AIC are highly similar  
# xm3c : LMA ~ Species + PI + Locality + combined_TCT + Pheno 

but since we took species-mean LMA, it is quite logical ... if we consider the relation to species is biased because we used species mean values, then the best model would be xm2b: LMA ~ Genus + PI + Locality + specimen_TCT + Pheno

anova(update(xm3c, . ~ ., method = "ML"), update(xm2b, . ~ ., method = "ML"))`

    Error in eigen(val, only.values = TRUE) : 
      infinite or missing values in 'x'

    anova(xm1b) 
    Denom. DF: 531 
                 numDF   F-value p-value
    Family           9 2137401.7  <.0001
    PI               1       0.1  0.7331
    Locality         1     118.4  <.0001
    specimen_TCT    10       5.2  <.0001
    Pheno            1       0.1  0.7080
    
    > anova(xm2b) # nothing but Species is significant... marginal corr of PI... 
    Denom. DF: 526 
                 numDF   F-value p-value
    Genus           14 247220.46  <.0001
    PI               1      2.76  0.0970
    Locality         1      0.10  0.7500
    specimen_TCT    10      0.56  0.8456
    Pheno            1      0.41  0.5226
    
    anova(xm3c) # nothing but Species is 
                # significant
    Denom. DF: 518 
                 numDF   F-value p-value
    Species         25 195602.70  <.0001
    PI               1      2.12  0.1463
    Locality         1      0.18  0.6691
    combined_TCT     7      0.75  0.6304
    Pheno            1      0.17  0.6784