We repeated the analysis for degree of polarisation and the raw RGB value sin the right handed polarisation filter excluding A. prasinus and the species belonging to the genus Xylonichus, Paraschizognathus because, unlike many scarabs, these species are known to have a bilayered mechanism that is not chiral.

Setting up

Libraries

Libraries sourced from an additional script.

source(here::here("Scripts/MacroEcol_1_Libraries.R"))

Data sets

Import polarisation data.

# elytra
PolBySpp <- read.csv(here::here("Data/FromCode/PolarzElytraBySpp.csv"))[-1] %>% 
  dplyr::arrange(phylogeny_name) %>%
  dplyr::filter(
    !str_detect(phylogeny_name, "Xylonichus_sp"),
    !str_detect(phylogeny_name, "Anoplognathus_prasinus"),
    !str_detect(phylogeny_name, "Paraschizognathus_ocularis"),
    !str_detect(phylogeny_name, "Paraschizognathus_prasinus"),
    !str_detect(phylogeny_name, "Paraschizognathus_olivaceous")
  )

PolByInd <- read.csv(here::here("Data/FromCode/PolarzElytraByInd.csv"))[-1] %>% 
  dplyr::filter(
    !str_detect(ind, "xyls"),
    !str_detect(ind, "prsi"),
    !str_detect(ind, "xyle"),
    !str_detect(ind, "pczo"),
    !str_detect(ind, "pczp"),
    !str_detect(ind, "pczv")
  )


# elytra and pronotum
#BeetlePol <- read.csv(here::here("Data/FromCode/PolarzElytraPronot.csv"))[-1]

All statistical analysis were done at the species level (polarisation data by species). Plots are produced with the polarisation data by individual.

And we imported the phylogeny from a .nwk file which contains 2000 phylogenetic trees. This allows us to test our hypothesis in multiple trees and not only on the MCC tree. This way, we account for any uncertainty in unresolved tips of our phylogeny.

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

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

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

To explore the correlations between ecological variables and polarisation

# PC Values
PCValuesF <- read.csv(here::here("Data/FromCode/PCsbySpp.csv"))[-1] %>% 
  dplyr::rename("PC1" = PC1_1, "PC2" = PC2_1) %>% 
  dplyr::arrange(phylogeny_name) %>%
  dplyr::filter(
    !str_detect(phylogeny_name, "Xylonichus_sp"),
    !str_detect(phylogeny_name, "Anoplognathus_prasinus"),
    !str_detect(phylogeny_name, "Paraschizognathus_ocularis"),
    !str_detect(phylogeny_name, "Paraschizognathus_prasinus"),
    !str_detect(phylogeny_name, "Paraschizognathus_olivaceous")
  )

# merge with polarisation
PolPCs <- merge(PCValuesF, PolBySpp)

# adapt to merge with phylogeny
PolPCs2 <- as.data.frame(PolPCs) # create new data frame
rownames(PolPCs2) <- PolPCs2[, 1] # make species the row names
PolPCsRows <- PolPCs2[, 2:length(PolPCs2)] # eliminate spp name (redundant)

To explore correlations with reflectivity

# Reflectivity data consolidated by species
ReflBySpp <- read.csv(here::here("Data/FromCode/ConsolidatedReflectivitySpp.csv"))[-1] %>% 
  dplyr::arrange(phylogeny_name) %>% 
  dplyr::filter(
    !str_detect(phylogeny_name, "Xylonichus_sp"),
    !str_detect(phylogeny_name, "Anoplognathus_prasinus"),
    !str_detect(phylogeny_name, "Paraschizognathus_ocularis"),
    !str_detect(phylogeny_name, "Paraschizognathus_prasinus"),
    !str_detect(phylogeny_name, "Paraschizognathus_olivaceous")
  )

# Merge with polarisation
PolRefl <- merge(ReflBySpp, PolBySpp)

# adapt to merge with phylogeny
PolRefl2 <- as.data.frame(PolRefl) # create new data frame
rownames(PolRefl2) <- PolRefl2[, 1] # make species the row names
PolReflRows <- PolRefl2[, 2:length(PolRefl2)] # eliminate spp name (redundant)

To explore correlations with absorptivity

AbsBySpp <- read.csv(here::here("Data/FromCode/ConsolidatedAbsorptivitySpp.csv"))[-1] %>%
  dplyr::arrange(phylogeny_name) %>%
  dplyr::filter(
    !str_detect(phylogeny_name, "Xylonichus_sp"),
    !str_detect(phylogeny_name, "Anoplognathus_prasinus"),
    !str_detect(phylogeny_name, "Paraschizognathus_ocularis"),
    !str_detect(phylogeny_name, "Paraschizognathus_prasinus"),
    !str_detect(phylogeny_name, "Paraschizognathus_olivaceous")
  )

# merge with polarisation
PolAbs <- merge(AbsBySpp, PolBySpp)

# adapt to merge with phylogeny
PolAbs2 <- as.data.frame(PolAbs) # create new data frame
rownames(PolAbs2) <- PolAbs2[, 1] # make species the row names
PolAbsRows <- PolAbs2[, 2:length(PolAbs2)] # eliminate spp name (redundant)

To explore correlations with reflectivity only in Right handed polarisation (as a proxy of the amount of melanin/pigments)

# Import equalised values for right handed polarisation
RHand0<-read.csv(here::here("Data/FromCode/PolarzElytra_LERGBValues.csv"))[-1]%>% 
  dplyr::filter(Type == "RCP") %>% 
  dplyr::select(-Type) %>% 
  dplyr::mutate(spp= substr(ind,1,4)) %>% 
  dplyr::group_by(spp) %>% 
  dplyr::summarise(
    meanRCP = mean(RGB_ave),
    sdRCP = sd(RGB_ave),
    )

# Join them with the phylogeny names
SppNames<-read.csv(here::here("Data/9_CodesAndSpecies.csv")) %>% 
  dplyr::mutate("spp" = substr(Ind,1,4)) %>% 
  dplyr::select(-Ind) %>% 
  distinct()


