In this script we use a combination of phylogenetically controlled models to evaluate if VIS and NIR absorptivity in the Christmas beetles’ elytra can be predicted by climate and size. Climate here is determined by the ecological variables: maximum temperature, number of days above 35°C, solar radiation, cloud cover, rain, water vapour and aridity summarised in two principal components.

Setting up data

Libraries

source(here::here("Scripts/MacroEcol_1_Libraries.R")) # we source the libraries from this file

Data Sets

Absorptivity data

Initially we import the absorptivity data obtained in the script “optical properties” and join it with the environmental data obtained in “Ecological variables”

Cons1aggA <- read.csv(here::here("Data/FromCode/ConsolidatedAbsorptivitySpp.csv"))[-1] %>%
  dplyr::arrange(phylogeny_name)
# PC Values
PCValuesF <- read.csv(here::here("Data/FromCode/PCsbySpp.csv"))[-1] %>%
  dplyr::arrange(phylogeny_name)

Cons1aggA <-
    PCValuesF %>% 
  filter(phylogeny_name %in% Cons1aggA$phylogeny_name) %>% 
  arrange(phylogeny_name) %>% 
  bind_cols(., Cons1aggA) %>%
  dplyr::select(1:3,5:10) %>% 
  dplyr::rename("phylogeny_name" = phylogeny_name...1, "PC1" = PC1_1, "PC2" = PC2_1)

Phylogenetic data

In our analyses we use a subset of the 1300 posterior sample trees to represent the phylogenetic information accounting for uncertainty in node ages and topology

# set subset limits
trees_subset_min<-1000
trees_subset_max<-2000

Phylogeny (multiple trees)

trees <- ape::read.tree(here::here("Data/XMAS_mat2b_bst2ef_set23nn2_pinct.nwk"))

MCCtree.raw <- 
  ape::read.nexus(here::here("Data/xmas_mat2b_bst2ef_set23nn2_pinct_med.tre"))

The MCC (Maximum clade credibility) tree used here is the BEAST MCC tree. We did not need to prun the tree. The equivalence between specie sin the tree and data frame was tested in previous steps (tab optical properties).

Merge

# Modify to make it compatible with tree tips
Cons1aggA <- as.data.frame(Cons1aggA) # convert to a data frame
rownames(Cons1aggA) <- Cons1aggA[, 1] # make species the row names 
ConsAggA <- Cons1aggA [,2:length(Cons1aggA)] # eliminate spp name (redundant)

The names between data and tree tips should match.

# Prune extra spp in the tree, not contain in the test sample
species.MCC <- as.data.frame(unique(Cons1aggA$phylogeny_name))

# Convert to "row names" (required for following steps)
row.names(species.MCC) <- species.MCC[, 1] 

# Make sure the names in data set and tree match
temp.MCC <- name.check(MCCtree.raw, species.MCC)

# This step is necessary because the tips are different.
MCCtree <- drop.tip(MCCtree.raw, temp.MCC$tree_not_data)

# Test if the species are the same
identical(
  length(name.check(MCCtree, Cons1aggA$phylogeny_name)$tree_not_data),
  length(Cons1aggA$phylogeny_name)
)

## [1] TRUE

Subsets

Create subsets for each spectral band. In these subsets the response variable is always called “response” so that we can use the same function to run various models.

In this script, the response variable is the reflectance on each spectral band or the NIR/VIS residuals.

ALLDataSetA <- 
  Cons1aggA %>% 
  dplyr::select (-VIS, -NIR, -FRS, -FRSP) %>% 
  dplyr::rename ("Response" = TOT) 

NIRDataSetA <- 
  Cons1aggA %>% 
  dplyr::select (-TOT, -VIS, -FRS, -FRSP) %>% 
  dplyr::rename ("Response" = NIR) 

VISDataSetA <- 
  Cons1aggA %>% 
  dplyr::select (-TOT, -NIR, -FRS, -FRSP) %>% 
  dplyr::rename ("Response" = VIS) 

NmoDataSetA <-
  Cons1aggA %>% 
  dplyr::select (-TOT, -FRS, -FRSP) %>% # keep VIS as predictor
  dplyr::rename ("Response" = NIR) # Raw NIR 

FRSDataSetA <-
  Cons1aggA %>% 
  dplyr::select (-TOT, -NIR, -VIS, -FRSP) %>% 
  dplyr::rename ("Response" = FRS) # PGLS Residuals

PRSDataSetA <-
  Cons1aggA %>% 
  dplyr::select (-TOT, -NIR, -VIS, -FRS) %>% 
  dplyr::rename ("Response" = FRSP) # PGLS Phyres residuals

Setting up models

The PGLS function has to be adapted to the data frame and model on each case.

Model 1

PGLS in the MCC

