Introduction

In this document we tested if the optical properties of Christmas Beetles can predict the heating of their elytra due to radiative heat in different stepectral bands (TOTAL: 400 to 1700nm, NIR: 700 to 1700 nm, VIS: 400 to 700 nm). We also tested if the presence of the elytra closed on top of the body can help reducing the heating rate of the beetle body.

We first processed the raw spectral data (reflectance and transmittance) and combined it with the irradiance spectra of the light source to obtain the values of reflectivity, transmissivity and absorptivity for each beetle in our sample (n=56 elytra representing 28 species of Christmas beetles)

Then we processed the heating raw data to obtain the total change in temperature after 5 minutes, and the maximum heating rate under each illumination (TOTAL, NIR and VIS).

For the single elytra experiments we used linear models and correlations to explore the relation between the optical properties and the two heating variables. For the beetle body experiments we used a bootstrap paired test to determine if the body of a recently dead beetle has a higher change in temperature when the elytra are open than when they are closed.

We provide details on data processing, statistical analysis and figures of the original manuscript submitted for peer review.

The following code is written in the R programming language.

The .Rmd file and Raw data can be found here: https://bit.ly/3OeBNU6 (This link contains information that may disclose author’s identity)

Libraries

Install packages

install.packages("ggplot2") # graphs

install.packages("wesanderson") # Colour palette

install.packages("gridExtra") # basic to combine multiple ggplot 

install.packages("rptR") # Repeatability analysis

install.packages("cowplot") # ggplot graphs alignment

install.packages("rsq") # Partial coeficients in the linear models

install.packages("pavo") # Processing spectral data

install.packages("dplyr") # Data wrangling

install.packages("jtools") # Reporting Results from Linear Models

install.packages("huxtable") # Reporting Results from Linear Models

install.packages("kableExtra") # To style tables in the HTML file

Load packages

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.1.3

library(wesanderson)

## Warning: package 'wesanderson' was built under R version 4.1.3

library(gridExtra)

## Warning: package 'gridExtra' was built under R version 4.1.3

library(rptR)

## Warning: package 'rptR' was built under R version 4.1.3

library(cowplot)

## Warning: package 'cowplot' was built under R version 4.1.3

library(rsq)

## Warning: package 'rsq' was built under R version 4.1.3

library(pavo)

## Warning: package 'pavo' was built under R version 4.1.3

library(tidyverse) 
library(dplyr) 
library(readr) 
library(jtools)

## Warning: package 'jtools' was built under R version 4.1.3

library(huxtable)

## Warning: package 'huxtable' was built under R version 4.1.3

library(kableExtra)

Data Sets

Raw Data

Spectral Data

# Reflectance
Reflect<-read.csv("../Data/1_Reflectancedata.csv") %>% 
  as.rspec(.) %>% 
  procspec(.,fixneg = "zero") 

# Transmittance
TopTrans<-read.csv("../Data/2_Transmittancetopdata.csv") %>% 
  filter (wl <= 1700 & wl >= 400) %>% 
  as.rspec(.) %>% 
  procspec(.,fixneg = "zero") 

BotTrans<-read.csv("../Data/3_Transmittanceventraldata.csv") %>% 
  filter (wl <= 1700 & wl >= 400) %>% 
  as.rspec(.) %>% 
  procspec(.,fixneg = "zero") 

# Solar Simulator and Filters
filters<-read.csv("../Data/4_SolarSimulatorAndFilters.csv") %>% 
  as.rspec(.) %>% 
  procspec(.,fixneg = "zero") 
  
# Solar Irradiance for the interval relevant for this manuscript
theorySun<-read.csv("../Data/5_TheoreticalIrradiance.csv") %>% 
  as.rspec(.) %>% 
  procspec(.,fixneg = "zero") 

# Standard Solar Irradiance from 280 to 4000 nm
# ASTM G173-03 Reference Spectra Derived from SMARTS v. 2.9.2 (AM1.5)
# Global tilt  W*m-2*nm-1
# from https://www.pveducation.org/pvcdrom/appendices/standard-solar-spectra 
theorySunALL<-read.csv("../Data/12_SolarIrradianceALLWavelengths.csv") %>%  
  as.rspec(., whichwl = "wl")

Heating Rates Data

Note: This data frame has been revised to ensure the alignment in times is correct for all species (columns)

TempRaw <- read.csv ("../Data/6_SingleElytronTempDataRaw.csv")

BodyRaw <- read.csv ("../Data/8_BeetleBodyTempDataRaw.csv")

Others

This data set includes information as species, batch (for repeatability or for the experiment), size and side of the platform

Other <- read.csv ("../Data/7_SizeAndSide.csv")

This data set is an example of the raw temperature data obtained by the thermometer:

Setup <- read.csv("../Data/9_ExampleSetup.csv")

Data Processing

Optical Properties

Reflectivity

Add filters

### One data frame for Total spectrum
forTOT<-Reflect
forTOT$irr<-filters$Full.f

### One data frame for NIR
forNIR<-Reflect
forNIR$irr<-filters$nir.f

### One data frame for vis
forVIS<-Reflect
forVIS$irr<-filters$vis.f

Multiplication function

# Define the following function to find the multiplication

Find.multiplication<-function(s){
  vector2<-rep("NA", length(s[ , 1]))
  for( i in 2 : length(s) - 1){
    for (y in 1 : length(s[ , 1])){
      vector2[y] <- s[y, i] * s[y, "irr"] # multiply each column by the last one
    }
    s[i] <- as.numeric(vector2) 
  }
  return(s)
}

Apply

# Total spectrum 
PrelValReflTOT <- 
  forTOT %>% 
  Find.multiplication(.) %>% # Apply function
  select (-wl, # the wavelengths were altered, eliminate
          -irr) %>% # irradiance not required anymore, eliminate
  mutate (wl = Reflect$wl) %>%  # add the correct wavelength column
  select(wl, everything()) %>% # set the wl column as column 1
  as.rspec(.) # convert to rspec object


# NIR

PrelValReflNIR <- 
  forNIR %>%   # Change data frame for NIR
  Find.multiplication(.) %>% 
  select (-wl, 
          -irr) %>% 
  mutate (wl = Reflect$wl) %>%  
  select(wl, everything()) %>% 
  as.rspec(.) 


# VIS

PrelValReflVIS <- 
  forVIS %>%  # Change data frame for VIS
  Find.multiplication(.) %>% 
  select (-wl, 
          -irr) %>% 
  mutate (wl = Reflect$wl) %>%  
  select(wl, everything()) %>% 
  as.rspec(.)

Calculate AUC

Area under the curve

# Calculate AUC for preliminary values of the samples

r.TOT <- summary(PrelValReflTOT)$B1 
r.NIR <- summary(PrelValReflNIR)$B1 
r.VIS <- summary(PrelValReflVIS)$B1

# Calculate AUC for the filters

AUCsun <- summary(filters[,c(1,5)])$B1 # Portal open -  full spectrum
AUCnir <- summary(filters[,c(1,7)])$B1 # When using filter to let pass NIR
AUCvis <- summary(filters[,c(1,6)])$B1 # When using filter to let pass vis

Standardise

Reflectivity<-data.frame("spp" = names(Reflect[2:57]),
                         "FUL" = r.TOT / AUCsun,
                         "NIR" = r.NIR / AUCnir,
                         "VIS" = r.VIS / AUCvis)
Reflectivity <-
  Reflectivity %>% 
  arrange (spp)

Transmissivity

Add filters

### One data frame for Total spectrum
TforTOT <- TopTrans
TforTOT$irr <- filters$Full.f

### One data frame for NIR
TforNIR <- TopTrans
TforNIR$irr <- filters$nir.f

### One data frame for vis
TforVIS <- TopTrans
TforVIS$irr <- filters$vis.f

Multiplication function

# Define the following function to find the multiplication

Find.multiplication<-function(s){
  vector2<-rep("NA", length(s[ , 1]))
  for( i in 2 : length(s) - 1){
    for (y in 1 : length(s[ , 1])){
      vector2[y] <- s[y, i] * s[y, "irr"] # multiply each column by the last one
    }
    s[i] <- as.numeric(vector2) 
  }
  return(s)
}

# This is the same function as the one in reflectivity

Apply

# Total spectrum 
PrelValTransTOT <- 
  TforTOT %>% 
  Find.multiplication(.) %>% # Apply function
  select (-wl, # the wavelengths were altered, eliminate
          -irr) %>% # irradiance not required anymore, eliminate
  mutate (wl = Reflect$wl) %>%  # add the correct wavelength column
  select(wl, everything()) %>% # set the wl column as column 1
  as.rspec(.) # convert to rspec object


# NIR

PrelValTransNIR <- 
  TforNIR %>%   # Change data frame for NIR
  Find.multiplication(.) %>% 
  select (-wl, 
          -irr) %>% 
  mutate (wl = Reflect$wl) %>%  
  select(wl, everything()) %>% 
  as.rspec(.) 


# VIS

PrelValTransVIS <- 
  TforVIS %>%  # Change data frame for VIS
  Find.multiplication(.) %>% 
  select (-wl, 
          -irr) %>% 
  mutate (wl = Reflect$wl) %>%  
  select(wl, everything()) %>% 
  as.rspec(.)

Calculate AUC

Area under the curve

# Calculate AUC for preliminary values of the samples

t.TOT <- summary(PrelValTransTOT)$B1 
t.NIR <- summary(PrelValTransNIR)$B1 
t.VIS <- summary(PrelValTransVIS)$B1

# Calculate AUC for the filters (same as reflectivity)

AUCsun <- summary(filters[,c(1,5)])$B1 # Portal open -  full spectrum
AUCnir <- summary(filters[,c(1,7)])$B1 # When using filter to let pass NIR
AUCvis <- summary(filters[,c(1,6)])$B1 # When using filter to let pass vis

Standardise

Transmissivity<-data.frame("spp" = names(TopTrans[2:57]),
                         "FUL" = t.TOT / AUCsun,
                         "NIR" = t.NIR / AUCnir,
                         "VIS" = t.VIS / AUCvis)
Transmissivity <-
  Transmissivity %>% 
  arrange (spp)

Ventral Transmissivity

Add filters

### One data frame for Total spectrum
TVforTOT <- BotTrans
TVforTOT$irr <- filters$Full.f

### One data frame for NIR
TVforNIR <- BotTrans
TVforNIR$irr <- filters$nir.f

### One data frame for vis
TVforVIS <- BotTrans
TVforVIS$irr <- filters$vis.f

Multiplication function

# Define the following function to find the multiplication

Find.multiplication<-function(s){
  vector2<-rep("NA", length(s[ , 1]))
  for( i in 2 : length(s) - 1){
    for (y in 1 : length(s[ , 1])){
      vector2[y] <- s[y, i] * s[y, "irr"] # multiply each column by the last one
    }
    s[i] <- as.numeric(vector2) 
  }
  return(s)
}

# This is the same function as the one in reflectivity

Apply

# Total spectrum 
PrelValTransVentTOT <- 
  TVforTOT %>% 
  Find.multiplication(.) %>% # Apply function
  select (-wl, # the wavelengths were altered, eliminate
          -irr) %>% # irradiance not required anymore, eliminate
  mutate (wl = Reflect$wl) %>%  # add the correct wavelength column
  select(wl, everything()) %>% # set the wl column as column 1
  as.rspec(.) # convert to rspec object


# NIR

PrelValTransVentNIR <- 
  TVforNIR %>%   # Change data frame for NIR
  Find.multiplication(.) %>% 
  select (-wl, 
          -irr) %>% 
  mutate (wl = Reflect$wl) %>%  
  select(wl, everything()) %>% 
  as.rspec(.) 


# VIS

