###################################################
# Melanitis leda phenotypic plasticity code
# Larval growth rate is not a major determinant of adult wing shape and eyespot size in the seasonally polyphenic butterfly Melanitis leda
# Freerk Molleman1, M. Elizabeth Moore2, Sridhar Halali3, Ullasa Kodandaramaiah4, Dheeraj Halali5, Erik van Bergen6, Paul M. Brakefield5, Vicencio Oostra7


# code written by Freerk Molleman

## figures in the MS and appendices are indicated when the figure is saved

# load libraries
library(readr)
library(readr)
library(Rmisc)
library(dplyr)
library(tidyr)
library(ggplot2)
library(cowplot)
library(ggfortify)
library(viridis)
library(plotly)
library(nlme)
library(lme4)
library(MuMIn)
library(scales)
library(lmerTest)
library(ggpubr)
library(gridExtra)
library(purrr)
library(RColorBrewer)
library(forcats)
library(broom)
library(scales)
library(MuMIn) # for model selection


theme_set(theme_classic())

packageVersion("MuMIn")


# import data where left and right were averaged, non-sexed individuals removed and parameters transformed.
PPc <- read_csv("MledaCleanR.xlsx")
PPc$expt <- as.factor(PPc$expt)
levels(PPc$expt)

## write.csv(as.data.frame(PPc), file = "MledaCleanR.csv")  # write PPc df to rescue csv file damaged by adding a sheet with readme before saving as excel file

# This PPc excludes a lot of data from the Temperature Humidity Experiment where sex was not recorded

# to get the non-sexed data for this experiment we can use
forePP <- read_csv("Mleda_PP_ForewingAveraged_CleanedIncNonSexed.csv")
THforePP <- subset(forePP, expt=="Temp_HR")

## First we needed to calculate relative size of forewing tip area (already present in current data file)
# PPc$relTAlg <- PPc$TAlg/Temp1PPstacks$WA
# hist(PPc$relTAlg)

## use square root to improve distribution
# PPc$sqrtRelTAlg <- sqrt(PPc$relTAlg)
# hist(PPc$sqrtRelTAlg)

# separate fore- from hind- wings 
forePPs <- subset(PPc, wingtype=="Forewing")
hindPPs <- subset(PPc, wingtype=="Hindwing")


# separate fore from hind wings even when sex is unknown
forePP <- subset(PP, wingtype=="Forewing")
forePP$expt <- as.factor(forePP$expt)

# select data for Tempertaure*Humdidity Experiment to make Reaction Norm graph
THforePP <- subset(forePP, expt=="Temp_HR")

#########################################################################################
############################# DENSITY PLOTS ############################################
#####################################################################################


## density plot of sqrt forewing eyespot size all experiments
# By population (India vs Ghana) and sex.

fore_ppPPesEx <- ggplot(forePPs, aes(x = sqrtES)) +
  geom_density(size = 1.5) +
  xlab("Sqrt relative eyespot size") +
  ylab("Density") +
  facet_grid(pop ~ sex, switch = "y", scale="free") +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    strip.text = element_text(color = "black", size = 16, face = "plain"),
    strip.placement = "outside",  # Place strip text outside
    strip.background = element_blank(),  # Remove strip background
    panel.border = element_rect(color = "white", fill = NA)
  )

fore_ppPPesEx

ggsave("Densityplot_EyespotSize_CountrySex.png", fore_ppPPesEx, height=10, width=15, units=c("cm"), dpi=600) ### Figure A2.1 ###

## density plot of wingshape all experiments
# By population (India vs Ghana) and sex.
fore_WS_dens <- ggplot(forePPs, aes(x = sqrtRelTAlg)) +
  geom_density(size = 1.5) +
  xlab("Wing shape") +
  ylab("Density") +
  facet_grid(pop ~ sex, switch = "y", scale="free") +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    strip.text = element_text(color = "black", size = 16, face = "plain"),
    strip.placement = "outside",  # Place strip text outside
    strip.background = element_blank(),  # Remove strip background
    panel.border = element_rect(color = "white", fill = NA)
  )

fore_WS_dens

ggsave("Densityplot_WinhShape_CountrySex.png", fore_WS_dens, height=10, width=15, units=c("cm"), dpi=600 ) ### Figure A2.5 ###

## Density Eyespot Size Temperature Experiment
dES_T1 <- ggplot(T1, aes(x = sqrtES)) +
  geom_density(size = 1.5) +
  xlab("Sqrt relative eyespot size") +
  ylab("Density") +
  facet_grid(Temp ~ sex, switch = "y", scale="free") +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    strip.text = element_text(color = "black", size = 16, face = "plain"),
    strip.placement = "outside",  # Place strip text outside
    strip.background = element_blank(),  # Remove strip background
    panel.border = element_rect(color = "white", fill = NA)
  )

dES_T1

ggsave("Densityplot_EyespotSize_T1.png", dES_T1, height=17, width=15, units=c("cm"), dpi=600) ### Figure A2.2 ###


## Density Eyespot Size Temperature Humidity Experiment
dES_TH <- ggplot(TH, aes(x = sqrtES)) +
  geom_density(size = 1.5) +
  xlab("Sqrt relative eyespot size") +
  ylab("Density") +
  facet_grid(Temp ~ sex, switch = "y", scale="free") +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    strip.text = element_text(color = "black", size = 16, face = "plain"),
    strip.placement = "outside",  # Place strip text outside
    strip.background = element_blank(),  # Remove strip background
    panel.border = element_rect(color = "white", fill = NA)
  )

dES_TH

ggsave("Densityplot_EyespotSize_TH.png", dES_TH, height=17, width=15, units=c("cm"), dpi=600) ### Figure A2.3 ###

## Density Eyespot Size Host Plant Experiment
dES_P <- ggplot(P4_8, aes(x = sqrtES)) +
  geom_density(size = 1.5) +
  xlab("Sqrt relative eyespot size") +
  ylab("Density") +
  facet_grid(Temp ~ sex, switch = "y", scale="free") +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    strip.text.x  = element_text(color = "black", size = 16, face = "plain"),
    strip.text.y = element_text(color = "black", size = 16, face = "italic"),
    strip.placement = "outside",  # Place strip text outside
    strip.background = element_blank(),  # Remove strip background
    panel.border = element_rect(color = "white", fill = NA)
  )

dES_P

ggsave("Densityplot_EyespotSize_Plant.png", dES_P, height=30, width=15, units=c("cm"), dpi=600)   ### Figure A2.4 ###


## Density Wing shape Temperature Experiment ############

dWS_T1 <- ggplot(T1, aes(x = sqrtRelTAlg)) +
  geom_density(size = 1.5) +
  xlab("Wing shape") +
  ylab("Density") +
  facet_grid(Temp ~ sex, switch = "y", scale="free") +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    strip.text = element_text(color = "black", size = 16, face = "plain"),
    strip.placement = "outside",  # Place strip text outside
    strip.background = element_blank(),  # Remove strip background
    panel.border = element_rect(color = "white", fill = NA)
  )

dWS_T1

ggsave("Densityplot_WingShape_T1.png", dWS_T1, height=17, width=15, units=c("cm"), dpi=600) ### Figure A2.6 ###


## Density Wing shape Temperature Humidity Experiment
dWS_TH <- ggplot(TH, aes(x = sqrtRelTAlg)) +
  geom_density(size = 1.5) +
  xlab("Wing shape") +
  ylab("Density") +
  facet_grid(Temp ~ sex, switch = "y", scale="free") +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    strip.text = element_text(color = "black", size = 16, face = "plain"),
    strip.placement = "outside",  # Place strip text outside
    strip.background = element_blank(),  # Remove strip background
    panel.border = element_rect(color = "white", fill = NA)
  )

dWS_TH

ggsave("Densityplot_WingShape_TH.png", dWS_TH, height=17, width=15, units=c("cm"), dpi=600) ### Figure A2.7 ###

## Density Wing shape Host Plant Experiment with N>7
dWS_P <- ggplot(P4_8, aes(x = sqrtRelTAlg)) +
  geom_density(size = 1.5) +
  xlab("Sqrt relative Wing shape") +
  ylab("Density") +
  facet_grid(Temp ~ sex, switch = "y", scale="free") +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    strip.text.x  = element_text(color = "black", size = 16, face = "plain"),
    strip.text.y = element_text(color = "black", size = 16, face = "italic"),
    strip.placement = "outside",  # Place strip text outside
    strip.background = element_blank(),  # Remove strip background
    panel.border = element_rect(color = "white", fill = NA)
  )

dWS_P

ggsave("Densityplot_WingShape_Plant.png", dWS_P, height=30, width=15, units=c("cm"), dpi=600) ### Figure A2.8 ###

str(P4_8)
######################################################
#### Data analysis per Experiment
#######################################################

# make sure all predictors are formatted to be categorical and thus have levels
forePPs$expt <- as.factor(forePPs$expt)
levels(forePPs$expt)
forePPs$Temp <- as.factor(forePPs$Temp)
levels(forePPs$Temp)
forePPs$plant <- as.factor(forePPs$plant)
levels(forePPs$plant)

# select data for Temperature Experiment 1 and remove unused levels
T1 <- subset(forePPs, expt=="Temp_1")
T1$Temp <- as.factor(T1$Temp)
# Remove unused levels
T1$Temp <- droplevels(T1$Temp)
levels(T1$Temp) # verify that there are 5 levels

# select data for Temperature Experiment 2 and remove unused levels
T2 <- subset(forePPs, expt=="Temp_2")
T2$Temp <- as.factor(T2$Temp)
T2$Temp <- droplevels(T2$Temp)
levels(T2$Temp) # verify that there are 3 levels

# select data for Temperature & Humidity Experiment  and remove unused levels
TH <- subset(forePPs, expt=="Temp_HR")
TH$Temp <- as.factor(TH$Temp)
TH$Temp <- droplevels(TH$Temp)
levels(TH$Temp) # verify that there are 4 levels

# select data for Plant experiment and remove unused levels
levels(forePPs$expt)
P4 <- subset(forePPs, expt=="Plant_4")
P4$Temp <- as.factor(P4$Temp)
P4$Temp <- droplevels(P4$Temp)
levels(P4$Temp) # verify that there is 1 level
P4$plant <- as.factor(P4$plant)
P4$Temp <- droplevels(P4$plant)
levels(P4$plant) # verify that there is 1 level


#### Data analysis per Experiment
# make sure all predictors are formatted to be levels
forePPs$expt <- as.factor(forePPs$expt)
levels(forePPs$expt)
forePPs$Temp <- as.factor(forePPs$Temp)
levels(forePPs$Temp)
forePPs$plant <- as.factor(forePPs$plant)
levels(forePPs$plant)

# select data for Temperature Experiment 1 and remove unused levels
T1 <- subset(forePPs, expt=="Temp_1")
T1$Temp <- as.factor(T1$Temp)
# Remove unused levels
T1$Temp <- droplevels(T1$Temp)
levels(T1$Temp) # verify that there are 5 levels

# select data for Temperature Experiment 2 and remove unused levels
T2 <- subset(forePPs, expt=="Temp_2")
T2$Temp <- as.factor(T2$Temp)
T2$Temp <- droplevels(T2$Temp)
levels(T2$Temp) # verify that there are 3 levels

# select data for Temperature & Humidity Experiment  and remove unused levels
TH <- subset(forePPs, expt=="Temp_HR")
TH$Temp <- as.factor(TH$Temp)
TH$Temp <- droplevels(TH$Temp)
levels(TH$Temp) # verify that there are 4 levels

# select data for Plant experiment and remove unused levels
levels(forePPs$expt)
P4 <- subset(forePPs, expt=="Plant_4")
P4$Temp <- as.factor(P4$Temp)
P4$Temp <- droplevels(P4$Temp)
levels(P4$Temp) # verify that there is 1 level
P4$plant <- as.factor(P4$plant)
P4$Temp <- droplevels(P4$plant)
levels(P4$plant) # verify that there is 1 level



####### REACTION NORMS #######


### Eyespot size #########

#######  REACTION NORM Temperature Experiment #######

### Temp 1 ###

RNT1b <- ggplot(T1, aes(x = Temp, y = sqrtES)) +
  geom_point(position = "jitter") +
  geom_violin(alpha = 0.1) +
  stat_summary(fun = "median", geom = "line", aes(group = 1), size=1.2) +
  facet_grid(cols = vars(sex)) +
    scale_y_continuous(limit = c(0, 0.3), breaks = c(0, 0.1, 0.2)) +
  # scale_fill_manual(values = c("darkgoldenrod", "lightgrey")) +
  theme_classic() +
  xlab("Temperature") +
  ylab("Sqrt relative eyespot size") +
  scale_x_discrete(labels = function(x) as.expression(bquote(.(as.character(x))))) +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    #   legend.position = c(0.15, 0.90),
    #   legend.direction = "horizontal",
    legend.background = element_rect(fill = "transparent"),
    strip.text.x = element_text(color = "black", size = 16, face = "bold"),
    legend.text = element_text(color = "black", size = 13, face = "plain"),
    legend.title = element_blank(),
    strip.placement = "outside",
    strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA)) 

RNT1b 

ggsave("ReactionNormTemp1EyeSpot.png", RNT1b , height=10, width=18, units=c("cm"), dpi=600)


## Temp 2 ### 

RNT2 <- ggplot(Temp2PPstacks_fore, aes(x = Temp, y = sqrtES)) +
  geom_point(position = "jitter") +
  geom_violin(alpha = 0.1) +
  stat_summary(fun = "median", geom = "line", aes(group = 1), size=1.2) +
  facet_grid(cols = vars(sex)) +
  #  scale_y_continuous(limit = c(0, 0.06), breaks = c(0, 0.02, 0.04, 0.06)) +
  # scale_fill_manual(values = c("darkgoldenrod", "lightgrey")) +
  theme_classic() +
  xlab("Temperature") +
  ylab("Sqrt relative eyespot size") +
  scale_x_discrete(labels = function(x) as.expression(bquote(.(as.character(x))))) +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    #   legend.position = c(0.15, 0.90),
    #   legend.direction = "horizontal",
    legend.background = element_rect(fill = "transparent"),
    strip.text.x = element_text(color = "black", size = 16, face = "bold"),
    legend.text = element_text(color = "black", size = 13, face = "plain"),
    legend.title = element_blank(),
    strip.placement = "outside",
    strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA)) 