comp_data <- comparative.data(
  phy = MCCtree, data = Cons1aggA,
  names.col = "phylogeny_name", vcv = TRUE,
  na.omit = FALSE, warn.dropped = TRUE
)

Model 2

PGLS Multiple Trees with 5 predictors + intercept

Source function

note that this function has to be adapted to the data frame and model on each case

source(here::here("Scripts/8_multiple_pgls_function_C.R")) # function C is for absorptivity as response

Define model

MuPGLSMod2A <- Response ~ PC1 + PC2 + size + PC1*size + PC2*size

Model 3

PGLS in multiple trees for NIR contains 6 predictors + intercept. The extra predictor here is the “VIS” absorptivity to account for the correlation of these two variables.

source(here::here("Scripts/9_multiple_pgls_function_D.R")) # function D: NIR explained by PCs and VIS

Define model

MuPGLSMod3A <- Response ~ PC1 + PC2 + size + PC1*size + PC2*size + VIS

Model 4

PGLS Multiple Trees with 3 predictors + intercept

We run simplified models with only one pc for each spectral band including only the predictors that were significant in the full model. (because the visible absorptivity seems to be more correlated to PC1 and the NIR to PC2)

Source function

note that this function has to be adapted to the data frame and model on each case

source(here::here("Scripts/10_multiple_pgls_function_E.R")) # function E: simplified. One PC and size

Define model

MuPGLSMod4A <- Response ~ PC + size + PC*size

Model 5

PGLS Multiple Trees with 4 predictors + intercept

We run simplified models with only one pc for each spectral band including only the predictors that were signifficant in the full model (because the visible reflectance seems to be more correlated to PC1 and the NIR to PC2). However, for NIR we need to include VIS reflectance as a predictor

Source function

note that this function has to be adapted to the data frame and model on each case

source(here::here("Scripts/11_multiple_pgls_function_F.R")) # F: One PC, size and VIS as predictors

Define model

MuPGLSMod5A <- Response ~ PC + size + PC*size + VIS

Results per spectral band

The PGLS model tests if the correlations we found in the previous step remain after correcting by phylogeny. This model does not consider instraspecific variation. We averaged both location and absorptivity and obtained only one value per species.

TOT

MCC Tree

pglsmodTOTA <- pgls(TOT ~ PC1 + PC2 + size + PC1*size + PC2*size,
  data = comp_data, param.CI = 0.95, lambda = "ML",
  bounds = list(lambda=c(1e-3,1)) # modify bounds (defaults are c(1e-6,1))
)

FinMccTotcA <- as.numeric(round(summary(pglsmodTOTA)$coefficients[,1],3))
FinMccTotpA <- as.numeric(round(summary(pglsmodTOTA)$coefficients[,4],3))

None of the variables has an effect in the total absorptivity.

multiple trees

runsTOTA<-lapply(trees[trees_subset_min:trees_subset_max],
                pgls_runC,
                model=MuPGLSMod2A,
                dataset=ALLDataSetA)

## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH

dfTOTA <- ldply(runsTOTA, data.frame)

write.csv(dfTOTA, here::here("Data/FromCode/MuTPglsResultsTOTAbsorptivity.csv"))

FinTotMA <- HPDinterval(as.mcmc(dfTOTA))

length(dfTOTA[,1])

## [1] 964

VIS

full model

MCC Tree

pglsmodVISA <- pgls(VIS ~ PC1 + PC2 + size + PC1*size + PC2*size,
  data = comp_data, param.CI = 0.95, lambda = "ML"
)

FinMccViscA <- as.numeric(round(summary(pglsmodVISA)$coefficients[,1],3)) 
FinMccVispA <- as.numeric(round(summary(pglsmodVISA)$coefficients[,4],3))

summary(pglsmodVISA)

## 
## Call:
## pgls(formula = VIS ~ PC1 + PC2 + size + PC1 * size + PC2 * size, 
##     data = comp_data, lambda = "ML", param.CI = 0.95)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.4374 -1.9677 -0.2706  2.2468  5.6877 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.963
##    lower bound : 0.000, p = 0.54376
##    upper bound : 1.000, p = 0.40483
##    95.0% CI   : (NA, NA)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##             Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)  85.4669    20.5219  4.1647 0.0004787 ***
## PC1           3.4397     7.3423  0.4685 0.6445140    
## PC2           5.6063     9.4504  0.5932 0.5596739    
## size         -1.6558     7.9155 -0.2092 0.8364250    
## PC1:size     -1.6442     3.3607 -0.4893 0.6299788    
## PC2:size     -2.3463     3.7321 -0.6287 0.5366743    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.966 on 20 degrees of freedom
## Multiple R-squared: 0.0246,  Adjusted R-squared: -0.2193 
## F-statistic: 0.1009 on 5 and 20 DF,  p-value: 0.9908

multiple trees