RHand <- merge(RHand0[,c(1,2)],SppNames) %>% 
 dplyr::arrange(phylogeny_name) %>%
  dplyr::filter(
    !str_detect(phylogeny_name, "Xylonichus_sp"),
    !str_detect(phylogeny_name, "Anoplognathus_prasinus"),
    !str_detect(phylogeny_name, "Paraschizognathus_ocularis"),
    !str_detect(phylogeny_name, "Paraschizognathus_prasinus"),
    !str_detect(phylogeny_name, "Paraschizognathus_olivaceous")
  )

# Name the right hand RGB "pol" to use the accessory scripts for multiple trees:
RHand <- 
  RHand %>% 
  dplyr::rename(Pol = meanRCP) %>% 
  dplyr::group_by(phylogeny_name) %>% 
  dplyr::summarise(Pol= mean(Pol))


# Combine with reflectivity

RHandRefl <- merge(ReflBySpp, RHand) 

RHandRefl2 <- as.data.frame(RHandRefl) # create new data frame
rownames(RHandRefl2) <- RHandRefl2[, 1] # make species the row names
RHandReflRows <- RHandRefl2[, 2:length(RHandRefl2)] # eliminate spp name 



# Combine with absorptivity

RHandAbs <- merge(AbsBySpp, RHand) 

RHandAbs2 <- as.data.frame(RHandAbs) # create new data frame
rownames(RHandAbs2) <- RHandAbs2[, 1] # make species the row names
RHandAbsRows <- RHandAbs2[, 2:length(RHandAbs2)] # eliminate spp name

We prune the species in the tree that are not present in the dataset

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


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

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

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

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

## [1] TRUE

identical(
  length(name.check(MCCtreeRefl, PolRefl2$phylogeny_name)$tree_not_data),
  length(PolRefl2$phylogeny_name)
)

## [1] TRUE

Optical properties

Pol Degree

This is an indication of “how much of the colour is circularly polarised to the left”. As a simplification, we assume that this is equivalent to most of the reflections produced by a chiral mechanism

Reflectivity

VIS

Question

Can polarisation predict visible reflectivity?

Correlation

cor.test(PolRefl$VIS, PolRefl$Pol) # by spp

## 
##  Pearson's product-moment correlation
## 
## data:  PolRefl$VIS and PolRefl$Pol
## t = 1.3748, df = 40, p-value = 0.1768
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.09782853  0.48503581
## sample estimates:
##       cor 
## 0.2124187

Linear model

ModL1 <- lm(PolRefl$VIS ~ PolRefl$Pol) # by spp
summary(ModL1)

## 
## Call:
## lm(formula = PolRefl$VIS ~ PolRefl$Pol)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.4974  -4.9589   0.4319   3.2515  18.8216 
## 
## Coefficients:
##             Estimate Std. Error t value       Pr(>|t|)    
## (Intercept)   13.374      1.308  10.228 0.000000000001 ***
## PolRefl$Pol   10.387      7.555   1.375          0.177    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.831 on 40 degrees of freedom
## Multiple R-squared:  0.04512,    Adjusted R-squared:  0.02125 
## F-statistic:  1.89 on 1 and 40 DF,  p-value: 0.1768

PGLS

Create the data frame of comparative data

comp_data <- comparative.data(
  phy = MCCtreeRefl,
  data = PolRefl, # aggregated data without the spp in the row names
  names.col = "phylogeny_name", # contains the column phylogeny name
  vcv = TRUE,
  na.omit = FALSE, warn.dropped = TRUE
)

Establish the model

pglsPDVIS <- pgls(VIS ~ Pol,
  data = comp_data, param.CI = 0.95, lambda = "ML"
)

summary(pglsPDVIS)

## 
## Call:
## pgls(formula = VIS ~ Pol, data = comp_data, lambda = "ML", param.CI = 0.95)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.7622 -1.6248 -0.1015  1.3323  7.4350 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.875
##    lower bound : 0.000, p = 0.15335
##    upper bound : 1.000, p = 0.32073
##    95.0% CI   : (NA, NA)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)   9.4574     5.9519  1.5890   0.1199
## Pol           9.3633    10.0150  0.9349   0.3554
## 
## Residual standard error: 2.514 on 40 degrees of freedom
## Multiple R-squared: 0.02138, Adjusted R-squared: -0.00308 
## F-statistic: 0.8741 on 1 and 40 DF,  p-value: 0.3554

PGLS multiple trees

Function

Note that the function has to be modified for the predictors and data frame that we are using

source(here::here("Scripts/6_multiple_pgls_function_A.R")) # script A is for Polarization as a response

Model

MuPGLSPol3 <- Response ~ Pol

Dataset

PolRefl3 <-
    PolReflRows %>% # phylogeny names in rows
    dplyr::select(Pol, VIS) %>% 
    dplyr::rename("Response" = VIS)

PolRefl3 <- as.data.frame(PolRefl3)

Apply

runsPol3 <- lapply(trees[trees_subset_min:trees_subset_max],
                   pgls_runA, model = MuPGLSPol3, dataset = PolRefl3)

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

dfPol3 <- ldply(runsPol3, data.frame)

length(dfPol3[,1])

## [1] 1000

Get results

HPDinterval(as.mcmc(dfPol3))

##                lower      upper
## estimateP1 3.1343929 13.3374537
## ts1        0.3411226  1.3949166
## ps1        0.1568928  0.7054531
## attr(,"Probability")
## [1] 0.95

Conclusion

There’s no correlation between visible reflectivity and polarisation at the species level.

NIR

Question

Is polarisation correlated with NIR reflectivity?

Correlation

Species level:

cor.test(PolRefl$NIR, PolRefl$Pol) # by spp

## 
##  Pearson's product-moment correlation
## 
## data:  PolRefl$NIR and PolRefl$Pol
## t = -0.98424, df = 40, p-value = 0.3309
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.4372662  0.1575225
## sample estimates:
##        cor 
## -0.1537708

Linear model

ModL2 <- lm(PolRefl$NIR ~ PolRefl$Pol + PolRefl$VIS) # by spp
summary(ModL2)