RNT2 

ggsave("ReactionNormTemp2EyeSpot.png", RNT2 , height=10, width=18, units=c("cm"), dpi=600) ### Figure A1.1 ###

# Reaction norms #
### Temp * HUMIDITY ###

RNTHall <- ggplot(THforePP, aes(x = Temp, y = sqrtES, fill = Humidity, color = Humidity)) +
  geom_point(position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.6), size = 2, aes(color = NULL)) +
  geom_violin(alpha = 0.0, size = 1.2, position = position_dodge(width = 0.6)) +
  stat_summary(fun = "median", geom = "line", aes(group = Humidity), size = 1.2, position = position_dodge(width = 0.6)) +
  scale_color_manual(values = c("black", "blue2")) +
  scale_fill_manual(values = c("black", "blue2"))+
  theme_classic() +
  xlab("Temperature") +
  ylab("Sqrt relative eyespot size") +
  scale_x_discrete(labels = function(x) as.expression(bquote(.(as.character(x))))) +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    legend.position = c(0.15, 0.80),
    #   legend.direction = "horizontal",
    legend.background = element_rect(fill = "transparent"),
    strip.text.x = element_text(color = "black", size = 16, face = "bold"),
    legend.text = element_text(color = "black", size = 13, face = "plain"),
    legend.title = element_blank(),
    strip.placement = "outside",
    strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA)) 

RNTHall

ggsave("ReactionNormTempHREyeSpotAllData.png", RNTHall , height=10, width=18, units=c("cm"), dpi=600)


# reaction norms
### For plants 4 with N > 7

# Filter the data to include only plants with at least 8 observations
P4_8 <- P4 %>%
  group_by(plant) %>%
  filter(n() >= 8) %>%
  ungroup()

# Calculate the average sqrtES for each plant to order plants by eyespot size of butterflies
#average_sqrtES <- P4_8 %>%
# group_by(plant) %>%
#summarise(average_sqrtES = mean(sqrtES, na.rm = TRUE)) %>%
#arrange(average_sqrtES) %>%
#pull(plant)

# Calculate the average growth rate for each plant to order plants by growth rate of butterflies
average_grwth <- P4_8 %>%
  group_by(plant) %>%
  summarise(average_grwth = mean(grwth_rate, na.rm = TRUE)) %>%
  arrange(average_grwth) %>%
  pull(plant)

# Create the plot with the filtered data
RNESplants <- ggplot(P4_8, aes(x = factor(plant, levels = average_grwth), y = sqrtES, color = sex)) +
  geom_jitter(position = position_dodge2(width = 0.6, preserve = "single"), size = 2) +
  geom_violin(alpha = 0.0, size = 1, position = position_dodge(width = 0.6)) +
  stat_summary(fun = "median", geom = "line", aes(group = plant), size = 1.2) +
  scale_color_manual(values = c("red", "cornflowerblue")) +
  scale_fill_manual(values = c("red", "cornflowerblue")) +
  scale_y_continuous(limit = c(0.15, 0.36), breaks = c(0.2, 0.3)) +
  theme_classic() +
  xlab("Plant species") +
  ylab("Eyespot size") +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 12, face = "italic"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    legend.position = c(0.20, 0.93),
    legend.direction = "horizontal",
    legend.background = element_blank(),
    legend.spacing.y = unit(-1.0, 'cm'),
    legend.spacing.x = unit(0.8, 'cm'),
    strip.text.x = element_blank(),
    legend.text = element_text(color = "black", size = 13, face = "italic"),
    legend.title = element_blank(),
    strip.placement = element_blank(),
    strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA)) +
  guides(fill = guide_legend(byrow = TRUE))


RNESplants

ggsave("ReactionNormPlant4EyeSpot8_orderedGrowthRate.png", RNESplants , height=10, width=22, units=c("cm"), dpi=600)

############### Combining the Reaction Norms for Eyespot size  ###
# Plant 4 RNESplants with N>7
# RNTHall Temperature*Humdity all data
# RNT1b Temperature experiment 1
# RNT2 Temperature experiment 2

# removing y-axis labels
RNT1r <- RNT1b + ylab("")
RNTHallr <- RNTHall + ylab("")
RNESplr <- RNESplants + ylab("")

# combine Temp1, T*HR and Plant4 eyespot reaction norms
RN_ES <- ggarrange(RNT1r, RNTHallr, RNESplr, ncol = 1, labels=c("a","b","c"), heights =c(10,10,10), font.label = list(size = 22), align = "v") + bgcolor("White")+ border(NA)

# give common y-axis name
RN_ESa <- annotate_figure(RN_ES , left = text_grob("Square root relative eyespot size", color = "black", face="bold", family = "sans", size=20, rot = 90))+ bgcolor("White") + border("White")
RN_ESa


ggsave("ReactionEyeSpotCombiNewest.png", RN_ESa , height=30, width=22, units=c("cm"), dpi=600)  ## Figure 4 ###

###################################################

### #########################  Wing Shape #########

#######  REACTION NORM Temperature Experiment #######

### Temp 1 ###

RNT1WS <- ggplot(T1, aes(x = Temp, y = sqrtRelTAlg)) +
  geom_point(position = "jitter") +
  geom_violin(alpha = 0.1) +
  stat_summary(fun = "median", geom = "line", aes(group = 1), size=1.2) +
  facet_grid(cols = vars(sex)) +
  #  scale_y_continuous(limit = c(0, 0.06), breaks = c(0, 0.02, 0.04, 0.06)) +
  # scale_fill_manual(values = c("darkgoldenrod", "lightgrey")) +
  theme_classic() +
  xlab("Temperature") +
  ylab("Wing shape") +
  scale_x_discrete(labels = function(x) as.expression(bquote(.(as.character(x))))) +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    #   legend.position = c(0.15, 0.90),
    #   legend.direction = "horizontal",
    legend.background = element_rect(fill = "transparent"),
    strip.text.x = element_text(color = "black", size = 16, face = "bold"),
    legend.text = element_text(color = "black", size = 13, face = "plain"),
    legend.title = element_blank(),
    strip.placement = "outside",
    strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA)) 

RNT1WS

ggsave("ReactionNormTemp1_WingShape.png", RNT1WS , height=10, width=18, units=c("cm"), dpi=600)


## Temp 2 ### 

RNT2WS <- ggplot(Temp2PPstacks_fore, aes(x = Temp, y = sqrtRelTAlg)) +
  geom_point(position = "jitter") +
  geom_violin(alpha = 0.1) +
  stat_summary(fun = "median", geom = "line", aes(group = 1), size=1.2) +
  facet_grid(cols = vars(sex)) +
  #  scale_y_continuous(limit = c(0, 0.06), breaks = c(0, 0.02, 0.04, 0.06)) +
   theme_classic() +
  xlab("Temperature") +
  ylab("Wing shape") +
  scale_x_discrete(labels = function(x) as.expression(bquote(.(as.character(x))))) +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    #   legend.position = c(0.15, 0.90),
    #   legend.direction = "horizontal",
    legend.background = element_rect(fill = "transparent"),
    strip.text.x = element_text(color = "black", size = 16, face = "bold"),
    legend.text = element_text(color = "black", size = 13, face = "plain"),
    legend.title = element_blank(),
    strip.placement = "outside",
    strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA)) 

RNT2WS 

ggsave("ReactionNormTemp2_WingShape.png", RNT2WS , height=10, width=18, units=c("cm"), dpi=600) ### Figure A1.2 ###

# Reaction norms #
### Temp * HUMIDITY ###

RNTHallWS <- ggplot(THforePP, aes(x = Temp, y = sqrtRelTAlg, fill = Humidity, color = Humidity)) +
  geom_point(position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.6), size = 2, aes(color = NULL)) +
  geom_violin(alpha = 0.0, size = 1.2, position = position_dodge(width = 0.6)) +
  stat_summary(fun = "median", geom = "line", aes(group = Humidity), size = 1.2, position = position_dodge(width = 0.6)) +
  scale_color_manual(values = c("black", "blue2")) +
  scale_fill_manual(values = c("black", "blue2"))+
  theme_classic() +
  xlab("Temperature") +
  ylab("Wing shape") +
  scale_x_discrete(labels = function(x) as.expression(bquote(.(as.character(x))))) +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    legend.position = c(0.85, 0.80),
    #   legend.direction = "horizontal",
    legend.background = element_rect(fill = "transparent"),
    strip.text.x = element_text(color = "black", size = 16, face = "bold"),
    legend.text = element_text(color = "black", size = 13, face = "plain"),
    legend.title = element_blank(),
    strip.placement = "outside",
    strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA)) 

RNTHallWS

ggsave("ReactionNorm_TempHR_WingShape.png", RNTHallWS , height=10, width=18, units=c("cm"), dpi=600)


# reaction norms
### For plants 4 with N > 7

# Filter the data to include only plants with at least 8 observations
P4_8 <- P4 %>%
  group_by(plant) %>%
  filter(n() >= 8) %>%
  ungroup()

# Calculate the average growth rate for each plant to order plants by growth rate of butterflies
average_grwth <- P4_8 %>%
  group_by(plant) %>%
  summarise(average_grwth = mean(grwth_rate, na.rm = TRUE)) %>%
  arrange(average_grwth) %>%
  pull(plant)

# Create the plot with the filtered data
RNESplantsWS <- ggplot(P4_8, aes(x = factor(plant, levels = average_grwth), y = sqrtRelTAlg, color = sex)) +
  geom_jitter(position = position_dodge2(width = 0.6, preserve = "single"), size = 2) +
  geom_violin(alpha = 0.0, size = 1, position = position_dodge(width = 0.6)) +
  stat_summary(fun = "median", geom = "line", aes(group = plant), size = 1.2) +
  scale_color_manual(values = c("red", "cornflowerblue")) +
  scale_fill_manual(values = c("red", "cornflowerblue")) +
  #scale_y_continuous(limit = c(0.15, 0.36), breaks = c(0.2, 0.3)) +
  theme_classic() +
  xlab("Plant species") +
  ylab("Wing shape") +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 12, face = "italic"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    legend.position = c(0.75, 0.89),
    legend.direction = "horizontal",
    legend.background = element_blank(),
    legend.spacing.y = unit(-1.0, 'cm'),
    legend.spacing.x = unit(0.8, 'cm'),
    strip.text.x = element_blank(),
    legend.text = element_text(color = "black", size = 13, face = "italic"),
    legend.title = element_blank(),
    strip.placement = element_blank(),
    strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA)) +
  guides(fill = guide_legend(byrow = TRUE))


RNESplantsWS

ggsave("ReactionNormPlant4_WingShape.png", RNESplantsWS , height=10, width=22, units=c("cm"), dpi=600)

############### Combining the Reaction Norms for Wing Shape  ###
# Plant 4 RNESplantsWS with N>7
# RNTHallWS Temperature*Humdity all data
# RNT1WS Temperature experiment 1
# RNT2WS Temperature experiment 2

# removing y-axis labels
RNT1WSr <- RNT1WS + ylab("")
RNTHallWSr <- RNTHallWS + ylab("")
RNESplWSr <- RNESplantsWS + ylab("")

# combine Temp1, T*HR and Plant4 eyespot reaction norms
RN_WS <- ggarrange(RNT1WSr, RNTHallWSr, RNESplWSr, ncol = 1, labels=c("a","b","c"), heights =c(10,10,10), font.label = list(size = 22), align = "v") + bgcolor("White")+ border(NA)

# give common y-axis name
RN_WSa <- annotate_figure(RN_WS , left = text_grob("           Wing shape", color = "black", face="bold", family = "sans", size=20, rot = 90))+ bgcolor("White") + border("White")
RN_WSa


ggsave("Reaction_WingShape_Combi.png", RN_WSa , height=30, width=22, units=c("cm"), dpi=600)  ### Figure 5 ###


###################################################################
### !!! ### growth rate  ######### !!! ###

#######  REACTION NORM Temperature Experiment #######

### Temp 1 ###

RNT1_GR <- ggplot(T1, aes(x = Temp, y = grwth_rate)) +
  geom_point(position = "jitter") +
  geom_violin(alpha = 0.1) +
  stat_summary(fun = "median", geom = "line", aes(group = 1), size=1.2) +
  facet_grid(cols = vars(sex)) +
  theme_classic() +
  xlab("Temperature") +
  ylab("Larval growth rate") +
  scale_x_discrete(labels = function(x) as.expression(bquote(.(as.character(x))))) +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    legend.background = element_rect(fill = "transparent"),
    strip.text.x = element_text(color = "black", size = 16, face = "bold"),
    legend.text = element_text(color = "black", size = 13, face = "plain"),
    legend.title = element_blank(),
    strip.placement = "outside",
    strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA)) 

RNT1_GR 

ggsave("ReactionNormTemp1Growth_rate.png", RNT1_GR , height=10, width=18, units=c("cm"), dpi=600)


## Temp 2 ### 

RNT2_GR <- ggplot(Temp2PPstacks_fore, aes(x = Temp, y = grwth_rate)) +
  geom_point(position = "jitter") +
  geom_violin(alpha = 0.1) +
  stat_summary(fun = "median", geom = "line", aes(group = 1), size=1.2) +
  facet_grid(cols = vars(sex)) +
  theme_classic() +
  xlab("Temperature") +
  ylab("Larval growth rate") +
  scale_x_discrete(labels = function(x) as.expression(bquote(.(as.character(x))))) +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    legend.background = element_rect(fill = "transparent"),
    strip.text.x = element_text(color = "black", size = 16, face = "bold"),
    legend.text = element_text(color = "black", size = 13, face = "plain"),
    legend.title = element_blank(),
    strip.placement = "outside",
    strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA)) 

RNT2_GR

ggsave("ReactionNormTemp2Growth_rate.png", RNT2_GR , height=10, width=18, units=c("cm"), dpi=600)

# Reaction norms #
### Temp * HUMIDITY ###