runsVISA<-lapply(trees[trees_subset_min:trees_subset_max],
                pgls_runC,
                model=MuPGLSMod2A,
                dataset=VISDataSetA)

## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH

dfVISA <- ldply (runsVISA, data.frame)

write.csv(dfVISA, here::here("Data/FromCode/MuTPglsResultsVISAbsorptivity.csv"))

FinVisMA <- HPDinterval(as.mcmc(dfVISA))

length(dfVISA[,1])

## [1] 965

simple model

MCC Tree

spglsmodVISA <- pgls(VIS ~ PC1 + size + PC1*size,
  data = comp_data, param.CI = 0.95, lambda = "ML"
)

summary(spglsmodVISA)

## 
## Call:
## pgls(formula = VIS ~ PC1 + size + PC1 * size, data = comp_data, 
##     lambda = "ML", param.CI = 0.95)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.5545 -1.7529 -0.0548  2.3300  5.8262 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.961
##    lower bound : 0.000, p = 0.55608
##    upper bound : 1.000, p = 0.3721
##    95.0% CI   : (NA, NA)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##             Estimate Std. Error t value   Pr(>|t|)    
## (Intercept)  78.4337    15.9803  4.9082 0.00006582 ***
## PC1           1.0890     5.3206  0.2047     0.8397    
## size          1.3355     5.8339  0.2289     0.8211    
## PC1:size     -0.5272     2.3801 -0.2215     0.8267    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.809 on 22 degrees of freedom
## Multiple R-squared: 0.004095,    Adjusted R-squared: -0.1317 
## F-statistic: 0.03016 on 3 and 22 DF,  p-value: 0.9927

sFinMccViscA <- as.numeric(round(summary(spglsmodVISA)$coefficients[,1],3)) 
sFinMccVispA <- as.numeric(round(summary(spglsmodVISA)$coefficients[,4],3))

multiple trees

First step is to modify the data frame.

SimpleVISA <- 
  VISDataSetA %>% 
  dplyr::select(phylogeny_name, Response, size, PC1) %>% 
  dplyr::rename("PC" = PC1)

SimpleVISA <- as.data.frame(SimpleVISA)

Now run the model

runsVISsimpleA<-lapply(trees[trees_subset_min:trees_subset_max],
                      pgls_runE,
                      model=MuPGLSMod4A,
                      dataset=SimpleVISA)

## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH

dfVISsimpleA <- ldply(runsVISsimpleA, data.frame)

FinVisSA <- HPDinterval(as.mcmc(dfVISsimpleA))

length(dfVISsimpleA[,1])

## [1] 970

NIR

Full model

MCC Tree

pglsmodNIRA <- pgls(NIR ~ PC1 + PC2 + size + PC1*size + PC2*size + VIS,
  data = comp_data, param.CI = 0.95, lambda = "ML"
)

FinMccNircA <- as.numeric(round(summary(pglsmodNIRA)$coefficients[,1],3)) 
FinMccNirpA <- as.numeric(round(summary(pglsmodNIRA)$coefficients[,4],3))

Results after correcting by phylogeny:

None of the variables has an effect in the NIR absorptivity.

multiple trees

runsNIRA<-lapply(trees[trees_subset_min:trees_subset_max],
                pgls_runD,
                model=MuPGLSMod3A,
                dataset=NmoDataSetA)

## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH

dfNIRA <- ldply(runsNIRA, data.frame)

write.csv(dfNIRA, here::here("Data/FromCode/MuTPglsResultsNmoAbsorptivity.csv"))

FinNirMA <- HPDinterval(as.mcmc(dfNIRA))

length(dfNIRA[,1])

## [1] 999

Simple model

MCC Tree

spglsmodNIRA <- pgls(NIR ~ PC2 + size + PC2*size + VIS,
  data = comp_data, param.CI = 0.95, lambda = "ML"
)

sFinMccNircA <- as.numeric(round(summary(spglsmodNIRA)$coefficients[,1],3)) 
sFinMccNirpA <- as.numeric(round(summary(spglsmodNIRA)$coefficients[,4],3))

multiple trees

First step is to modify the data frame.

SimpleNIRA <- 
  NmoDataSetA %>% 
  dplyr::select(phylogeny_name, Response, size, VIS, PC2) %>% 
  dplyr::rename("PC" = PC2)

SimpleNIRA <- as.data.frame(SimpleNIRA)
head(SimpleNIRA)

phylogeny_name	Response	size	VIS	PC
Anomala_antigua	37.5	1.77	88.5	-1.9
Anoplognathus_aureus	9.56	1.71	55.9	-0.412
Anoplognathus_brunnipennis	23.6	2.33	68.9	0.461
Anoplognathus_concolor	14.6	2.17	65.4	1.43
Anoplognathus_hirsutus	35.7	2.32	79.5	1.05
Anoplognathus_macleayi	36	3.03	75.9	-3.88