PrelValTransVentVIS <- 
  TVforVIS %>%  # Change data frame for VIS
  Find.multiplication(.) %>% 
  select (-wl, 
          -irr) %>% 
  mutate (wl = Reflect$wl) %>%  
  select(wl, everything()) %>% 
  as.rspec(.)

Calculate AUC

Area under the curve

# Calculate AUC for preliminary values of the samples

t.v.TOT <- summary(PrelValTransVentTOT)$B1 
t.v.NIR <- summary(PrelValTransVentNIR)$B1 
t.v.VIS <- summary(PrelValTransVentVIS)$B1

# Calculate AUC for the filters (same as reflectivity)

AUCsun <- summary(filters[,c(1,5)])$B1 # Portal open -  full spectrum
AUCnir <- summary(filters[,c(1,7)])$B1 # When using filter to let pass NIR
AUCvis <- summary(filters[,c(1,6)])$B1 # When using filter to let pass vis

Standardise

Transmiss.Ventral<-data.frame("spp" = names(BotTrans[2:57]),
                         "FUL" = t.v.TOT / AUCsun,
                         "NIR" = t.v.NIR / AUCnir,
                         "VIS" = t.v.VIS / AUCvis)
Transmiss.Ventral <-
  Transmiss.Ventral %>% 
  arrange (spp)

Consolidated

In this step we calculated absorptivity as 100 - (Reflectivity + transmissivity)

Optical <-
  data.frame (Reflectivity,
              Transmissivity, 
              Transmiss.Ventral) %>% 
  select (-spp.1, -spp.2) 

names(Optical) <- c( "spp", 
                     "Rf_TOT", "Rf_NIR", "Rf_VIS", # Reflectivity
                     "Td_TOT", "Td_NIR", "Td_VIS", # Dorsal Transmissivity
                     "Tv_TOT", "Tv_NIR", "Tv_VIS") # Ventral Transmissivity

Optical <-
  Optical %>% 
  mutate (DifferenceTransTOT = (Tv_TOT - Td_TOT) / Td_TOT * 100,
          DifferenceTransNIR = (Tv_NIR - Td_NIR) / Td_NIR * 100,
          DifferenceTransVIS = (Tv_VIS - Td_VIS) / Td_VIS * 100
          ) %>% # Percentaages Difference between dorsal and ventral
  mutate (Ab_TOT = 100 - (Rf_TOT + Td_TOT),   
          Ab_NIR = 100 - (Rf_NIR + Td_NIR),
          Ab_VIS = 100 - (Rf_VIS + Td_VIS)) # Calculate Absorptivity

Heating

Single Elytra

Delta temp after 5 min

ind <- (names(TempRaw)[3:96])

t0TOT <- 
  TempRaw %>% 
  filter (time == 280) %>%  # t0 time
  select (-light, -time) %>% 
  as.numeric (.)

tfTOT <- 
  TempRaw %>% 
  filter (time == 580) %>% # end time of the first illumination 
  select (-light, -time) %>% 
  as.numeric (.)

t0NIR <- 
  TempRaw %>% 
  filter (time == 880) %>%  # t0 time
  select (-light, -time) %>% 
  as.numeric (.)

tfNIR <- 
  TempRaw %>% 
  filter (time == 1180) %>% # end time of the first illumination 
  select (-light, -time) %>% 
  as.numeric (.)

t0VIS <- 
  TempRaw %>% 
  filter (time == 1480) %>%  # t0 time
  select (-light, -time) %>% 
  as.numeric (.)

tfVIS <- 
  TempRaw %>% 
  filter (time == 1780) %>% # end time of the first illumination 
  select (-light, -time) %>% 
  as.numeric (.)


TimePoints <- data.frame( ind, t0TOT, tfTOT, t0NIR, tfNIR, t0VIS, tfVIS )

DeltaT5 <- 
  TimePoints %>% 
  mutate (
    DT5_TOT = tfTOT - t0TOT,
    DT5_NIR = tfNIR - t0NIR,
    DT5_VIS = tfVIS - t0VIS ) %>% 
  select (1,8,9,10) %>% 
  arrange (ind)

Max Heating Rates

Calculated as the max. slope in degrees celsius / sec. Max. change in a 20 sec interval divided by 20

# TOTA LIGHT (nir + vis)

# Prepare subset
SlopesTOT2 <-
  TempRaw %>% 
  select (-light) %>%  # remove the column "light"
  filter (time >= 300 & time <= 600) %>%  # keep the interval with the light on
  select (-time) %>% 
  select (sort(colnames(.)))

SlopesTOT <- SlopesTOT2 # rename because it will be overwritten by the function

# Run this function to find the slope for each 20 sec

for (i in 1: length (SlopesTOT)) { # for each column
  for (ii in 2: length (SlopesTOT[,i])){ # take each value (start in the second)
    SlopesTOT [ii-1,i] <- # overwrite the previous value with the result
      (abs( (SlopesTOT[ii-1,i]) - (SlopesTOT[(ii),i])))/ 20 # divide by 20
    }
  }

SlopesTOT <- SlopesTOT[1:15,] # Remove last value. it does not contain result

# Run this function to find the max slope

maxHRTOT <- rep("NA",length(SlopesTOT)) # empty vector
for (i in 1: length (SlopesTOT)) { # for each column in the data frame
  maxHRTOT[i]<-round(max(as.numeric(SlopesTOT[,i])),3) # find the max value
}


# NIR

# Prepare subset
SlopesNIR2 <-
  TempRaw %>% 
  select (-light) %>%  # remove the column "light"
  filter (time >= 900 & time <= 1200) %>%  # keep the interval with NIR
  select (-time) %>% 
  select (sort(colnames(.)))

SlopesNIR <- SlopesNIR2 # rename 

# Run this function to find the slope for each 20 sec

for (i in 1: length (SlopesNIR)) { 
  for (ii in 2: length (SlopesNIR[,i])){ 
    SlopesNIR [ii-1,i] <- 
      (abs( (SlopesNIR[ii-1,i]) - (SlopesNIR[(ii),i])))/ 20 
    }
  }

SlopesNIR <- SlopesNIR[1:15,] # Remove last value. 

# Run this function to find the max slope

maxHRNIR <- rep("NA",length(SlopesNIR)) 
for (i in 1: length (SlopesNIR)) { 
  maxHRNIR[i]<-round(max(as.numeric(SlopesNIR[,i])),3) 
}


# VIS

# Prepare subset
SlopesVIS2 <-
  TempRaw %>% 
  select (-light) %>%  # remove the column "light"
  filter (time >= 1500 & time <= 1800) %>%  # keep the interval with VIS
  select (-time) %>% 
  select (sort(colnames(.)))

SlopesVIS <- SlopesVIS2 # rename 

# Run this function to find the slope for each 20 sec

for (i in 1: length (SlopesVIS)) { 
  for (ii in 2: length (SlopesVIS[,i])){ 
    SlopesVIS [ii-1,i] <- 
      (abs( (SlopesVIS[ii-1,i]) - (SlopesVIS[(ii),i])))/ 20 
    }
  }

SlopesVIS <- SlopesVIS[1:15,] # Remove last value. 

# Run this function to find the max slope

maxHRVIS <- rep("NA",length(SlopesVIS)) 
for (i in 1: length (SlopesVIS)) { 
  maxHRVIS[i]<-round(max(as.numeric(SlopesVIS[,i])),3) 
}


# COMBINE 

maxHR <- data.frame (
  "ind" = names(SlopesTOT2), 
  maxHRTOT = as.numeric(maxHRTOT), 
  maxHRNIR = as.numeric(maxHRNIR), 
  maxHRVIS = as.numeric(maxHRVIS))

Consolidated

HeatData <- data.frame (Other, DeltaT5, maxHR) %>% 
  select(-ind.1, -ind.2) %>% 
  select(ind, spp, genus, 
         batch, Thermocouple, side, 
         size, everything())
  
Repeatability <- 
  HeatData %>%  
  filter (batch == "repeatability" | batch == "both") %>% 
  arrange (ind)

SingleElytra <- 
  HeatData %>% 
  filter (batch == "experiment" | batch == "both") %>% 
  arrange (ind)

Beetle Body

Delta temp after 5 min

In this case, the time to cool down between treatments was double due to a higer mass

codexp <- (names(BodyRaw)[3:length(BodyRaw)])

B_t0TOT <- 
  BodyRaw %>% 
  filter (time == 280) %>%  # t0 time
  select (-light, -time) %>% 
  as.numeric (.)

B_tfTOT <- 
  BodyRaw %>% 
  filter (time == 580) %>% # end time of the first illumination 
  select (-light, -time) %>% 
  as.numeric (.)

B_t0NIR <- 
  BodyRaw %>% 
  filter (time == 1180) %>%  # t0 time
  select (-light, -time) %>% 
  as.numeric (.)

B_tfNIR <- 
  BodyRaw %>% 
  filter (time == 1480) %>% # end time of the first illumination 
  select (-light, -time) %>% 
  as.numeric (.)

B_t0VIS <- 
  BodyRaw %>% 
  filter (time == 2080) %>%  # t0 time
  select (-light, -time) %>% 
  as.numeric (.)

B_tfVIS <- 
  BodyRaw %>% 
  filter (time == 2380) %>% # end time of the first illumination 
  select (-light, -time) %>% 
  as.numeric (.)


B_TimePoints <- data.frame( codexp, 
                          B_t0TOT, B_tfTOT, 
                          B_t0NIR, B_tfNIR, 
                          B_t0VIS, B_tfVIS )

BodyDeltaT5 <- 
  B_TimePoints %>% 
  mutate (
    B_DT5_TOT = B_tfTOT - B_t0TOT,
    B_DT5_NIR = B_tfNIR - B_t0NIR,
    B_DT5_VIS = B_tfVIS - B_t0VIS ) %>% 
  select (1,8,9,10) %>% 
  arrange (codexp) %>% 
  separate (codexp, sep = 5, c("Beetle","conditions")) %>% 
  separate (conditions, sep = 1, c("Thermocouple", "ElytraPosition"))

Consolidated

Reorder the data frame to find the difference between elytra open and closed

# "e" = Thermocouple under the elytra 
# "i" = Thermocouple inside the body 
# "c" =  elytra closed
# "o" =  elytra open

# Thermocouple outside the body, right under the elytra:
DT5_ThermoExt <-
  BodyDeltaT5 %>% 
  filter (Thermocouple == "e") %>% # Keep only: thermocouple external
  select (-Thermocouple) %>%  # column not needed anymore
  gather (SpectralBand, DT5, starts_with("B_")) %>% # column for spectral band
  unite (Var1, SpectralBand, ElytraPosition) %>% # fuse spectral band + Elytra
  spread (Var1, DT5) %>%  # obtain one col for open and close at each band
  mutate (
    PairDifTOT = round (( B_DT5_TOT_c - B_DT5_TOT_o),4), 
    PairDifNIR = round (( B_DT5_NIR_c - B_DT5_NIR_o),4),
    PairDifVIS = round (( B_DT5_VIS_c - B_DT5_VIS_o),4)
  ) %>% # Find paired difference for each band. Elytra closed - open
  select (1, 8:10) # Keep only the differences


DT5_ThermoInt <-
  BodyDeltaT5 %>% 
  filter (Thermocouple == "i") %>%  # Keep only: thermocouple internal
  select (-Thermocouple) %>% 
  gather (SpectralBand, DT5, starts_with("B_")) %>% # column for spectral band
  unite (Var1, SpectralBand, ElytraPosition) %>% # fuse spectral band + Elytra
  spread (Var1, DT5) %>%  # obtain one col for open and close at each band
  mutate (
    PairDifTOT = round (( B_DT5_TOT_c - B_DT5_TOT_o),4), 
    PairDifNIR = round (( B_DT5_NIR_c - B_DT5_NIR_o),4),
    PairDifVIS = round (( B_DT5_VIS_c - B_DT5_VIS_o),4)
  ) %>% # Find paired difference for each band. Elytra closed - open
  select (1, 8:10) # Keep only the differences