RNTHall_GR <- ggplot(THforePP, aes(x = Temp, y = grwth_rate, fill = Humidity, color = Humidity)) +
  geom_point(position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.6), size = 2, aes(color = NULL)) +
  geom_violin(alpha = 0.0, size = 1.2, position = position_dodge(width = 0.6)) +
  stat_summary(fun = "median", geom = "point", aes(group = Humidity), size = 4, shape = 5, stroke = 2, position = position_dodge(width = 0.6)) +
    scale_color_manual(values = c("black", "blue2")) +
  scale_fill_manual(values = c("black", "blue2"))+
  theme_classic() +
  xlab("Temperature") +
  ylab("Larval growth rate") +
  scale_x_discrete(labels = function(x) as.expression(bquote(.(as.character(x))))) +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    legend.position = c(0.15, 0.80),
    #   legend.direction = "horizontal",
    legend.background = element_rect(fill = "transparent"),
    strip.text.x = element_text(color = "black", size = 16, face = "bold"),
    legend.text = element_text(color = "black", size = 13, face = "plain"),
    legend.title = element_blank(),
    strip.placement = "outside",
    strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA)) 

RNTHall_GR

ggsave("ReactionNormTempHRGrowthRateAllData.png", RNTHall_GR , height=10, width=18, units=c("cm"), dpi=600)


# reaction norms
### For plants 4 with N>7

# Data were before filtered to get plants swith N>7 thus Filter the data to include only plants with at least 8 observations



# Her its assumed that you already calculated the average growth rate for each plant to order plants by growth rate of butterflies for the reaction norms for eyespot size


# Create the plot with the filtered data
RNESplantsGR <- ggplot(P4_8, aes(x = factor(plant, levels = average_grwth), y = grwth_rate, color = sex)) +
  geom_jitter(position = position_dodge2(width = 0.6, preserve = "single"), size = 2) +
  geom_violin(alpha = 0.0, size = 1, position = position_dodge(width = 0.6)) +
  stat_summary(fun = "median", geom = "line", aes(group = plant), size = 1.2) +
  scale_color_manual(values = c("red", "cornflowerblue")) +
  scale_fill_manual(values = c("red", "cornflowerblue")) +
 # scale_y_continuous(limit = c(0.15, 0.36), breaks = c(0.2, 0.3)) +
  theme_classic() +
  xlab("Plant species") +
  ylab("Larval growth rate") +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 12, face = "italic"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    legend.position = c(0.12, 0.75),
    legend.direction = "vertical",
    legend.background = element_blank(),
    legend.spacing.y = unit(-1.0, 'cm'),
    legend.spacing.x = unit(0.8, 'cm'),
    strip.text.x = element_blank(),
    legend.text = element_text(color = "black", size = 13, face = "italic"),
    legend.title = element_blank(),
    strip.placement = element_blank(),
    strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA)) +
  guides(fill = guide_legend(byrow = TRUE))


RNESplantsGR

ggsave("ReactionNormPlant4_GrowthRate_8_orderedGrowthRate.png", RNESplantsGR , height=10, width=22, units=c("cm"), dpi=600)

############### Combining the Reaction Norms for Eyespot size
# Plant 4 RNESplants with N>7
# RNTHall Temperature*Humdity all data
# RNT1b Temperature experiment 1
# RNT2 Temperature experiment 2

# removing y-axis labels
RNT1GRr <- RNT1_GR + ylab("")
RNTHallGRr <- RNTHall_GR + ylab("")
RNESplGRr <- RNESplantsGR + ylab("")

# combine Temp1, T*HR and Plant4 eyespot reaction norms
RN_GR <- ggarrange(RNT1GRr, RNTHallGRr, RNESplGRr, ncol = 1, labels=c("a","b","c"), heights =c(10,10,10), font.label = list(size = 22), align = "v") + bgcolor("White")+ border(NA)

# give common y-axis name
RN_GRa <- annotate_figure(RN_GR , left = text_grob("Larval growth rate", color = "black", face="bold", family = "sans", size=20, rot = 90))+ bgcolor("White") + border("White")
RN_GRa


ggsave("ReactionNorm_GrowthRate.png", RN_GRa , height=30, width=22, units=c("cm"), dpi=600) ## Figure 3 ###


#############################
#### Questions of size  ######

### How are pupal mass and wing area related to each other?

## Wing are vs pupal mass

### Temp 1 ###
str(T1)
ra <-range(T1$pup_mss, na.rm = TRUE)
ra

Wa_PM <-lmer(WAd  ~ pup_mss + Temp + sex + Temp*sex + (1|sleeve), data=T1)
summary(Wa_PM)
anova(Wa_PM)

  
  
WA_PM_T1 <- ggplot(T1, aes(x = pup_mss, y = WAd, color = Temp)) +
  geom_point() +
  geom_smooth(method='lm',se=FALSE, na.rm = TRUE) +
    facet_grid(cols = vars(sex)) +  # , rows = vars(Temp)
  scale_x_continuous(limit = c(0.25, 0.7), breaks = c(0.3, 0.4, 0.5, 0.6)) +
  scale_color_manual(values=c("blue", "#0066FF", "orange", "red","darkred"))+
  theme_classic() +
  xlab("Pupal mass (grams)") +
  ylab("Wing area index") +
    theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    #   legend.position = c(0.15, 0.90),
    #   legend.direction = "horizontal",
    legend.background = element_rect(fill = "transparent"),
    strip.text.x = element_text(color = "black", size = 16, face = "bold"),
    legend.text = element_text(color = "black", size = 13, face = "plain"),
    legend.title = element_blank(),
    strip.placement = "outside",
    strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA)) 

WA_PM_T1

ggsave("WingArea_pupal_massT1_byTemp.png", WA_PM_T1 , height=10, width=18, units=c("cm"), dpi=600) ### Figure A3.1. ###

#######  REACTION NORM for WING AREA  #######

### Temp 1: in Temperature Experiment   ###

RNT1bWA <- ggplot(T1, aes(x = Temp, y = WA)) +
  geom_point(position = "jitter") +
  geom_violin(alpha = 0.1) +
  stat_summary(fun = "median", geom = "line", aes(group = 1), size=1.2) +
  facet_grid(cols = vars(sex)) +
  #  scale_y_continuous(limit = c(0, 0.06), breaks = c(0, 0.02, 0.04, 0.06)) +
  # scale_fill_manual(values = c("darkgoldenrod", "lightgrey")) +
  theme_classic() +
  xlab("Temperature") +
  ylab("Wing area index") +
  scale_x_discrete(labels = function(x) as.expression(bquote(.(as.character(x))))) +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    #   legend.position = c(0.15, 0.90),
    #   legend.direction = "horizontal",
    legend.background = element_rect(fill = "transparent"),
    strip.text.x = element_text(color = "black", size = 16, face = "bold"),
    legend.text = element_text(color = "black", size = 13, face = "plain"),
    legend.title = element_blank(),
    strip.placement = "outside",
    strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA)) 

RNT1bWA 

ggsave("ReactionNormTemp1WingArea.png", RNT1bWA , height=10, width=18, units=c("cm"), dpi=600)


## Temp 2 ### 

RNT2WA <- ggplot(Temp2PPstacks_fore, aes(x = Temp, y = WA)) +
  geom_point(position = "jitter") +
  geom_violin(alpha = 0.1) +
  stat_summary(fun = "median", geom = "line", aes(group = 1), size=1.2) +
  facet_grid(cols = vars(sex)) +
  #  scale_y_continuous(limit = c(0, 0.06), breaks = c(0, 0.02, 0.04, 0.06)) +
  # scale_fill_manual(values = c("darkgoldenrod", "lightgrey")) +
  theme_classic() +
  xlab("Temperature") +
  ylab("Wing area index") +
  scale_x_discrete(labels = function(x) as.expression(bquote(.(as.character(x))))) +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    #   legend.position = c(0.15, 0.90),
    #   legend.direction = "horizontal",
    legend.background = element_rect(fill = "transparent"),
    strip.text.x = element_text(color = "black", size = 16, face = "bold"),
    legend.text = element_text(color = "black", size = 13, face = "plain"),
    legend.title = element_blank(),
    strip.placement = "outside",
    strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA)) 

RNT2WA 

ggsave("ReactionNormTemp2WingArea.png", RNT2WA , height=10, width=18, units=c("cm"), dpi=600)

# Reaction norms #
### Temp * HUMIDITY ###

RNTHallWA <- ggplot(THforePP, aes(x = Temp, y = WA, fill = Humidity, color = Humidity)) +
  geom_point(position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.6), size = 2, aes(color = NULL)) +
  geom_violin(alpha = 0.0, size = 1.2, position = position_dodge(width = 0.6)) +
  stat_summary(fun = "median", geom = "line", aes(group = Humidity), size = 1.2, position = position_dodge(width = 0.6)) +
  scale_color_manual(values = c("black", "blue2")) +
  scale_fill_manual(values = c("black", "blue2"))+
  theme_classic() +
  xlab("Temperature") +
  ylab("Sqrt relative eyespot size") +
  scale_x_discrete(labels = function(x) as.expression(bquote(.(as.character(x))))) +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    legend.position = c(0.85, 0.80),
    #   legend.direction = "horizontal",
    legend.background = element_rect(fill = "transparent"),
    strip.text.x = element_text(color = "black", size = 16, face = "bold"),
    legend.text = element_text(color = "black", size = 13, face = "plain"),
    legend.title = element_blank(),
    strip.placement = "outside",
    strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA)) 

RNTHallWA

ggsave("ReactionNormTempHRWingArea.png", RNTHallWA , height=10, width=18, units=c("cm"), dpi=600)


# reaction norms Wing Area
### For plants 4 with N > 7

# Create the plot with the filtered data
RNESplantsWA <- ggplot(P4_8, aes(x = factor(plant, levels = average_grwth), y = WA, color = sex)) +
  geom_jitter(position = position_dodge2(width = 0.6, preserve = "single"), size = 2) +
  geom_violin(alpha = 0.0, size = 1, position = position_dodge(width = 0.6)) +
  stat_summary(fun = "median", geom = "line", aes(group = plant), size = 1.2) +
  scale_color_manual(values = c("red", "cornflowerblue")) +
  scale_fill_manual(values = c("red", "cornflowerblue")) +
  #scale_y_continuous(limit = c(0.15, 0.36), breaks = c(0.2, 0.3)) +
  theme_classic() +
  xlab("Plant species") +
  ylab("Wing area index") +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 12, face = "italic"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    legend.position = c(0.20, 0.93),
    legend.direction = "horizontal",
    legend.background = element_blank(),
    legend.spacing.y = unit(-1.0, 'cm'),
    legend.spacing.x = unit(0.8, 'cm'),
    strip.text.x = element_blank(),
    legend.text = element_text(color = "black", size = 13, face = "italic"),
    legend.title = element_blank(),
    strip.placement = element_blank(),
    strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA)) +
  guides(fill = guide_legend(byrow = TRUE))


RNESplantsWA

ggsave("ReactionNormPlant4WingAREA8_orderedGrowthRate.png", RNESplantsWA , height=10, width=22, units=c("cm"), dpi=600)

############### Combining the Reaction Norms for Wing Area  ###
# Plant 4 RNESplants with N>7
# RNTHall Temperature*Humdity all data
# RNT1b Temperature experiment 1
# RNT2 Temperature experiment 2

# removing y-axis labels
RNT1rWA <- RNT1bWA + ylab("")
RNTHallrWA <- RNTHallWA + ylab("")
RNESplrWA <- RNESplantsWA + ylab("")

# combine Temp1, T*HR and Plant4 eyespot reaction norms
RN_WA<- ggarrange(RNT1rWA, RNTHallrWA, RNESplrWA, ncol = 1, labels=c("a","b","c"), heights =c(10,10,10), font.label = list(size = 22), align = "v") + bgcolor("White")+ border(NA)

# give common y-axis name
RN_WAa <- annotate_figure(RN_WA , left = text_grob("Wing area index", color = "black", face="bold", family = "sans", size=20, rot = 90))+ bgcolor("White") + border("White")
RN_WAa


ggsave("ReactionWINGAREACombi.png", RN_WAa , height=30, width=22, units=c("cm"), dpi=600) ### Figure A3.2 ###



########## END of REACTION NORMS #############



########### Statistical tests   ###############

# Temperature 1 


#### using Temp as a fixed factor AND SLEEVE as a random effect
## with growth rate
# eyespot size

# Subset data to complete cases for the predictors used in the model
complete_cases <- complete.cases(T1$sqrtES, T1$grwth_rate, T1$Temp, T1$sex, T1$sleeve)
T1c <- T1[complete_cases, ]  # Filter data to include only complete cases

# create column GR0 (growth rate average zero) where the average is zero by subtracting the average from each growth rate observation
# Calculate the mean of the growth rate column
mean_grwth_rate <- mean(T1c$grwth_rate)

# Create the new column GR0 by subtracting the mean from each growth rate value
T1c$GR0 <- T1c$grwth_rate - mean_grwth_rate

options(na.action = "na.fail")
T1Fr<-lmer(sqrtES ~ GR0 + Temp + sex + Temp*sex + GR0*sex + (1|sleeve), data=T1c) ## full model
summary(T1Fr)
anova(T1Fr)

dT1Fr <- dredge(T1Fr)   # Model selection based on cAIC
dT1Fr
write.csv(as.data.frame(dT1Fr), file = "modelsES_T1_0.csv")  # write results of dredge into CSV file

# top model
T1FrESt<-lmer(sqrtES ~ GR0 + Temp + (1|sleeve), data=T1c)
summary(T1FrESt)
anova(T1FrESt)

## wing shape

# Subset data to complete cases for the predictors used in the model
complete_cases <- complete.cases(T1$sqrtRelTAlg, T1$grwth_rate, T1$Temp, T1$sex, T1$sleeve)
T1c <- T1[complete_cases, ]  # Filter data to include only complete cases

T1FrWS<-lmer(sqrtRelTAlg  ~ GR0 + Temp + sex  +Temp*sex + GR0*sex + (1|sleeve), data=T1c)
summary(T1FrWS)
anova(T1FrWS)

options(na.action = "na.fail")


T1FrWSr <- dredge(T1FrWS)   # Model selection based on cAIC
T1FrWSr
write.csv(as.data.frame(T1FrWSr), file = "modelsWS_T1_WA_0.csv")  # write results of dredge into CSV file

# top model
T1FrWSt<-lmer(sqrtRelTAlg ~ GR0 + Temp + sex + GR0*sex + (1|sleeve), data=T1c)
summary(T1FrWSt)
anova(T1FrWSt)

## with larval development time and pupal mass
# eyespot size
complete_cases <- complete.cases(T1$sqrtES, T1$pup_mss, T1$larv_time, T1$Temp, T1$sex, T1$sleeve)
T1c <- T1[complete_cases, ]  # Filter data to include only complete cases