Now run the model

runsNIRsimpleA<-lapply(trees[trees_subset_min:trees_subset_max],
                      pgls_runF,
                      model=MuPGLSMod5A,
                      dataset=SimpleNIRA)

## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH

dfNIRsimpleA <- ldply(runsNIRsimpleA, data.frame)

FinNirSA <- HPDinterval(as.mcmc(dfNIRsimpleA))

length(dfNIRsimpleA[,1])

## [1] 1000

FRS

Residuals between NIR ~ VIS corrected by phylogeny. Extracted as model$residuals

Full model

MCC Tree

pglsmodFRSA <- pgls(FRS ~ PC1 + PC2 + size + PC1*size + PC2*size,
  data = comp_data, param.CI = 0.95, lambda = "ML"
)

FinMccFRScA <- as.numeric(round(summary(pglsmodFRSA)$coefficients[,1],3)) 
FinMccFRSpA <- as.numeric(round(summary(pglsmodFRSA)$coefficients[,4],3))

multiple trees

Residuals between NIR ~ VIS corrected by phylogeny. Extracted as model$residuals

runsFRSA<-lapply(trees[trees_subset_min:trees_subset_max],
                pgls_runC,
                model=MuPGLSMod2A,
                dataset=FRSDataSetA)

## ERROR : Problem with optim:52ERROR: ABNORMAL_TERMINATION_IN_LNSRCH

dfFRSA <- ldply(runsFRSA, data.frame)

write.csv(dfFRSA, here::here("Data/FromCode/MuTPglsFRSultsFRSAbsorptivity.csv"))

FinFRSMA <- HPDinterval(as.mcmc(dfFRSA))

length(dfFRSA[,1])

## [1] 1000

Simple model

MCC Tree

spglsmodFRSA <- pgls(FRS ~ PC2 + size + PC2*size,
  data = comp_data, param.CI = 0.95, lambda = "ML"
)

sFinMccFRScA <- as.numeric(round(summary(spglsmodFRSA)$coefficients[,1],3)) 
sFinMccFRSpA <- as.numeric(round(summary(spglsmodFRSA)$coefficients[,4],3))

multiple trees

First step is to modify the data frame.

SimpleFRSA <- 
  FRSDataSetA %>% 
  dplyr::select(phylogeny_name, Response, size, PC2) %>% 
  dplyr::rename(PC = PC2)

SimpleFRSA <- as.data.frame(SimpleFRSA)

Now run the model

runsFRSsimpleA <- lapply(trees[trees_subset_min:trees_subset_max],
                      pgls_runE,
                      model=MuPGLSMod4A,
                      dataset=SimpleFRSA) 

dfFRSsimpleA <- ldply(runsFRSsimpleA, data.frame)

FinFRSSA <-HPDinterval(as.mcmc(dfFRSsimpleA))

length(dfFRSsimpleA[,1])

## [1] 1001

PRS

Residuals between NIR ~ VIS corrected by phylogeny. Phylogenetic signal eliminated by the PGLS correction. Extracted as model$phyres.

Full model

MCC Tree

pglsmodPRSA <- pgls(FRSP ~ PC1 + PC2 + size + PC1*size + PC2*size,
  data = comp_data, param.CI = 0.95, lambda = "ML"
)

FinMccPRScA <- as.numeric(round(summary(pglsmodPRSA)$coefficients[,1],3)) 
FinMccPRSpA <- as.numeric(round(summary(pglsmodPRSA)$coefficients[,4],3))

multiple trees

Residuals between NIR ~ VIS corrected by phylogeny. Phylogenetic signal eliminated by the PGLS correction. Extracted as model$phyres.

runsPRSA<-lapply(trees[trees_subset_min:trees_subset_max],
                pgls_runC,
                model=MuPGLSMod2A,
                dataset=PRSDataSetA) 

dfPRSA <- ldply(runsPRSA, data.frame)

write.csv(dfPRSA, here::here("Data/FromCode/MuTPglsFRSultsPRSAbsorptivity.csv"))

FinPRSMA <- HPDinterval(as.mcmc(dfPRSA))

length(dfPRSA[,1])

## [1] 1001

Simple model

MCC Tree

spglsmodPRSA <- pgls(FRSP ~ PC2 + size + PC2*size,
  data = comp_data, param.CI = 0.95, lambda = "ML"
)

sFinMccPRScA <- as.numeric(round(summary(spglsmodPRSA)$coefficients[,1],3)) 
sFinMccPRSpA <- as.numeric(round(summary(spglsmodPRSA)$coefficients[,4],3))

multiple trees

First step is to modify the data frame.

SimplePRSA <- 
  PRSDataSetA %>% 
  dplyr::select(phylogeny_name, Response, size, PC2) %>% 
  dplyr::rename(PC = PC2)

SimplePRSA <- as.data.frame(SimplePRSA)