BodyData <- data.frame (DT5_ThermoExt, DT5_ThermoInt)

Experimental Set up

Figure 1

setup<-Setup[Setup$time<1900,]

par(mgp=c(2.1,0.9,0))
plot(setup$time, setup$T1,
     yaxt = "n", xaxt = "n", type = "l",
     lwd = 2, col = "gray", 
     ylim = c(18,26), xlab = "Time(s)", ylab = "Temperature (°C)")
axis(2, cex.axis = 0.9)
axis(1, at = c(0,600,1200,1800), label = rep("", 4))
axis(1, at = c(0,600,1200,1800),  
     line = -0.5, lwd = 0, 
     cex.axis = 0.9)
lines(setup$time, setup$T2,  lwd = 2, col = "black")
lines(setup$time, setup$tch, lwd = 2, col = "#5BBCD6")
lines(setup$time, setup$twb, lwd = 2, col = "#00A08A")
abline(v = 300,  col = "gray", lty = 3)
abline(v = 600,  col = "gray", lty = 3)
abline(v = 900,  col = "gray", lty = 3)
abline(v = 1200, col = "gray", lty = 3)
abline(v = 1500, col = "gray", lty = 3)
abline(v = 1800, col = "gray", lty = 3)

Single Elytra

First, fuse the heating and spectral data

Single <- 
  data.frame ( Optical, # Optical properties
               SingleElytra ) %>%  # Heating Data
  rename (Beetle = spp) %>%  
  select (1:10, 14:16, 21, 22:29)  # Keep only the relevant columns

Stats

Linear Models

Heat ~ Absorptivity + Size + Side

options("jtools-digits" = 4)
set_summ_defaults(digits = 4)

Delta T5

# TOTAL
mod1B <- lm(DT5_TOT ~ Ab_TOT + size + side, data = Single)

mod1BPer <- c (
  anova(mod1B)[1,2]/sum(anova(mod1B)[,2])*100, 
  anova(mod1B)[2,2]/sum(anova(mod1B)[,2])*100, 
  anova(mod1B)[3,2]/sum(anova(mod1B)[,2])*100)

# NIR
mod2B <- lm(DT5_NIR ~ Ab_NIR + size + side, data = Single)

mod2BPer <- c ( 
  anova(mod2B)[1,2]/sum(anova(mod2B)[,2])*100, 
  anova(mod2B)[2,2]/sum(anova(mod2B)[,2])*100, 
  anova(mod2B)[3,2]/sum(anova(mod2B)[,2])*100)

# VIS
mod3B <- lm(DT5_VIS ~ Ab_VIS + size + side, data = Single)

mod3BPer <- c ( 
  anova(mod3B)[1,2]/sum(anova(mod3B)[,2])*100,
  anova(mod3B)[2,2]/sum(anova(mod3B)[,2])*100,
  anova(mod3B)[3,2]/sum(anova(mod3B)[,2])*100)


export_summs(mod1B, mod2B, mod3B, 
             error_format = "SE = {std.error} p = {p.value}",
             model.names = c("TOT", "NIR", "VIS"))

	TOT	NIR	VIS
(Intercept)	2.0240 ***	1.1302 ***	0.2011
	SE = 0.3389 p = 0.0000	SE = 0.1922 p = 0.0000	SE = 0.2350 p = 0.3960
Ab_TOT	0.0284 ***
	SE = 0.0043 p = 0.0000
size	0.2737 *	0.0688	0.2219 ***
	SE = 0.1238 p = 0.0315	SE = 0.0833 p = 0.4126	SE = 0.0530 p = 0.0001
sideB	-0.1451	-0.0461	-0.1215 **
	SE = 0.0994 p = 0.1504	SE = 0.0660 p = 0.4885	SE = 0.0433 p = 0.0070
Ab_NIR		0.0101 ***
		SE = 0.0023 p = 0.0001
Ab_VIS			0.0078 **
			SE = 0.0025 p = 0.0033
N	56	56	56
R2	0.5365	0.3217	0.4343
* p < 0.001; p < 0.01; * p < 0.05.

data.frame( "Var" = c("OpticalProp", "Size", "Side"),
            "TOT" = mod1BPer,
            "NIR" = mod2BPer,
            "VIS" = mod3BPer)  %>% 
  kbl() %>%
  kable_paper("hover", full_width = F)

Var	TOT	NIR	VIS
OpticalProp	47.081185	30.5720888	14.08036
Size	4.670822	0.9654150	20.75990
Side	1.899194	0.6347567	8.58811

maxHR

Heating Rates : degrees per second (max slope)

# TOTAL
mod1C<-lm(maxHRTOT ~ Ab_TOT + size + side, data = Single)

mod1CPer <- c (
anova(mod1C)[1,2]/sum(anova(mod1C)[,2])*100,
anova(mod1C)[2,2]/sum(anova(mod1C)[,2])*100,
anova(mod1C)[3,2]/sum(anova(mod1C)[,2])*100)

# NIR
mod2C<-lm(maxHRNIR ~ Ab_NIR + size + side, data = Single)

mod2CPer <- c (
anova(mod2C)[1,2]/sum(anova(mod2C)[,2])*100,
anova(mod2C)[2,2]/sum(anova(mod2C)[,2])*100,
anova(mod2C)[3,2]/sum(anova(mod2C)[,2])*100)

# VIS
mod3C<-lm(maxHRVIS ~ Ab_VIS + size + side, data = Single)

mod3CPer <- c (
anova(mod3C)[1,2]/sum(anova(mod3C)[,2])*100,
anova(mod3C)[2,2]/sum(anova(mod3C)[,2])*100,
anova(mod3C)[3,2]/sum(anova(mod3C)[,2])*100)


export_summs(mod1C, mod2C, mod3C, 
             error_format = "SE = {std.error} p = {p.value} ", 
             model.names = c("TOT", "NIR", "VIS"))

	TOT	NIR	VIS
(Intercept)	0.0322 *	0.0263 ***	0.0026
	SE = 0.0125 p = 0.0127	SE = 0.0043 p = 0.0000	SE = 0.0097 p = 0.7884
Ab_TOT	0.0006 ***
	SE = 0.0002 p = 0.0006
size	0.0036	-0.0035	0.0007
	SE = 0.0046 p = 0.4321	SE = 0.0019 p = 0.0675	SE = 0.0022 p = 0.7630
sideB	-0.0079 *	-0.0025	-0.0052 **
	SE = 0.0037 p = 0.0355	SE = 0.0015 p = 0.0966	SE = 0.0018 p = 0.0054
Ab_NIR		0.0002 ***
		SE = 0.0001 p = 0.0001
Ab_VIS			0.0003 *
			SE = 0.0001 p = 0.0121
N	56	56	56
R2	0.2991	0.3033	0.2427
* p < 0.001; p < 0.01; * p < 0.05.

data.frame( "Var" = c("OpticalProp", "Size", "Side"),
            "TOT" = mod1CPer,
            "NIR" = mod2CPer,
            "VIS" = mod3CPer)  %>% 
  kbl() %>%
  kable_paper("hover", full_width = F)

Var	TOT	NIR	VIS
OpticalProp	22.527798	22.214424	11.6158168
Size	1.102470	4.283228	0.3434188
Side	6.283053	3.837064	12.3083013

Correlations

Involving Optical properties

What is the correlation between Trasmissivity and Reflectivity on each spectral band?

Cor1TOT <- cor.test (Single$Rf_TOT, Single$Td_TOT) # Total

Cor1NIR <- cor.test (Single$Rf_NIR, Single$Td_NIR) # NIR

Cor1VIS <- cor.test (Single$Rf_VIS, Single$Td_VIS) #VIS


# Extract parameters
Cor1Est <- round (
  c(parse_number(toString(Cor1TOT[4])), 
    parse_number(toString(Cor1NIR[4])), 
    parse_number(toString(Cor1VIS[4]))), 4)

Cor1pval <- round (
  c(parse_number(toString(Cor1TOT[3])), 
    parse_number(toString(Cor1NIR[3])), 
    parse_number(toString(Cor1VIS[3]))), 4)

Cor1CInum <-  round (
  c((Cor1TOT[9]$conf.int[1:2]),
    (Cor1NIR[9]$conf.int[1:2]),
    (Cor1VIS[9]$conf.int[1:2])), 4)

Cor1CI <- c(
  paste(Cor1CInum[1], Cor1CInum[2], sep = " : " ),
  paste(Cor1CInum[3], Cor1CInum[4], sep = " : " ),
  paste(Cor1CInum[5], Cor1CInum[6], sep = " : " )
)            


data.frame (
  "SpectralBand" = c("TOT", "NIR", "VIS"),
  "PearsonR" = Cor1Est,
  "ConfInterval" = Cor1CI ,
  "pVal" = Cor1pval
   ) %>% 
  kbl() %>%
  kable_paper("hover", full_width = F)

SpectralBand	PearsonR	ConfInterval	pVal
TOT	-0.3253	-0.5419 : -0.0682	0.0144
NIR	-0.4700	-0.6523 : -0.2363	0.0003
VIS	-0.0272	-0.288 : 0.2374	0.8422

Involving Spectral Bands

Cor2Reflec <- cor.test (Single$Rf_NIR, Single$Rf_VIS) # Reflectivity

Cor2Transm <- cor.test (Single$Td_NIR, Single$Td_VIS) # Transmissivity

Cor2Absorp <- cor.test (Single$Ab_NIR, Single$Ab_VIS) # Absorptivity


# Extract parameters
Cor2Est <- round (
  c(parse_number(toString(Cor2Reflec[4])), 
    parse_number(toString(Cor2Transm[4])), 
    parse_number(toString(Cor2Absorp[4]))), 4)

Cor2pval <- round (
  c(parse_number(toString(Cor2Reflec[3])), 
    parse_number(toString(Cor2Transm[3])), 
    parse_number(toString(Cor2Absorp[3]))), 4)

Cor2CInum <-  round (
  c((Cor2Reflec[9]$conf.int[1:2]),
    (Cor2Transm[9]$conf.int[1:2]),
    (Cor2Absorp[9]$conf.int[1:2])), 4)

Cor2CI <- c(
  paste(Cor2CInum[1], Cor2CInum[2], sep = " : " ),
  paste(Cor2CInum[3], Cor2CInum[4], sep = " : " ),
  paste(Cor2CInum[5], Cor2CInum[6], sep = " : " )
)            


data.frame (
  "OpticalProp" = c("Reflectivity", "Transmissivity", "Absorptivity"),
  "PearsonR" = Cor2Est,
  "ConfInterval" = Cor2CI ,
  "pVal" = Cor2pval
  ) %>% 
  kbl() %>%
  kable_paper("hover", full_width = F)

OpticalProp	PearsonR	ConfInterval
Reflectivity	0.5851	0.3808 : 0.735
Transmissivity	0.7502	0.607 : 0.8462
Absorptivity	0.7854	0.6585 : 0.8689

Involving Size

Reflectivity

Cor3RfTOT <- cor.test (Single$Rf_TOT, Single$size) 
Cor3RfNIR <- cor.test (Single$Rf_NIR, Single$size) 
Cor3RfVIS <- cor.test (Single$Rf_VIS, Single$size) 

# Extract parameters
Cor3Est <- round (
  c(parse_number(toString(Cor3RfTOT[4])), 
    parse_number(toString(Cor3RfNIR[4])), 
    parse_number(toString(Cor3RfVIS[4]))), 4)

Cor3pval <- round (
  c(parse_number(toString(Cor3RfTOT[3])), 
    parse_number(toString(Cor3RfNIR[3])), 
    parse_number(toString(Cor3RfVIS[3]))), 4)