options(na.action = "na.fail")  # to not allow missing values for dredge
# no random effect bcs of singularity
T1Frc<-lm(sqrtES ~ larv_time + pup_mss + Temp + sex + larv_time*pup_mss + Temp*sex, data=T1c)
summary(T1Frc)
anova(T1Frc)

T1Frcr <- dredge(T1Frc)   # Model selection based on cAIC
T1Frcr
write.csv(as.data.frame(T1Frcr), file = "modelsES_T1comp.csv")  # write results of dredge into CSV file

# top model
# in this case, full model is top model

# wing shape
complete_cases <- complete.cases(T1$sqrtRelTAlg, T1$pup_mss, T1$larv_time, T1$Temp, T1$sex, T1$sleeve)
T1c <- T1[complete_cases, ]  # Filter data to include only complete cases

options(na.action = "na.fail")  # to not allow missing values for dredge
T1FrWSc<-lmer(sqrtRelTAlg  ~ larv_time + pup_mss + Temp + sex + larv_time*pup_mss + Temp*sex + (1|sleeve), data=T1c)
summary(T1FrWSc)
anova(T1FrWSc)


T1FrWScr <- dredge(T1FrWSc)   # Model selection based on cAIC
T1FrWScr
write.csv(as.data.frame(T1FrWScr), file = "modelsWS_T1comp.csv")  # write results of dredge into CSV file

# top model (no interactions)
T1FrWSct<-lmer(sqrtRelTAlg  ~ larv_time + pup_mss + Temp + sex  + (1|sleeve), data=T1c)
summary(T1FrWSct)
anova(T1FrWSct)

##########################################

######## Plots of growth rate vs phenotype ########

### growth rate vs eyespot size


# Temperature Experiment 1 

EST1a <- ggplot(T1, aes(x = grwth_rate, y = sqrtES, color=Temp)) +
  geom_point(shape=19, size=3) +
  scale_color_manual(values=c("blue", "#0066FF", "orange", "red","darkred"))+
  geom_smooth(method='lm',se=FALSE) +
  facet_grid(cols = vars(sex)) +
  scale_y_continuous(limit = c(0.0, 0.30), breaks = c(0.0, 0.1, 0.2, 0.3)) +
  scale_x_continuous(limit = c(0.015, 0.055), breaks = c(0, 0.02, 0.04)) +
  theme_classic() +
  xlab("Larval growth rate") +
  ylab("Square root relative eyespot size") +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    legend.position = c(0.4, 0.95),
    legend.direction = "horizontal",
    legend.background = element_blank(),
    strip.text.x = element_text(color = "black", size = 16, face = "bold"),
    legend.text = element_text(color = "black", size = 13, face = "plain"),
    legend.title = element_blank(),
    strip.placement = "outside",
    strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA)
  )

EST1a

ggsave("EyespotSizeGrowthRateRegsTemperature1.png", EST1a, height=15, width=20, units=c("cm"), dpi=600)

###############################

## Temperature * Humidity

#### using Temp as a fixed factor AND SLEEVE as a random effect
## with growth rate
# eyespot size

complete_cases <- complete.cases(TH$sqrtES, TH$grwth_rate, TH$Temp, TH$sex, TH$sleeve)
THc <- TH[complete_cases, ]  # Filter data to include only complete cases

# create column GR0 (growth rate average zero) where the average is zero by subtracting the average from each growth rate observation
# Calculate the mean of the growth rate column
mean_grwth_rate <- mean(THc$grwth_rate)

# Create the new column GR0 by subtracting the mean from each growth rate value
THc$GR0 <- THc$grwth_rate - mean_grwth_rate

THFr<-lmer(sqrtES ~ GR0 + Temp + Humidity + sex + Temp*sex + GR0*sex + (1|sleeve), data=THc)
summary(THFr)
anova(THFr)

# This complete model is rather over-parametrised compared to the observations available, so we use this simpler Limited model after eliminating interactions and sex that were p > 0.5
THFrL<-lmer(sqrtES ~ GR0 + Temp + (1|sleeve), data=THc)
summary(THFrL)
anova(THFrL)

options(na.action = "na.fail")  # to not allow missing values for dredge

THFrr <- dredge(THFr)   # Model selection based on cAIC
THFrr
write.csv(as.data.frame(THFrr), file = "modelsES_TH_0.csv")  # write results of dredge into CSV file

# top model
THFrt<-lmer(sqrtES ~ GR0 + (1|sleeve), data=THc)
summary(THFrt)
anova(THFrt)

THFrTT<-lmer(sqrtES ~ GR0 + Temp +(1|sleeve), data=THc)
anova(THFrTT)


### wing shape for TH
complete_cases <- complete.cases(TH$sqrtRelTAlg, TH$grwth_rate, TH$Temp, TH$sex, TH$sleeve)
THc <- TH[complete_cases, ]  # Filter data to include only complete cases

# create column GR0 (growth rate average zero) where the average is zero by subtracting the average from each growth rate observation
# Calculate the mean of the growth rate column
mean_grwth_rate <- mean(THc$grwth_rate)

# Create the new column GR0 by subtracting the mean from each growth rate value
THc$GR0 <- THc$grwth_rate - mean_grwth_rate

THFrWS<-lmer(sqrtRelTAlg  ~ GR0 + Temp + Humidity + sex +Temp*sex + GR0*sex + (1|sleeve), data=THc)
summary(THFrWS)
anova(THFrWS)

# This complete model is rather over parametrised compared to the number of observations available, so we use this simpler Limited model after eliminaing interactions and sex that were p > 0.5
THFrWSL<-lmer(sqrtRelTAlg ~ GR0 + Temp + (1|sleeve), data=THc)
summary(THFrWSL)
anova(THFrWSL)

options(na.action = "na.fail")  # to not allow missing values for dredge

THFrWSr <- dredge(THFrWS)   # Model selection based on cAIC
THFrWSr
write.csv(as.data.frame(THFrWSr), file = "modelsWS_TH_0.csv")  # write results of dredge into CSV file

# top model
THFrWSrt<-lmer(sqrtRelTAlg ~ GR0 + (1|sleeve), data=THc)
summary(THFrWSrt)
anova(THFrWSrt)

## with larval development time and pupal mass
# eyespot size

complete_cases <- complete.cases(TH$sqrtES, TH$pup_mss, TH$larv_time, TH$Temp, TH$sex, TH$sleeve)
THc <- TH[complete_cases, ]  # Filter data to include only complete cases

THFrc<-lmer(sqrtES ~ larv_time + pup_mss + Temp + Humidity + sex + larv_time*pup_mss + Temp*sex + (1|sleeve), data=THc)
summary(THFrc)
anova(THFrc)

options(na.action = "na.fail")  # to not allow missing values for dredge

THFrcr <- dredge(THFrc)   # Model selection based on cAIC
THFrcr
write.csv(as.data.frame(THFrcr), file = "modelsES_THcomp.csv")  # write results of dredge into CSV file

# top model
THFrcrt<-lmer(sqrtES ~ larv_time + (1|sleeve), data=THc)
summary(THFrcrt)
anova(THFrcrt)

# wing shape
complete_cases <- complete.cases(TH$sqrtRelTAlg, TH$pup_mss, TH$larv_time, TH$Temp, TH$sex, TH$sleeve)
THc <- TH[complete_cases, ]  # Filter data to include only complete cases

THFrWSc<-lmer(sqrtRelTAlg  ~ larv_time + pup_mss + Temp + Humidity + sex + larv_time*pup_mss + Temp*sex + (1|sleeve), data=THc)
summary(THFrWSc)
anova(THFrWSc)

THFrWScr <- dredge(THFrWSc)   # Model selection based on cAIC
THFrWScr
write.csv(as.data.frame(THFrWScr), file = "modelsWS_THcomp.csv")  # write results of dredge into CSV file

# top model
THFrWScrt<-lmer(sqrtES ~ larv_time + pup_mss + (1|sleeve), data=THc)
summary(THFrWScrt)
anova(THFrWScrt)

## graph growth rate vs eyespot size Temperature and Humidity Experiment

ESTHa <- ggplot(TH, aes(x = grwth_rate, y = sqrtES, color=Temp, shape=Humidity)) +
  geom_point(size=3) +
  scale_color_manual(values=c("blue","orange","red","darkred"))+
  geom_smooth(method='lm',se=FALSE) +
  facet_grid(cols = vars(sex)) +
  scale_y_continuous(limit = c(0.00, 0.31), breaks = c(0.0, 0.1, 0.2, 0.3)) +
  scale_x_continuous(limit = c(0.015, 0.055), breaks = c(0, 0.02, 0.04)) +
  theme_classic() +
  xlab("Larval growth rate") +
  ylab("Square root relative eyespot size") +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    legend.position = c(0.3, 0.95),
    legend.direction = "horizontal",
    legend.spacing = unit(-0.7, "cm"),
    legend.background = element_blank(),
    strip.text.x =element_blank(),
    legend.text = element_text(color = "black", size = 13, face = "plain"),
    legend.title = element_blank(),
    strip.placement = "outside",
    #strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.text.y =  element_blank(),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA)
  ) 

ESTHa


ggsave("EyespotSizeGrowthRateRegsTempHRexp.png", ESTHa, height=15, width=20, units=c("cm"), dpi=600)


#### Plants 4  regressions

#### using Temp as a fixed factor AND SLEEVE as a random effect
## with growth rate
# eyespot size

complete_cases <- complete.cases(P4$sqrtES, P4$grwth_rate, P4$plant, P4$sex, P4$sleeve, P4$WA)
P4c <- P4[complete_cases, ]  # Filter data to include only complete cases

# create column GR0 (growth rate average zero) where the average is zero by subtracting the average from each growth rate observation
# Calculate the mean of the growth rate column
mean_grwth_rate <- mean(P4c$grwth_rate)

# Create the new column GR0 by subtracting the mean from each growth rate value
P4c$GR0 <-P4c$grwth_rate - mean_grwth_rate


P4Fr<-lmer(sqrtES ~ GR0 + plant + sex +  GR0*sex + (1|sleeve), data=P4c)
summary(P4Fr)
anova(P4Fr)

options(na.action = "na.fail")  # to not allow missing values for dredge

P4Frr <- dredge(P4Fr)   # Model selection based on cAIC
P4Frr
write.csv(as.data.frame(P4Frr), file = "modelsES_P4comp_0.csv")  # write results of dredge into CSV file

# top model
P4Frrt <-lmer(sqrtES ~ GR0 + sex + GR0*sex +(1|sleeve), data=P4c)
summary(P4Frrt)
anova(P4Frrt)

## wing shape ###
complete_cases <- complete.cases(P4$sqrtRelTAlg, P4$grwth_rate, P4$plant, P4$sex, P4$sleeve, P4$WA)
P4c <- P4[complete_cases, ]  # Filter data to include only complete cases

# create column GR0 (growth rate average zero) where the average is zero by subtracting the average from each growth rate observation
# Calculate the mean of the growth rate column
mean_grwth_rate <- mean(P4c$grwth_rate)

# Create the new column GR0 by subtracting the mean from each growth rate value
P4c$GR0 <-P4c$grwth_rate - mean_grwth_rate


P4FrWS<-lmer(sqrtRelTAlg  ~ GR0 + plant + sex +  GR0*sex + (1|sleeve), data=P4c)
summary(P4FrWS)
anova(P4FrWS)

P4FrWSr <- dredge(P4FrWS)   # Model selection based on AICc
P4FrWSr

write.csv(as.data.frame(P4FrWSr), file = "modelsWS_P4_0.csv")  # write results of dredge into CSV file

# top model
P4FrWSt<-lmer(sqrtRelTAlg  ~ GR0 + sex + GR0*sex + (1|sleeve), data=P4c)
summary(P4FrWSt)
anova(P4FrWSt)


## with larval development time and pupal mass
# eyespot size
complete_cases <- complete.cases(P4$sqrtES, P4$pup_mss, P4$larv_time, P4$plant, P4$sex, P4$sleeve)
P4c <- P4[complete_cases, ]  # Filter data to include only complete cases

P4Frc<-lmer(sqrtES ~ larv_time + pup_mss + plant + sex + larv_time*pup_mss  + (1|sleeve), data=P4c)
summary(P4Frc)
anova(P4Frc)

P4Frcr <- dredge(P4Frc)   # Model selection based on cAIC
P4Frcr

write.csv(as.data.frame(P4Frcr), file = "modelsES_P4compo.csv")  # write results of dredge into CSV file

# top model
P4Frct<-lmer(sqrtES  ~ larv_time + sex + (1|sleeve), data=P4c)
summary(P4Frct)
anova(P4Frct)

## wing shape
complete_cases <- complete.cases(P4$sqrtRelTAlg , P4$pup_mss, P4$larv_time, P4$plant, P4$sex, P4$sleeve)
P4c <- P4[complete_cases, ]  # Filter data to include only complete cases

P4FrWSc<-lmer(sqrtRelTAlg  ~ larv_time + pup_mss + plant + sex + larv_time*pup_mss + (1|sleeve), data=P4c)
summary(P4FrWSc)
anova(P4FrWSc)


P4FrWScr <- dredge(P4FrWSc)   # Model selection based on cAIC
P4FrWScr

write.csv(as.data.frame(P4FrWScr), file = "modelsWS_P4comp.csv")  # write results of dredge into CSV file

# top model
P4FrWSct<-lmer(sqrtRelTAlg  ~ larv_time + pup_mss +  sex  + (1|sleeve), data=P4c)
summary(P4FrWSct)
anova(P4FrWSct)


# graph

PlantES <- ggplot(P4, aes(x = grwth_rate, y = sqrtES, color=plant)) +
  geom_point(shape=1, size=3) +
  scale_fill_brewer(palette="Set2")+
  #scale_color_manual(values=c("blue", "#0066FF", "orange", "red","darkred"))+
  geom_smooth(method='lm',se=FALSE)  +
  facet_grid(cols = vars(sex)) +
  scale_y_continuous(limit = c(0.15, 0.30), breaks = c(0.2, 0.3)) +
  scale_x_continuous(limit = c(0.015, 0.055), breaks = c(0, 0.02, 0.04)) +
  theme_classic() +
  xlab("Larval growth rate") +
  ylab("Sqrt rel. eyespot size") +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    legend.text = element_text(color = "black", size = 13, face = "italic"),
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.background = element_blank(),
    legend.spacing.y = unit(-1.0, 'cm'),
    legend.spacing.x = unit(0.8, 'cm'),
    strip.text.x = element_blank(),
    legend.title = element_blank(),
    strip.placement = "outside",
    strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA))  +
  guides(fill = guide_legend(byrow = TRUE))