## 
## Call:
## lm(formula = PolRefl$NIR ~ PolRefl$Pol + PolRefl$VIS)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.9085 -2.8564  0.3102  3.1581 11.5848 
## 
## Coefficients:
##             Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)   24.854      1.588  15.653 < 0.0000000000000002 ***
## PolRefl$Pol  -20.238      4.938  -4.099             0.000204 ***
## PolRefl$VIS    1.063      0.101  10.531     0.00000000000058 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.363 on 39 degrees of freedom
## Multiple R-squared:  0.746,  Adjusted R-squared:  0.733 
## F-statistic: 57.26 on 2 and 39 DF,  p-value: 0.000000000002483

PGLS

The data frame of comparative data is the same as for visible light

pglsPDNIR <- pgls(NIR ~ Pol + VIS,
  data = comp_data, param.CI = 0.95, lambda = "ML"
)

summary(pglsPDNIR)

## 
## Call:
## pgls(formula = NIR ~ Pol + VIS, data = comp_data, lambda = "ML", 
##     param.CI = 0.95)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.1396 -0.6860  0.0745  0.7585  2.7823 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.000
##    lower bound : 0.000, p = 1    
##    upper bound : 1.000, p = 0.00023001
##    95.0% CI   : (NA, 0.832)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##              Estimate Std. Error t value              Pr(>|t|)    
## (Intercept)  24.85385    1.58780 15.6530 < 0.00000000000000022 ***
## Pol         -20.23845    4.93766 -4.0988             0.0002036 ***
## VIS           1.06337    0.10098 10.5309    0.0000000000005795 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.048 on 39 degrees of freedom
## Multiple R-squared: 0.746,   Adjusted R-squared: 0.733 
## F-statistic: 57.26 on 2 and 39 DF,  p-value: 0.000000000002483

PGLS multiple trees

Function

Note that the function has to be modified for the predictors and data frame that we are using

source(here::here("Scripts/13_multiple_pgls_function_H.R")) # script H is for Polarization as a predictor

Model

MuPGLSPol4 <- Response ~ Pol + VIS

Dataset

PolRefl4 <-
    PolReflRows %>% # phylogeny names in rows
    dplyr::select(Pol, VIS, NIR) %>% 
    dplyr::rename("Response" = NIR)

PolRefl4 <- as.data.frame(PolRefl4)
length(PolRefl4$Pol)

## [1] 42

Apply

runsPol4 <- lapply(trees[trees_subset_min:trees_subset_max],
                   pgls_runH, model = MuPGLSPol4, dataset = PolRefl4)

dfPol4 <- ldply(runsPol4, data.frame)

length(dfPol4[,1])

## [1] 1001

Get results

HPDinterval(as.mcmc(dfPol4))

##                              lower                   upper
## estimateP1 -20.2384529413432687761 -19.0080070342819738016
## ts1         -4.0987982503145756752  -3.3330352524410287174
## ps1          0.0002035776028979086   0.0018902982560042858
## estimateP2   1.0446033217726025111   1.0633688517341663982
## tP2         10.3421952599614090929  10.5308729889663279522
## pP2          0.0000000000005795364   0.0000000000009783285
## attr(,"Probability")
## [1] 0.95005

FRS

Question

Is polarisation correlated with the NIR residuals corrected by phylogeny?

Correlation

Species level:

cor.test(PolRefl$FRS, PolRefl$Pol) # by spp

## 
##  Pearson's product-moment correlation
## 
## data:  PolRefl$FRS and PolRefl$Pol
## t = -4.2365, df = 40, p-value = 0.0001296
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.7359742 -0.3040243
## sample estimates:
##        cor 
## -0.5565308

Linear model

ModL3 <- lm(PolRefl$FRS ~ PolRefl$Pol) # by spp
summary(ModL3)

## 
## Call:
## lm(formula = PolRefl$FRS ~ PolRefl$Pol)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.9032 -2.7627  0.3213  3.1546 11.6020 
## 
## Coefficients:
##             Estimate Std. Error t value     Pr(>|t|)    
## (Intercept)  -5.4607     0.8246  -6.622 0.0000000632 ***
## PolRefl$Pol -20.1846     4.7644  -4.237      0.00013 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.308 on 40 degrees of freedom
## Multiple R-squared:  0.3097, Adjusted R-squared:  0.2925 
## F-statistic: 17.95 on 1 and 40 DF,  p-value: 0.0001296

PGLS

The data frame of comparative data is the same as for visible light

pglsPDFRS <- pgls(FRS ~ Pol,
  data = comp_data, param.CI = 0.95, lambda = "ML"
)

summary(pglsPDFRS)

## 
## Call:
## pgls(formula = FRS ~ Pol, data = comp_data, lambda = "ML", param.CI = 0.95)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.13828 -0.66352  0.07717  0.75763  2.78645 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.000
##    lower bound : 0.000, p = 1    
##    upper bound : 1.000, p = 0.00021248
##    95.0% CI   : (NA, 0.826)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##              Estimate Std. Error t value      Pr(>|t|)    
## (Intercept)  -5.46072    0.82461 -6.6222 0.00000006321 ***
## Pol         -20.18460    4.76444 -4.2365     0.0001296 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.035 on 40 degrees of freedom
## Multiple R-squared: 0.3097,  Adjusted R-squared: 0.2925 
## F-statistic: 17.95 on 1 and 40 DF,  p-value: 0.0001296

PGLS multiple trees

Dataset

PolRefl5 <-
    PolReflRows %>% # phylogeny names in rows
    dplyr::select(Pol, FRS) %>% 
    dplyr::rename("Response" = FRS)

PolRefl5 <- as.data.frame(PolRefl5)

Apply

runsPol5 <- lapply(trees[trees_subset_min:trees_subset_max],
                   pgls_runA, model = MuPGLSPol3, dataset = PolRefl5)

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

dfPol5 <- ldply(runsPol5, data.frame)

length(dfPol5[,1])

## [1] 1000

Get results

HPDinterval(as.mcmc(dfPol5))

##                     lower        upper
## estimateP1 -20.1846083423 -19.44409171
## ts1         -4.2365125038  -3.54355037
## ps1          0.0001295554   0.00102159
## attr(,"Probability")
## [1] 0.95