Cor3CInum <-  round (
  c((Cor3RfTOT[9]$conf.int[1:2]),
    (Cor3RfNIR[9]$conf.int[1:2]),
    (Cor3RfVIS[9]$conf.int[1:2])), 4)

Cor3CI <- c(
  paste(Cor3CInum[1], Cor3CInum[2], sep = " : " ),
  paste(Cor3CInum[3], Cor3CInum[4], sep = " : " ),
  paste(Cor3CInum[5], Cor3CInum[6], sep = " : " )
)            

data.frame (
  "SpectralBand" = c("TOT", "NIR", "VIS"),
  "PearsonR" = Cor3Est,
  "ConfInterval" = Cor3CI ,
  "pVal" = Cor3pval
   ) %>% 
  kbl() %>%
  kable_paper("hover", full_width = F)

SpectralBand	PearsonR	ConfInterval	pVal
TOT	0.1243	-0.1433 : 0.375	0.3613
NIR	0.1301	-0.1375 : 0.38	0.3393
VIS	0.0618	-0.2044 : 0.3195	0.6507

Transmissivity

Cor4TdTOT <- cor.test (Single$Td_TOT, Single$size) 
Cor4TdNIR <- cor.test (Single$Td_NIR, Single$size) 
Cor4TdVIS <- cor.test (Single$Td_VIS, Single$size) 

# Extract parameters
Cor4Est <- round (
  c(parse_number(toString(Cor4TdTOT[4])), 
    parse_number(toString(Cor4TdNIR[4])), 
    parse_number(toString(Cor4TdVIS[4]))), 4)

Cor4pval <- round (
  c(parse_number(toString(Cor4TdTOT[3])), 
    parse_number(toString(Cor4TdNIR[3])), 
    parse_number(toString(Cor4TdVIS[3]))), 4)

Cor4CInum <-  round (
  c((Cor4TdTOT[9]$conf.int[1:2]),
    (Cor4TdNIR[9]$conf.int[1:2]),
    (Cor4TdVIS[9]$conf.int[1:2])), 4)

Cor4CI <- c(
  paste(Cor4CInum[1], Cor4CInum[2], sep = " : " ),
  paste(Cor4CInum[3], Cor4CInum[4], sep = " : " ),
  paste(Cor4CInum[5], Cor4CInum[6], sep = " : " )
)            

data.frame (
  "SpectralBand" = c("TOT", "NIR", "VIS"),
  "PearsonR" = Cor4Est,
  "ConfInterval" = Cor4CI ,
  "pVal" = Cor4pval
   ) %>% 
  kbl() %>%
  kable_paper("hover", full_width = F)

SpectralBand	PearsonR	ConfInterval	pVal
TOT	-0.3429	-0.5557 : -0.0879	0.0097
NIR	-0.3661	-0.5738 : -0.1142	0.0055
VIS	-0.2041	-0.4432 : 0.0622	0.1314

Absorptivity

Cor5AbTOT <- cor.test (Single$Ab_TOT, Single$size) 
Cor5AbNIR <- cor.test (Single$Ab_NIR, Single$size) 
Cor5AbVIS <- cor.test (Single$Ab_VIS, Single$size) 

# Extract parameters
Cor5Est <- round (
  c(parse_number(toString(Cor5AbTOT[4])), 
    parse_number(toString(Cor5AbNIR[4])), 
    parse_number(toString(Cor5AbVIS[4]))), 4)

Cor5pval <- round (
  c(parse_number(toString(Cor5AbTOT[3])), 
    parse_number(toString(Cor5AbNIR[3])), 
    parse_number(toString(Cor5AbVIS[3]))), 4)

Cor5CInum <-  round (
  c((Cor5AbTOT[9]$conf.int[1:2]),
    (Cor5AbNIR[9]$conf.int[1:2]),
    (Cor5AbVIS[9]$conf.int[1:2])), 4)

Cor5CI <- c(
  paste(Cor5CInum[1], Cor5CInum[2], sep = " : " ),
  paste(Cor5CInum[3], Cor5CInum[4], sep = " : " ),
  paste(Cor5CInum[5], Cor5CInum[6], sep = " : " )
)            

data.frame (
  "SpectralBand" = c("TOT", "NIR", "VIS"),
  "PearsonR" = Cor5Est,
  "ConfInterval" = Cor5CI ,
  "pVal" = Cor5pval
   ) %>% 
  kbl() %>%
  kable_paper("hover", full_width = F)

SpectralBand	PearsonR	ConfInterval	pVal
TOT	0.1975	-0.0689 : 0.4377	0.1445
NIR	0.2497	-0.0141 : 0.481	0.0634
VIS	0.0760	-0.1907 : 0.3322	0.5779

Figure 2

Data wrangling

Figure 2 and 3 need data frames in a special order:

Stacked1 <- 
  Single %>% 
  select (1:7) %>% 
  gather (set, Percentage, (2:7)) %>% # combine columns
  mutate (band = substr(set,4,6),
          optic = substr(set,1,2)) %>% # create columns for band and properties
  select (-set) %>% 
  spread (band,Percentage) %>%  # separate according to band 
  arrange (optic, Beetle)

# TOT index
indexT <- 
  Stacked1 %>% 
  filter (optic == "Rf") %>% 
  arrange (TOT) %>% 
  mutate (TOT_index = seq (1,56,1)) %>% 
  arrange (Beetle)

Stacked1$TOT_factor <- rep (indexT$TOT_index,2)

# NIR Index
indexN <- 
  Stacked1 %>% 
  filter (optic == "Rf") %>% 
  arrange (NIR) %>% 
  mutate (NIR_index = seq (1,56,1)) %>% 
  arrange (Beetle)

Stacked1$NIR_factor <- rep (indexN$NIR_index,2)


# VIS Index
indexV <- 
  Stacked1 %>% 
  filter (optic == "Rf") %>% 
  arrange (VIS) %>% 
  mutate (VIS_index = seq (1,56,1)) %>% 
  arrange (Beetle)

Stacked1$VIS_factor <- rep (indexV$VIS_index,2)

# Finally, also include absorptivity in TOT for figure 2

Stacked1$AbsTOT <- rep(Single$Ab_TOT,2)

Plot

Single <-
  Single %>% 
  mutate (genus = SingleElytra$genus)

## Plot Delta Temp 5 min ~ Absorptivity
dummy <- lm(DT5_TOT ~ Ab_TOT, data = Single)

ggplot(Single,aes(Ab_TOT, DT5_TOT, label = Beetle, group = genus))+
  geom_point(aes(fill=genus),alpha=0.5,cex=2,shape=21)+
  theme_light() +
  theme(axis.text.x = element_text(size = 15))+
  theme(axis.text.y = element_text(size = 15))+
  theme(axis.title.x = element_text(size = 17))+
  theme(axis.title.y = element_text(size = 17))+
  theme(axis.title.y = element_text(margin = 
                                      margin(t = 0, r = 20, b = 0, l = 0)))+
  labs(x="Absorptivity (%)",y="Delta T [5] ")+
  geom_text(size=3, hjust=-0.05, vjust=-0.4,col="gray")+
  scale_fill_manual(values = wes_palette("Darjeeling1", 
                                         n = 8,type = "continuous"))+
  scale_x_continuous(position="top", limits = c(20,80), expand = c(0, 0))+
  scale_y_continuous(limits=c(3, 5.2), expand = c(0.05, 0.05)) +
  theme(legend.position = "none")+
  geom_abline(intercept = dummy$coefficients[1], 
              slope = dummy$coefficients[2], col="black",size=0.6)

# Plot Delta Temp 5 min ~ Size
dummy2 <- lm(DT5_TOT ~ size , data = Single)

ggplot(Single,aes(size, DT5_TOT, label = Beetle, group = genus))+
  geom_point(aes(fill = genus), alpha = 0.5, cex = 2, shape = 21)+
  theme_light() +
  theme(axis.text.x = element_text(size = 15))+
  theme(axis.text.y = element_text(size = 15))+
  theme(axis.title.x = element_text(size = 17))+
  theme(axis.title.y = element_text(size = 17))+
  theme(axis.title.y = element_text(margin = 
                                      margin(t = 0, r = 20, b = 0, l = 0)))+
  labs(x="Size (cm)",y="Delta T [5]")+
  geom_text(size=3, hjust=-0.05, vjust=-0.4,col="gray")+
  scale_fill_manual(values = wes_palette("Darjeeling1", 
                                         n = 8,type = "continuous"))+
  scale_x_continuous(position="top", expand = c(0.1, 0.1))+
  scale_y_continuous(limits=c(3, 5.2), expand = c(0.05, 0.05)) +
  geom_abline(intercept = dummy2$coefficients[1], 
              slope = dummy2$coefficients[2], col="black",size=0.6)+
  theme(legend.position="bottom")+
  theme(legend.title=element_blank())+
  theme(legend.background = element_rect(colour="black"))+
  theme(legend.key.height = unit(0.4, 'cm'),
        legend.key.width= unit(0.6, 'cm'))+
  guides(fill=guide_legend(nrow=3,byrow=TRUE))+
  theme(legend.text = element_text(size=9))

# Barplot 

ggplot(Stacked1, aes(AbsTOT, y = TOT))+
  geom_col(aes(fill = optic), width = 0.5, 
           position = position_stack(reverse=TRUE))+
  theme_light() +
  scale_x_continuous(limits=c(20, 80), expand = c(0, 0)) +
  scale_y_continuous(limits=c(0, 100), expand = c(0, 0)) +
  theme(legend.position= c(0.85,0.2))+
  theme(legend.title=element_blank())+
  scale_y_reverse()+
  theme(axis.title.y = element_text(margin = 
                                      margin(t = 0, r = 10, b = 0, l = 0)))+
  theme(axis.text.y = element_text(size = 14))+
  theme(axis.title.x = element_text(margin = 
                                      margin(t = 10, r = 0, b = 0, l = 0)))+
  scale_fill_manual(values = c("black","gray"))+
  labs(x=" ",y=" Non-absorbed Light (%) ")+
  theme(legend.key.height = unit(0.3, 'cm'),
        legend.key.width= unit(0.3, 'cm'))

Figure 3

Barplots

# NIR

bp2<-ggplot(Stacked1, aes(x = reorder (Beetle, NIR_factor, FUN = sum ),
                          y = NIR))+
  geom_col(aes(fill = optic), width = 0.7,
           position = position_stack(reverse=TRUE))+
  coord_flip()+
  theme_minimal()+
  ylim(0,100)+
  scale_fill_manual(values = c("black","gray"))+
  theme(legend.position="none")+
  labs(x=" ",y=" NIR Light (%)")+
  theme(axis.text.y = element_text(size=7))+
  theme(axis.text.x = element_text(size = 10))+
  theme(axis.title.x = element_text(margin = 
                                    margin(t = 10, r = 0, b = 0, l = 0)))

# VIS

bp3<-ggplot(Stacked1, aes(x = reorder (Beetle, VIS_factor, FUN = sum ), 
                          y = VIS))+
  geom_col(aes(fill = optic), width = 0.7,
           position = position_stack(reverse=TRUE))+
  coord_flip()+
  theme_minimal()+
  ylim(0,50)+
  scale_fill_manual(values = c("black","gray"))+
  theme(legend.position="none")+
  labs(x=" ",y=" VIS Light (%)")+
  theme(axis.text.y = element_text(size=7))+
  theme(axis.text.x = element_text(size = 10))+
  theme(axis.title.x = element_text(margin = 
                                      margin(t = 10, r = 0, b = 0, l = 0)))

# Join them

grid.arrange(bp2,bp3,ncol=2)

Examples of the diversity in optical properties