Now run the model

runsPRSsimpleA <- lapply(trees[trees_subset_min:trees_subset_max],
                      pgls_runE,
                      model=MuPGLSMod4A,
                      dataset=SimplePRSA) 

dfPRSsimpleA <- ldply(runsPRSsimpleA, data.frame)

FinPRSSA <-HPDinterval(as.mcmc(dfPRSsimpleA))

length(dfPRSsimpleA[,1])

## [1] 1001

Final results

Full models

MCC Tree

PredictorA <- c("PC1","PC1 p-val" , 
               "PC2", "PC2 p-val", 
               "Size", "Size p-val",
               "VIS", "VIS p-val",
               "PC1:size", "PC1:size p-val", 
               "PC2:size", "PC2:size p-val")

# Arrange the vectors to build the data frame
# New vectors alternate coefficient and p-values to match "predictors"

TotiA <- as.character(as.numeric(c(rbind((c(
            FinMccTotcA[2:4],"NA",FinMccTotcA[5:6])),
                (c(FinMccTotpA[2:4],"NA",FinMccTotpA[5:6]))))))

VisiA <- as.character(as.numeric(c(rbind((c(
            FinMccViscA[2:4],"NA",FinMccViscA[5:6])),
                (c(FinMccVispA[2:4],"NA",FinMccVispA[5:6]))))))

NiriA <- as.character(as.numeric(c(rbind(
            FinMccNircA[2:7], FinMccNirpA[2:7]))))

FRSiA <- as.character(as.numeric(c(rbind((c(
            FinMccFRScA[2:4],"NA",FinMccFRScA[5:6])),
                (c(FinMccFRSpA[2:4],"NA",FinMccFRSpA[5:6]))))))

PRSiA <- as.character(as.numeric(c(rbind((c(
            FinMccPRScA[2:4],"NA",FinMccPRScA[5:6])),
                (c(FinMccPRSpA[2:4],"NA",FinMccPRSpA[5:6]))))))
  
# This section colours the significant p-values and their coefficients

#Visi[1] <- cell_spec(Visi[1],bold = TRUE, color="#D40481")
#Visi[2] <- cell_spec(Visi[2],bold = TRUE)

#Visi[5] <- cell_spec(Visi[5],bold = TRUE, color="#D40481")
#Visi[6] <- cell_spec(Visi[6],bold = TRUE)

#Visi[9] <- cell_spec(Visi[9],bold = TRUE, color="#D40481")
#Visi[10] <- cell_spec(Visi[10],bold = TRUE)

#Niri[3] <- cell_spec(Niri[3],bold = TRUE, color="#D40481")
#Niri[4] <- cell_spec(Niri[4],bold = TRUE)

#Niri[7] <- cell_spec(Niri[7],bold = TRUE, color="#D40481")
#Niri[8] <- cell_spec(Niri[8],bold = TRUE)

#Niri[11] <- cell_spec(Niri[11],bold = TRUE, color="#D40481")
#Niri[12] <- cell_spec(Niri[12],bold = TRUE)

#FRSi[3] <- cell_spec(FRSi[3],bold = TRUE, color="#D40481")
#FRSi[4] <- cell_spec(FRSi[4],bold = TRUE)

#FRSi[11] <- cell_spec(FRSi[11],bold = TRUE, color="#D40481")
#FRSi[12] <- cell_spec(FRSi[12],bold = TRUE)

# Assemble the table

ResultspglsA <- data.frame("Predictor" = PredictorA,  
                          "TOT" = TotiA,
                          "VIS" = VisiA,
                          "NIR" = NiriA,
                          "FRS" = FRSiA,
                          "PRS" = PRSiA) 

ResultspglsA %>% 
  kbl(align ="c", escape = FALSE) %>%
  kable_classic() %>% 
  add_indent(c(1, 3, 5, 7, 9))

Predictor	TOT	VIS	NIR	FRS	PRS
PC1	6.729	3.44	-1.563	-1.532	-0.336
PC1 p-val	0.516	0.645	0.732	0.728	0.83
PC2	-3.478	5.606	-11.891	-11.861	-3.073
PC2 p-val	0.767	0.56	0.047	0.041	0.095
Size	4.856	-1.656	8.861	8.854	2.537
Size p-val	0.607	0.836	0.067	0.061	0.087
VIS	NA	NA	1.354	NA	NA
VIS p-val	NA	NA	0	NA	NA
PC1:size	-3.29	-1.644	0.616	0.601	0.131
PC1:size p-val	0.476	0.63	0.768	0.765	0.851
PC2:size	0.657	-2.346	4.334	4.321	1.083
PC2:size p-val	0.889	0.537	0.067	0.059	0.14

multiple trees