Absorptivity

VIS

Question

Can polarisation predict visible absorptivity?

Correlation

cor.test(PolAbs$VIS, PolAbs$Pol) # by spp

## 
##  Pearson's product-moment correlation
## 
## data:  PolAbs$VIS and PolAbs$Pol
## t = -0.012444, df = 19, p-value = 0.9902
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.4340068  0.4293611
## sample estimates:
##          cor 
## -0.002854879

Linear model

ModL1Abs <- lm(PolAbs$VIS ~ PolAbs$Pol) # by spp
summary(ModL1Abs)

## 
## Call:
## lm(formula = PolAbs$VIS ~ PolAbs$Pol)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -23.6026  -3.7674  -0.0883   8.9576  14.7073 
## 
## Coefficients:
##             Estimate Std. Error t value            Pr(>|t|)    
## (Intercept)   79.571      3.206  24.816 0.00000000000000061 ***
## PolAbs$Pol    -0.190     15.270  -0.012                0.99    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.05 on 19 degrees of freedom
## Multiple R-squared:  8.15e-06,   Adjusted R-squared:  -0.05262 
## F-statistic: 0.0001549 on 1 and 19 DF,  p-value: 0.9902

PGLS

Create the data frame of comparative data

CompDataAbs <- comparative.data(
  phy = MCCtreeAbs,
  data = PolAbs, # aggregated data without the spp in the row names
  names.col = "phylogeny_name", # contains the column phylogeny name
  vcv = TRUE,
  na.omit = FALSE, warn.dropped = TRUE
)

Establish the model

pglsPDVISAbs <- pgls(VIS ~ Pol,
  data = CompDataAbs, param.CI = 0.95, lambda = "ML"
)

summary(pglsPDVISAbs)

## 
## Call:
## pgls(formula = VIS ~ Pol, data = CompDataAbs, lambda = "ML", 
##     param.CI = 0.95)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.1700 -1.5808  0.0347  2.7116  6.2288 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.973
##    lower bound : 0.000, p = 0.14709
##    upper bound : 1.000, p = 0.55963
##    95.0% CI   : (NA, NA)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##             Estimate Std. Error t value     Pr(>|t|)    
## (Intercept)  80.0709    10.8307  7.3929 0.0000005301 ***
## Pol           8.5603    19.3596  0.4422       0.6634    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.134 on 19 degrees of freedom
## Multiple R-squared: 0.01019, Adjusted R-squared: -0.04191 
## F-statistic: 0.1955 on 1 and 19 DF,  p-value: 0.6634

PGLS multiple trees

Dataset

PolAbs3 <-
    PolAbsRows %>% # phylogeny names in rows
    dplyr::select(Pol, VIS) %>% 
    dplyr::rename("Response" = VIS)

PolAbs3 <- as.data.frame(PolAbs3)

Apply

runsPol6 <- lapply(trees[trees_subset_min:trees_subset_max],
                   pgls_runA, model = MuPGLSPol3, dataset = PolAbs3)

## 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 
## 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

dfPol6 <- ldply(runsPol6, data.frame)

length(dfPol6[,1])

## [1] 960

Get results

HPDinterval(as.mcmc(dfPol6))

##                  lower      upper
## estimateP1 -0.19004396 17.5596701
## ts1        -0.01244528  0.9539945
## ps1         0.35206704  0.9902008
## attr(,"Probability")
## [1] 0.95

Conclusion

There’s no correlation between visible absorptivity and polarisation at the species level.

NIR

Question

Is polarisation correlated with NIR absorptivity?

Correlation

Species level:

cor.test(PolAbs$NIR, PolAbs$Pol) # by spp

## 
##  Pearson's product-moment correlation
## 
## data:  PolAbs$NIR and PolAbs$Pol
## t = 0.33542, df = 19, p-value = 0.741
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.3671231  0.4921109
## sample estimates:
##        cor 
## 0.07672295

Linear model

ModL2Abs <- lm(PolAbs$NIR ~ PolAbs$Pol + PolAbs$VIS) # by spp
summary(ModL2Abs)

## 
## Call:
## lm(formula = PolAbs$NIR ~ PolAbs$Pol + PolAbs$VIS)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13.0173  -1.0129   0.6361   2.4295  11.9533 
## 
## Coefficients:
##             Estimate Std. Error t value      Pr(>|t|)    
## (Intercept) -69.4722    10.2698  -6.765 0.00000244642 ***
## PolAbs$Pol    7.6601     8.4612   0.905         0.377    
## PolAbs$VIS    1.3414     0.1271  10.553 0.00000000388 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.123 on 18 degrees of freedom
## Multiple R-squared:  0.8617, Adjusted R-squared:  0.8463 
## F-statistic: 56.06 on 2 and 18 DF,  p-value: 0.00000001855

PGLS

The data frame of comparative data is the same as for visible light

pglsPDNIRAbs <- pgls(NIR ~ Pol + VIS,
  data = CompDataAbs, param.CI = 0.95, lambda = "ML"
)

summary(pglsPDNIRAbs)

## 
## Call:
## pgls(formula = NIR ~ Pol + VIS, data = CompDataAbs, lambda = "ML", 
##     param.CI = 0.95)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.1264 -0.2433  0.1528  0.5835  2.8708 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.000
##    lower bound : 0.000, p = 1    
##    upper bound : 1.000, p = 0.025945
##    95.0% CI   : (NA, 0.932)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##              Estimate Std. Error t value       Pr(>|t|)    
## (Intercept) -69.47219   10.26976 -6.7647 0.000002446418 ***
## Pol           7.66009    8.46123  0.9053         0.3773    
## VIS           1.34142    0.12712 10.5526 0.000000003879 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.47 on 18 degrees of freedom
## Multiple R-squared: 0.8617,  Adjusted R-squared: 0.8463 
## F-statistic: 56.06 on 2 and 18 DF,  p-value: 0.00000001855

PGLS multiple trees

Dataset

PolAbs4 <-
    PolAbsRows %>% # phylogeny names in rows
    dplyr::select(Pol, VIS, NIR) %>% 
    dplyr::rename("Response" = NIR)