PlotT1B<-ggplot() +
  geom_ribbon(data=Reflect,
              mapping = aes(x = wl, y = 100 - anom2,
                            ymin=TopTrans$anom2t, ymax=100-anom2), 
              fill="#b6ada6", alpha=0.7)+
  
  geom_ribbon(data=Reflect,
              mapping = aes(x = wl, y = 100 - anom2,
                            ymin=100,ymax=100-anom2), 
              colour="black", fill="#ffd781", alpha=0.3)+
  
  geom_area(data = TopTrans,
            aes(x = wl, y = anom2t),
            fill="#00a08a", colour="black", alpha=0.5) +
  
  scale_x_continuous(limits=c(400, 1700),expand = c(0, 0))+
  scale_y_continuous(limits=c(0, 100),breaks = c(50,100),
                     expand = c(0, 0),
                     sec.axis = sec_axis(~100-.,breaks = c(100,50))) +
  theme_classic()+
  ylab( " " )+
  theme(plot.margin = unit(c(0.1,0.1,0.05,0.1), "cm"))+
  theme(axis.text.x = element_blank())+
  theme(axis.title.x = element_blank())


PlotT2B<-ggplot() +
  
  geom_ribbon(data=Reflect,
              mapping = aes(x=wl,y=100-fchi1,
                            ymin=TopTrans$fchi1t,ymax=100-fchi1),
              fill="#b6ada6", alpha=0.7)+
  
  geom_ribbon(data=Reflect,
              mapping = aes(x=wl,y=100-fchi1,
                            ymin=100,ymax=100-fchi1), 
              colour="black", fill="#ffd781", alpha=0.3)+
  
  geom_area(data=TopTrans,
            aes(x=wl, y=fchi1t),
            fill="#00a08a", colour="black", alpha=0.5) +
  
  scale_x_continuous(limits=c(400, 1700),
                     expand = c(0, 0))+
  scale_y_continuous(limits=c(0, 100),breaks = c(50,100), 
                     expand = c(0, 0),
                     sec.axis = sec_axis(~100-.,breaks = c(100,50))) +
  theme_classic()+
  ylab( " " )+
  theme(plot.margin = unit(c(0.1,0.1,0.05,0.1), "cm"))+
  theme(axis.text.x = element_blank())+
  theme(axis.title.x = element_blank())


PlotT3B<-ggplot() +
  
  geom_ribbon(data=Reflect,
              mapping =aes(x=wl,y=100-ocul2,
                           ymin=TopTrans$ocul2t,ymax=100-ocul2),
              fill="#b6ada6", alpha=0.7)+
  
  geom_ribbon(data=Reflect,
              mapping = aes(x=wl,y=100-ocul2,
                            ymin=100,ymax=100-ocul2), 
              colour="black", fill="#ffd781", alpha=0.3)+
  
  geom_area(data=TopTrans, 
            aes(x=wl, y=ocul2t),
            fill="#00a08a", colour="black", alpha=0.5) +
  
  scale_x_continuous(limits=c(400, 1700),
                     expand = c(0, 0))+
  scale_y_continuous(limits=c(0, 100),breaks = c(50,100), 
                     expand = c(0, 0),
                     sec.axis = sec_axis(~100-.,breaks = c(100,50))) +
  theme_classic()+
  ylab( " " )+
  theme(plot.margin = unit(c(0.1,0.1,0.05,0.1), "cm"))+
  theme(axis.text.x = element_blank())+
  theme(axis.title.x = element_blank())

PlotT4B<-ggplot() +
  geom_ribbon(data=Reflect,
              mapping = aes(x=wl,y=100-aurs1,
                            ymin=TopTrans$aurs1t,ymax=100-aurs1),
              fill="#b6ada6", alpha=0.7)+
  
  geom_ribbon(data=Reflect,
              mapping = aes(x=wl,y=100-aurs1,
                            ymin=100,ymax=100-aurs1), 
              colour="black", fill="#ffd781", alpha=0.3)+
  
  geom_area(data=TopTrans,
            aes(x=wl, y=aurs1t),
            fill="#00a08a", colour="black", alpha=0.5) +
  
  scale_x_continuous(limits=c(400, 1700),
                     expand = c(0, 0))+
  scale_y_continuous(limits=c(0, 100), breaks = c(50,100), 
                     expand = c(0, 0),
                     sec.axis = sec_axis(~100-.,breaks = c(100,50))) +
  theme_classic()+
  theme(axis.title.x = element_blank())+
  ylab( " " )+
  theme(axis.text.x = element_blank())+
  theme(plot.margin = unit(c(0.1,0.1,0.05,0.1), "cm"))

PlotT5B<-ggplot() +
  geom_ribbon(data=Reflect,
              mapping = aes(x=wl,y=100-clor1,
                            ymin=TopTrans$clor1t,ymax=100-clor1), 
              fill="#b6ada6", alpha=0.7)+
  
  geom_ribbon(data=Reflect,
              mapping = aes(x=wl,y=100-clor1,
                            ymin=100,ymax=100-clor1), 
              colour="black", fill="#ffd781", alpha=0.3)+
  
  geom_area(data=TopTrans, aes(x=wl, y=clor1t),
            fill="#00a08a", colour="black", alpha=0.5) +
  
  scale_x_continuous(limits=c(400, 1700),
                     expand = c(0, 0))+
  scale_y_continuous(limits=c(0, 100), breaks = c(50,100), 
                     expand = c(0, 0),
                     sec.axis = sec_axis(~100-.,breaks = c(100,50))) +
  theme_classic()+
  theme(axis.title.x = element_blank())+
  ylab( " " )+
  theme(axis.text.x = element_blank())+
  theme(plot.margin = unit(c(0.1,0.1,0.04,0.01), "cm"))


plot_grid(PlotT4B, PlotT3B, PlotT5B, PlotT1B, PlotT2B,
          align = "v", nrow=5, ncol=1, rel_widths = c(1, 1))

Beetle Body

Stats

We used bootstrap due to limited n.

Body

Thermocouple placed inside the body of the beetle

DT5_ThermoInt

Beetle	PairDifTOT	PairDifNIR	PairDifVIS
RDB01	0.2	-0.1	0.3
RDB02	-0.4	-0.2	0
RDB03	-0.4	-0.2	0.1
RDB04	0	0	0
RDB05	-0.3	-0.2	0.1
RDB06	-0.4	-0.3	0.1
RDB07	-0.4	-0.4	-0.1
RDB08	0.1	-0.2	0.1
RDB09	-0.4	-0.1	0.1
RDB10	-0.4	-0.2	0.1
RDB11	-0.9	-0.6	-0.2

# TOT
out1<-rep(NA, 1000) 
for (ii in 1:1000){ 
  samp1<-sample(DT5_ThermoInt$PairDifTOT, 
                size = length(DT5_ThermoInt$PairDifTOT), 
                replace = TRUE) 
  out1[ii]<-mean(samp1) 
}

out1.sorted <- sort(out1)
ci1 <- round(
    quantile(out1.sorted,c(0.025,0.975), type=7), # Confidence interval
    4) 
ci1r <- paste(ci1[1],ci1[2],sep = " : ") 


# NIR
out2<-rep(NA, 1000) 
for (ii in 1:1000){ 
  samp2<-sample(DT5_ThermoInt$PairDifNIR, 
                size = length(DT5_ThermoInt$PairDifNIR), 
                replace = TRUE) 
  out2[ii]<-mean(samp2) 
}

out2.sorted <- sort(out2)
ci2 <- round(
    quantile(out2.sorted,c(0.025,0.975), type=7),# Confidence interval
    4) 
ci2r <- paste(ci2[1],ci2[2],sep = " : ") 

# VIS
out3<-rep(NA, 1000) 
for (ii in 1:1000){ 
  samp3<-sample(DT5_ThermoInt$PairDifVIS, 
                size = length(DT5_ThermoInt$PairDifVIS), 
                replace = TRUE) 
  out3[ii]<-mean(samp3) 
}

out3.sorted <- sort(out3)
ci3 <- round(
    quantile(out3.sorted,c(0.025,0.975), type=7), # Confidence interval
    4)  
ci3r <- paste(ci3[1],ci3[2],sep = " : ") 

data.frame(
  "Band" = c("TOT", "NIR", "VIS"),
  "Mean" = round(c(mean(out1), mean(out2), mean(out3)),4),
  "C.I." = as.character(c(ci1r, ci2r, ci3r))
    ) %>% 
  kbl() %>%
  kable_paper("hover", full_width = F)

Band	Mean	C.I.
TOT	-0.2991	-0.4727 : -0.1361
NIR	-0.2269	-0.3273 : -0.1452
VIS	0.0546	-0.0091 : 0.1273

Under the elytra

Thermocouple placed under the elytra

Control

DT5_ThermoExt

Beetle	PairDifTOT	PairDifNIR	PairDifVIS
RDB01	0.2	0	0.3
RDB02	-0.1	-0.3	-0.2
RDB03	-0.7	-0.4	-0.2
RDB04	0.4	0.2	0.1
RDB05	0.4	-0.1	0.3
RDB06	0.4	0.2	0.1
RDB07	0.1	0	-0.1
RDB08	0.4	-0.2	0
RDB09	-0.1	0.1	-0.1
RDB10	-0.1	0	0
RDB11	-0.7	-0.5	-0.3

# TOT
out5<-rep(NA, 1000) 
for (ii in 1:1000){ 
  samp5<-sample(DT5_ThermoExt$PairDifTOT, 
                size = length(DT5_ThermoExt$PairDifTOT), 
                replace = TRUE) 
  out5[ii]<-mean(samp5) 
}

out5.sorted<-sort(out5)
ci5 <- round(
    quantile(out5.sorted,c(0.025,0.975), type=7), # Confidence interval
    4) 
ci5r <- paste(ci5[1],ci5[2],sep = " : ") 

# NIR
out6<-rep(NA, 1000) 
for (ii in 1:1000){ 
  samp6<-sample(DT5_ThermoExt$PairDifNIR, 
                size = length(DT5_ThermoExt$PairDifNIR), 
                replace = TRUE) 
  out6[ii]<-mean(samp6) 
}

out6.sorted<-sort(out6)
ci6 <- round(
    quantile(out6.sorted,c(0.025,0.975), type=7), # Confidence interval
    4) 
ci6r <- paste(ci6[1],ci6[2],sep = " : ") 


# VIS
out7<-rep(NA, 1000) 
for (ii in 1:1000){ 
  samp7<-sample(DT5_ThermoExt$PairDifVIS, 
                size = length(DT5_ThermoExt$PairDifVIS), 
                replace = TRUE) 
  out7[ii]<-mean(samp7) 
}

out7.sorted<-sort(out7)
ci7 <- round(
    quantile(out7.sorted,c(0.025,0.975), type=7), # Confidence interval
    4) 
ci7r <- paste(ci7[1],ci7[2],sep = " : ") 

data.frame(
  "Band" = c("TOT", "NIR", "VIS"),
  "Mean" = round(c(mean(out1), mean(out2), mean(out3)),4),
  "C.I." = as.character(c(ci5r, ci6r, ci7r))
    ) %>% 
  kbl() %>%
  kable_paper("hover", full_width = F)

Band	Mean	C.I.
TOT	-0.2991	-0.2184 : 0.2273
NIR	-0.2269	-0.2364 : 0.0273
VIS	0.0546	-0.1182 : 0.1

Figure 4

# Prepare Data Sets

BodyDiff <-
  DT5_ThermoInt %>% 
  gather (key = B, Difference, -Beetle) %>% 
  mutate (Band = substr(B,8,10)) %>% 
  select (-B) %>% 
  mutate (Band = factor (Band, 
                         levels = c("TOT", "NIR", "VIS")) )
  

BodyBoots <- 
  data.frame (out1, out2, out3) %>% # bootstrap results 
  rename (TOT = out1,
          NIR = out2,
          VIS = out3) %>% 
  gather (key = Band, Difference) %>% 
  mutate (Band = factor (Band, 
                         levels = c("TOT", "NIR", "VIS")) )
  