rvTOTlA <- as.character(
            round(as.numeric(c(
              FinTotMA[1,1],FinTotMA[3,1],FinTotMA[4,1],FinTotMA[6,1],
              FinTotMA[7,1],FinTotMA[9,1],"NA","NA",
              FinTotMA[10,1],FinTotMA[12,1],
              FinTotMA[13,1],FinTotMA[15,1])),3))

rvTOTuA <- as.character(
            format(round(as.numeric(c(
              FinTotMA[1,2],FinTotMA[3,2],FinTotMA[4,2],FinTotMA[6,2],
              FinTotMA[7,2],FinTotMA[9,2],"NA","NA",
              FinTotMA[10,2],FinTotMA[12,2],
              FinTotMA[13,2],FinTotMA[15,2])),3),nsmall=3))


rvVISlA <- as.character(
          format(round(as.numeric(c(
              FinVisMA[1,1],FinVisMA[3,1],FinVisMA[4,1],FinVisMA[6,1],
              FinVisMA[7,1],FinVisMA[9,1],"NA","NA",
              FinVisMA[10,1],FinVisMA[12,1],
              FinVisMA[13,1],FinVisMA[15,1])),3),nsmall=3))

rvVISuA <- as.character(
            format(round(as.numeric(c(
              FinVisMA[1,2],FinVisMA[3,2],FinVisMA[4,2],FinVisMA[6,2],
              FinVisMA[7,2],FinVisMA[9,2],"NA","NA",
              FinVisMA[10,2],FinVisMA[12,2],
              FinVisMA[13,2],FinVisMA[15,2])),3),nsmall=3))


rvNIRlA <- as.character(
            format(round(as.numeric(c(
              FinNirMA[1,1],FinNirMA[3,1],FinNirMA[4,1],FinNirMA[6,1],
              FinNirMA[7,1],FinNirMA[9,1],FinNirMA[10,1],FinNirMA[12,1],
              FinNirMA[13,1],FinNirMA[15,1], FinNirMA[16,1],FinNirMA[18,1]
              )),3),nsmall=3))

rvNIRuA <- as.character(
           format(round(as.numeric(c(
              FinNirMA[1,2],FinNirMA[3,2],FinNirMA[4,2],FinNirMA[6,2],
              FinNirMA[7,2],FinNirMA[9,2],FinNirMA[10,2],FinNirMA[12,2],
              FinNirMA[13,2],FinNirMA[15,2], FinNirMA[16,2],FinNirMA[18,2]
              )),3),nsmall=3))

rvFRSlA <- as.character(
            format(round(as.numeric(c(
              FinFRSMA[1,1],FinFRSMA[3,1],FinFRSMA[4,1],FinFRSMA[6,1],
              FinFRSMA[7,1],FinFRSMA[9,1],"NA","NA",
              FinFRSMA[10,1],FinFRSMA[12,1],
              FinFRSMA[13,1],FinFRSMA[15,1])),3),nsmall=3))

rvFRSuA <- as.character(
            format(round(as.numeric(c(
              FinFRSMA[1,2],FinFRSMA[3,2],FinFRSMA[4,2],FinFRSMA[6,2],
              FinFRSMA[7,2],FinFRSMA[9,2],"NA","NA",
              FinFRSMA[10,2],FinFRSMA[12,2],
              FinFRSMA[13,2],FinFRSMA[15,2])),3),nsmall=3))

rvPRSlA <- as.character(
            format(round(as.numeric(c(
              FinPRSMA[1,1],FinPRSMA[3,1],FinPRSMA[4,1],FinPRSMA[6,1],
              FinPRSMA[7,1],FinPRSMA[9,1],"NA","NA",
              FinPRSMA[10,1],FinPRSMA[12,1],
              FinPRSMA[13,1],FinPRSMA[15,1])),3),nsmall=3))

rvPRSuA <- as.character(
            format(round(as.numeric(c(
              FinPRSMA[1,2],FinPRSMA[3,2],FinPRSMA[4,2],FinPRSMA[6,2],
              FinPRSMA[7,2],FinPRSMA[9,2],"NA","NA",
              FinPRSMA[10,2],FinPRSMA[12,2],
              FinPRSMA[13,2],FinPRSMA[15,2])),3),nsmall=3))


rvTOTiA<-paste(rvTOTlA,rvTOTuA, sep = " ; ")
rvVISiA<-paste(rvVISlA,rvVISuA, sep = " ; ")
rvNIRiA<-paste(rvNIRlA,rvNIRuA, sep = " ; ")
rvFRSiA<-paste(rvFRSlA,rvFRSuA, sep = " ; ")
rvPRSiA<-paste(rvPRSlA,rvPRSuA, sep = " ; ")

#rvVISi[1] <- cell_spec(rvVISi[1],bold = TRUE, color="#367D91")
#rvVISi[2] <- cell_spec(rvVISi[2],bold = TRUE, color="#B2B9BF")