PlantES


ggsave("EyespotSizeGrowthRateRegsPLANTs.png", PlantES, height=15, width=20, units=c("cm"), dpi=600)

##### Combining the PLOTS for eyespot size ##########################

# removing y-axis labels
EST1r <- EST1a + ylab("") + xlab("")
ESTHnewr <- ESTHa + ylab("") + xlab("")
PlantESr <- PlantES + ylab("")

# combine Temp1, T*HR and Plant4 eyespot size vs growth rate
Growth_ES <- ggarrange(EST1r, ESTHnewr, PlantESr, ncol = 1, labels=c("a","b","c"), heights =c(10,10,17), font.label = list(size = 22), align = "v") + bgcolor("White") + border(NA)
Growth_ES

# give common y-axis name
Growth_ESa <- annotate_figure(Growth_ES , left = text_grob("                Square root relative eyespot size", color = "black", face="bold", family = "sans", size=20, rot = 90))+ bgcolor("White") + border("White")
Growth_ESa


ggsave("GrowthRateEyeSpotCombiNEW.png", Growth_ESa , height=34, width=24, units=c("cm"), dpi=600) ### Figure 8 ###





##########################################
## !!!!!!!!!!!!!!!!!!!!!!!!!!! #######################

######## Plots of Eyespot size vs Wing Shape  ###################

# statistical model model
ModES_WS_T1 <-lmer(sqrtRelTAlg  ~ sqrtES + Temp  +  sex  + (1|sleeve), data=T1)
summary(ModES_WS_T1)
anova(ModES_WS_T1)


# Temperature Experiment 1 

ES_WST1 <- ggplot(T1, aes(x = sqrtES, y = sqrtRelTAlg, color=Temp)) +
  geom_point(shape=19, size=3) +
  scale_color_manual(values=c("blue", "#0066FF", "orange", "red","darkred"))+
  geom_smooth(method='lm',se=FALSE) +
  facet_grid(cols = vars(sex)) +
  scale_y_continuous(limit = c(0.0, 0.30), breaks = c(0.0, 0.1, 0.2, 0.3)) +
  #scale_x_continuous(limit = c(0.015, 0.055), breaks = c(0, 0.02, 0.04)) +
  theme_classic() +
  xlab("Square root relative eyespot size") +
  ylab("Wing shape") +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    legend.position = c(0.4, 0.95),
    legend.direction = "horizontal",
    legend.background = element_blank(),
    strip.text.x = element_text(color = "black", size = 16, face = "bold"),
    legend.text = element_text(color = "black", size = 13, face = "plain"),
    legend.title = element_blank(),
    strip.placement = "outside",
    strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA)
  )

ES_WST1

ggsave("EyespotSize_WingShape_Temperature1.png", ES_WST1, height=15, width=20, units=c("cm"), dpi=600)  

###############################

## Temperature * Humidity

# stats model
ModES_WS_TH <-lmer(sqrtRelTAlg  ~ sqrtES + Temp  +  sex  + (1|sleeve), data=TH)
summary(ModES_WS_TH)
anova(ModES_WS_TH)

## graph eyespot size vs wing shape, Temperature and Humidity Experiment

ES_WSTH <- ggplot(TH, aes(x = sqrtES, y = sqrtRelTAlg, color=Temp, shape=Humidity)) +
  geom_point(size=3) +
  scale_color_manual(values=c("blue","orange","red","darkred"))+
  geom_smooth(method='lm',se=FALSE) +
  facet_grid(cols = vars(sex)) +
  scale_y_continuous(limit = c(0.00, 0.27), breaks = c(0.0, 0.1, 0.2)) +
  #scale_x_continuous(limit = c(0.015, 0.055), breaks = c(0, 0.02, 0.04)) +
  theme_classic() +
  xlab("Square root relative eyespot size") +
  ylab("Wing shape") +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    legend.position = c(0.3, 0.95),
    legend.direction = "horizontal",
    legend.spacing = unit(-0.7, "cm"),
    legend.background = element_blank(),
    strip.text.x =element_blank(),
    legend.text = element_text(color = "black", size = 13, face = "plain"),
    legend.title = element_blank(),
    strip.placement = "outside",
    #strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.text.y =  element_blank(),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA)
  ) 

ES_WSTH


ggsave("EyespotSize_WingShape_TempHR.png", ES_WSTH, height=15, width=20, units=c("cm"), dpi=600)


#### Plants 4  regressions


ModES_WS_Plant <-lmer(sqrtRelTAlg  ~ sqrtES + plant  +  sex  + (1|sleeve), data=P4)
summary(ModES_WS_Plant)
anova(ModES_WS_Plant)
# graph
PlantES_WS <- ggplot(P4, aes(x = sqrtES, y = sqrtRelTAlg, color=plant)) +
  geom_point(shape=1, size=3) +
  scale_fill_brewer(palette="Set2")+
  #scale_color_manual(values=c("blue", "#0066FF", "orange", "red","darkred"))+
  geom_smooth(method='lm',se=FALSE)  +
  facet_grid(cols = vars(sex)) +
  scale_y_continuous(limit = c(0.05, 0.23), breaks = c(0.1, 0.2)) +
  scale_x_continuous(limit = c(0.15, 0.38), breaks = c(0.2, 0.3)) +
  theme_classic() +
  xlab("Square root relative eyespot size") +
  ylab("Wing shape") +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    legend.text = element_text(color = "black", size = 13, face = "italic"),
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.background = element_blank(),
    legend.spacing.y = unit(-1.0, 'cm'),
    legend.spacing.x = unit(0.8, 'cm'),
    strip.text.x = element_blank(),
    legend.title = element_blank(),
    strip.placement = "outside",
    strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA))  +
  guides(fill = guide_legend(byrow = TRUE))

PlantES_WS


ggsave("EyespotSize_WingShape_PLANT.png", PlantES_WS, height=15, width=20, units=c("cm"), dpi=600)

##### Combining the PLOTS for eyespot size vs wing shape ##########################

# removing y-axis labels
ES_WST1r<- ES_WST1 + ylab("") + xlab("")
ES_WSTHr <- ES_WSTH + ylab("") + xlab("")
PlantES_WSr <- PlantES_WS + ylab("")

# combine Temp1, T*HR and Plant4 eyespot size vs wing shape
ES_WS <- ggarrange(ES_WST1r, ES_WSTHr, PlantES_WSr, ncol = 1, labels=c("a","b","c"), heights =c(10,10,17), font.label = list(size = 22), align = "v") + bgcolor("White") + border(NA)
ES_WS

# give common y-axis name
ES_WSa <- annotate_figure(ES_WS , left = text_grob("                Wing shape", color = "black", face="bold", family = "sans", size=20, rot = 90))+ bgcolor("White") + border("White")
ES_WSa


ggsave("EyeSpot_WingShapeCombi.png", ES_WSa , height=34, width=24, units=c("cm"), dpi=600)  #### Figure 6 ###


######## Plots of growth rate vs phenotype ########

### growth rate vs wing shape

## regressions for wing shape  ######################

### Temperature 1 wing shape vs growth rate
# Temp1PPstacks_fore$Temp <- factor(T1$Temp, levels = rev(levels(Temp1PPstacks_fore$Temp))) # to make sure Temp treatments are in correct order

# Temperature 1 plot
T1RegWS <- ggplot(T1, aes(x = grwth_rate, y = sqrtRelTAlg, color=Temp)) +
  geom_point(shape=19, size=3) +
  scale_color_manual(values=c("blue", "#0066FF", "orange", "red","darkred"))+
  geom_smooth(method='lm',se=FALSE) +
  facet_grid(cols = vars(sex)) +
  scale_y_continuous(limit = c(0.0, 0.27), breaks = c(0.0, 0.1, 0.2)) +
  scale_x_continuous(limit = c(0.015, 0.055), breaks = c(0, 0.02, 0.04)) +
  theme_classic() +
  xlab("Larval growth rate") +
  ylab("Forewing-tip shape") +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    legend.position = c(0.4, 0.95),
    legend.direction = "horizontal",
    legend.background = element_blank(),
    strip.text.x = element_text(color = "black", size = 16, face = "bold"),
    legend.text = element_text(color = "black", size = 13, face = "plain"),
    legend.title = element_blank(),
    strip.placement = "outside",
    strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA)
  )

T1RegWS

ggsave("GrowthRate_WingShape_Temperature1.png", T1RegWS, height=15, width=20, units=c("cm"), dpi=600)

# Temperature and Humidity plot
THregWS <- ggplot(TH, aes(x = grwth_rate, y = sqrtRelTAlg, color=Temp, shape=Humidity)) +
  geom_point(size=3) +
  scale_color_manual(values=c("blue","orange","red","darkred"))+
  geom_smooth(method='lm',se=FALSE) +
  facet_grid(cols = vars(sex)) +
  scale_y_continuous(limit = c(0.05, 0.26), breaks = c(0.1, 0.2)) +
  scale_x_continuous(limit = c(0.015, 0.055), breaks = c(0, 0.02, 0.04)) +
  theme_classic() +
  xlab("Larval growth rate") +
  ylab("Forewing-tip shape") +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    legend.position = c(0.3, 0.95),
    legend.direction = "horizontal",
    legend.spacing = unit(-0.7, "cm"),
    legend.background = element_blank(),
    strip.text.x =element_blank(),
    legend.text = element_text(color = "black", size = 13, face = "plain"),
    legend.title = element_blank(),
    strip.placement = "outside",
    #strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.text.y =  element_blank(),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA)
  ) 

THregWS

ggsave("GrowthRate_WingShape_TemperatureHR.png", THregWS, height=15, width=20, units=c("cm"), dpi=600)

# Host Plant Experiment
PregWS <- ggplot(P4, aes(x = grwth_rate, y = sqrtRelTAlg, color=plant)) +
  geom_point(shape=1, size=3) +
  scale_fill_brewer(palette="Set2")+
  #scale_color_manual(values=c("blue", "#0066FF", "orange", "red","darkred"))+
  geom_smooth(method='lm',se=FALSE)  +
  facet_grid(cols = vars(sex)) +
  scale_y_continuous(limit = c(0.05, 0.19), breaks = c(0.05,0.1,0.15)) +
  scale_x_continuous(limit = c(0.015, 0.055), breaks = c(0, 0.02, 0.04)) +
  theme_classic() +
  xlab("Larval growth rate") +
  ylab("Forewing-tip shape") +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    legend.text = element_text(color = "black", size = 13, face = "italic"),
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.background = element_blank(),
    legend.spacing.y = unit(-1.0, 'cm'),
    legend.spacing.x = unit(0.8, 'cm'),
    strip.text.x = element_blank(),
    legend.title = element_blank(),
    strip.placement = "outside",
    strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA))  +
  guides(fill = guide_legend(byrow = TRUE))

PregWS

# removing y-axis labels
T1regrWS <- T1RegWS  + ylab("") + xlab("")
THregrWS <- THregWS + ylab("") + xlab("")
PregrWS <- PregWS + ylab("")

# combine Temp1, T*HR and Plant4 Wing area vs wing shape
GR_WS <- ggarrange(T1regrWS, THregrWS, PregrWS, ncol = 1, labels=c("a","b","c"), heights =c(12,10,17), font.label = list(size = 22), align = "v") + bgcolor("White") +border(NA)
GR_WS

# give common y-axis name
GR_WSa <- annotate_figure(GR_WS , left = text_grob("                 Wing shape", color = "black", face="bold", family = "sans", size=20, rot = 90))+ bgcolor("White") + border("White")
GR_WSa


ggsave("GrowthRate_Wingshape_Combi.png", GR_WSa , height=34, width=24, units=c("cm"), dpi=600)  ### Figure 9 #####

## using Temp as a fixed factor AND SLEEVE as a random effect
T1$WAd <- ((T1$WA)/2000)  # to make scale of Wing area similar to that of grwth_rate

T1FWS<-lmer(sqrtRelTAlg ~ grwth_rate + Temp + WAd + sex  + (1|sleeve), data=T1)
summary(T1FWS)
anova(T1FWS)


##############################################################################
##### WING AREA vs WING SHAPE ###############################################

## plot wing area against tip shape
T1RegWA <- ggplot(T1, aes(x = WA, y = sqrtRelTAlg, color=Temp)) +
  geom_point(shape=19, size=3) +
  scale_color_manual(values=c("blue", "#0066FF", "orange", "red","darkred"))+
  geom_smooth(method='lm',se=FALSE) +
  facet_grid(cols = vars(sex)) +
  scale_y_continuous(limit = c(0.0, 0.25), breaks = c(0.0, 0.1, 0.2)) +
  #scale_x_continuous(limit = c(0.015, 0.055), breaks = c(0, 0.02, 0.04)) +
  theme_classic() +
  xlab("Wing Area Index") +
  ylab("Wing shape") +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    legend.position = c(0.4, 0.95),
    legend.direction = "horizontal",
    legend.background = element_blank(),
    strip.text.x = element_text(color = "black", size = 16, face = "bold"),
    legend.text = element_text(color = "black", size = 13, face = "plain"),
    legend.title = element_blank(),
    strip.placement = "outside",
    strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA)
  )

T1RegWA

ggsave("ForeWingShapeTemperature1GrowthRateNEW.png", T1Reg, height=15, width=20, units=c("cm"), dpi=600)


### Growth rate vs wing shape Temperature* Humidity

TH$WAd <- ((TH$WA)/2000)  # to make scale of Wing area similar to that of grwth_rate

THFWS<-lmer(sqrtRelTAlg ~ grwth_rate + Temp + WAd + sex + sex*Temp + (1|sleeve), data=TH)
summary(THFWS)
anova(THFWS)