# Package for half violins
source("https://raw.githubusercontent.com/datavizpyr/data/master/half_flat_violinplot.R")  

# plot

ggplot(NULL, aes(x = Band, Difference)) +
  geom_flat_violin(data = BodyBoots, 
                   fill="#f80000", alpha= 0.5) +
  geom_dotplot(data = BodyDiff,
               binaxis='y', stackdir='center', dotsize=0.5,
               colour="black",fill="black",alpha=0.5)+
  ylim(-1,0.3)+
  theme_bw()+
  geom_hline(yintercept=0, linetype="dashed", color = "#02401b")+
  theme(axis.title.y = 
          element_text(size=17,
                       margin = margin(t = 0, r = 10, b = 0, l = 0)))+
  theme(axis.title.x =
          element_text(size=15,
                       margin = margin(t = 15, r = 0, b = 0, l = 0)))+
  theme(axis.text.x=element_text(face="bold",size=14))+
  theme(axis.text.y=element_text(size=14, angle=90, hjust=0.6))+
  labs(y = paste ("Difference ΔT (°C/5 min)"), x = "Spectral Band")

Supplementary Materials

Light Source

Figure S2

FiltersN <- 
  filters %>% 
  select (wl, nirfilterTransmittance, visfilterTransmittance) %>% 
  gather (key = filters, value = Transm, -wl) %>%  
  mutate (Band = substr(filters, 1, 3 ) ) %>% 
  select (-filters)


plotS1D <- ggplot() +
  geom_area(data= theorySunALL[theorySunALL$wl<400,],
            mapping = aes(x=wl, y=gobalIrrad),
            fill="#8D53FC", alpha = 0.3) +
  geom_area(data= theorySunALL[theorySunALL$wl>400 & 
                                 theorySunALL$wl<700,],
            mapping = aes(x=wl, y=gobalIrrad),
            fill="#00A08A", alpha = 0.3)+
  geom_area(data= theorySunALL[theorySunALL$wl>700 &
                               theorySunALL$wl<1700,],
            mapping = aes(x=wl, y=gobalIrrad),
            fill="#FF0000", alpha = 0.3) +
  geom_line(data=theorySunALL,
            mapping = aes(x=wl, y=gobalIrrad))+
  theme_bw()+
  theme(axis.title.x = element_blank())+
  ylab(expression("   Irradiance \n (W/m^2/nm)"))+ 
  theme(axis.title.y = element_text(size=7),
        axis.text.y = element_text(size=7),
        axis.text.x = element_text(size=7),
        axis.title.x = element_blank()) +
  theme(plot.margin = unit(c(0.4, 0.4, 0, 0.4), "cm"))+
  xlim(250,3000)

plotS1A <- ggplot() +
  geom_area(aes(x=filters$wl, y=filters$halfsun),
            fill="#F98400", alpha = 0.3) +
  geom_line(aes(x=filters$wl, y=filters$halfsun),
            colour="#E6B200") +
  geom_line(data=theorySun,mapping = aes(x=wl,
                                         y=irradiance))+
  theme_bw()+
  theme(axis.ticks.x = element_blank(),
        axis.text.x = element_blank(),
        axis.title.x = element_blank())+
  ylab(expression("   Irradiance \n (W/m^2/nm)"))+
  theme(axis.title.y = element_text(size=7),
        axis.text.y = element_text(size=7)) +
  scale_x_continuous(breaks=seq(300,1700,200)) +
  theme(plot.margin = unit(c(0.4, 0.4, 0, 0.4), "cm"))

plotS1B <- ggplot() +
  geom_area(aes(x=FiltersN$wl[FiltersN$Band == "vis"], 
                y=FiltersN$Transm[FiltersN$Band == "vis"]),
            fill="#00A08A", alpha = 0.3) +
  geom_area(aes(x=FiltersN$wl[FiltersN$Band =="nir"], 
                y=FiltersN$Transm[FiltersN$Band =="nir"]),
            fill="#FF0000", alpha = 0.3)+
  theme_bw()+
  xlab(" ")+
  ylab("Transmittance (%)")+
  theme(axis.title.y = element_text(size=7),
        axis.text.y = element_text(size=7)) +
  theme(axis.ticks.x = element_blank(),
        axis.text.x = element_blank(),
        axis.title.x = element_blank())+
  scale_x_continuous(breaks=seq(300,1700,200)) +
  theme(plot.margin = unit(c(0.2, 0.4, 0.3, 0.4), "cm"))

plotS1C <- ggplot() +
  geom_line(data=filters,mapping = aes(x=wl,
                                       y=Full.f),
            colour="#E6B200")+
  geom_area(aes(x=filters$wl, y=filters$vis.f),
            fill="#00A08A", alpha = 0.3) +
  geom_line(aes(x=filters$wl, y=filters$vis.f),
            colour="#00A08A", alpha = 0.9) +
  geom_area(aes(x=filters$wl, y=filters$nir.f),
            fill="#FF0000", alpha = 0.3)+
  geom_line(aes(x=filters$wl, y=filters$nir.f),
            colour="#FF0000", alpha = 0.9)+
  geom_line(data=filters,mapping = aes(x=wl,
                                       y=Full.f),
            colour="#E6B200")+
  theme_bw()+
  xlab("Wavelength (nm)")+
  ylab(expression("   Irradiance \n (W/m^2/nm)"))+
  theme(axis.title.y = element_text(size=7),
        axis.title.x = element_text(size=7),
        axis.text.y = element_text(size=7),
        axis.text.x = element_text(size=7)) +
  scale_x_continuous(breaks=seq(300,1700,200)) +
  theme(plot.margin = unit(c(0, 0.4, 0, 0.4), "cm"))

plot_grid(plotS1D, plotS1A, plotS1B,plotS1C, 
          align = "v", axis = "bt", rel_widths = c(1, 1),
          ncol=1)

## Warning: Removed 1000 row(s) containing missing values (geom_path).

Percentages in the different bands:

# Area under the curve for each band

# Sun from Standard Solar Spectrum
# ASTM G173-03 Reference Spectra Derived from SMARTS v. 2.9.2 (AM1.5)
# Global tilt  W*m-2*nm-1

SolarAUCTot <- summary(theorySunALL)$B1
SolarAUCUV <- summary(theorySunALL[theorySunALL$wl<400,])$B1

## Warning: cannot calculate red chroma; wavelength range not between 605 and 700
## nm

## Warning: cannot calculate yellow chroma; wavelength range not between 550 and
## 625 nm

## Warning: cannot calculate green chroma; wavelength range not between 510 and 605
## nm

## Warning: cannot calculate blue chroma; wavelength range not between 400 and 510
## nm

## Warning: cannot calculate UV chroma; wavelength range not below 400 nm

## Warning: cannot calculate violet chroma; wavelength below 415 nm

SolarAUCVIS <- summary(theorySunALL[theorySunALL$wl>=400 & theorySunALL$wl<=700,])$B1
SolarAUCNIR <- summary(theorySunALL[theorySunALL$wl>700 & theorySunALL$wl<=1700,])$B1

## Warning: cannot calculate red chroma; wavelength range not between 605 and 700
## nm

## Warning: cannot calculate yellow chroma; wavelength range not between 550 and
## 625 nm

## Warning: cannot calculate green chroma; wavelength range not between 510 and 605
## nm

## Warning: cannot calculate blue chroma; wavelength range not between 400 and 510
## nm

## Warning: cannot calculate UV chroma; wavelength range not below 400 nm

## Warning: cannot calculate violet chroma; wavelength below 415 nm

SolarAUCNIRlong <- summary(theorySunALL[theorySunALL$wl>1700,])$B1

## Warning: cannot calculate red chroma; wavelength range not between 605 and 700
## nm

## Warning: cannot calculate yellow chroma; wavelength range not between 550 and
## 625 nm

## Warning: cannot calculate green chroma; wavelength range not between 510 and 605
## nm

## Warning: cannot calculate blue chroma; wavelength range not between 400 and 510
## nm

## Warning: cannot calculate UV chroma; wavelength range not below 400 nm

## Warning: cannot calculate violet chroma; wavelength below 415 nm

# Percentages per spectral band
SolarAUCUV/SolarAUCTot*100 # UV spectral band

## [1] 4.551361

SolarAUCVIS/SolarAUCTot*100 # VIS spectral band

## [1] 43.08467

SolarAUCNIR/SolarAUCTot*100 # NIR spectral band

## [1] 46.89958

SolarAUCNIRlong/SolarAUCTot*100 # NIR not considered

## [1] 5.464388

# Ratio between the full spectrum and the solar simulator in the relevant range

(AUCsun)/(SolarAUCVIS + SolarAUCNIR)

## [1] 1.414769

# Ratio per each band: 

(AUCnir)/SolarAUCNIR

## [1] 1.435048

(AUCvis)/SolarAUCVIS

## [1] 1.10968

Repeatability

Table

# Delta Temp after 5 min

RepDTdata<-
  Repeatability %>% 
  select (spp, DT5_TOT, DT5_NIR, DT5_VIS) 

RModelDT_Full <- rpt(DT5_TOT ~ (1 | spp), "spp",
                   data = RepDTdata,
                   nboot = 1000, npermut = 0,
                   datatype = "Gaussian")

## Bootstrap Progress:

summary(RModelDT_Full)

## 
## Repeatability estimation using the lmm method
## 
## Call = rpt(formula = DT5_TOT ~ (1 | spp), grname = "spp", data = RepDTdata, datatype = "Gaussian", nboot = 1000, npermut = 0)
## 
## Data: 44 observations
## ----------------------------------------
## 
## spp (6 groups)
## 
## Repeatability estimation overview: 
##       R     SE   2.5%  97.5% P_permut  LRT_P
##    0.91  0.101  0.601  0.967       NA      0
## 
## Bootstrapping and Permutation test: 
##             N   Mean Median   2.5%  97.5%
## boot     1000   0.87  0.899  0.601  0.967
## permut      1     NA     NA     NA     NA
## 
## Likelihood ratio test: 
## logLik full model = -11.36359
## logLik red. model = -47.50959
## D  = 72.3, df = 1, P = 9.28e-18
## 
## ----------------------------------------

RModelDT_nir <- rpt(DT5_NIR ~ (1 | spp), "spp",
                  data = RepDTdata,
                  nboot = 1000, npermut = 0,
                  datatype = "Gaussian")

## Bootstrap Progress:

summary(RModelDT_nir)

## 
## Repeatability estimation using the lmm method
## 
## Call = rpt(formula = DT5_NIR ~ (1 | spp), grname = "spp", data = RepDTdata, datatype = "Gaussian", nboot = 1000, npermut = 0)
## 
## Data: 44 observations
## ----------------------------------------
## 
## spp (6 groups)
## 
## Repeatability estimation overview: 
##       R     SE   2.5%  97.5% P_permut  LRT_P
##   0.795  0.145  0.361  0.919       NA      0
## 
## Bootstrapping and Permutation test: 
##             N   Mean Median   2.5%  97.5%
## boot     1000  0.739  0.779  0.361  0.919
## permut      1     NA     NA     NA     NA
## 
## Likelihood ratio test: 
## logLik full model = 0.7551644
## logLik red. model = -20.97594
## D  = 43.5, df = 1, P = 2.16e-11
## 
## ----------------------------------------

RModelDT_vis <- rpt(DT5_VIS ~ (1 | spp), "spp",
                  data = RepDTdata,
                  nboot = 1000, npermut = 0,
                  datatype = "Gaussian")

## Bootstrap Progress:

summary(RModelDT_vis)