PolAbs4 <- as.data.frame(PolAbs4)

Apply

runsPol7 <- lapply(trees[trees_subset_min:trees_subset_max],
                   pgls_runH, model = MuPGLSPol4, dataset = PolAbs4)

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

dfPol7 <- ldply(runsPol7, data.frame)

length(dfPol7[,1])

## [1] 1000

Get results

HPDinterval(as.mcmc(dfPol7))

##                         lower              upper
## estimateP1  7.660094914901437  7.660099409000678
## ts1         0.905317536902053  0.905318037132289
## ps1         0.377255765964919  0.377256023817372
## estimateP2  1.341424192136824  1.341424236910118
## tP2        10.552631804039306 10.552633288832789
## pP2         0.000000003878582  0.000000003878591
## attr(,"Probability")
## [1] 0.95

FRS

Question

Is polarisation correlated with the NIR residuals corrected by phylogeny?

Correlation

Species level:

cor.test(PolAbs$FRS, PolAbs$Pol) # by spp

## 
##  Pearson's product-moment correlation
## 
## data:  PolAbs$FRS and PolAbs$Pol
## t = 0.93021, df = 19, p-value = 0.3639
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.2450606  0.5874641
## sample estimates:
##      cor 
## 0.208705

Linear model

ModL3Abs <- lm(PolAbs$FRS ~ PolAbs$Pol) # by spp
summary(ModL3Abs)

## 
## Call:
## lm(formula = PolAbs$FRS ~ PolAbs$Pol)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13.0703  -0.8889   0.6560   2.5562  11.8883 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    8.352      1.729   4.829 0.000116 ***
## PolAbs$Pol     7.661      8.236   0.930 0.363932    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.96 on 19 degrees of freedom
## Multiple R-squared:  0.04356,    Adjusted R-squared:  -0.006781 
## F-statistic: 0.8653 on 1 and 19 DF,  p-value: 0.3639

PGLS

The data frame of comparative data is the same as for visible light

pglsPDFRSAbs <- pgls(FRS ~ Pol,
  data = CompDataAbs, param.CI = 0.95, lambda = "ML"
)

summary(pglsPDFRSAbs)

## 
## Call:
## pgls(formula = FRS ~ Pol, data = CompDataAbs, lambda = "ML", 
##     param.CI = 0.95)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.1391 -0.2135  0.1575  0.6139  2.8552 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.000
##    lower bound : 0.000, p = 1    
##    upper bound : 1.000, p = 0.024803
##    95.0% CI   : (NA, 0.931)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##             Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)   8.3518     1.7294  4.8294 0.0001165 ***
## Pol           7.6612     8.2360  0.9302 0.3639322    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.431 on 19 degrees of freedom
## Multiple R-squared: 0.04356, Adjusted R-squared: -0.006781 
## F-statistic: 0.8653 on 1 and 19 DF,  p-value: 0.3639

PGLS multiple trees

Dataset

PolAbs5 <-
    PolAbsRows %>% # phylogeny names in rows
    dplyr::select(Pol, FRS) %>% 
    dplyr::rename("Response" = FRS)

PolAbs5 <- as.data.frame(PolAbs5)

Apply

runsPol8 <- lapply(trees[trees_subset_min:trees_subset_max],
                   pgls_runA, model = MuPGLSPol3, dataset = PolAbs5)

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

dfPol8 <- ldply(runsPol8, data.frame)

length(dfPol8[,1])

## [1] 1000

Get results

HPDinterval(as.mcmc(dfPol8))

##                lower     upper
## estimateP1 7.6612072 7.6612117
## ts1        0.9302082 0.9302088
## ps1        0.3639320 0.3639323
## attr(,"Probability")
## [1] 0.95

Right handed Pol

Reflectivity

VIS

Question

Can the darkness in the underlaying pigment affect the visible reflectivity?

Correlation

cor.test(RHandRefl$VIS, RHandRefl$Pol) # by spp

## 
##  Pearson's product-moment correlation
## 
## data:  RHandRefl$VIS and RHandRefl$Pol
## t = 4.3441, df = 40, p-value = 0.00009303
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.3167468 0.7423607
## sample estimates:
##       cor 
## 0.5661731

Linear model

ModL1RH <- lm(RHandRefl$VIS ~ RHandRefl$Pol) # by spp
summary(ModL1RH)

## 
## Call:
## lm(formula = RHandRefl$VIS ~ RHandRefl$Pol)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -7.168 -3.737 -1.561  1.309 21.874 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     2.5867     2.8693   0.901    0.373    
## RHandRefl$Pol   2.6162     0.6022   4.344 0.000093 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.762 on 40 degrees of freedom
## Multiple R-squared:  0.3206, Adjusted R-squared:  0.3036 
## F-statistic: 18.87 on 1 and 40 DF,  p-value: 0.00009303

PGLS

Create the data frame of comparative data

comp_dataRH <- comparative.data(
  phy = MCCtreeRefl,
  data = RHandRefl, # aggregated data without the spp in the row names
  names.col = "phylogeny_name", # contains the column phylogeny name
  vcv = TRUE,
  na.omit = FALSE, warn.dropped = TRUE
)

Establish the model

pglsPDVISRH <- pgls(VIS ~ Pol,
  data = comp_dataRH, param.CI = 0.95, lambda = "ML"
)

summary(pglsPDVISRH)

## 
## Call:
## pgls(formula = VIS ~ Pol, data = comp_dataRH, lambda = "ML", 
##     param.CI = 0.95)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.9627 -1.4411 -0.1444  0.8486  7.8598 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.948
##    lower bound : 0.000, p = 0.0007746
##    upper bound : 1.000, p = 0.14876
##    95.0% CI   : (0.719, NA)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##             Estimate Std. Error t value    Pr(>|t|)    
## (Intercept)  1.06687    5.25844  0.2029      0.8403    
## Pol          2.57717    0.45862  5.6194 0.000001619 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.112 on 40 degrees of freedom
## Multiple R-squared: 0.4412,  Adjusted R-squared: 0.4272 
## F-statistic: 31.58 on 1 and 40 DF,  p-value: 0.000001619