#rvVISi[5] <- cell_spec(rvVISi[5],bold = TRUE, color="#D40481")
#rvVISi[6] <- cell_spec(rvVISi[6],bold = TRUE)

#rvVISi[9] <- cell_spec(rvVISi[9],bold = TRUE, color="#367D91")
#rvVISi[10] <- cell_spec(rvVISi[10],bold = TRUE,  color="#B2B9BF")

#rvNIRi[3] <- cell_spec(rvNIRi[3],bold = TRUE, color="#D40481")
#rvNIRi[4] <- cell_spec(rvNIRi[4],bold = TRUE)

#rvNIRi[7] <- cell_spec(rvNIRi[7],bold = TRUE, color="#D40481")
#rvNIRi[8] <- cell_spec(rvNIRi[8],bold = TRUE)

#rvNIRi[11] <- cell_spec(rvNIRi[11],bold = TRUE, color="#D40481")
#rvNIRi[12] <- cell_spec(rvNIRi[12],bold = TRUE)

#rvFRSi[3] <- cell_spec(rvFRSi[3],bold = TRUE, color="#D40481")
#rvFRSi[4] <- cell_spec(rvFRSi[4],bold = TRUE)

#rvFRSi[11] <- cell_spec(rvFRSi[11],bold = TRUE, color="#D40481")
#rvFRSi[12] <- cell_spec(rvFRSi[12],bold = TRUE)

ResultspglsMultipleA <- data.frame("Predictor" = PredictorA,  
                          "TOT" = rvTOTiA,
                          "VIS" = rvVISiA,
                          "NIR" = rvNIRiA,
                          "FRS" = rvFRSiA,
                          "PRS" = rvPRSiA) 

ResultspglsMultipleA %>% 
  kbl(align ="c", escape = FALSE) %>%
  kable_classic() %>% 
  add_indent(c(1, 3, 5, 7, 9))

Predictor	TOT	VIS	NIR	FRS	PRS
PC1	1.375 ; 6.728	2.048 ; 6.212	-2.165 ; -0.619	-2.209 ; -0.685	-0.336 ; -0.336
PC1 p-val	0.491 ; 0.885	0.426 ; 0.784	0.626 ; 0.888	0.610 ; 0.874	0.830 ; 0.830
PC2	-4.602 ; 1.473	1.332 ; 6.757	-13.665 ; -10.700	-13.836 ; -10.758	-3.073 ; -3.073
PC2 p-val	0.705 ; 1.000	0.475 ; 0.890	0.019 ; 0.077	0.016 ; 0.069	0.095 ; 0.095
Size	0.117 ; 4.856	-3.854 ; 0.345	7.631 ; 9.380	7.712 ; 9.378	2.537 ; 2.537
Size p-val	0.607 ; 0.975	0.614 ; 0.964	0.049 ; 0.108	0.042 ; 0.096	0.087 ; 0.087
VIS	NA ; NA	NA ; NA	1.300 ; 1.365	NA ; NA	NA ; NA
VIS p-val	NA ; NA	NA ; NA	0.000 ; 0.000	NA ; NA	NA ; NA
PC1:size	-3.289 ; -0.802	-2.948 ; -0.957	0.193 ; 0.904	0.201 ; 0.919	0.131 ; 0.131
PC1:size p-val	0.461 ; 0.856	0.391 ; 0.756	0.645 ; 0.918	0.627 ; 0.909	0.851 ; 0.851
PC2:size	-0.987 ; 1.237	-2.890 ; -0.907	3.856 ; 4.960	3.844 ; 5.001	1.083 ; 1.083
PC2:size p-val	0.794 ; 0.999	0.437 ; 0.815	0.030 ; 0.109	0.025 ; 0.098	0.140 ; 0.140

Plots

Plot PC1

toplotAPC1A <- 
  Cons1aggA %>% 
  dplyr::select (2, 4:6, 8) %>% 
  gather (key = spectralband, 
          value = Absorptivity, -PC1, -size) %>%    # prepare for ggplot
  mutate(PC1n = PC1 * -1) %>% 
  mutate(spectralband = replace(spectralband, spectralband == "FRS", "Residuals"))

# Arrange facets in custom order
toplotAPC1A$spectralband <- factor(toplotAPC1A$spectralband, levels = c("TOT", "VIS", "Residuals"))

ggplot(toplotAPC1A, aes(x = PC1n, y = Absorptivity, colour = spectralband)) +
  geom_point(size=2,alpha=0.4) +
  geom_point(size=2, pch=21, colour="black", alpha=0.7) +
  theme_bw() +
  facet_wrap(~ spectralband) +
  scale_colour_manual(values = c("TOT" = "#FFB000", "VIS" =  "#648FFF", "Residuals" = "#FF2C85")) +
  theme(legend.position = "none")

??scale_color_manual