## Plot Wing area to wing shape Temp*Hum
THregWA <- ggplot(TH, aes(x = WA, y = sqrtRelTAlg, color=Temp, shape=Humidity)) +
  geom_point(size=3) +
  scale_color_manual(values=c("blue","orange","red","darkred"))+
  geom_smooth(method='lm',se=FALSE) +
  facet_grid(cols = vars(sex)) +
  scale_y_continuous(limit = c(0.05, 0.25), breaks = c(0.1, 0.2)) +
  #scale_x_continuous(limit = c(0.015, 0.055), breaks = c(0, 0.02, 0.04)) +
  theme_classic() +
  xlab("Wing area index") +
  ylab("Wing shape") +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    legend.position = c(0.3, 0.95),
    legend.direction = "horizontal",
    legend.spacing = unit(-0.7, "cm"),
    legend.background = element_blank(),
    strip.text.x =element_blank(),
    legend.text = element_text(color = "black", size = 13, face = "plain"),
    legend.title = element_blank(),
    strip.placement = "outside",
    #strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.text.y =  element_blank(),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA)
  ) 

THregWA

ggsave("WINGAREA_ForeWingShapeTemperatureHR.png", THregWA, height=15, width=20, units=c("cm"), dpi=600) 

## Plant 4 growth rate vs wing shape

P4$WAd <- ((P4$WA)/2000)  # to make scale of Wing area similar to that of grwth_rate

P4FWS<-lmer(sqrtRelTAlg ~ grwth_rate + Temp + WAd + sex + (1|sleeve), data=TH)
summary(P4FWS)
anova(P4FWS)




# plot wing area against wing shape for plants data
PregWA <- ggplot(P4, aes(x = WA, y = sqrtRelTAlg, color=plant)) +
  geom_point(shape=1, size=3) +
  scale_fill_brewer(palette="Set2")+
  #scale_color_manual(values=c("blue", "#0066FF", "orange", "red","darkred"))+
  geom_smooth(method='lm',se=FALSE)  +
  facet_grid(cols = vars(sex)) +
  scale_y_continuous(limit = c(0.05, 0.19), breaks = c(0.05,0.1,0.15)) +
  #scale_x_continuous(limit = c(0.015, 0.055), breaks = c(0, 0.02, 0.04)) +
  theme_classic() +
  xlab("Wing area index") +
  ylab("Forewing-tip shape") +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    legend.text = element_text(color = "black", size = 13, face = "italic"),
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.background = element_blank(),
    legend.spacing.y = unit(-1.0, 'cm'),
    legend.spacing.x = unit(0.8, 'cm'),
    strip.text.x = element_blank(),
    legend.title = element_blank(),
    strip.placement = "outside",
    strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA))  +
  guides(fill = guide_legend(byrow = TRUE))

PregWA



# removing y-axis labels
T1regrWA <- T1RegWA  + ylab("") + xlab("")
THregrWA <- THregWA + ylab("") + xlab("")
PregrWA <- PregWA + ylab("")

# combine Temp1, T*HR and Plant4 Wing area vs wing shape
WA_WS <- ggarrange(T1regrWA, THregrWA, PregrWA, ncol = 1, labels=c("a","b","c"), heights =c(12,10,17), font.label = list(size = 22), align = "v") + bgcolor("White") +border(NA)
WA_WS

# give common y-axis name
WA_WSa <- annotate_figure(WA_WS , left = text_grob("                 Wing shape", color = "black", face="bold", family = "sans", size=20, rot = 90))+ bgcolor("White") + border("White")
WA_WSa


ggsave("WINGAREAWingshapeCombi.png", WA_WSa , height=34, width=24, units=c("cm"), dpi=600)  ### Figure A3.3 ###


##############################################################################
##### WING AREA vs EYESPOT SIZE ###############################################

T1WA_WS<-lmer(sqrtES ~  Temp + WAd + sex + (WAd*sex) + (1|sleeve), data=T1)
summary(T1WA_WS)
anova(T1WA_WS)

## plot wing area against tip shape
T1RegWA_ES <- ggplot(T1, aes(x = WA, y = sqrtES, color=Temp)) +
  geom_point(shape=19, size=3) +
  scale_color_manual(values=c("blue", "#0066FF", "orange", "red","darkred"))+
  geom_smooth(method='lm',se=FALSE) +
  facet_grid(cols = vars(sex)) +
  scale_y_continuous(limit = c(0.0, 0.27), breaks = c(0.1, 0.2)) +
  #scale_x_continuous(limit = c(0.015, 0.055), breaks = c(0, 0.02, 0.04)) +
  theme_classic() +
  xlab("Wing area") +
  ylab("Square root relative eyespot size") +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    legend.position = c(0.4, 0.95),
    legend.direction = "horizontal",
    legend.background = element_blank(),
    strip.text.x = element_text(color = "black", size = 16, face = "bold"),
    legend.text = element_text(color = "black", size = 13, face = "plain"),
    legend.title = element_blank(),
    strip.placement = "outside",
    strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA)
  )

T1RegWA_ES

ggsave("WingSize_WingShapeTemperature1.png", T1RegWA_ES, height=15, width=20, units=c("cm"), dpi=600)


### Growth rate vs Square root relative eyespot size Temperature* Humidity


THWA_WS<-lmer(sqrtES ~ grwth_rate + Temp + WAd + sex + sex*Temp + (1|sleeve), data=TH)
summary(THWA_WS)
anova(THWA_WS)



## Plot Wing area to Square root relative eyespot size Temp*Hum
THregWA <- ggplot(TH, aes(x = WA, y = sqrtES, color=Temp, shape=Humidity)) +
  geom_point(size=3) +
  scale_color_manual(values=c("blue","orange","red","darkred"))+
  geom_smooth(method='lm',se=FALSE) +
  facet_grid(cols = vars(sex)) +
  scale_y_continuous(limit = c(0.05, 0.25), breaks = c(0.1, 0.2)) +
  #scale_x_continuous(limit = c(0.015, 0.055), breaks = c(0, 0.02, 0.04)) +
  theme_classic() +
  xlab("Wing area index") +
  ylab("Square root relative eyespot size") +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    legend.position = c(0.3, 0.95),
    legend.direction = "horizontal",
    legend.spacing = unit(-0.7, "cm"),
    legend.background = element_blank(),
    strip.text.x =element_blank(),
    legend.text = element_text(color = "black", size = 13, face = "plain"),
    legend.title = element_blank(),
    strip.placement = "outside",
    #strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.text.y =  element_blank(),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA)
  ) 

THregWA

ggsave("WINGAREA_ForeWingShapeTemperatureHR.png", THregWA, height=15, width=20, units=c("cm"), dpi=600)

## Plant 4 growth rate vs Square root relative eyespot size

P4$WAd <- ((P4$WA)/2000)  # to make scale of Wing area similar to that of grwth_rate

P4FWS<-lmer(sqrtES ~ grwth_rate + Temp + WAd + sex + (1|sleeve), data=TH)
summary(P4FWS)
anova(P4FWS)




# plot wing area against Square root relative eyespot size for plants data
PregWA <- ggplot(P4, aes(x = WA, y = sqrtES, color=plant)) +
  geom_point(shape=1, size=3) +
  scale_fill_brewer(palette="Set2")+
  #scale_color_manual(values=c("blue", "#0066FF", "orange", "red","darkred"))+
  geom_smooth(method='lm',se=FALSE)  +
  facet_grid(cols = vars(sex)) +
  scale_y_continuous(limit = c(0.05, 0.19), breaks = c(0.05,0.1,0.15)) +
  #scale_x_continuous(limit = c(0.015, 0.055), breaks = c(0, 0.02, 0.04)) +
  theme_classic() +
  xlab("Wing area index") +
  ylab("Forewing-tip shape") +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    legend.text = element_text(color = "black", size = 13, face = "italic"),
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.background = element_blank(),
    legend.spacing.y = unit(-1.0, 'cm'),
    legend.spacing.x = unit(0.8, 'cm'),
    strip.text.x = element_blank(),
    legend.title = element_blank(),
    strip.placement = "outside",
    strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA))  +
  guides(fill = guide_legend(byrow = TRUE))

PregWA



# removing y-axis labels
T1regrWA <- T1RegWA  + ylab("") + xlab("")
THregrWA <- THregWA + ylab("") + xlab("")
PregrWA <- PregWA + ylab("")

# combine Temp1, T*HR and Plant4 Wing area vs Square root relative eyespot size
WA_WS <- ggarrange(T1regrWA, THregrWA, PregrWA, ncol = 1, labels=c("a","b","c"), heights =c(12,10,17), font.label = list(size = 22), align = "v") + bgcolor("White") +border(NA)
WA_WS

# give common y-axis name
WA_WSa <- annotate_figure(WA_WS , left = text_grob("                 Forewing-tip shape", color = "black", face="bold", family = "sans", size=20, rot = 90))+ bgcolor("White") + border("White")
WA_WSa


ggsave("WINGAREAWingshapeCombi.png", WA_WSa , height=34, width=24, units=c("cm"), dpi=600)  ### Figure A3.3 ### ALTernative


######################################################################################################################################
#### Potting pupal mass (y), development time (x) and wing shape (color gradient).
########################################################################################################################################

## Dev time by size reaction norm, Plotting pupal mass (y), development time (x) and wing shape (color gradient) for Temperature experiment shape=19, 

PMDTT1 <- ggplot(T1, aes(y = pup_mss, x = larv_time, color = sqrtRelTAlg, linetype = Temp)) +
  geom_point(aes(shape = Temp), size = 3, fill = "black") + # Filled black points
  geom_smooth(aes(group = Temp), method = 'lm', se = FALSE) + # Group by Temp for line mapping
  scale_linetype_manual(values = c("solid", "dashed", "dotted", "dotdash", "longdash")) + # Specify linetypes for Temp levels
  scale_shape_manual(values = c(17, 16, 15, 18, 6)) + # Filled shapes: circle, triangle (up and down facing), square, and diamond
  scale_color_gradient(low = "lightblue", high = "black") + # Increase color contrast from light blue to black
  facet_grid(cols = vars(sex)) +
  theme_classic() +
  xlab("Larval development time (days)") +
  ylab("Pupal mass (grams)") +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    legend.direction = "vertical",
    legend.background = element_blank(),
    strip.text.x = element_text(color = "black", size = 16, face = "bold"),
    legend.text = element_text(color = "black", size = 13, face = "plain"),
    legend.title = element_text(color = "black", size = 13, face = "plain"),
    strip.placement = "outside",
    strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA)
  )+
  labs(color = "Wing shape")

PMDTT1


ggsave("Devtime_Pupmass_ForeWingShape_Temperature1.png", PMDTT1, height=18, width=40, units=c("cm"), dpi=600) #### Figure A3.4. ###

## Make plots for each treatment separately, faceted by sex
T1_19 <- subset(T1, Temp == '19 C')
T1_22 <- subset(T1, Temp == '22 C')
T1_25 <- subset(T1, Temp == '25 C')
T1_25 <- subset(T1, Temp == '25 C')
T1_31 <- subset(T1, Temp == '31 C')

# Temp is 19 C (for methods loess span = 3,)

range_sqrtRelTAlg <- range(T1_19$sqrtRelTAlg) # find the range to use in scale_fill_gradient

DP19 <- ggplot(T1_19, aes(y = pup_mss, x = larv_time, color = sqrtRelTAlg)) +
  geom_point(size = 3) + # Points
  geom_smooth(method = 'lm',  se = FALSE) + # Line mapping
  stat_summary_2d(
    aes(fill = ..value.., z = sqrtRelTAlg), 
    fun = "mean", bins = 5, geom = "tile", 
    color = NA, alpha = 0.5) + # Heatmap based on sqrtRelTAlg values
  scale_color_gradient(low = "lightblue", high = "black") + # Increase color contrast from light blue to black
  scale_fill_gradient(low = "lightblue", high = "black", guide = FALSE, limits = c(0.13, 0.25)) + # Adjust fill scale limits
  scale_alpha_continuous(range = c(0.3, 1), guide = FALSE) + # Set alpha range
  scale_x_continuous(limit = c(35, 55), breaks = c(40, 50)) + 
  scale_y_continuous(limit = c(0.3, 0.7), breaks = c(0.3, 0.4, 0.5, 0.6, 0.7)) +
  facet_grid(cols = vars(sex)) +
  theme_classic() +
  xlab("Larval development time (days)") +
  ylab("Pupal mass (grams)") +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    legend.direction = "vertical",
    legend.background = element_blank(),
    strip.text.x = element_text(color = "black", size = 16, face = "bold"),
    legend.text = element_text(color = "black", size = 13, face = "plain"),
    legend.title = element_text(color = "black", size = 13, face = "plain"),
    strip.placement = "outside",
    strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA)
  ) +
  labs(color = "Wing shape")

DP19


# Temp is 22 C (for methods loess span = 3,)
range_sqrtRelTAlg <- range(T1_22$sqrtRelTAlg) # find the range to use in scale_fill_gradient

DP22 <- ggplot(T1_22, aes(y = pup_mss, x = larv_time, color = sqrtRelTAlg)) +
  geom_point(size = 3) + # Points
  geom_smooth(method = 'lm',  se = FALSE) + # Line mapping
  stat_summary_2d(
    aes(fill = ..value.., z = sqrtRelTAlg), 
    fun = "mean", bins = 5, geom = "tile", 
    color = NA, alpha = 0.5) + # Heatmap based on sqrtRelTAlg values
  scale_color_gradient(low = "lightblue", high = "black") + # Increase color contrast from light blue to black
  scale_fill_gradient(low = "lightblue", high = "black", guide = FALSE, limits = c(0.03, 0.24)) + # Adjust fill scale limits
  scale_alpha_continuous(range = c(0.3, 1), guide = FALSE) + # Set alpha range

  facet_grid(cols = vars(sex)) +
  theme_classic() +
  xlab("Larval development time (days)") +
  ylab("Pupal mass (grams)") +
  theme(
    axis.title.x = element_text(color = "white", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    legend.direction = "vertical",
    legend.background = element_blank(),
    strip.text.x = element_text(color = "white", size = 16, face = "bold"),
    legend.text = element_text(color = "black", size = 13, face = "plain"),
    legend.title = element_text(color = "black", size = 13, face = "plain"),
    strip.placement = "outside",
    strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA)
  ) +
  labs(color = "Wing shape")

DP22

# Temp is 25 C 
range_sqrtRelTAlg <- range(T1_25$sqrtRelTAlg) # find the range to use in scale_fill_gradient