## 
## Repeatability estimation using the lmm method
## 
## Call = rpt(formula = DT5_VIS ~ (1 | spp), grname = "spp", data = RepDTdata, datatype = "Gaussian", nboot = 1000, npermut = 0)
## 
## Data: 44 observations
## ----------------------------------------
## 
## spp (6 groups)
## 
## Repeatability estimation overview: 
##       R     SE   2.5%  97.5% P_permut  LRT_P
##   0.781  0.163  0.275  0.914       NA      0
## 
## Bootstrapping and Permutation test: 
##             N   Mean Median   2.5%  97.5%
## boot     1000  0.714  0.758  0.275  0.914
## permut      1     NA     NA     NA     NA
## 
## Likelihood ratio test: 
## logLik full model = 8.440956
## logLik red. model = -12.25749
## D  = 41.4, df = 1, P = 6.21e-11
## 
## ----------------------------------------

# Heating Rates

RepHRdata<-
  Repeatability %>% 
  select (spp, maxHRTOT, maxHRNIR, maxHRVIS)

RModelHR_Full <- rpt(maxHRTOT ~ (1 | spp), "spp",
                   data = RepHRdata,
                   nboot = 1000, npermut = 0,
                   datatype = "Gaussian")

## Bootstrap Progress:

summary(RModelHR_Full)

## 
## Repeatability estimation using the lmm method
## 
## Call = rpt(formula = maxHRTOT ~ (1 | spp), grname = "spp", data = RepHRdata, datatype = "Gaussian", nboot = 1000, npermut = 0)
## 
## Data: 44 observations
## ----------------------------------------
## 
## spp (6 groups)
## 
## Repeatability estimation overview: 
##       R     SE   2.5%  97.5% P_permut  LRT_P
##   0.424  0.187      0  0.708       NA  0.001
## 
## Bootstrapping and Permutation test: 
##             N   Mean Median   2.5%  97.5%
## boot     1000  0.379  0.383      0  0.708
## permut      1     NA     NA     NA     NA
## 
## Likelihood ratio test: 
## logLik full model = 126.561
## logLik red. model = 122.0428
## D  = 9.04, df = 1, P = 0.00132
## 
## ----------------------------------------

RModelHR_nir <- rpt(maxHRNIR ~ (1 | spp), "spp",
                   data = RepHRdata,
                   nboot = 1000, npermut = 0,
                   datatype = "Gaussian")

## Bootstrap Progress:

summary(RModelHR_nir)

## 
## Repeatability estimation using the lmm method
## 
## Call = rpt(formula = maxHRNIR ~ (1 | spp), grname = "spp", data = RepHRdata, datatype = "Gaussian", nboot = 1000, npermut = 0)
## 
## Data: 44 observations
## ----------------------------------------
## 
## spp (6 groups)
## 
## Repeatability estimation overview: 
##       R     SE   2.5%  97.5% P_permut  LRT_P
##   0.637   0.17  0.179  0.838       NA      0
## 
## Bootstrapping and Permutation test: 
##             N   Mean Median   2.5%  97.5%
## boot     1000  0.588  0.618  0.179  0.838
## permut      1     NA     NA     NA     NA
## 
## Likelihood ratio test: 
## logLik full model = 162.2648
## logLik red. model = 151.2175
## D  = 22.1, df = 1, P = 1.3e-06
## 
## ----------------------------------------

RModelHR_vis <- rpt(maxHRVIS ~ (1 | spp), "spp",
                  data = RepHRdata,
                  nboot = 1000, npermut = 0,
                  datatype = "Gaussian")

## Bootstrap Progress:

summary(RModelHR_vis)

## 
## Repeatability estimation using the lmm method
## 
## Call = rpt(formula = maxHRVIS ~ (1 | spp), grname = "spp", data = RepHRdata, datatype = "Gaussian", nboot = 1000, npermut = 0)
## 
## Data: 44 observations
## ----------------------------------------
## 
## spp (6 groups)
## 
## Repeatability estimation overview: 
##       R     SE   2.5%  97.5% P_permut  LRT_P
##    0.22  0.151      0  0.518       NA   0.06
## 
## Bootstrapping and Permutation test: 
##             N   Mean Median   2.5%  97.5%
## boot     1000  0.204  0.189      0  0.518
## permut      1     NA     NA     NA     NA
## 
## Likelihood ratio test: 
## logLik full model = 158.329
## logLik red. model = 157.1193
## D  = 2.42, df = 1, P = 0.0599
## 
## ----------------------------------------

Figures

# Delta Temp after 5 min

PlotDTFull<-
  ggplot( aes(x=spp, y=DT5_TOT), data=RepDTdata) +
  geom_boxplot(aes( fill=spp, alpha = 0.5)) +
  geom_dotplot(binaxis='y', stackdir='center', dotsize=0.7,
               colour="black",fill="gray",alpha=0.7)+
  theme_minimal() +
  theme(
    legend.position="none",
    plot.title = element_text(size=11)
  ) +
  ggtitle("Delta Temp WIDE") +
  xlab("") +
  scale_fill_manual(values = wes_palette(
    "Darjeeling1", n = 6,type = "continuous"))+
  ylab("ΔT [5] (°C/5 min)")

PlotDTnir<-
  ggplot( aes(x=spp, y=DT5_NIR), data=RepDTdata) +
  geom_boxplot(aes( fill=spp, alpha = 0.5)) +
  geom_dotplot(binaxis='y', stackdir='center', dotsize=0.7,
               colour="black",fill="gray",alpha=0.7)+
  theme_minimal() +
  theme(
    legend.position="none",
    plot.title = element_text(size=11)
  ) +
  ggtitle("Delta Temp NIR") +
  xlab("") +
  scale_fill_manual(values = wes_palette("Darjeeling1", 
                                         n = 6,type = "continuous"))+
  ylab("ΔT [5] (°C/5 min)")

PlotDTvis<-
  ggplot( aes(x=spp, y=DT5_VIS), data=RepDTdata) +
  geom_boxplot(aes( fill=spp, alpha = 0.5)) +
  geom_dotplot(binaxis='y', stackdir='center', dotsize=0.7,
               colour="black",fill="gray",alpha=0.7)+
  theme_minimal() +
  theme(
    legend.position="none",
    plot.title = element_text(size=11)
  ) +
  ggtitle("Delta Temp VIS") +
  xlab("") +
  scale_fill_manual(values = wes_palette("Darjeeling1", 
                                         n = 6,type = "continuous"))+
  ylab("ΔT [5] (°C/5 min)")

plot_grid(PlotDTFull, PlotDTnir, PlotDTvis,align = "v",
          nrow=1,ncol=3, rel_widths = c(1, 1))

# Heating Rates

PlotHRFull<-
  ggplot( aes(x=spp, y=maxHRTOT), data=RepHRdata) +
  geom_boxplot(aes( fill=spp, alpha = 0.5)) +
  geom_dotplot(binaxis='y', stackdir='center', dotsize=0.7,
               colour="black",fill="gray",alpha=0.7)+
  theme_minimal() +
  theme(
    legend.position="none",
    plot.title = element_text(size=11)
  ) +
  ggtitle("max Heating Rate WIDE") +
  xlab("") +
  scale_fill_manual(values = wes_palette("Darjeeling1", 
                                         n = 6,type = "continuous"))

PlotHRnir<-
  ggplot( aes(x=spp, y=maxHRNIR), data=RepHRdata) +
  geom_boxplot(aes( fill=spp, alpha = 0.5)) +
  geom_dotplot(binaxis='y', stackdir='center', dotsize=0.7,
               colour="black",fill="gray",alpha=0.7)+
  theme_minimal() +
  theme(
    legend.position="none",
    plot.title = element_text(size=11)
  ) +
  ggtitle("max Heating Rate NIR") +
  xlab("") +
  scale_fill_manual(values = wes_palette("Darjeeling1", 
                                         n = 6,type = "continuous"))

PlotHRvis<-
  ggplot( aes(x=spp, y=maxHRVIS), data=RepHRdata) +
  geom_boxplot(aes( fill=spp, alpha = 0.5)) +
  geom_dotplot(binaxis='y', stackdir='center', dotsize=0.7,
               colour="black",fill="gray",alpha=0.7)+
  theme_minimal() +
  theme(
    legend.position="none",
    plot.title = element_text(size=11)
  ) +
  ggtitle("max Heating Rate VIS") +
  xlab("") +
  scale_fill_manual(values = wes_palette("Darjeeling1", 
                                         n = 6,type = "continuous"))

plot_grid(PlotHRFull, PlotHRnir, PlotHRvis,align = "v",
          nrow=1,ncol=3, rel_widths = c(1, 1))

Reflectivity

Delta T5

# TOTAL
mod1D <- lm(DT5_TOT ~ Rf_TOT + size + side, data = Single)

mod1DPer <- c (
  anova(mod1D)[1,2]/sum(anova(mod1D)[,2])*100, 
  anova(mod1D)[2,2]/sum(anova(mod1D)[,2])*100, 
  anova(mod1D)[3,2]/sum(anova(mod1D)[,2])*100)

# NIR
mod2D <- lm(DT5_NIR ~ Rf_NIR + size + side, data = Single)

mod2DPer <- c ( 
  anova(mod2D)[1,2]/sum(anova(mod2D)[,2])*100, 
  anova(mod2D)[2,2]/sum(anova(mod2D)[,2])*100, 
  anova(mod2D)[3,2]/sum(anova(mod2D)[,2])*100)

# VIS
mod3D <- lm(DT5_VIS ~ Rf_VIS + size + side, data = Single)

mod3DPer <- c ( 
  anova(mod3D)[1,2]/sum(anova(mod3D)[,2])*100,
  anova(mod3D)[2,2]/sum(anova(mod3D)[,2])*100,
  anova(mod3D)[3,2]/sum(anova(mod3D)[,2])*100)


export_summs(mod1D, mod2D, mod3D, 
             error_format = "SE = {std.error} p = {p.value}",
             model.names = c("TOT", "NIR", "VIS"))

	TOT	NIR	VIS
(Intercept)	3.7782 ***	1.5649 ***	0.8795 ***
	SE = 0.3530 p = 0.0000	SE = 0.2110 p = 0.0000	SE = 0.1341 p = 0.0000
Rf_TOT	-0.0264 ***
	SE = 0.0059 p = 0.0000
size	0.5159 ***	0.2062 *	0.2416 ***
	SE = 0.1410 p = 0.0006	SE = 0.0838 p = 0.0172	SE = 0.0554 p = 0.0001
sideB	-0.1079	-0.0227	-0.1191 *
	SE = 0.1156 p = 0.3548	SE = 0.0686 p = 0.7422	SE = 0.0454 p = 0.0114
Rf_NIR		-0.0100 ***
		SE = 0.0025 p = 0.0002
Rf_VIS			-0.0068 *
			SE = 0.0033 p = 0.0440
N	56	56	56
R2	0.3876	0.2843	0.3817
* p < 0.001; p < 0.01; * p < 0.05.

data.frame( "Var" = c("OpticalProp", "Size", "Side"),
            "TOT" = mod1DPer,
            "NIR" = mod2DPer,
            "VIS" = mod3DPer)  %>% 
  kbl() %>%
  kable_paper("hover", full_width = F)

Var	TOT	NIR	VIS
OpticalProp	21.11500	19.6704960	5.200977
Size	16.62151	8.6084900	24.787868
Side	1.02668	0.1505468	8.181039

maxHR

# TOTAL
mod1E<-lm(maxHRTOT ~ Rf_TOT + size + side, data = Single)

mod1EPer <- c (
anova(mod1E)[1,2]/sum(anova(mod1E)[,2])*100,
anova(mod1E)[2,2]/sum(anova(mod1E)[,2])*100,
anova(mod1E)[3,2]/sum(anova(mod1E)[,2])*100)

# NIR
mod2E<-lm(maxHRNIR ~ Rf_NIR + size + side, data = Single)

