Polarisation analysis
Laura Ospina-Rozo
In this script we study the correlations between the degree of polarisation and the ecological variables (PCs) and other optical properties, i.e. reflectivity and absorptivity. In addition, we used the raw RGB value sin the RCP filter as a coarse indicator of the presence of pigments (melanin). Here it is called “Right handed Pol”. We tested for correlations between the Right handed Polarisation and reflectivity and absorptivity.
The part called “Mechanisms” contains only exploratory analysis. The formal equivalents can be found in “Quality Control (Polarisation)”
Setting up
Libraries
Libraries sourced from an additional script.
Data sets
Import polarisation data.
# elytra
PolBySpp <- read.csv(here::here("Data/FromCode/PolarzElytraBySpp.csv"))[-1]
PolByInd <- read.csv(here::here("Data/FromCode/PolarzElytraByInd.csv"))[-1]
# 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)
# 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)
# 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)
# 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 transmissivity
TraBySpp <- read.csv(here::here("Data/FromCode/ConsolidatedTransmbySpp.csv"
)) %>%
dplyr::rename(phylogeny_name = X) %>%
dplyr::arrange(phylogeny_name)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)
# 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
# Combine with Transmissivity
RHandTra <- merge(TraBySpp, RHand)
RHandTra2 <- as.data.frame(RHandTra) # create new data frame
rownames(RHandTra2) <- RHandTra2[, 1] # make species the row names
RHandTraRows <- RHandTra2[, 2:length(RHandTra2)] # 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))
# Convert to "row names" (required for following steps)
row.names(PolAbs.species.MCC) <- PolAbs.species.MCC[, 1]
# Make sure the names in data set and tree match
PolAbs.temp.MCC <- name.check(MCCtree, PolAbs.species.MCC)
# This step is necessary because the tips are different.
MCCtreeAbs <- drop.tip(MCCtree, PolAbs.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] TRUEBody part
It is possible that different selective pressures apply to the elytra than to other parts of the beetle body. There are various species that seem to have stronger structural coloration in the pronotum. They may have lost the polarization in their elytra.
Are there differences in the polarisation between the elytra and pronotum of the same individual?
Statistics
Conduct a paired t-test to compare the mean of the difference to 0. A negative difference would mean that in general the pronotum is more polarized than the elytra.
##
## One Sample t-test
##
## data: BeetlePol$Difference
## t = -1.2468, df = 173, p-value = 0.2142
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.019894362 0.004490845
## sample estimates:
## mean of x
## -0.007701759The mean difference in polarisation between the two body parts is not significantly different than 0.
Plots
hist(BeetlePol$Difference,
breaks = 20,
xlab = "Difference in polarization (Elytra-Pronot)",
ylab = "Frequency"
)
These beetles with higher polarisation in their elytra are:
## [1] "atki01" "atki02" "brun01" "brun04" "brun05" "ecry01" "fchi01" "fchi02" "hirs03" "lats01" "lats03" "pvul02" "pvul03" "pvul04" "rayn01" "rayn03" "smgr01"These beetles with higher polarisation in their pronotum are:
## [1] "aurs01" "aurs02" "clor09" "meli02" "meli03" "neus01" "oliv03" "opal01" "pvul01" "rubi01" "rubi02" "vrdi01" "vrdi03"Climate
Question
Is polarization correlated to the extracted climate variables?
Statistics
Pagel’s lambda
pagelPDPol <- PolPCsRows$Pol # Define which trait we want to test
names(pagelPDPol) <- rownames(PolPCsRows) # Row names = tree tips
phylosig(MCCtree, pagelPDPol, method = "lambda", test = TRUE, nsim = 999)##
## Phylogenetic signal lambda : 0.955198
## logL(lambda) : 40.0202
## LR(lambda=0) : 28.1436
## P-value (based on LR test) : 0.000000112637Phylogenetic signal is high.
PGLS in the MCC
Create the data frame of comparative data
# extract sizes
SizesDf <- PolRefl %>%
dplyr::select("phylogeny_name", "size")
# merge size into ecological variables dataset
PolPCsWithSizes <- merge(SizesDf, PolPCs)
comp_data2 <- comparative.data(
phy = MCCtree,
data = PolPCsWithSizes, # 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
pglsPDPol <- pgls(Pol ~ PC1 + PC2 + size + PC1 * size + PC2 * size,
data = comp_data2, param.CI = 0.95, lambda = "ML"
)
summary(pglsPDPol)##
## Call:
## pgls(formula = Pol ~ PC1 + PC2 + size + PC1 * size + PC2 * size,
## data = comp_data2, lambda = "ML", param.CI = 0.95)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.073617 -0.016438 -0.005165 0.018150 0.076266
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.908
## lower bound : 0.000, p = 0.0000078697
## upper bound : 1.000, p = 0.0037317
## 95.0% CI : (0.710, 0.983)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0245801 0.1416766 0.1735 0.8631
## PC1 -0.0037788 0.0508632 -0.0743 0.9411
## PC2 -0.1073466 0.0641932 -1.6722 0.1021
## size 0.0107276 0.0493041 0.2176 0.8288
## PC1:size -0.0078721 0.0218726 -0.3599 0.7208
## PC2:size 0.0358297 0.0257598 1.3909 0.1718
##
## Residual standard error: 0.03551 on 41 degrees of freedom
## Multiple R-squared: 0.2871, Adjusted R-squared: 0.2001
## F-statistic: 3.302 on 5 and 41 DF, p-value: 0.01348
Conclusion
Polarization does not directly correlate with these ecological variables
PGLS in 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/7_multiple_pgls_function_B.R")) # script B is for Polarization as a responseModel
Dataset
# Add the phylogeny name again because in this function it is needed:
PDPol2 <-
PolPCsWithSizes %>%
dplyr::select(phylogeny_name, Pol, size, PC1, PC2) %>%
dplyr::rename("Polarization" = Pol)
PDPol2 <- as.data.frame(PDPol2)
rownames(PDPol2) <- PDPol2[, 1]Apply
runsPol2 <- lapply(trees[trees_subset_min:trees_subset_max],
pgls_runB, model = MuPGLSPol2, dataset = PDPol2)## 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_LNSRCHGet results
## lower upper
## estimateP1 -0.019961932 0.012065708
## ts1 -0.398236915 0.218698747
## ps1 0.731524167 0.999803405
## estimateP2 -0.150983585 -0.075873842
## tP2 -2.341620113 -1.284685625
## pP2 0.017177465 0.173785078
## estimateP3 -0.004033164 0.037508114
## tP3 -0.081294101 0.739204491
## pP3 0.472676043 0.999673203
## estimateP4 -0.014311118 -0.001079027
## tP4 -0.659690073 -0.050508450
## pP4 0.513141587 0.939636609
## estimateP5 0.024501773 0.053619651
## tP5 0.974480549 2.027298118
## pP5 0.035180472 0.282547057
## attr(,"Probability")
## [1] 0.950455Optical 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
##
## Pearson's product-moment correlation
##
## data: PolRefl$VIS and PolRefl$Pol
## t = 1.5434, df = 45, p-value = 0.1297
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.06728212 0.48044869
## sample estimates:
## cor
## 0.2242168
Linear model
##
## Call:
## lm(formula = PolRefl$VIS ~ PolRefl$Pol)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.4480 -4.9439 0.6226 3.3241 18.8607
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13.288 1.131 11.744 0.00000000000000268 ***
## PolRefl$Pol 10.617 6.879 1.543 0.13
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.57 on 45 degrees of freedom
## Multiple R-squared: 0.05027, Adjusted R-squared: 0.02917
## F-statistic: 2.382 on 1 and 45 DF, p-value: 0.1297
PGLS
Create the data frame of comparative data
comp_data <- comparative.data(
phy = MCCtree,
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
##
## Call:
## pgls(formula = VIS ~ Pol, data = comp_data, lambda = "ML", param.CI = 0.95)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.463 -1.473 0.033 1.128 6.399
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.821
## lower bound : 0.000, p = 0.59187
## upper bound : 1.000, p = 0.15519
## 95.0% CI : (NA, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.4628 5.1935 2.2072 0.03245 *
## Pol 8.5383 8.8468 0.9651 0.33964
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.204 on 45 degrees of freedom
## Multiple R-squared: 0.02028, Adjusted R-squared: -0.001492
## F-statistic: 0.9315 on 1 and 45 DF, p-value: 0.3396
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 responseModel
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
## 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## [1] 980Get results
## lower upper
## estimateP1 1.5389948 12.2133812
## ts1 0.2060000 1.5457472
## ps1 0.1200562 0.8031361
## attr(,"Probability")
## [1] 0.95Conclusion
There’s no correlation between visible reflectivity and polarisation at the species level.
NIR
Question
Is polarisation correlated with NIR reflectivity?
Correlation
Species level:
##
## Pearson's product-moment correlation
##
## data: PolRefl$NIR and PolRefl$Pol
## t = -1.7419, df = 45, p-value = 0.08836
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.50225033 0.03862105
## sample estimates:
## cor
## -0.2513332
Linear model
##
## Call:
## lm(formula = PolRefl$NIR ~ PolRefl$Pol + PolRefl$VIS)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.7760 -3.7231 0.3623 1.8078 26.2201
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 26.4224 2.2774 11.602 0.00000000000000563 ***
## PolRefl$Pol -29.5900 7.0475 -4.199 0.000129 ***
## PolRefl$VIS 1.1102 0.1488 7.460 0.00000000243789919 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.56 on 44 degrees of freedom
## Multiple R-squared: 0.5863, Adjusted R-squared: 0.5675
## F-statistic: 31.18 on 2 and 44 DF, p-value: 0.000000003685
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
## -3.3801 -1.0779 -0.2568 0.5219 7.7747
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.602
## lower bound : 0.000, p = 0.34197
## upper bound : 1.000, p = 0.021873
## 95.0% CI : (NA, 0.985)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 30.4608 4.3375 7.0227 0.000000010590 ***
## Pol -28.4643 8.2204 -3.4627 0.001203 **
## VIS 1.1146 0.1471 7.5774 0.000000001644 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.824 on 44 degrees of freedom
## Multiple R-squared: 0.5847, Adjusted R-squared: 0.5658
## F-statistic: 30.97 on 2 and 44 DF, p-value: 0.000000004027
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 predictorModel
Dataset
PolRefl4 <-
PolReflRows %>% # phylogeny names in rows
dplyr::select(Pol, VIS, NIR) %>%
dplyr::rename("Response" = NIR)
PolRefl4 <- as.data.frame(PolRefl4)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] 1001Get results
## lower upper
## estimateP1 -30.5787807936124665 -26.984834060583118
## ts1 -4.1986326287349671 -3.208252902704267
## ps1 0.0001285882928512 0.002492216158802
## estimateP2 1.0936560546755665 1.133847417482905
## tP2 7.3334922240313078 7.843723842514416
## pP2 0.0000000006540746 0.000000003672992
## attr(,"Probability")
## [1] 0.95005FRS
Question
Is polarisation correlated with the NIR residuals corrected by phylogeny?
Correlation
Species level:
##
## Pearson's product-moment correlation
##
## data: PolRefl$FRS and PolRefl$Pol
## t = -4.2698, df = 45, p-value = 0.00009974
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.7140250 -0.2953307
## sample estimates:
## cor
## -0.5369585
Linear model
##
## Call:
## lm(formula = PolRefl$FRS ~ PolRefl$Pol)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.7187 -3.4801 0.3577 2.0168 26.3987
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.270 1.119 -2.924 0.00539 **
## PolRefl$Pol -29.038 6.801 -4.270 0.0000997 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.495 on 45 degrees of freedom
## Multiple R-squared: 0.2883, Adjusted R-squared: 0.2725
## F-statistic: 18.23 on 1 and 45 DF, p-value: 0.00009974
PGLS
The data frame of comparative data is the same as for visible light
##
## Call:
## pgls(formula = FRS ~ Pol, data = comp_data, lambda = "ML", param.CI = 0.95)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.3650 -1.0477 -0.2635 0.5654 7.8514
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.603
## lower bound : 0.000, p = 0.34922
## upper bound : 1.000, p = 0.023464
## 95.0% CI : (NA, 0.986)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.74353 3.95005 0.1882 0.851541
## Pol -27.88818 8.01185 -3.4809 0.001123 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.809 on 45 degrees of freedom
## Multiple R-squared: 0.2121, Adjusted R-squared: 0.1946
## F-statistic: 12.12 on 1 and 45 DF, p-value: 0.001123
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)
dfPol5 <- ldply(runsPol5, data.frame)
length(dfPol5[,1])## [1] 1001Get results
## lower upper
## estimateP1 -29.78511659011 -26.361994828
## ts1 -4.26978422766 -3.205786973
## ps1 0.00009973628 0.002479313
## attr(,"Probability")
## [1] 0.95005Absorptivity
VIS
Question
Can polarisation predict visible absorptivity?
Correlation
##
## Pearson's product-moment correlation
##
## data: PolAbs$VIS and PolAbs$Pol
## t = 0.041288, df = 24, p-value = 0.9674
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.3801654 0.3944917
## sample estimates:
## cor
## 0.008427517
Linear model
##
## Call:
## lm(formula = PolAbs$VIS ~ PolAbs$Pol)
##
## Residuals:
## Min 1Q Median 3Q Max
## -23.4765 -4.4632 0.5597 7.5187 14.9412
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 79.3015 2.3443 33.828 <0.0000000000000002 ***
## PolAbs$Pol 0.5093 12.3345 0.041 0.967
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.99 on 24 degrees of freedom
## Multiple R-squared: 7.102e-05, Adjusted R-squared: -0.04159
## F-statistic: 0.001705 on 1 and 24 DF, p-value: 0.9674
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.0851 -1.5359 0.0797 2.6208 6.0661
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.969
## lower bound : 0.000, p = 0.29867
## upper bound : 1.000, p = 0.4437
## 95.0% CI : (NA, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 80.7282 9.2640 8.7142 0.000000006724 ***
## Pol 5.9897 15.8311 0.3784 0.7085
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.695 on 24 degrees of freedom
## Multiple R-squared: 0.005929, Adjusted R-squared: -0.03549
## F-statistic: 0.1432 on 1 and 24 DF, p-value: 0.7085
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## [1] 963Get results
## lower upper
## estimateP1 0.45439180 11.9958349
## ts1 0.04128687 0.8167620
## ps1 0.43312728 0.9674088
## attr(,"Probability")
## [1] 0.9501558Conclusion
There’s no correlation between visible absorptivity and polarisation at the species level.
NIR
Question
Is polarisation correlated with NIR absorptivity?
Correlation
Species level:
##
## Pearson's product-moment correlation
##
## data: PolAbs$NIR and PolAbs$Pol
## t = 1.0748, df = 24, p-value = 0.2932
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1887276 0.5555297
## sample estimates:
## cor
## 0.2142898
Linear model
##
## Call:
## lm(formula = PolAbs$NIR ~ PolAbs$Pol + PolAbs$VIS)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.1275 -4.5971 0.2726 4.3498 11.4918
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -73.006 11.617 -6.285 0.00000206065 ***
## PolAbs$Pol 19.381 8.761 2.212 0.0372 *
## PolAbs$VIS 1.345 0.145 9.274 0.00000000311 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.095 on 23 degrees of freedom
## Multiple R-squared: 0.7987, Adjusted R-squared: 0.7812
## F-statistic: 45.62 on 2 and 23 DF, p-value: 0.000000009877
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.6008 -0.9192 0.5430 1.3462 3.6475
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.698
## lower bound : 0.000, p = 0.047281
## upper bound : 1.000, p = 0.09444
## 95.0% CI : (0.024, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -79.24713 11.51810 -6.8802 0.0000005146022 ***
## Pol 18.59591 9.55406 1.9464 0.06392 .
## VIS 1.34625 0.13193 10.2041 0.0000000005214 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.878 on 23 degrees of freedom
## Multiple R-squared: 0.8256, Adjusted R-squared: 0.8104
## F-statistic: 54.42 on 2 and 23 DF, p-value: 0.000000001902
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)
dfPol7 <- ldply(runsPol7, data.frame)
length(dfPol7[,1])## [1] 1001Get results
## lower upper
## estimateP1 17.1499526859974480 21.154935391163267
## ts1 1.8199538150269474 2.345445447530623
## ps1 0.0263609336822412 0.079977595961537
## estimateP2 1.3184074324238098 1.360615900689254
## tP2 9.6488193841426106 10.803000242408297
## pP2 0.0000000001229536 0.000000001323613
## attr(,"Probability")
## [1] 0.95005FRS
Question
Is polarisation correlated with the NIR residuals corrected by phylogeny?
Correlation
Species level:
##
## Pearson's product-moment correlation
##
## data: PolAbs$FRS and PolAbs$Pol
## t = 2.2597, df = 24, p-value = 0.03319
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.03759808 0.69366079
## sample estimates:
## cor
## 0.41885
Linear model
##
## Call:
## lm(formula = PolAbs$FRS ~ PolAbs$Pol)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.1317 -4.5898 0.3022 4.3367 11.4591
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.061 1.630 3.105 0.00483 **
## PolAbs$Pol 19.379 8.576 2.260 0.03319 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.946 on 24 degrees of freedom
## Multiple R-squared: 0.1754, Adjusted R-squared: 0.1411
## F-statistic: 5.106 on 1 and 24 DF, p-value: 0.03319
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.6041 -0.9185 0.5429 1.3454 3.6435
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.697
## lower bound : 0.000, p = 0.047271
## upper bound : 1.000, p = 0.087831
## 95.0% CI : (0.024, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.0393 4.2966 -0.2419 0.81092
## Pol 18.5950 9.3502 1.9887 0.05825 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.839 on 24 degrees of freedom
## Multiple R-squared: 0.1415, Adjusted R-squared: 0.1057
## F-statistic: 3.955 on 1 and 24 DF, p-value: 0.05825
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)
dfPol8 <- ldply(runsPol8, data.frame)
length(dfPol8[,1])## [1] 1001Get results
## lower upper
## estimateP1 17.15984924 21.18197601
## ts1 1.85696827 2.39465696
## ps1 0.02383071 0.07389166
## attr(,"Probability")
## [1] 0.95005Right handed Pol
Reflectivity
VIS
Question
Can the darkness in the underlaying pigment affect the visible reflectivity?
Correlation
##
## Pearson's product-moment correlation
##
## data: RHandRefl$VIS and RHandRefl$Pol
## t = 4.5919, df = 45, p-value = 0.00003531
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.3314552 0.7331079
## sample estimates:
## cor
## 0.5648573
Linear model
##
## Call:
## lm(formula = RHandRefl$VIS ~ RHandRefl$Pol)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.183 -3.552 -1.193 1.611 21.904
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.0506 2.5633 1.190 0.24
## RHandRefl$Pol 2.4853 0.5412 4.592 0.0000353 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.563 on 45 degrees of freedom
## Multiple R-squared: 0.3191, Adjusted R-squared: 0.3039
## F-statistic: 21.09 on 1 and 45 DF, p-value: 0.00003531
PGLS
Create the data frame of comparative data
comp_dataRH <- comparative.data(
phy = MCCtree,
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.7136 -1.4749 -0.2582 0.7616 7.7002
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.944
## lower bound : 0.000, p = 0.0018141
## upper bound : 1.000, p = 0.10278
## 95.0% CI : (0.699, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.71490 5.10664 0.5316 0.5976
## Pol 2.48798 0.42974 5.7896 0.0000006404 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.007 on 45 degrees of freedom
## Multiple R-squared: 0.4269, Adjusted R-squared: 0.4142
## F-statistic: 33.52 on 1 and 45 DF, p-value: 0.0000006404
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 responseModel
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## [1] 992Get results
## lower upper
## estimateP1 2.32471130950857 2.675270900988
## ts1 5.32476946879407 6.392733318809
## ps1 0.00000001819721 0.000002109167
## attr(,"Probability")
## [1] 0.9495968Conclusion
PENDING.
NIR
Question
Can the darkness in the underlaying pigment affect the NIR reflectivity?
Correlation
Species level:
##
## Pearson's product-moment correlation
##
## data: RHandRefl$NIR and RHandRefl$Pol
## t = 4.8707, df = 45, p-value = 0.00001413
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.3613412 0.7484313
## sample estimates:
## cor
## 0.5875417
Linear model
##
## Call:
## lm(formula = RHandRefl$NIR ~ RHandRefl$Pol + RHandRefl$VIS)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.9226 -4.3421 -0.2965 2.4271 24.2257
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 20.1240 3.4000 5.919 0.000000444 ***
## RHandRefl$Pol 2.1381 0.8566 2.496 0.016379 *
## RHandRefl$VIS 0.6956 0.1947 3.573 0.000871 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.266 on 44 degrees of freedom
## Multiple R-squared: 0.4925, Adjusted R-squared: 0.4694
## F-statistic: 21.35 on 2 and 44 DF, p-value: 0.0000003315
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
## -3.6004 -1.8025 -0.3117 0.8898 8.8684
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.827
## lower bound : 0.000, p = 0.00854
## upper bound : 1.000, p = 0.10273
## 95.0% CI : (0.340, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 26.50712 5.84639 4.5339 0.00004418 ***
## Pol 1.97254 0.80127 2.4618 0.0178146 *
## VIS 0.73731 0.19706 3.7416 0.0005261 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.286 on 44 degrees of freedom
## Multiple R-squared: 0.548, Adjusted R-squared: 0.5274
## F-statistic: 26.67 on 2 and 44 DF, p-value: 0.00000002592
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 responseModel
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] 1001Get results
## lower upper
## estimateP1 1.568982596475 2.343066714
## ts1 2.227658400835 2.925192813
## ps1 0.004720836798 0.028982730
## estimateP2 0.662245251570 0.813247260
## tP2 3.317189090327 4.402433444
## pP2 0.000006930049 0.001525176
## attr(,"Probability")
## [1] 0.95005FRS
Question
Can the darkness in the underlaying pigment affect the NIR reflectivity?
Correlation
Species level:
##
## Pearson's product-moment correlation
##
## data: RHandRefl$FRS and RHandRefl$Pol
## t = 1.7038, df = 45, p-value = 0.09531
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.04410901 0.49812836
## sample estimates:
## cor
## 0.246176
Linear model
##
## Call:
## lm(formula = RHandRefl$FRS ~ RHandRefl$Pol)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.3853 -5.1370 -0.4918 3.1890 26.0037
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -11.366 3.438 -3.306 0.00187 **
## RHandRefl$Pol 1.237 0.726 1.704 0.09531 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.462 on 45 degrees of freedom
## Multiple R-squared: 0.0606, Adjusted R-squared: 0.03973
## F-statistic: 2.903 on 1 and 45 DF, p-value: 0.09531
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
## -4.0152 -1.7088 -0.6232 0.9233 9.9783
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.875
## lower bound : 0.000, p = 0.005182
## upper bound : 1.000, p = 0.2351
## 95.0% CI : (0.400, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.26850 6.28085 -0.6796 0.50024
## Pol 1.09243 0.61722 1.7699 0.08352 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.464 on 45 degrees of freedom
## Multiple R-squared: 0.06508, Adjusted R-squared: 0.04431
## F-statistic: 3.133 on 1 and 45 DF, p-value: 0.08352
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] 1001Get results
## lower upper
## estimateP1 0.86062030 1.3845539
## ts1 1.50152500 2.3154223
## ps1 0.02085044 0.1321538
## attr(,"Probability")
## [1] 0.95005Absorptivity
VIS
Question
Can polarisation predict visible absorptivity?
Correlation
##
## Pearson's product-moment correlation
##
## data: RHandAbs$VIS and RHandAbs$Pol
## t = -0.82271, df = 24, p-value = 0.4188
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.5196329 0.2369354
## sample estimates:
## cor
## -0.1656164
Linear model
##
## Call:
## lm(formula = RHandAbs$VIS ~ RHandAbs$Pol)
##
## Residuals:
## Min 1Q Median 3Q Max
## -23.7907 -5.6429 0.3487 7.6899 13.0918
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 83.909 5.863 14.312 0.000000000000301 ***
## RHandAbs$Pol -1.110 1.349 -0.823 0.419
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.852 on 24 degrees of freedom
## Multiple R-squared: 0.02743, Adjusted R-squared: -0.0131
## F-statistic: 0.6769 on 1 and 24 DF, p-value: 0.4188
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
## -7.9689 -1.5659 0.0845 2.3200 5.5711
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.958
## lower bound : 0.000, p = 0.52354
## upper bound : 1.000, p = 0.34748
## 95.0% CI : (NA, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 83.17557 9.94708 8.3618 0.0000000143 ***
## Pol -0.52023 1.18225 -0.4400 0.6638
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.619 on 24 degrees of freedom
## Multiple R-squared: 0.008003, Adjusted R-squared: -0.03333
## F-statistic: 0.1936 on 1 and 24 DF, p-value: 0.6638
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## [1] 971Get results
## lower upper
## estimateP1 -1.1101595 -0.09741585
## ts1 -0.8227129 -0.08824638
## ps1 0.4187710 0.91363493
## attr(,"Probability")
## [1] 0.9495366Conclusion
There’s no correlation between visible absorptivity and polarisation at the species level.
NIR
Question
Is polarisation correlated with NIR absorptivity?
Correlation
Species level:
##
## Pearson's product-moment correlation
##
## data: RHandAbs$NIR and RHandAbs$Pol
## t = -0.47705, df = 24, p-value = 0.6376
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.4667489 0.3017610
## sample estimates:
## cor
## -0.09691968
Linear model
##
## Call:
## lm(formula = RHandAbs$NIR ~ RHandAbs$Pol + RHandAbs$VIS)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.142 -4.009 1.767 4.109 16.372
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -74.2454 14.2904 -5.195 0.000028762 ***
## RHandAbs$Pol 0.5027 1.0801 0.465 0.646
## RHandAbs$VIS 1.3596 0.1611 8.438 0.000000017 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.777 on 23 degrees of freedom
## Multiple R-squared: 0.7581, Adjusted R-squared: 0.7371
## F-statistic: 36.04 on 2 and 23 DF, p-value: 0.00000008153
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
## -3.6619 -1.2537 0.6996 1.4723 5.9650
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.882
## lower bound : 0.000, p = 0.019671
## upper bound : 1.000, p = 0.46582
## 95.0% CI : (0.239, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -81.03156 13.82652 -5.8606 0.000005671739 ***
## Pol 0.51156 0.89074 0.5743 0.5713
## VIS 1.35538 0.14503 9.3458 0.000000002697 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.342 on 23 degrees of freedom
## Multiple R-squared: 0.7922, Adjusted R-squared: 0.7741
## F-statistic: 43.85 on 2 and 23 DF, p-value: 0.00000001421
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)
dfRHPol7 <- ldply(runsRHPol7, data.frame)
length(dfRHPol7[,1])## [1] 1001Get results
## lower upper
## estimateP1 0.173351686790757 0.706918816313237
## ts1 0.174959664134327 0.885683520350157
## ps1 0.384954689714818 0.862642397828346
## estimateP2 1.309132205945973 1.376597122439301
## tP2 8.654408621330273 10.025365076043126
## pP2 0.000000000393297 0.000000008174573
## attr(,"Probability")
## [1] 0.95005FRS
Question
Is polarisation correlated with the NIR residuals corrected by phylogeny?
Correlation
Species level:
##
## Pearson's product-moment correlation
##
## data: RHandAbs$FRS and RHandAbs$Pol
## t = 0.46896, df = 24, p-value = 0.6433
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.3032555 0.4654613
## sample estimates:
## cor
## 0.09528982
Linear model
##
## Call:
## lm(formula = RHandAbs$FRS ~ RHandAbs$Pol)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.128 -4.098 1.842 4.087 16.505
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.0768 4.5312 1.120 0.274
## RHandAbs$Pol 0.4891 1.0429 0.469 0.643
##
## Residual standard error: 7.614 on 24 degrees of freedom
## Multiple R-squared: 0.00908, Adjusted R-squared: -0.03221
## F-statistic: 0.2199 on 1 and 24 DF, p-value: 0.6433
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
## -3.6342 -1.2773 0.6938 1.4857 5.9980
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.883
## lower bound : 0.000, p = 0.019636
## upper bound : 1.000, p = 0.45962
## 95.0% CI : (0.239, NA)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.07307 6.41864 -0.3230 0.7495
## Pol 0.50711 0.86494 0.5863 0.5631
##
## Residual standard error: 2.296 on 24 degrees of freedom
## Multiple R-squared: 0.01412, Adjusted R-squared: -0.02696
## F-statistic: 0.3437 on 1 and 24 DF, p-value: 0.5631
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)
dfRHPol8 <- ldply(runsRHPol8, data.frame)
length(dfRHPol8[,1])## [1] 1001Get results
## lower upper
## estimateP1 0.1747109 0.6713303
## ts1 0.1923224 0.8844235
## ps1 0.3841219 0.8476692
## attr(,"Probability")
## [1] 0.95005Transmissivity
VIS
Question
Can polarisation predict visible transmissivity?
Correlation
##
## Pearson's product-moment correlation
##
## data: RHandTra2$VIS and RHandTra2$Pol
## t = -0.1629, df = 24, p-value = 0.872
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.4152401 0.3587364
## sample estimates:
## cor
## -0.03323364
Plot:
par(mfrow=c(1,2))
plot(RHandTra2$Pol, RHandTra2$VIS,
ylab=("Transmissivity in VIS (%)"),
xlab = ("Reflectance in RCP filter"))
plot(RHandAbs2$Pol, RHandAbs2$VIS,
ylab=("Absorptivity in VIS (%)"),
xlab = ("Reflectance in RCP filter"))
Neither absorptivity nor transmissivity can be explained by the reflectance in RCP filter in the subset of 20 species. To further understand these correlations, a larger n is needed.
Mechanisms
We tested the correlations between optical properties, i.e. reflectivity and absorptivity, and polarisation in different subsets of species since these correlations may depend on the underlaying mechanisms.
Subsets
We will compare three scenarios: 1) all species included 2) excluding the species with a white underlay since they clearly have a different mechanism that doe snot involve polarisation. 3) some species may be combining additional structures or disorder with chiral structures. This would change the correlations between polarisation and reflectivity, thus in this scenario we also exclude these cases (clustering apart from the other species).
Subset no underlay:
# reflectivity
NoUnderlayRefl <-
PolRefl %>%
dplyr::filter(phylogeny_name != "Paraschizognathus_ocularis" &
phylogeny_name != "Paraschizognathus_prasinus" &
phylogeny_name != "Paraschizognathus_olivaceous" &
phylogeny_name != "Anoplognathus_prasinus" &
phylogeny_name != "Xylonichus_sp")
# absorptivity
NoUnderlayAbs <-
PolAbs %>%
dplyr::filter(phylogeny_name != "Paraschizognathus_ocularis" &
phylogeny_name != "Paraschizognathus_prasinus" &
phylogeny_name != "Paraschizognathus_olivaceous" &
phylogeny_name != "Anoplognathus_prasinus" &
phylogeny_name != "Xylonichus_sp")Subset no special cases:
# reflectivity
NoSpecialCasesRefl <-
PolRefl %>%
dplyr::filter(phylogeny_name != "Paraschizognathus_ocularis" &
phylogeny_name != "Paraschizognathus_prasinus" &
phylogeny_name != "Paraschizognathus_olivaceous" &
phylogeny_name != "Anoplognathus_prasinus" &
phylogeny_name != "Xylonichus_sp" &
phylogeny_name != "Epichrysus_lamprimoides" &
phylogeny_name != "Anoplostethus_laetus" &
phylogeny_name != "Anoplostethus_opalinus" &
phylogeny_name != "Anoplostethus_roseus" &
phylogeny_name != "Anoplognathus_parvulus" &
phylogeny_name != "Anoplognathus_aureus")
# absorptivity
NoSpecialCasesAbs <-
PolAbs %>%
dplyr::filter(phylogeny_name != "Paraschizognathus_ocularis" &
phylogeny_name != "Paraschizognathus_prasinus" &
phylogeny_name != "Paraschizognathus_olivaceous" &
phylogeny_name != "Anoplognathus_prasinus" &
phylogeny_name != "Xylonichus_sp" &
phylogeny_name != "Epichrysus_lamprimoides" &
phylogeny_name != "Anoplostethus_laetus" &
phylogeny_name != "Anoplostethus_opalinus" &
phylogeny_name != "Anoplostethus_roseus" &
phylogeny_name != "Anoplognathus_parvulus" &
phylogeny_name != "Anoplognathus_aureus")Subset no chiral composites:
# reflectivity
NoCCRefl <-
PolRefl %>%
dplyr::filter(phylogeny_name != "Epichrysus_lamprimoides" &
phylogeny_name != "Anoplostethus_laetus" &
phylogeny_name != "Anoplostethus_opalinus" &
phylogeny_name != "Anoplostethus_roseus" &
phylogeny_name != "Anoplognathus_parvulus" &
phylogeny_name != "Anoplognathus_aureus")
# absorptivity
NoCCAbs <-
PolAbs %>%
dplyr::filter(phylogeny_name != "Epichrysus_lamprimoides" &
phylogeny_name != "Anoplostethus_laetus" &
phylogeny_name != "Anoplostethus_opalinus" &
phylogeny_name != "Anoplostethus_roseus" &
phylogeny_name != "Anoplognathus_parvulus" &
phylogeny_name != "Anoplognathus_aureus")Reflectivity
VIS
All species
##
## Pearson's product-moment correlation
##
## data: PolRefl$VIS and PolRefl$Pol
## t = 1.5434, df = 45, p-value = 0.1297
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.06728212 0.48044869
## sample estimates:
## cor
## 0.2242168
No underlay
##
## Pearson's product-moment correlation
##
## data: NoUnderlayRefl$VIS and NoUnderlayRefl$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
No special cases
##
## Pearson's product-moment correlation
##
## data: NoSpecialCasesRefl$VIS and NoSpecialCasesRefl$Pol
## t = -1.5685, df = 34, p-value = 0.126
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.54203774 0.07519306
## sample estimates:
## cor
## -0.2597597
No chiral composites
##
## Pearson's product-moment correlation
##
## data: NoCCRefl$VIS and NoCCRefl$Pol
## t = -1.4969, df = 39, p-value = 0.1425
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.50456471 0.08031186
## sample estimates:
## cor
## -0.233098NIR
All species
##
## Pearson's product-moment correlation
##
## data: PolRefl$NIR and PolRefl$Pol
## t = -1.7419, df = 45, p-value = 0.08836
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.50225033 0.03862105
## sample estimates:
## cor
## -0.2513332
No underlay
##
## Pearson's product-moment correlation
##
## data: NoUnderlayRefl$NIR and NoUnderlayRefl$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
No special cases
##
## Pearson's product-moment correlation
##
## data: NoSpecialCasesRefl$NIR and NoSpecialCasesRefl$Pol
## t = -3.4904, df = 34, p-value = 0.001356
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.7205635 -0.2226514
## sample estimates:
## cor
## -0.5136167
No chiral composites
##
## Pearson's product-moment correlation
##
## data: NoCCRefl$NIR and NoCCRefl$Pol
## t = -3.7868, df = 39, p-value = 0.0005151
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.7124920 -0.2508605
## sample estimates:
## cor
## -0.5184949FRS
All species
##
## Pearson's product-moment correlation
##
## data: PolRefl$FRS and PolRefl$Pol
## t = -4.2698, df = 45, p-value = 0.00009974
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.7140250 -0.2953307
## sample estimates:
## cor
## -0.5369585
No underlay
##
## Pearson's product-moment correlation
##
## data: NoUnderlayRefl$FRS and NoUnderlayRefl$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
No special cases
##
## Pearson's product-moment correlation
##
## data: NoSpecialCasesRefl$FRS and NoSpecialCasesRefl$Pol
## t = -4.0927, df = 34, p-value = 0.0002481
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.7596570 -0.3031927
## sample estimates:
## cor
## -0.5745024
No chiral composites
##
## Pearson's product-moment correlation
##
## data: NoCCRefl$FRS and NoCCRefl$Pol
## t = -3.8721, df = 39, p-value = 0.0004006
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.7181806 -0.2617452
## sample estimates:
## cor
## -0.5269616Absorptivity
VIS
All species
##
## Pearson's product-moment correlation
##
## data: PolAbs$VIS and PolAbs$Pol
## t = 0.041288, df = 24, p-value = 0.9674
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.3801654 0.3944917
## sample estimates:
## cor
## 0.008427517
No underlay
##
## Pearson's product-moment correlation
##
## data: NoUnderlayAbs$VIS and NoUnderlayAbs$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
No special cases
##
## Pearson's product-moment correlation
##
## data: NoSpecialCasesAbs$VIS and NoSpecialCasesAbs$Pol
## t = 2.8367, df = 14, p-value = 0.01319
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1547983 0.8463786
## sample estimates:
## cor
## 0.6041449
No chiral composites
##
## Pearson's product-moment correlation
##
## data: NoCCAbs$VIS and NoCCAbs$Pol
## t = 3.2942, df = 19, p-value = 0.003815
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.2314862 0.8209405
## sample estimates:
## cor
## 0.6029232NIR
All species
##
## Pearson's product-moment correlation
##
## data: PolAbs$NIR and PolAbs$Pol
## t = 1.0748, df = 24, p-value = 0.2932
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1887276 0.5555297
## sample estimates:
## cor
## 0.2142898
No underlay
##
## Pearson's product-moment correlation
##
## data: NoUnderlayAbs$NIR and NoUnderlayAbs$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
No special cases
##
## Pearson's product-moment correlation
##
## data: NoSpecialCasesAbs$NIR and NoSpecialCasesAbs$Pol
## t = 2.6934, df = 14, p-value = 0.01748
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1245892 0.8374076
## sample estimates:
## cor
## 0.5842138
No chiral composites
##
## Pearson's product-moment correlation
##
## data: NoCCAbs$NIR and NoCCAbs$Pol
## t = 3.6547, df = 19, p-value = 0.001685
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.2917107 0.8409388
## sample estimates:
## cor
## 0.6424902FRS
All species
##
## Pearson's product-moment correlation
##
## data: PolAbs$FRS and PolAbs$Pol
## t = 2.2597, df = 24, p-value = 0.03319
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.03759808 0.69366079
## sample estimates:
## cor
## 0.41885
No underlay
##
## Pearson's product-moment correlation
##
## data: NoUnderlayAbs$FRS and NoUnderlayAbs$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
No special cases
##
## Pearson's product-moment correlation
##
## data: NoSpecialCasesAbs$FRS and NoSpecialCasesAbs$Pol
## t = 0.68241, df = 14, p-value = 0.5061
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.3471587 0.6199864
## sample estimates:
## cor
## 0.1794233
No chiral composites
##
## Pearson's product-moment correlation
##
## data: NoCCAbs$FRS and NoCCAbs$Pol
## t = 1.7531, df = 19, p-value = 0.0957
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.06978245 0.69317408
## sample estimates:
## cor
## 0.373145Plots
Degree of polarisation
ReflByInd <- read.csv(here::here("Data/FromCode/ConsolidatedReflectivityInd.csv"))[-1]
AbsByInd0 <- read.csv(here::here("Data/FromCode/ConsolidatedAbsorptivityInd.csv"))[-1]
AbsByInd0$A_Res <- residuals ( lm (AbsByInd0$A_NIR~AbsByInd0$A_VIS))
# Mean and sd per species Absorptivity
AbsByInd <-
AbsByInd0 %>%
dplyr::mutate("spp" = substr(ind, 1, 4)) %>%
dplyr::select(spp, A_VIS, A_NIR, A_Res) %>%
dplyr::group_by(spp) %>%
dplyr::summarise(
meanVISA = mean(A_VIS),
meanNIRA = mean(A_NIR),
meanResA = mean(A_Res),
sdVISA = sd(A_VIS),
sdNIRA = sd(A_NIR),
sdResA = sd(A_Res),
)
AbsByInd <- as.data.frame(AbsByInd)
#mean and sd of transmittance per species
TraByInd <-
AbsByInd0 %>%
dplyr::mutate("spp" = substr(ind, 1, 4)) %>%
dplyr::select(spp, T_VIS) %>%
dplyr::group_by(spp) %>%
dplyr::summarise(
meanVIST = mean(T_VIS),
sdVIST = sd(T_VIS)
)
TraByInd <- as.data.frame(TraByInd)
# Mean and sd per species Reflectivity including pol
ToPlotPolRefl <-
merge(ReflByInd[,c(1,3,4,5)], PolByInd, by = "ind") %>%
dplyr::select(-phylogeny_name) %>%
dplyr::group_by(spp) %>%
dplyr::summarise(
meanVIS = mean(R_VIS),
meanNIR = mean(R_NIR),
meanRes = mean(Res),
meanPol = mean(Polarization),
sdVIS = sd(R_VIS),
sdNIR = sd(R_NIR),
sdRes = sd(Res),
sdPol = sd(Polarization)
)
ToPlotPolRefl <- as.data.frame(ToPlotPolRefl)
rownames(ToPlotPolRefl) <- ToPlotPolRefl$spp
ToPlotPolAbs <- merge(
ToPlotPolRefl[,c("spp", "meanPol", "sdPol")], AbsByInd)
rownames(ToPlotPolAbs) <- ToPlotPolAbs$spp
ToPlotPolTra <- merge(
ToPlotPolRefl[,c("spp", "meanPol", "sdPol")], TraByInd)
rownames(ToPlotPolTra) <- ToPlotPolTra$sppSubset for plotting
UnderlayRefl <- ToPlotPolRefl %>%
dplyr::filter(spp == "xyls" |
spp == "xyle" |
spp == "pczo" |
spp == "pczp" |
spp == "pczv" |
spp == "prsi")
UnderlayAbs <- ToPlotPolAbs %>%
dplyr::filter(spp == "xyls" |
spp == "xyle" |
spp == "pczo" |
spp == "pczp" |
spp == "pczv" |
spp == "prsi")
NoUnderlayRefl <- ToPlotPolRefl %>%
dplyr::filter(spp != "xyls" |
spp != "xyle" |
spp != "pczo" |
spp != "pczp" |
spp != "pczv" |
spp != "prsi")
NoUnderlayAbs <- ToPlotPolAbs %>%
dplyr::filter(spp != "xyls" |
spp != "xyle" |
spp != "pczo" |
spp != "pczp" |
spp != "pczv" |
spp != "prsi")Predictions from PGLS for the NIR reflectivity and absorptivity
## set values for predictors:
PolDeg <- seq(range(ToPlotPolRefl$meanPol)[1],
range(ToPlotPolRefl$meanPol)[2],0.01) # A range of Pol values
new1<-data.frame("Pol" = PolDeg) # data frame
RefDeg<-predict(pglsPDFRS, newdata = new1,
type = "response") # Expected reflectivity
trend1<-data.frame(PolDeg,RefDeg) # Pol values and expected reflectivity joint
AbsDeg<-predict(pglsPDFRSAbs, newdata = new1,
type = "response") # Expected absorptivity
trend2<-data.frame(PolDeg,AbsDeg) # Pol values and expected reflectivity jointPredictions from PGLS for the NIR absorptivity
## set values for predictors:
PolDeg <- seq(range(ToPlotPolRefl$meanPol)[1],
range(ToPlotPolRefl$meanPol)[2],0.01) # A range of Pol values
new1<-data.frame("Pol" = PolDeg) # data frame
RefDeg<-predict(pglsPDFRS, newdata = new1,
type = "response") # Expected reflectivity
trend1<-data.frame(PolDeg,RefDeg) # Pol values and expected reflectivity jointNudgeXVal = 0.02
NudgeYVal = 0.9
PVRPlot <- ggplot(data = ToPlotPolRefl, aes(x = meanPol, y = meanVIS),
size = 2.2, alpha = 0.6) +
geom_text(
label=rownames(ToPlotPolRefl),
nudge_x = NudgeXVal, nudge_y = NudgeYVal,
col="gray", size=3
) +
# error bars
geom_errorbar(aes(
ymin = meanVIS - sdVIS,
ymax = meanVIS + sdVIS),
col = "#cecec2"
) +
geom_errorbarh(aes(
xmin = meanPol - sdPol,
xmax = meanPol + sdPol
),
col = "#cecec2"
) +
geom_point(data = NoUnderlayRefl, size = 2.2, pch = 21, fill = "#648fff", color = "black", alpha = 0.7) +
geom_point(data = UnderlayRefl, size = 2.4, pch = 22, fill = "#F7E949", color = "black", alpha = 0.7) +
theme_minimal() +
theme(legend.position = "none") +
xlab("Degree of Polarisation (L - R)/(L + R)")+
ylab("Visible reflectivity (%)") +
ylim(0, 50)
PNRPlot <- ggplot(data = ToPlotPolRefl, aes(x = meanPol, y = meanRes),
size = 2.2, alpha = 0.6) +
geom_text(
label=rownames(ToPlotPolRefl),
nudge_x = NudgeXVal, nudge_y = NudgeYVal,
col="gray", size=3
) +
# error bars
geom_errorbar(aes(
ymin = meanRes - sdRes,
ymax = meanRes + sdRes),
col = "#cecec2"
) +
geom_errorbarh(aes(
xmin = meanPol - sdPol,
xmax = meanPol + sdPol
),
col = "#cecec2"
) +
geom_point(data = NoUnderlayRefl, size = 2.2, pch = 21, fill = "#648fff", color = "black", alpha = 0.7) +
geom_point(data = UnderlayRefl, size = 2.4, pch = 22, fill = "#F7E949", color = "black", alpha = 0.7) +
geom_line(data=trend1, aes(x= PolDeg, y= RefDeg), color = "black") +
theme_minimal() +
theme(legend.position = "none") +
xlab("Degree of Polarisation (L - R)/(L + R)")+
ylab("NIR (Res) reflectivity (%)")
PVAPlot <- ggplot(data = ToPlotPolAbs, aes(x = meanPol, y = meanVISA),
size = 2.2, alpha = 0.6) +
geom_text(
label=rownames(ToPlotPolAbs),
nudge_x = NudgeXVal, nudge_y = NudgeYVal,
col="gray", size=3
) +
# error bars
geom_errorbar(aes(
ymin = meanVISA - sdVISA,
ymax = meanVISA + sdVISA),
col = "#cecec2"
) +
geom_errorbarh(aes(
xmin = meanPol - sdPol,
xmax = meanPol + sdPol
),
col = "#cecec2"
) +
geom_point(data = NoUnderlayAbs, size = 2.2, pch = 21, fill = "#648fff", color = "black", alpha = 0.7) +
geom_point(data = UnderlayAbs, size = 2.4, pch = 22, fill = "#F7E949", color = "black", alpha = 0.7) +
theme_minimal() +
theme(legend.position = "none") +
xlab("Degree of Polarisation (L - R)/(L + R)")+
ylab("Visible absorptivity (%)") +
ylim(0, 100)
PNAPlot <- ggplot(data = ToPlotPolAbs, aes(x = meanPol, y = meanResA),
size = 2.2, alpha = 0.6) +
geom_text(
label=rownames(ToPlotPolAbs),
nudge_x = NudgeXVal, nudge_y = NudgeYVal,
col="gray", size=3
) +
# error bars
geom_errorbar(aes(
ymin = meanResA - sdResA,
ymax = meanResA + sdResA),
col = "#cecec2"
) +
geom_errorbarh(aes(
xmin = meanPol - sdPol,
xmax = meanPol + sdPol
),
col = "#cecec2"
) +
geom_point(data = NoUnderlayAbs, size = 2.2, pch = 21, fill = "#648fff", color = "black", alpha = 0.7) +
geom_point(data = UnderlayAbs, size = 2.4, pch = 22, fill = "#F7E949", color = "black", alpha = 0.7) +
geom_line(data=trend2, aes(x= PolDeg, y= AbsDeg), color = "black") +
theme_minimal() +
theme(legend.position = "none") +
xlab("Degree of Polarisation (L - R)/(L + R)")+
ylab("NIR (Res) absorptivity (%)")
grid.arrange(PVRPlot, PNRPlot, PVAPlot, PNAPlot, nrow = 2)## Warning: Removed 3 rows containing missing values (`geom_errorbarh()`).
## Removed 3 rows containing missing values (`geom_errorbarh()`).
Right handed polarisation
Setting up data frames by individual for plots
# Import equalised values for right handed polarisation
RHand0I<-read.csv(here::here("Data/FromCode/PolarzElytra_LERGBValues.csv"))[-1]%>%
dplyr::filter(Type == "RCP") %>%
dplyr::select(-Type) %>%
dplyr::rename("RHPol" = RGB_ave)
# Load degree of polarisation per individual
PolDeg0I <- read.csv(here::here("Data/FromCode/PolarzElytraByInd.csv"))[-1] %>%
dplyr::select(ind, Polarization)
# Merge RCP with degree of polarisation
RHandWPolDeg <- merge(RHand0I, PolDeg0I)
# Join them with the phylogeny names
SppNamesI<-read.csv(here::here("Data/9_CodesAndSpecies.csv")) %>%
dplyr::rename("ind" = Ind)
RHandWPolDegByInd <- merge(RHandWPolDeg, SppNamesI)
# Add 4-letter spp code
RHandWPolDegByInd <-
RHandWPolDegByInd %>%
dplyr::mutate(spp = substr(ind, 1, 4))
# Combine with visible reflectivity
RHandPolDegReflByInd <- merge(RHandWPolDegByInd, ReflByInd) %>%
dplyr::select(-R_ALL, -R_NIR, -Res, -size)And re order them to plot
ToPlotRHPolRefl <-
RHandPolDegReflByInd %>%
dplyr::select(-phylogeny_name) %>%
dplyr::group_by(spp) %>%
dplyr::summarise(
meanVIS = mean(R_VIS),
meanRHPol = mean(RHPol),
meanPol = mean(Polarization),
sdVIS = sd(R_VIS),
sdRHPol = sd(RHPol),
sdPol = sd(Polarization)
)
ToPlotRHPolRefl <- as.data.frame(ToPlotRHPolRefl)
rownames(ToPlotRHPolRefl) <- ToPlotRHPolRefl$sppPredictions from PGLS for the VIS reflectivity
## set values for predictors:
RHPolSeq <- seq(range(ToPlotRHPolRefl$meanRHPol)[1],
range(ToPlotRHPolRefl$meanRHPol)[2],0.01) # A range of Pol values
new1RH<-data.frame("Pol" = RHPolSeq) # data frame
RefRHPolPred <- predict(pglsPDVISRH, newdata = new1RH,
type = "response") # Expected reflectivity
trend1RH<-data.frame(RHPolSeq, RefRHPolPred) # Pol values and expected reflectivity jointPlots
NudgeXVal = 0.3
NudgeYVal = 0.7
RCPVRPlot <- ggplot(data = ToPlotRHPolRefl, aes(x = meanRHPol, y = meanVIS),
size = 2.2, alpha = 0.6) +
geom_text(
label=rownames(ToPlotRHPolRefl),
nudge_x = NudgeXVal, nudge_y = NudgeYVal,
col="gray", size=3
) +
# error bars
geom_errorbar(aes(
ymin = meanVIS - sdVIS,
ymax = meanVIS + sdVIS),
col = "#cecec2"
) +
geom_errorbarh(aes(
xmin = meanRHPol - sdRHPol,
xmax = meanRHPol + sdRHPol
),
col = "#cecec2"
) +
geom_point(data = ToPlotRHPolRefl, size = 2.2, pch = 21,
fill = "#648fff", color = "black", alpha = 0.7) +
geom_line(data=trend1RH, aes(x= RHPolSeq, y= RefRHPolPred), color = "black") +
theme_minimal() +
theme(legend.position = "none") +
xlab("Right Handed Polarisation (RGB mean)")+
ylab("Visible reflectivity (%)") +
ylim(0, 40)
RCPVRPlot ## Warning: Removed 3 rows containing missing values (`geom_errorbarh()`).