DP25 <- ggplot(T1_25, aes(y = pup_mss, x = larv_time, color = sqrtRelTAlg)) +
  geom_point(size = 3) + # Points
  geom_smooth(method = 'lm',  se = FALSE) + # Line mapping
  stat_summary_2d(
    aes(fill = ..value.., z = sqrtRelTAlg), 
    fun = "mean", bins = 5, geom = "tile", 
    color = NA, alpha = 0.5) + # Heatmap based on sqrtRelTAlg values
  scale_color_gradient(low = "lightblue", high = "black") + # Increase color contrast from light blue to black
  scale_fill_gradient(low = "lightblue", high = "black", guide = FALSE, limits = c(0.1, 0.25)) + # Adjust fill scale limits
  scale_alpha_continuous(range = c(0.3, 1), guide = FALSE) + # Set alpha range
  scale_x_continuous(limit = c(20, 40), breaks = c(25, 35)) + 
  scale_y_continuous(limit = c(0.31, 0.65), breaks = c(0.4, 0.5, 0.6)) +
  facet_grid(cols = vars(sex)) +
  theme_classic() +
  xlab("Larval development time (days)") +
  ylab("Pupal mass (grams)") +
  theme(
    axis.title.x = element_text(color = "white", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    legend.direction = "vertical",
    legend.background = element_blank(),
    strip.text.x = element_text(color = "white", size = 16, face = "bold"),
    legend.text = element_text(color = "black", size = 13, face = "plain"),
    legend.title = element_text(color = "black", size = 13, face = "plain"),
    strip.placement = "outside",
    strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA)
  ) +
  labs(color = "Wing shape")

DP25

# Temp is 28 C 
range_sqrtRelTAlg <- range(T1_28$sqrtRelTAlg) # find the range to use in scale_fill_gradient

## code that controls the guide better !!!
DP28 <- ggplot(T1_28, aes(y = pup_mss, x = larv_time, color = sqrtRelTAlg, fill = sqrtES)) +
  geom_point(size = 3) + # Points
  geom_smooth(method = 'lm',  se = FALSE) + # Line mapping
  stat_summary_2d(
    aes(fill = ..value.., z = sqrtRelTAlg), 
    fun = "mean", bins = 5, geom = "tile", 
    color = NA, alpha = 0.5) + # Heatmap based on sqrtES values
  scale_color_gradientn(colors = c("lightblue", "black"), values = scales::rescale(c(0, 1)), guide = FALSE) + # Set color gradient without legend
  scale_fill_gradientn(colors = c("lightblue", "black"), values = scales::rescale(c(0, 1)), 
                       breaks = c(0.05, 0.1, 0.15, 0.20, 0.25), labels = c("0.05", "0.10", "0.15", "0.20", "0.25"),
                       guide = guide_colorbar(title = element_blank())) + # Set fill gradient with legend
  scale_alpha_continuous(range = c(0.3, 1), guide = FALSE) + # Set alpha range
  scale_x_continuous(limit = c(13, 23), breaks = c(15, 20)) + 
  scale_y_continuous(limit = c(0.3, 0.56), breaks = c(0.3, 0.4, 0.5)) + 
  facet_grid(cols = vars(sex)) +
  theme_classic() +
  xlab("Larval development time (days)") +
  ylab("Pupal mass (grams)") +
  theme(
    axis.title.x = element_text(color = "white", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    legend.direction = "vertical",
    legend.background = element_blank(),
    strip.text.x = element_text(color = "white", size = 16, face = "bold"),
    legend.text = element_text(color = "black", size = 13, face = "plain"),
    legend.title = element_text(color = "black", size = 13, face = "plain"),
    strip.placement = "outside",
    strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA)
  )

DP28

# Temp is 31 C 
range_sqrtRelTAlg <- range(T1_31$sqrtRelTAlg) # find the range to use in scale_fill_gradient

DP31 <- ggplot(T1_31, aes(y = pup_mss, x = larv_time, color = sqrtRelTAlg, fill = sqrtES)) +
  geom_point(size = 3) + # Points
  geom_smooth(method = 'lm',  se = FALSE) + # Line mapping
  stat_summary_2d(
    aes(fill = ..value.., z = sqrtRelTAlg), 
    fun = "mean", bins = 5, geom = "tile", 
    color = NA, alpha = 0.5) + # Heatmap based on sqrtES values
  scale_color_gradientn(colors = c("lightblue", "black"), values = scales::rescale(c(0, 1)), guide = FALSE) + # Set color gradient without legend
  scale_fill_gradientn(colors = c("lightblue", "black"), values = scales::rescale(c(0, 1)), 
                       breaks = c(0.05, 0.1, 0.15, 0.20, 0.25), labels = c("0.05", "0.10", "0.15", "0.20", "0.25"),
                       guide = guide_colorbar(title = element_blank())) + # Set fill gradient with legend
  scale_alpha_continuous(range = c(0.3, 1), guide = FALSE) + # Set alpha range
  scale_x_continuous(limit = c(13, 19), breaks = c(14, 16,18)) + 
  scale_y_continuous(limit = c(0.3, 0.56), breaks = c(0.3, 0.4, 0.5)) + 
  facet_grid(cols = vars(sex)) +
  theme_classic() +
  xlab("Larval development time (days)") +
  ylab("Pupal mass (grams)") +
  theme(
    axis.title.x = element_text(color = "black", size = 20, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    legend.direction = "vertical",
    legend.background = element_blank(),
    strip.text.x = element_text(color = "white", size = 16, face = "bold"),
    legend.text = element_text(color = "black", size = 13, face = "plain"),
    legend.title = element_text(color = "black", size = 13, face = "plain"),
    strip.placement = "outside",
    strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA)
  )

DP31




# combine the reaction norms for size at maturity for the Temperature experiment
# removing y-axis labels
DP19r <- DP19 + ylab("") + xlab("")
DP22r <- DP22 + ylab("") + xlab("") + labs(color = "") #+ theme(strip.text.x =  element_blank())
DP25r <- DP25 + ylab("") + xlab("") + labs(color = "") #+ theme(strip.text.x =  element_blank())
DP28r <- DP28 + ylab("") + xlab("") + labs(color = "") #+ theme(strip.text.x = element_blank())
DP31r <- DP31 + ylab("") + labs(color = "") #+ theme(strip.text.x = element_blank())

DP_T1 <- ggarrange(DP19r, DP22r, DP25r, DP28r, DP31r, ncol = 1, labels=c("19 C","22 C","25 C","28 C","31 C"), heights =c(11,11,11,11,11), font.label = list(size = 20), align = "v") + bgcolor("White")+ border(NA)
DP_T1

# give common y-axis name
DP_T1a <- annotate_figure(DP_T1 , left = text_grob("Pupal mass (grams)", color = "black", face="bold", family = "sans", size=20, rot = 90))+ bgcolor("White") + border("White")
DP_T1a


ggsave("SizeMaturirtRN_wingshape.png", DP_T1a , height=43, width=22, units=c("cm"), dpi=600) ### Figure A3.4b ###


#####################################

#### EYESPOT SIZE

##########################################


### growth rate vs eyespot size

## regressions for eyespot size  ######################

### Temperature 1 wing shape vs growth rate
# Temp1PPstacks_fore$Temp <- factor(T1$Temp, levels = rev(levels(Temp1PPstacks_fore$Temp))) # to make sure Temp treatments are in correct order

## using Temp as a fixed factor AND SLEEVE as a random effect
#T1$WAd <- ((T1$WA)/2000)  # to make scale of Wing area similar to that of grwth_rate

T1F_ES<-lmer(sqrtES ~ grwth_rate + Temp + WAd + sex  + (1|sleeve), data=T1)
summary(T1F_ES)
anova(T1F_ES)

T1mod<-lmer(sqrtES ~ larv_time + pup_mss + larv_time*pup_mss + Temp + sex  + (1|sleeve), data=T1)
summary(T1mod)
anova(T1mod)


## plot wing area against eyespot size
T1_WA_ES <- ggplot(T1, aes(x = WAd, y = sqrtES, color=Temp)) +
  geom_point(shape=19, size=3) +
  scale_color_manual(values=c("blue", "#0066FF", "orange", "red","darkred"))+
  geom_smooth(method='lm',se=FALSE) +
  facet_grid(cols = vars(sex)) +
  scale_y_continuous(limit = c(0.0, 0.27), breaks = c(0.0, 0.1, 0.2)) +
  #scale_x_continuous(limit = c(0.015, 0.055), breaks = c(0, 0.02, 0.04)) +
  theme_classic() +
  xlab("Wing Area Index") +
  ylab("Sqrt relative eyespot size") +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    legend.position = c(0.4, 0.95),
    legend.direction = "horizontal",
    legend.background = element_blank(),
    strip.text.x = element_text(color = "black", size = 16, face = "bold"),
    legend.text = element_text(color = "black", size = 13, face = "plain"),
    legend.title = element_blank(),
    strip.placement = "outside",
    strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA)
  )

T1_WA_ES

ggsave("EyespotsizeTemperature1GrowthRate.png", T1_WA_ES, height=15, width=20, units=c("cm"), dpi=600)


### Growth rate vs wing shape Temperature* Humidity

# TH$WAd <- ((TH$WA)/2000)  # to make scale of Wing area similar to that of grwth_rate

THF_ES<-lmer(sqrtES ~ grwth_rate + Temp + WAd + sex  + (1|sleeve), data=TH)
summary(THF_ES)
anova(THF_ES)

THreg <- ggplot(TH, aes(x = grwth_rate, y = sqrtRelTAlg, color=Temp, shape=Humidity)) +
  geom_point(size=3) +
  scale_color_manual(values=c("blue","orange","red","darkred"))+
  geom_smooth(method='lm',se=FALSE) +
  facet_grid(cols = vars(sex)) +
  scale_y_continuous(limit = c(0.05, 0.25), breaks = c(0.1, 0.2)) +
  scale_x_continuous(limit = c(0.015, 0.055), breaks = c(0, 0.02, 0.04)) +
  theme_classic() +
  xlab("Larval growth rate") +
  ylab("Forewing-tip shape") +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    legend.position = c(0.3, 0.95),
    legend.direction = "horizontal",
    legend.spacing = unit(-0.7, "cm"),
    legend.background = element_blank(),
    strip.text.x =element_blank(),
    legend.text = element_text(color = "black", size = 13, face = "plain"),
    legend.title = element_blank(),
    strip.placement = "outside",
    #strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.text.y =  element_blank(),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA)
  ) 

THreg


ggsave("GrowthRateForeWingShapeTemperatureHR_new.png", THreg, height=15, width=20, units=c("cm"), dpi=600)


## Plot Wing area to wing shape
THregWA <- ggplot(TH, aes(x = WA, y = sqrtRelTAlg, color=Temp, shape=Humidity)) +
  geom_point(size=3) +
  scale_color_manual(values=c("blue","orange","red","darkred"))+
  geom_smooth(method='lm',se=FALSE) +
  facet_grid(cols = vars(sex)) +
  scale_y_continuous(limit = c(0.05, 0.25), breaks = c(0.1, 0.2)) +
  #scale_x_continuous(limit = c(0.015, 0.055), breaks = c(0, 0.02, 0.04)) +
  theme_classic() +
  xlab("Wing area index") +
  ylab("Wing shape") +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    legend.position = c(0.3, 0.95),
    legend.direction = "horizontal",
    legend.spacing = unit(-0.7, "cm"),
    legend.background = element_blank(),
    strip.text.x =element_blank(),
    legend.text = element_text(color = "black", size = 13, face = "plain"),
    legend.title = element_blank(),
    strip.placement = "outside",
    #strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.text.y =  element_blank(),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA)
  ) 

THregWA

ggsave("EINGAREA_ForeWingShapeTemperatureHR.png", THregWA, height=15, width=20, units=c("cm"), dpi=600)

## Plant 4 growth rate vs wing shape

P4$WAd <- ((P4$WA)/2000)  # to make scale of Wing area similar to that of grwth_rate

P4FWS<-lmer(sqrtRelTAlg ~ grwth_rate + Temp + WAd + sex + (1|sleeve), data=TH)
summary(P4FWS)
anova(P4FWS)

Preg <- ggplot(P4, aes(x = grwth_rate, y = sqrtRelTAlg, color=plant)) +
  geom_point(shape=1, size=3) +
  scale_fill_brewer(palette="Set2")+
  #scale_color_manual(values=c("blue", "#0066FF", "orange", "red","darkred"))+
  geom_smooth(method='lm',se=FALSE)  +
  facet_grid(cols = vars(sex)) +
  scale_y_continuous(limit = c(0.05, 0.19), breaks = c(0.05,0.1,0.15)) +
  scale_x_continuous(limit = c(0.015, 0.055), breaks = c(0, 0.02, 0.04)) +
  theme_classic() +
  xlab("Larval growth rate") +
  ylab("Forewing-tip shape") +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    legend.text = element_text(color = "black", size = 13, face = "italic"),
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.background = element_blank(),
    legend.spacing.y = unit(-1.0, 'cm'),
    legend.spacing.x = unit(0.8, 'cm'),
    strip.text.x = element_blank(),
    legend.title = element_blank(),
    strip.placement = "outside",
    strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA))  +
  guides(fill = guide_legend(byrow = TRUE))

Preg


# plot wing area against wing shape for plants data
PregWA <- ggplot(P4, aes(x = WA, y = sqrtRelTAlg, color=plant)) +
  geom_point(shape=1, size=3) +
  scale_fill_brewer(palette="Set2")+
  #scale_color_manual(values=c("blue", "#0066FF", "orange", "red","darkred"))+
  geom_smooth(method='lm',se=FALSE)  +
  facet_grid(cols = vars(sex)) +
  scale_y_continuous(limit = c(0.05, 0.19), breaks = c(0.05,0.1,0.15)) +
  #scale_x_continuous(limit = c(0.015, 0.055), breaks = c(0, 0.02, 0.04)) +
  theme_classic() +
  xlab("Wing area index") +
  ylab("Forewing-tip shape") +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    legend.text = element_text(color = "black", size = 13, face = "italic"),
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.background = element_blank(),
    legend.spacing.y = unit(-1.0, 'cm'),
    legend.spacing.x = unit(0.8, 'cm'),
    strip.text.x = element_blank(),
    legend.title = element_blank(),
    strip.placement = "outside",
    strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA))  +
  guides(fill = guide_legend(byrow = TRUE))

PregWA