mod2EPer <- c (
anova(mod2E)[1,2]/sum(anova(mod2E)[,2])*100,
anova(mod2E)[2,2]/sum(anova(mod2E)[,2])*100,
anova(mod2E)[3,2]/sum(anova(mod2E)[,2])*100)

# VIS
mod3E<-lm(maxHRVIS ~ Rf_VIS + size + side, data = Single)

mod3EPer <- c (
anova(mod3E)[1,2]/sum(anova(mod3E)[,2])*100,
anova(mod3E)[2,2]/sum(anova(mod3E)[,2])*100,
anova(mod3E)[3,2]/sum(anova(mod3E)[,2])*100)


export_summs(mod1E, mod2E, mod3E, 
             error_format = "SE = {std.error} p = {p.value} ", 
             model.names = c("TOT", "NIR", "VIS"))

	TOT	NIR	VIS
(Intercept)	0.0724 ***	0.0372 ***	0.0267 ***
	SE = 0.0112 p = 0.0000	SE = 0.0043 p = 0.0000	SE = 0.0054 p = 0.0000
Rf_TOT	-0.0007 ***
	SE = 0.0002 p = 0.0003
size	0.0092 *	-0.0003	0.0014
	SE = 0.0045 p = 0.0439	SE = 0.0017 p = 0.8776	SE = 0.0022 p = 0.5319
sideB	-0.0065	-0.0018	-0.0050 **
	SE = 0.0037 p = 0.0792	SE = 0.0014 p = 0.2148	SE = 0.0018 p = 0.0081
Rf_NIR		-0.0003 ***
		SE = 0.0001 p = 0.0000
Rf_VIS			-0.0003 *
			SE = 0.0001 p = 0.0375
N	56	56	56
R2	0.3173	0.3777	0.2132
* p < 0.001; p < 0.01; * p < 0.05.

data.frame( "Var" = c("OpticalProp", "Size", "Side"),
            "TOT" = mod1EPer,
            "NIR" = mod2EPer,
            "VIS" = mod3EPer)  %>% 
  kbl() %>%
  kable_paper("hover", full_width = F)

Var	TOT	NIR	VIS
OpticalProp	20.987713	35.8854246	8.779891
Size	6.534403	0.0023255	1.057666
Side	4.209703	1.8868104	11.487366

Transmissivity

Delta T5

# TOTAL
mod1F <- lm(DT5_TOT ~ Td_TOT + size + side, data = Single)

mod1FPer <- c (
  anova(mod1F)[1,2]/sum(anova(mod1F)[,2])*100, 
  anova(mod1F)[2,2]/sum(anova(mod1F)[,2])*100, 
  anova(mod1F)[3,2]/sum(anova(mod1F)[,2])*100)

# NIR
mod2F <- lm(DT5_NIR ~ Td_NIR + size + side, data = Single)

mod2FPer <- c ( 
  anova(mod2F)[1,2]/sum(anova(mod2F)[,2])*100, 
  anova(mod2F)[2,2]/sum(anova(mod2F)[,2])*100, 
  anova(mod2F)[3,2]/sum(anova(mod2F)[,2])*100)

# VIS
mod3F <- lm(DT5_VIS ~ Td_VIS + size + side, data = Single)

mod3FPer <- c ( 
  anova(mod3F)[1,2]/sum(anova(mod3F)[,2])*100,
  anova(mod3F)[2,2]/sum(anova(mod3F)[,2])*100,
  anova(mod3F)[3,2]/sum(anova(mod3F)[,2])*100)


export_summs(mod1F, mod2F, mod3F, 
             error_format = "SE = {std.error} p = {p.value}",
             model.names = c("TOT", "NIR", "VIS"))

	TOT	NIR	VIS
(Intercept)	3.7136 ***	1.3302 ***	0.9257 ***
	SE = 0.4458 p = 0.0000	SE = 0.2720 p = 0.0000	SE = 0.1399 p = 0.0000
Td_TOT	-0.0137 *
	SE = 0.0067 p = 0.0456
size	0.3103	0.1394	0.2087 ***
	SE = 0.1684 p = 0.0710	SE = 0.1016 p = 0.1761	SE = 0.0564 p = 0.0005
sideB	-0.2134	-0.0696	-0.1342 **
	SE = 0.1293 p = 0.1048	SE = 0.0773 p = 0.3724	SE = 0.0451 p = 0.0044
Td_NIR		-0.0014
		SE = 0.0028 p = 0.6059
Td_VIS			-0.0093 *
			SE = 0.0044 p = 0.0371
N	56	56	56
R2	0.2174	0.0724	0.3852
* p < 0.001; p < 0.01; * p < 0.05.

data.frame( "Var" = c("OpticalProp", "Size", "Side"),
            "TOT" = mod1FPer,
            "NIR" = mod2FPer,
            "VIS" = mod3FPer)  %>% 
  kbl() %>%
  kable_paper("hover", full_width = F)

Var	TOT	NIR	VIS
OpticalProp	11.614523	1.976260	9.686435
Size	6.023813	3.814659	18.345616
Side	4.101043	1.444339	10.484757

maxHR

# TOTAL
mod1G<-lm(maxHRTOT ~ Td_TOT + size + side, data = Single)

mod1GPer <- c (
anova(mod1G)[1,2]/sum(anova(mod1G)[,2])*100,
anova(mod1G)[2,2]/sum(anova(mod1G)[,2])*100,
anova(mod1G)[3,2]/sum(anova(mod1G)[,2])*100)

# NIR
mod2G<-lm(maxHRNIR ~ Td_NIR + size + side, data = Single)

mod2GPer <- c (
anova(mod2G)[1,2]/sum(anova(mod2G)[,2])*100,
anova(mod2G)[2,2]/sum(anova(mod2G)[,2])*100,
anova(mod2G)[3,2]/sum(anova(mod2G)[,2])*100)

# VIS
mod3G<-lm(maxHRVIS ~ Td_VIS + size + side, data = Single)

mod3GPer <- c (
anova(mod3G)[1,2]/sum(anova(mod3G)[,2])*100,
anova(mod3G)[2,2]/sum(anova(mod3G)[,2])*100,
anova(mod3G)[3,2]/sum(anova(mod3G)[,2])*100)


export_summs(mod1G, mod2G, mod3G, 
             error_format = "SE = {std.error} p = {p.value} ", 
             model.names = c("TOT", "NIR", "VIS"))

	TOT	NIR	VIS
(Intercept)	0.0602 ***	0.0279 ***	0.0269 ***
	SE = 0.0141 p = 0.0001	SE = 0.0060 p = 0.0000	SE = 0.0058 p = 0.0000
Td_TOT	-0.0001
	SE = 0.0002 p = 0.6498
size	0.0060	-0.0013	0.0004
	SE = 0.0053 p = 0.2676	SE = 0.0022 p = 0.5673	SE = 0.0023 p = 0.8574
sideB	-0.0090 *	-0.0029	-0.0056 **
	SE = 0.0041 p = 0.0322	SE = 0.0017 p = 0.0978	SE = 0.0019 p = 0.0041
Td_NIR		0.0000
		SE = 0.0001 p = 0.7858
Td_VIS			-0.0002
			SE = 0.0002 p = 0.1826
N	56	56	56
R2	0.1219	0.0613	0.1733
* p < 0.001; p < 0.01; * p < 0.05.

data.frame( "Var" = c("OpticalProp", "Size", "Side"),
            "TOT" = mod1GPer,
            "NIR" = mod2GPer,
            "VIS" = mod3GPer)  %>% 
  kbl() %>%
  kable_paper("hover", full_width = F)

Var	TOT	NIR	VIS
OpticalProp	1.039351	0.6815375	2.762681
Size	2.964358	0.3180801	0.265086
Side	8.182405	5.1295563	14.299062

Dorsal/Ventral Transm

Stats

Difference was calculated as ventral - dorsal / dorsal * 100. So, if the difference is positive, the ventral transmittance would be larger than dorsal.

Tram <-
  Optical %>% 
  select (1, 11:13 )


ttestTOT <- t.test(Tram$DifferenceTransTOT,mu=0)
ttestNIR <- t.test(Tram$DifferenceTransNIR,mu=0)
ttestVIS <- t.test(Tram$DifferenceTransVIS,mu=0)

tt1mean <- c(mean(Tram$DifferenceTransTOT),
             mean(Tram$DifferenceTransNIR),
             mean(Tram$DifferenceTransVIS))

tt1tstat <- round (
  c(parse_number(toString(ttestTOT[1])), 
    parse_number(toString(ttestNIR[1])), 
    parse_number(toString(ttestVIS[1]))), 4)


tt1pval <- round (
  c(parse_number(toString(ttestTOT[3])), 
    parse_number(toString(ttestNIR[3])), 
    parse_number(toString(ttestVIS[3]))), 4)

tt1CInum <-  round (
  c((ttestTOT[4]$conf.int[1:2]),
    (ttestNIR[4]$conf.int[1:2]),
    (ttestVIS[4]$conf.int[1:2])), 4)

ttCI <- c(
  paste(Cor1CInum[1], Cor1CInum[2], sep = " : " ),
  paste(Cor1CInum[3], Cor1CInum[4], sep = " : " ),
  paste(Cor1CInum[5], Cor1CInum[6], sep = " : " )
)            


data.frame (
  "SpectralBand" = c("TOT", "NIR", "VIS"),
  "Mean" = tt1mean,
  "t-stat" = tt1tstat,
  "C.I." = ttCI,
  "p-value" = tt1pval 
    ) %>% 
  kbl() %>%
  kable_paper("hover", full_width = F)

SpectralBand	Mean	t.stat	C.I.	p.value
TOT	10.412575	2.8217	-0.5419 : -0.0682	0.0066
NIR	9.758875	2.7390	-0.6523 : -0.2363	0.0083
VIS	17.979794	3.1409	-0.288 : 0.2374	0.0027

Plot

ggplot(Tram, aes(x = reorder (spp, DifferenceTransTOT),
                          y = DifferenceTransTOT))+
  geom_col(aes(), width = 0.7,
           position = position_stack(reverse=TRUE))+
  coord_flip()+
  theme_minimal()+
  ylim(-30,130)+
  theme(legend.position="none")+
  labs(x=" ",y=" Difference in Transmittance (%)")+
  theme(axis.text.y = element_text(size=6))+
  theme(axis.text.x = element_text(size = 10))+
  theme(axis.title.x = element_text(margin = 
                                    margin(t = 8, r = 0, b = 0, l = 0)))

HeatingRatesFullCode

Authors removed for double-blind review

12/01/2022

Introduction

Libraries

Install packages

Load packages

Data Sets

Raw Data

Spectral Data

Heating Rates Data

Others

Data Processing

Optical Properties

Reflectivity

Add filters

Multiplication function

Apply

Calculate AUC

Standardise

Transmissivity

Add filters

Multiplication function

Apply

Calculate AUC

Standardise

Ventral Transmissivity

Add filters

Multiplication function

Apply

Calculate AUC

Standardise

Consolidated

Heating

Single Elytra

Delta temp after 5 min

Max Heating Rates

Consolidated

Beetle Body

Delta temp after 5 min

Consolidated

Experimental Set up

Figure 1

Single Elytra

Stats

Linear Models

Delta T5

maxHR

Correlations

Involving Optical properties

Involving Spectral Bands

Involving Size

Reflectivity

Transmissivity

Absorptivity

Figure 2

Data wrangling

Plot

Figure 3

Beetle Body

Stats

Body

Under the elytra

Figure 4

Supplementary Materials

Light Source

Figure S2

Repeatability

Table

Figures

Reflectivity

Delta T5

maxHR

Transmissivity

Delta T5

maxHR

Dorsal/Ventral Transm

Stats

Plot