PGLS multiple trees

Function

Note that the function has to be modified for the predictors and data frame that we are using

source(here::here("Scripts/6_multiple_pgls_function_A.R")) # script A is for Polarization as a response

Model

MuPGLSPol3 <- Response ~ Pol

Dataset

RHandRefl3 <-
    RHandReflRows %>% # phylogeny names in rows
    dplyr::select(Pol, VIS) %>% 
    dplyr::rename("Response" = VIS)

RHandRefl3 <- as.data.frame(RHandRefl3)

Apply

runsRHPol3 <- lapply(trees[trees_subset_min:trees_subset_max],
                   pgls_runA, model = MuPGLSPol3, dataset = RHandRefl3)

## 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

dfRHPol3 <- ldply(runsRHPol3, data.frame)

length(dfRHPol3[,1])

## [1] 989

Get results

HPDinterval(as.mcmc(dfRHPol3))

##                       lower          upper
## estimateP1 2.38593639691057 2.815684962779
## ts1        5.18844277798936 6.317166520323
## ps1        0.00000003188983 0.000005007898
## attr(,"Probability")
## [1] 0.950455

Conclusion

PENDING.

NIR

Question

Can the darkness in the underlaying pigment affect the NIR reflectivity?

Correlation

Species level:

cor.test(RHandRefl$NIR, RHandRefl$Pol) # by spp

## 
##  Pearson's product-moment correlation
## 
## data:  RHandRefl$NIR and RHandRefl$Pol
## t = 5.5789, df = 40, p-value = 0.000001846
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4475682 0.8038311
## sample estimates:
##       cor 
## 0.6615138

Linear model

ModL2RH <- lm(RHandRefl$NIR ~ RHandRefl$Pol + RHandRefl$VIS) # by spp
summary(ModL2RH)

## 
## Call:
## lm(formula = RHandRefl$NIR ~ RHandRefl$Pol + RHandRefl$VIS)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.1608 -3.3970  0.2871  3.4140  8.9805 
## 
## Coefficients:
##               Estimate Std. Error t value       Pr(>|t|)    
## (Intercept)    19.2340     2.3793   8.084 0.000000000732 ***
## RHandRefl$Pol   1.7444     0.5998   2.908        0.00597 ** 
## RHandRefl$VIS   0.7617     0.1298   5.868 0.000000789023 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.73 on 39 degrees of freedom
## Multiple R-squared:  0.7013, Adjusted R-squared:  0.686 
## F-statistic: 45.79 on 2 and 39 DF,  p-value: 0.00000000005837

PGLS

The data frame of comparative data is the same as for visible light

pglsPDNIRRH <- pgls(NIR ~ Pol + VIS,
  data = comp_dataRH, param.CI = 0.95, lambda = "ML"
)

summary(pglsPDNIRRH)

## 
## Call:
## pgls(formula = NIR ~ Pol + VIS, data = comp_dataRH, lambda = "ML", 
##     param.CI = 0.95)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.4817 -1.0627 -0.1391  0.8433  2.9004 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.547
##    lower bound : 0.000, p = 0.32623
##    upper bound : 1.000, p = 0.0027707
##    95.0% CI   : (NA, 0.933)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##             Estimate Std. Error t value     Pr(>|t|)    
## (Intercept) 20.98072    3.29190  6.3734 0.0000001570 ***
## Pol          1.23339    0.60744  2.0305      0.04917 *  
## VIS          0.80997    0.13593  5.9589 0.0000005907 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.356 on 39 degrees of freedom
## Multiple R-squared: 0.6922,  Adjusted R-squared: 0.6764 
## F-statistic: 43.85 on 2 and 39 DF,  p-value: 0.000000000105

PGLS multiple trees

Function

Note that the function has to be modified for the predictors and data frame that we are using

source(here::here("Scripts/13_multiple_pgls_function_H.R")) # script A is for Polarization as a response

Model

MuPGLSPol4 <- Response ~ Pol + VIS

Dataset

RHandRefl4 <-
    RHandReflRows %>% # phylogeny names in rows
    dplyr::select(Pol, VIS, NIR) %>% 
    dplyr::rename("Response" = NIR)

RHandRefl4 <- as.data.frame(RHandRefl4)

Apply

runsRHPol4 <- lapply(trees[trees_subset_min:trees_subset_max],
                   pgls_runH, model = MuPGLSPol4, dataset = RHandRefl4)

dfRHPol4 <- ldply(runsRHPol4, data.frame)

length(dfRHPol4[,1])

## [1] 1001

Get results

HPDinterval(as.mcmc(dfRHPol4))

##                     lower          upper
## estimateP1 1.094772337238 1.399996461732
## ts1        1.787603918459 2.288182630668
## ps1        0.024712922809 0.075434993034
## estimateP2 0.772205679303 0.835260771457
## tP2        5.603581786223 6.172203061114
## pP2        0.000000199416 0.000001658501
## attr(,"Probability")
## [1] 0.95005

FRS

Question

Can the darkness in the underlaying pigment affect the NIR reflectivity?

Correlation

Species level:

cor.test(RHandRefl$FRS, RHandRefl$Pol) # by spp

## 
##  Pearson's product-moment correlation
## 
## data:  RHandRefl$FRS and RHandRefl$Pol
## t = 1.8638, df = 40, p-value = 0.0697
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.02325426  0.54019604
## sample estimates:
##      cor 
## 0.282675

Linear model

ModL3RH <- lm(RHandRefl$FRS ~ RHandRefl$Pol) # by spp
summary(ModL3RH)

## 
## Call:
## lm(formula = RHandRefl$FRS ~ RHandRefl$Pol)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.0141 -4.1551  0.1555  4.2004  9.2638 
## 
## Coefficients:
##               Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)   -11.9167     2.4766  -4.812 0.0000215 ***
## RHandRefl$Pol   0.9688     0.5198   1.864    0.0697 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.973 on 40 degrees of freedom
## Multiple R-squared:  0.07991,    Adjusted R-squared:  0.0569 
## F-statistic: 3.474 on 1 and 40 DF,  p-value: 0.0697