# removing y-axis labels
T1regrWA <- T1RegWA  + ylab("") + xlab("")
THregrWA <- THregWA + ylab("") + xlab("")
PregrWA <- PregWA + ylab("")

# combine Temp1, T*HR and Plant4 growth rate vs wing shape
WA_WS <- ggarrange(T1regrWA, THregrWA, PregrWA, ncol = 1, labels=c("a","b","c"), heights =c(12,10,17), font.label = list(size = 22), align = "v") + bgcolor("White") +border(NA)
WA_WS

# give common y-axis name
WA_WSa <- annotate_figure(WA_WS , left = text_grob("                 Forewing-tip shape", color = "black", face="bold", family = "sans", size=20, rot = 90))+ bgcolor("White") + border("White")
WA_WSa


ggsave("WINGAREAWingshapeCombi.png", WA_WSa , height=34, width=24, units=c("cm"), dpi=600) ### Figure A3.3. ###

######################################################################################################################################
#### Potting development time (x) pupal mass (y), and eyespot size (color gradient).
########################################################################################################################################

## Dev time by size reaction norm, Plotting pupal mass (y), development time (x) and eyespot size (color gradient) for Temperature experiment shape=19, 

ESPMDTT1 <- ggplot(T1, aes(y = pup_mss, x = larv_time, color = sqrtES, linetype = Temp)) +
  geom_point(aes(shape = Temp), size = 3, fill = "black") + # Filled black points
  geom_smooth(aes(group = Temp), method = 'lm', se = FALSE) + # Group by Temp for line mapping
  scale_linetype_manual(values = c("solid", "dashed", "dotted", "dotdash", "longdash")) + # Specify linetypes for Temp levels
  scale_shape_manual(values = c(17, 16, 15, 18, 6)) + # Filled shapes: circle, triangle (up and down facing), square, and diamond
  scale_color_gradient(low = "lightblue", high = "black") + # Increase color contrast from light blue to black
  facet_grid(cols = vars(sex)) +
  theme_classic() +
  xlab("Larval development time (days)") +
  ylab("Pupal mass (grams)") +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    legend.direction = "vertical",
    legend.background = element_blank(),
    strip.text.x = element_text(color = "black", size = 16, face = "bold"),
    legend.text = element_text(color = "black", size = 13, face = "plain"),
    legend.title = element_text(color = "black", size = 13, face = "plain"),
    strip.placement = "outside",
    strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA)
  )+
  labs(color = "Eyespot size")

ESPMDTT1


ggsave("Devtime_Pupmass_EysepotSize_Temperature1.png", ESPMDTT1, height=18, width=40, units=c("cm"), dpi=600)

##%%%### 
#    Make plots for each treatment separately, faceted by sex
## data were subsetted before

# Temp is 19 C

range_sqrtES <- range(T1_19$sqrtES) # find the range to use in scale_fill_gradient

DP19ES <- ggplot(T1_19, aes(y = pup_mss, x = larv_time, color = sqrtES)) +
  geom_point(size = 3) + # Points
  geom_smooth(method = 'lm',  se = FALSE) + # Line mapping
  stat_summary_2d(
    aes(fill = ..value.., z = sqrtES), 
    fun = "mean", bins = 5, geom = "tile", 
    color = NA, alpha = 0.5) + # Heatmap based on sqrtRelTAlg values
  scale_color_gradient(low = "lightblue", high = "black") + # Increase color contrast from light blue to black
  scale_fill_gradient(low = "lightblue", high = "black", guide = FALSE, limits = c(0.026, 0.081)) + # Adjust fill scale limits
  scale_alpha_continuous(range = c(0.3, 1), guide = FALSE) + # Set alpha range
  scale_x_continuous(limit = c(35, 55), breaks = c(40, 50)) + 
  facet_grid(cols = vars(sex)) +
  theme_classic() +
  xlab("Larval development time (days)") +
  ylab("Pupal mass (grams)") +
  theme(
    axis.title.x = element_text(color = "black", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    legend.direction = "vertical",
    legend.background = element_blank(),
    strip.text.x = element_text(color = "black", size = 16, face = "bold"),
    legend.text = element_text(color = "black", size = 13, face = "plain"),
    legend.title = element_text(color = "black", size = 13, face = "plain"),
    strip.placement = "outside",
    strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA)
  ) +
  labs(color = "Eyespot size")

DP19ES


# Temp is 22 C (for methods loess span = 3,)
range_sqrtES <- range(T1_22$sqrtES) # find the range to use in scale_fill_gradient

DP22ES <- ggplot(T1_22, aes(y = pup_mss, x = larv_time, color = sqrtES)) +
  geom_point(size = 3) + # Points
  geom_smooth(method = 'lm',  se = FALSE) + # Line mapping
  stat_summary_2d(
    aes(fill = ..value.., z = sqrtES), 
    fun = "mean", bins = 5, geom = "tile", 
    color = NA, alpha = 0.5) + # Heatmap based on sqrtES values
  scale_color_gradient(low = "lightblue", high = "black") + # Increase color contrast from light blue to black
  scale_fill_gradient(low = "lightblue", high = "black", guide = FALSE, limits = c(0.095, 0.261)) + # Adjust fill scale limits
  scale_alpha_continuous(range = c(0.3, 1), guide = FALSE) + # Set alpha range
  facet_grid(cols = vars(sex)) +
  theme_classic() +
  xlab("Larval development time (days)") +
  ylab("Pupal mass (grams)") +
  theme(
    axis.title.x = element_text(color = "white", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    legend.direction = "vertical",
    legend.background = element_blank(),
    strip.text.x = element_text(color = "white", size = 16, face = "bold"),
    legend.text = element_text(color = "black", size = 13, face = "plain"),
    legend.title = element_text(color = "black", size = 13, face = "plain"),
    strip.placement = "outside",
    strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA)
  ) +
  labs(color = "Eyespot size")

DP22ES

# Temp is 25 C 
range_sqrtES <- range(T1_25$sqrtES) # find the range to use in scale_fill_gradient

DP25ES <- ggplot(T1_25, aes(y = pup_mss, x = larv_time, color = sqrtES)) +
  geom_point(size = 3) + # Points
  geom_smooth(method = 'lm',  se = FALSE) + # Line mapping
  stat_summary_2d(
    aes(fill = ..value.., z = sqrtES), 
    fun = "mean", bins = 5, geom = "tile", 
    color = NA, alpha = 0.5) + # Heatmap based on sqrtES values
  scale_color_gradient(low = "lightblue", high = "black") + # Increase color contrast from light blue to black
  scale_fill_gradient(low = "lightblue", high = "black", guide = FALSE, limits = c(0.036, 0.224)) + # Adjust fill scale limits
  scale_alpha_continuous(range = c(0.3, 1), guide = FALSE) + # Set alpha range
  scale_x_continuous(limit = c(20, 37), breaks = c(25, 35)) + 
  facet_grid(cols = vars(sex)) +
  theme_classic() +
  xlab("Larval development time (days)") +
  ylab("Pupal mass (grams)") +
  theme(
    axis.title.x = element_text(color = "white", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    legend.direction = "vertical",
    legend.background = element_blank(),
    strip.text.x = element_text(color = "white", size = 16, face = "bold"),
    legend.text = element_text(color = "black", size = 13, face = "plain"),
    legend.title = element_text(color = "black", size = 13, face = "plain"),
    strip.placement = "outside",
    strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA)
  ) +
  labs(color = "Eyespot size")

DP25ES

# Temp is 28 C 
range_sqrtES <- range(T1_28$sqrtES) # find the range to use in scale_fill_gradient

## code that controls the guide better !!!
DP28ES <- ggplot(T1_28, aes(y = pup_mss, x = larv_time, color = sqrtES, fill = sqrtES)) +
  geom_point(size = 3) + # Points
  geom_smooth(method = 'lm',  se = FALSE) + # Line mapping
  stat_summary_2d(
    aes(fill = ..value.., z = sqrtES), 
    fun = "mean", bins = 5, geom = "tile", 
    color = NA, alpha = 0.5) + # Heatmap based on sqrtES values
  scale_color_gradientn(colors = c("lightblue", "black"), values = scales::rescale(c(0, 1)), guide = FALSE) + # Set color gradient without legend
  scale_fill_gradientn(colors = c("lightblue", "black"), values = scales::rescale(c(0, 1)), 
                       breaks = c(0.15, 0.20, 0.25), labels = c("0.15", "0.20", "0.25"),
                       guide = guide_colorbar(title = element_blank())) + # Set fill gradient with legend
  scale_alpha_continuous(range = c(0.3, 1), guide = FALSE) + # Set alpha range
  scale_x_continuous(limit = c(13, 23), breaks = c(15, 20)) + 
  scale_y_continuous(limit = c(0.25, 0.55), breaks = c(0.30, 0.40, 0.50)) + 
  facet_grid(cols = vars(sex)) +
  theme_classic() +
  xlab("Larval development time (days)") +
  ylab("Pupal mass (grams)") +
  theme(
    axis.title.x = element_text(color = "white", size = 16, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    legend.direction = "vertical",
    legend.background = element_blank(),
    strip.text.x = element_text(color = "white", size = 16, face = "bold"),
    legend.text = element_text(color = "black", size = 13, face = "plain"),
    legend.title = element_text(color = "black", size = 13, face = "plain"),
    strip.placement = "outside",
    strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA)
  )

DP28ES


# Temp is 31 C 
range_sqrtES <- range(T1_31$sqrtES) # find the range to use in scale_fill_gradient

DP31ES <- ggplot(T1_31, aes(y = pup_mss, x = larv_time, color = sqrtES)) +
  geom_point(size = 3) + # Points
  geom_smooth(method = 'lm',  se = FALSE) + # Line mapping
  stat_summary_2d(
    aes(fill = ..value.., z = sqrtES), 
    fun = "mean", bins = 5, geom = "tile", 
    color = NA, alpha = 0.5) + # Heatmap based on sqrtES values
  scale_color_gradient(low = "lightblue", high = "black") + # Increase color contrast from light blue to black
  scale_fill_gradient(low = "lightblue", high = "black", guide = FALSE, limits = c(0.133, 0.242)) + # Adjust fill scale limits
  scale_alpha_continuous(range = c(0.3, 1), guide = FALSE) + # Set alpha range
  scale_x_continuous(limit = c(13, 21), breaks = c(15, 20)) + 
  scale_y_continuous(limit = c(0.25, 0.55), breaks = c(0.30, 0.40, 0.50)) +
  facet_grid(cols = vars(sex)) +
  theme_classic() +
  xlab("Larval development time (days)") +
  ylab("Pupal mass (grams)") +
  theme(
    axis.title.x = element_text(color = "black", size = 20, face = "bold"),
    axis.title.y = element_text(color = "black", size = 16, face = "bold"),
    axis.text.x = element_text(color = "black", size = 13, face = "plain"),
    axis.text.y = element_text(color = "black", size = 13, face = "plain"),
    legend.key.size = unit(3, "line"),
    legend.direction = "vertical",
    legend.background = element_blank(),
    strip.text.x = element_text(color = "white", size = 16, face = "bold"),
    legend.text = element_text(color = "black", size = 13, face = "plain"),
    legend.title = element_text(color = "black", size = 13, face = "plain"),
    strip.placement = "outside",
    strip.text.y = element_text(size = 16, color = "black", face = "bold"),
    strip.background = element_rect(fill = NA, colour = NA),
    panel.border = element_rect(color = "white", fill = NA)
  ) +
  labs(color = "Eyespot size")

DP31ES



# combine these reaction norms for size at maturity for Temperature experiment
# removing y-axis labels
DP19ESr <- DP19ES + ylab("") + xlab("")
DP22ESr <- DP22ES + ylab("") + xlab("") + labs(color = "") #+ theme(strip.text.x =  element_blank())
DP25ESr <- DP25ES + ylab("") + xlab("") + labs(color = "") #+ theme(strip.text.x =  element_blank())
DP28ESr <- DP28ES + ylab("") + xlab("")  # + theme(strip.text.x = element_blank())
DP31ESr <- DP31ES + ylab("") + labs(color = "") #+ theme(strip.text.x = element_blank())

DP_ES_T1 <- ggarrange(DP19ESr, DP22ESr, DP25ESr, DP28ESr, DP31ESr, ncol = 1, labels=c("19 C","22 C","25 C","28 C","31 C"), heights =c(11,11,11,11,11), font.label = list(size = 20), align = "v") + bgcolor("White")+ border(NA)
DP_ES_T1

# give common y-axis name
DP_ES_T1a <- annotate_figure(DP_ES_T1 , left = text_grob("Pupal mass (grams)", color = "black", face="bold", family = "sans", size=20, rot = 90))+ bgcolor("White") + border("White")
DP_ES_T1a


ggsave("SizeMaturirtRN_EysespotSize.png", DP_ES_T1a , height=43, width=22, units=c("cm"), dpi=600) ### Figure A3.3a ###


### Combine eyespot size and wing shape for Age and Size at Maturity
DP_ES_T1ar <- DP_ES_T1a  + xlab("")
DP_T1r <- DP_T1 + ylab("") + xlab("") 

ASM<- ggarrange(DP_ES_T1a, DP_T1, ncol = 2, heights =c(31,31), font.label = list(size = 20), align = "h") + bgcolor("White")+ border(NA)
ASM

# give common x-axis name
ASMa <- annotate_figure(ASM , bottom = text_grob("Larval development time (days)", color = "black", face="bold", family = "sans", size=20, rot = 0))+ bgcolor("White") + border("White")
ASMa

ggsave("AgeSizeMaturirtRN_EysespotSize&WingShape.png", ASM , height=43, width=30, units=c("cm"), dpi=600) #### Figure A3.4. ###



#### Counting sleeves per treatment #####


# Group by Temp and sex, and count unique sleeves
sleeves_counts <- P4c %>%
  group_by(plant, sex) %>%
  summarise(Number_of_Sleeves = n_distinct(sleeve))

# Display the result
sleeves_counts

## counting surviving butterflies per sleeve

# Group by sleeve and count records of grwth_rate
grwth_rate_counts <- T1c %>%
  group_by(sleeve) %>%
  summarise(Number_of_records = n())

# Display the result
grwth_rate_counts

# Calculate the average number of records per sleeve
average_records_per_sleeve <- mean(grwth_rate_counts$Number_of_records)

# Display the result
average_records_per_sleeve