Plot PC2

toplotAPC2A <- 
  Cons1aggA %>% 
  dplyr::select (3:6, 8) %>% 
  gather (key = spectralband, 
          value = Absorptivity, -PC2, -size)  %>%  # prepare for ggplot
  mutate(PC2n = PC2 * -1) %>% 
  mutate(spectralband = replace(spectralband, spectralband == "FRS", "Residuals"))

# Arrange facets in custom order
toplotAPC2A$spectralband <- factor(toplotAPC2A$spectralband, levels = c("TOT", "VIS", "Residuals"))

ggplot(toplotAPC2A, aes(x = PC2n, y = Absorptivity, colour = spectralband))+
  geom_point(size=2,alpha=0.4)+
  geom_point(size=2, pch=21, colour="black", alpha=0.7)+
  theme_bw()+
  facet_wrap(~ spectralband) +
  scale_colour_manual(values = c("TOT" = "#FFB000", "VIS" =  "#648FFF", "Residuals" = "#FF2C85")) +
  theme(legend.position = "none")

Pagel’s lambda

pagelTOTAbs <- Cons1aggA$TOT # Define which trait you want to test
names(pagelTOTAbs) <- rownames(Cons1aggA) # Row names = tree tips
phylosig(MCCtree, pagelTOTAbs, method = "lambda", test = TRUE, nsim = 999)

## 
## Phylogenetic signal lambda : 0.971361 
## logL(lambda) : -100.778 
## LR(lambda=0) : 1.56805 
## P-value (based on LR test) : 0.21049

# nsim = 999 means testing with 999 randomisations

pagelVISAbs <- Cons1aggA$VIS # Define which trait you want to test
names(pagelVISAbs) <- rownames(Cons1aggA) # Row names = tree tips
phylosig(MCCtree, pagelVISAbs, method = "lambda", test = TRUE, nsim = 999)

## 
## Phylogenetic signal lambda : 0.964661 
## logL(lambda) : -95.2276 
## LR(lambda=0) : 0.93322 
## P-value (based on LR test) : 0.334028

# nsim = 999 means testing with 999 randomisations

pagelNIRAbs <- Cons1aggA$NIR # Define which trait you want to test
names(pagelNIRAbs) <- rownames(Cons1aggA) # Row names = tree tips
phylosig(MCCtree, pagelNIRAbs, method = "lambda", test = TRUE, nsim = 999)

## 
## Phylogenetic signal lambda : 0.973529 
## logL(lambda) : -105.985 
## LR(lambda=0) : 2.193 
## P-value (based on LR test) : 0.138639

# nsim = 999 means testing with 999 randomisations

pagelFRSAbs <- Cons1aggA$FRS # Define which trait you want to test
names(pagelFRSAbs) <- rownames(Cons1aggA) # Row names = tree tips
phylosig(MCCtree, pagelFRSAbs, method = "lambda", test = TRUE, nsim = 999)

## 
## Phylogenetic signal lambda : 0.842634 
## logL(lambda) : -86.0733 
## LR(lambda=0) : 5.35705 
## P-value (based on LR test) : 0.0206386

# nsim = 999 means testing with 999 randomisations

pagelPRSAbs <- Cons1aggA$FRSP # Define which trait you want to test
names(pagelPRSAbs) <- rownames(Cons1aggA) # Row names = tree tips
phylosig(MCCtree, pagelPRSAbs, method = "lambda", test = TRUE, nsim = 999)

## 
## Phylogenetic signal lambda : 0.0000760803 
## logL(lambda) : -56.1662 
## LR(lambda=0) : -0.00103419 
## P-value (based on LR test) : 1

# nsim = 999 means testing with 999 randomisations

Conclusions

For Visible reflectivity the interaction between Size and PC1 was signifficant. Lower PC1 = higher humidity (vapour, rain and clouds) lower aridity. Smaller beetles reflect more visible light in humid environments

For NIR reflectivity the interaction between Size and PC2 was signifficant. Lower PC2 = higher solar radiation, higher max temp, more days above 35 and more aridity Smaller beetles reflect less NIR light in hot/arid environments.

Absorptivity PGLS

Setting up data

Libraries

Data Sets

Setting up models

Model 1

Model 2

Model 3

Model 4

Model 5

Results per spectral band

TOT

MCC Tree

multiple trees

VIS

full model

MCC Tree

multiple trees

simple model

MCC Tree

multiple trees

NIR

Full model

MCC Tree

multiple trees

Simple model

MCC Tree

multiple trees

FRS

Full model

MCC Tree

multiple trees

Simple model

MCC Tree

multiple trees

PRS

Full model

MCC Tree

multiple trees

Simple model

MCC Tree

multiple trees

Final results

Full models

MCC Tree

multiple trees

Plots

Plot PC1

Plot PC2

Pagel’s lambda

Conclusions