PGLS

The data frame of comparative data is the same as for visible light

pglsPDFRSRH <- pgls(FRS ~ Pol,
  data = comp_dataRH, param.CI = 0.95, lambda = "ML"
)

summary(pglsPDFRSRH)

## 
## Call:
## pgls(formula = FRS ~ Pol, data = comp_dataRH, lambda = "ML", 
##     param.CI = 0.95)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.8756 -1.0530 -0.2348  1.1716  2.9551 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.693
##    lower bound : 0.000, p = 0.081072
##    upper bound : 1.000, p = 0.0086855
##    95.0% CI   : (NA, 0.967)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -9.11035    3.75985 -2.4231  0.02001 *
## Pol          0.44772    0.48340  0.9262  0.35990  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.525 on 40 degrees of freedom
## Multiple R-squared: 0.021,   Adjusted R-squared: -0.00348 
## F-statistic: 0.8578 on 1 and 40 DF,  p-value: 0.3599

PGLS multiple trees

Dataset

RHandRefl5 <-
    RHandReflRows %>% # phylogeny names in rows
    dplyr::select(Pol, FRS) %>% 
    dplyr::rename("Response" = FRS)

RHandRefl5 <- as.data.frame(RHandRefl5)

Apply

runsRHPol5 <- lapply(trees[trees_subset_min:trees_subset_max],
                   pgls_runA, model = MuPGLSPol3, dataset = RHandRefl5)

dfRHPol5 <- ldply(runsRHPol5, data.frame)

length(dfRHPol5[,1])

## [1] 1001

Get results

HPDinterval(as.mcmc(dfRHPol5))

##                lower     upper
## estimateP1 0.3552530 0.5722316
## ts1        0.7702818 1.1862122
## ps1        0.2425368 0.4456585
## attr(,"Probability")
## [1] 0.95005

Absorptivity

VIS

Question

Can polarisation predict visible absorptivity?

Correlation

cor.test(RHandAbs$VIS, RHandAbs$Pol) # by spp

## 
##  Pearson's product-moment correlation
## 
## data:  RHandAbs$VIS and RHandAbs$Pol
## t = -0.92958, df = 19, p-value = 0.3643
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.5873711  0.2451940
## sample estimates:
##        cor 
## -0.2085692

Linear model

ModL1AbsRH <- lm(RHandAbs$VIS ~ RHandAbs$Pol) # by spp
summary(ModL1AbsRH)

## 
## Call:
## lm(formula = RHandAbs$VIS ~ RHandAbs$Pol)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -24.1117  -5.4171   0.0063   8.3618  12.1358 
## 
## Coefficients:
##              Estimate Std. Error t value      Pr(>|t|)    
## (Intercept)    86.056      7.390   11.64 0.00000000043 ***
## RHandAbs$Pol   -1.594      1.715   -0.93         0.364    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.81 on 19 degrees of freedom
## Multiple R-squared:  0.0435, Adjusted R-squared:  -0.006841 
## F-statistic: 0.8641 on 1 and 19 DF,  p-value: 0.3643

PGLS

Create the data frame of comparative data

CompDataAbsRH <- comparative.data(
  phy = MCCtreeAbs,
  data = RHandAbs, # aggregated data without the spp in the row names
  names.col = "phylogeny_name", # contains the column phylogeny name
  vcv = TRUE,
  na.omit = FALSE, warn.dropped = TRUE
)

Establish the model

pglsPDVISAbsRH <- pgls(VIS ~ Pol,
  data = CompDataAbsRH, param.CI = 0.95, lambda = "ML"
)

summary(pglsPDVISAbsRH)

## 
## Call:
## pgls(formula = VIS ~ Pol, data = CompDataAbsRH, lambda = "ML", 
##     param.CI = 0.95)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.0943 -1.3240  0.2026  3.3121  5.6247 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.960
##    lower bound : 0.000, p = 0.30126
##    upper bound : 1.000, p = 0.44529
##    95.0% CI   : (NA, NA)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##             Estimate Std. Error t value     Pr(>|t|)    
## (Intercept) 83.23122   11.24337  7.4027 0.0000005202 ***
## Pol         -0.44213    1.38708 -0.3188       0.7534    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.061 on 19 degrees of freedom
## Multiple R-squared: 0.005319,    Adjusted R-squared: -0.04703 
## F-statistic: 0.1016 on 1 and 19 DF,  p-value: 0.7534

PGLS multiple trees

Dataset

RHandAbs3 <-
    RHandAbsRows %>% # phylogeny names in rows
    dplyr::select(Pol, VIS) %>% 
    dplyr::rename("Response" = VIS)

RHandAbs3 <- as.data.frame(RHandAbs3)

Apply

runsRHPol6 <- lapply(trees[trees_subset_min:trees_subset_max],
                   pgls_runA, model = MuPGLSPol3, dataset = RHandAbs3)

## 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

dfRHPol6 <- ldply(runsRHPol6, data.frame)

length(dfRHPol6[,1])

## [1] 965

Get results

HPDinterval(as.mcmc(dfRHPol6))

##                 lower      upper
## estimateP1 -1.5940794 0.06334685
## ts1        -0.9295755 0.05191069
## ps1         0.3642516 0.95914174
## attr(,"Probability")
## [1] 0.9502591

Conclusion

There’s no correlation between visible absorptivity and polarisation at the species level.

NIR

Question

Is polarisation correlated with NIR absorptivity?

Correlation

Species level:

cor.test(RHandAbs$NIR, RHandAbs$Pol) # by spp

## 
##  Pearson's product-moment correlation
## 
## data:  RHandAbs$NIR and RHandAbs$Pol
## t = -0.69617, df = 19, p-value = 0.4948
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.5518296  0.2939895
## sample estimates:
##        cor 
## -0.1577128

Linear model

ModL2AbsRH <- lm(RHandAbs$NIR ~ RHandAbs$Pol + RHandAbs$VIS) # by spp
summary(ModL2AbsRH)

## 
## Call:
## lm(formula = RHandAbs$NIR ~ RHandAbs$Pol + RHandAbs$VIS)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.4573  -1.9427   0.5778   2.7432  14.0929 
## 
## Coefficients:
##              Estimate Std. Error t value      Pr(>|t|)    
## (Intercept)  -70.9359    12.1572  -5.835 0.00001583165 ***
## RHandAbs$Pol   0.4077     1.0112   0.403         0.692    
## RHandAbs$VIS   1.3522     0.1323  10.220 0.00000000638 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.233 on 18 degrees of freedom
## Multiple R-squared:  0.8567, Adjusted R-squared:  0.8407 
## F-statistic: 53.79 on 2 and 18 DF,  p-value: 0.00000002554

PGLS

The data frame of comparative data is the same as for visible light

pglsPDNIRAbsRH <- pgls(NIR ~ Pol + VIS,
  data = CompDataAbsRH, param.CI = 0.95, lambda = "ML"
)

summary(pglsPDNIRAbsRH)

## 
## Call:
## pgls(formula = NIR ~ Pol + VIS, data = CompDataAbsRH, lambda = "ML", 
##     param.CI = 0.95)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.9919 -0.4666  0.1388  0.6588  3.3847 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.000
##    lower bound : 0.000, p = 1    
##    upper bound : 1.000, p = 0.060342
##    95.0% CI   : (NA, NA)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##              Estimate Std. Error t value       Pr(>|t|)    
## (Intercept) -70.93587   12.15716 -5.8349 0.000015831622 ***
## Pol           0.40768    1.01120  0.4032         0.6916    
## VIS           1.35222    0.13231 10.2204 0.000000006378 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.497 on 18 degrees of freedom
## Multiple R-squared: 0.8567,  Adjusted R-squared: 0.8407 
## F-statistic: 53.79 on 2 and 18 DF,  p-value: 0.00000002554

PGLS multiple trees

Dataset

RHandAbs4 <-
    RHandAbsRows %>% # phylogeny names in rows
    dplyr::select(Pol, VIS, NIR) %>% 
    dplyr::rename("Response" = NIR)

RHandAbs4 <- as.data.frame(RHandAbs4)

Apply

runsRHPol7 <- lapply(trees[trees_subset_min:trees_subset_max],
                   pgls_runH, model = MuPGLSPol4, dataset = RHandAbs4)

## 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

dfRHPol7 <- ldply(runsRHPol7, data.frame)

length(dfRHPol7[,1])

## [1] 994

Get results

HPDinterval(as.mcmc(dfRHPol7))

##                         lower              upper
## estimateP1  0.407675987360903  0.407676342898483
## ts1         0.403159649045050  0.403159984140517
## ps1         0.691579037712050  0.691579279810262
## estimateP2  1.352220773403229  1.352220816847392
## tP2        10.220434063137793 10.220435396596766
## pP2         0.000000006378462  0.000000006378475
## attr(,"Probability")
## [1] 0.9496982

FRS

Question

Is polarisation correlated with the NIR residuals corrected by phylogeny?

Correlation

Species level:

cor.test(RHandAbs$FRS, RHandAbs$Pol) # by spp

## 
##  Pearson's product-moment correlation
## 
## data:  RHandAbs$FRS and RHandAbs$Pol
## t = 0.41542, df = 19, p-value = 0.6825
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.3511961  0.5058437
## sample estimates:
##        cor 
## 0.09487424

Linear model

ModL3AbsRH <- lm(RHandAbs$FRS ~ RHandAbs$Pol) # by spp
summary(ModL3AbsRH)

## 
## Call:
## lm(formula = RHandAbs$FRS ~ RHandAbs$Pol)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.4229  -1.9535   0.4778   2.7771  14.1443 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept)    7.7789     4.1484   1.875   0.0762 .
## RHandAbs$Pol   0.3999     0.9626   0.415   0.6825  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.067 on 19 degrees of freedom
## Multiple R-squared:  0.009001,   Adjusted R-squared:  -0.04316 
## F-statistic: 0.1726 on 1 and 19 DF,  p-value: 0.6825

PGLS

The data frame of comparative data is the same as for visible light

pglsPDFRSAbsRH <- pgls(FRS ~ Pol,
  data = CompDataAbsRH, param.CI = 0.95, lambda = "ML"
)

summary(pglsPDFRSAbsRH)

## 
## Call:
## pgls(formula = FRS ~ Pol, data = CompDataAbsRH, lambda = "ML", 
##     param.CI = 0.95)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.9836 -0.4692  0.1147  0.6670  3.3970 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.000
##    lower bound : 0.000, p = 1    
##    upper bound : 1.000, p = 0.059965
##    95.0% CI   : (NA, NA)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  7.77892    4.14840  1.8752  0.07623 .
## Pol          0.39989    0.96262  0.4154  0.68249  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.457 on 19 degrees of freedom
## Multiple R-squared: 0.009001,    Adjusted R-squared: -0.04316 
## F-statistic: 0.1726 on 1 and 19 DF,  p-value: 0.6825

PGLS multiple trees

Dataset

RHandAbs5 <-
    RHandAbsRows %>% # phylogeny names in rows
    dplyr::select(Pol, FRS) %>% 
    dplyr::rename("Response" = FRS)

RHandAbs5 <- as.data.frame(RHandAbs5)

Apply

runsRHPol8 <- lapply(trees[trees_subset_min:trees_subset_max],
                   pgls_runA, model = MuPGLSPol3, dataset = RHandAbs5)

## 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

dfRHPol8 <- ldply(runsRHPol8, data.frame)

length(dfRHPol8[,1])

## [1] 994

Get results

HPDinterval(as.mcmc(dfRHPol8))

##                lower     upper
## estimateP1 0.3998936 0.3998940
## ts1        0.4154209 0.4154212
## ps1        0.6824874 0.6824877
## attr(,"Probability")
## [1] 0.9496982

Polarisation analysis without outliers

Laura Ospina-Rozo

Setting up

Libraries

Data sets

Optical properties

Pol Degree

Reflectivity

VIS

NIR

FRS

Absorptivity

VIS

NIR

FRS

Right handed Pol

Reflectivity

VIS

NIR

FRS

Absorptivity

VIS

NIR

FRS