---
title: "Identification of proximal cue(s) for pupal polymorphism in Chlosyne lacinia"
output: html_document
date: "2024-03-12"
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```



load dataset 

```{r}
cc.df <- read.csv("pupal color master sheet w charcter vectors.csv") 
cc.df
attach(cc.df)
```


graphs for proporation of blackness and eclosion time


```{r} 
library(ggplot2)

#proportion of blackness graph only treatment and sex
ggplot(data = cc.df, aes(x = sex, y = X..black))+ 
  geom_boxplot(fill = c("pink","lightblue","pink","lightblue"))+ 
  theme_bw()+ 
  facet_wrap(Treatment ~ .)+
  xlab("Treatment Color")+ 
  ylab("Percentage of Blackness (%)")+ 
  ggtitle("The proportion of blackness in experimentally manipulated C. lacinia pupae") 

#eclosion time vs treatment and sex
ggplot(data = cc.df, aes(x = sex, y = pupal.duration))+ 
  geom_boxplot(fill = c("pink","lightblue","pink","lightblue"))+ 
  theme_bw()+ 
  facet_wrap(Treatment ~ .)+
  xlab("Treatment Color")+ 
  ylab("Time to Eclosion (days)")+ 
  ggtitle("Eclosion time in black and white C. lacinia pupae") 

#larval color
ggplot(data = cc.df, aes(x = larval.color, y = pupal.duration))+ 
  geom_boxplot(fill=c("maroon","darkgray","red"))+ 
  theme_bw()+ 
  xlab("larval morph")+ 
  ylab("Pupal Duration (days)")+ 
  ggtitle("The time to eclosion in C. lacinia pupae of different larval morphs") 

```

data analysis 

```{r}
#model selection for proportion of blackness in both treatments 
library(AICcmodavg) 

blm1 <- lm(X..black~Treatment) 
summary(blm1) 
#pvalue = 5.07e-16 

blm2 <- lm(X..black ~ sex) 
summary(blm2)
#pvalue = 0.7911 

blm3 <- lm(X..black~Treatment*sex) 
summary(blm3) 
#pvalue = 3.295e-14 

blm4 <- lm(X..black~pupal.duration) 
summary(blm4) 
#pvalue = 0.02842 


blm5 <- lm(X..black~Treatment*sex*pupal.duration) 
summary(blm5) 
#pvalue = 3.212e-10

blm6 <- lm(X..black~sex*pupal.duration) 
summary(blm6)
#pvalue = 0.1776

blm7 <- lm(X..black~larval.color) 
summary(blm7) 
#pvalue = 0.6003

blm8 <- lm(X..black~Treatment*sex*pupal.duration*larval.color) 
summary(blm8) 
#pvalue = 2.23e-05

blm9 <- lm(X..black~sex*pupal.duration*larval.color) 
summary(blm9) 
#pvalue = 0.5969

blm10 <- lm(X..black~pupal.duration*larval.color)
summary(blm10) 
#pvalue = 0.2074

#define list of models
b.models <- list(blm1,blm2,blm3,blm4,blm5,blm6,blm7,blm8,blm9,blm10)

#specify model names
b.mod.names <- c('1','2','3','4','5','6','7','8','9','10') 

#calculate AIC of each model
bAIC <- aictab(cand.set = b.models, modnames = b.mod.names) 
bAIC 

summary(bAIC) 

################# 
# model average # 
################# 
library(MuMIn)

#models blm1 & blm3    

b.avg <- model.avg(blm1,blm3, revised.var=TRUE)
summary(b.avg) 

plot(b.avg)

#separating out black only too see relationship between blackness and eclosion time
library(dplyr) 

c.black.df <- cc.df[grep("black", cc.df$Treatment), ] #character sex 
c.white.df <- cc.df[grep("white", cc.df$Treatment), ]

#mean percent of melanization for each treatment
black.mean <- mean(c.black.df$X..black) 
#87.95748
white.mean <- mean(c.white.df$X..black)
#17.49526 

#separating sexes in black pupae to see mean of melanization 
f.b.df <- c.black.df[grep("F", c.black.df$sex), ] 
m.b.df <- c.black.df[grep("M", c.black.df$sex), ] 

f.b.mean <- mean(f.b.df$X..black) 
#84.46077
m.b.mean <- mean(m.b.df$X..black)
#93.20255


b.ttest <- t.test(c.black.df$X..black ~ c.black.df$sex) 
b.ttest 
#pvalue = 0.03911 

b.o.lm1 <- lm(c.black.df$pupal.duration ~ c.black.df$sex) 
summary(b.o.lm1) 
#pvalue = 0.3708 

b.o.lm2 <- lm(c.black.df$pupal.duration ~ c.black.df$X..black) 
summary(b.o.lm2) 
#pvalue = 0.8723 

b.o.lm3 <- lm(c.black.df$pupal.duration ~ c.black.df$larval.color) 
summary(b.o.lm3) 
#pvalue = 0.2662 

b.o.lm4 <- lm(c.black.df$pupal.duration ~ c.black.df$sex*c.black.df$X..black) 
summary(b.o.lm4) 
#pvalue = 0.6203 

b.o.lm5 <- lm(c.black.df$pupal.duration ~ c.black.df$sex*c.black.df$larval.color) 
summary(b.o.lm5) 
#pvalue = 0.5979 

b.o.lm6 <- lm(c.black.df$pupal.duration ~ c.black.df$larval.color*c.black.df$X..black) 
summary(b.o.lm6) 
#pvalue = 0.6179 

b.o.lm7 <- lm(c.black.df$pupal.duration ~ c.black.df$sex*c.black.df$X..black*c.black.df$larval.color) 
summary(b.o.lm7) 
#pvalue = 0.2819

w.ttest <- t.test(c.white.df$X..black ~ c.white.df$sex) 
w.ttest 
#pvalue = 0.5061

w.o.lm1 <- lm(c.white.df$pupal.duration ~ c.white.df$sex) 
summary(w.o.lm1) 
#pvalue = 0.639

w.o.lm2 <- lm(c.white.df$pupal.duration ~ c.white.df$X..black) 
summary(w.o.lm2) 
#pvalue = 0.2474

w.o.lm3 <- lm(c.white.df$pupal.duration ~ c.white.df$larval.color) 
summary(w.o.lm3) 
#pvalue = 0.6224 

w.o.lm4 <- lm(c.white.df$pupal.duration ~ c.white.df$sex*c.white.df$X..black) 
summary(w.o.lm4) 
#pvalue = 0.5847 

w.o.lm5 <- lm(c.white.df$pupal.duration ~ c.white.df$sex*c.white.df$larval.color) 
summary(w.o.lm5) 
#pvalue = 0.3873 

w.o.lm6 <- lm(c.white.df$pupal.duration ~ c.white.df$larval.color*c.white.df$X..black) 
summary(w.o.lm6) 
#pvalue = 0.8418

w.o.lm7 <- lm(c.white.df$pupal.duration ~ c.white.df$sex*c.white.df$X..black*c.white.df$larval.color) 
summary(w.o.lm7) 
#pvalue = 0.4007

eclosionaov <- aov(pupal.duration ~ Treatment) 
summary(eclosionaov) 
#pvalue = 0.0198 

eclosionlm <- lm(pupal.duration ~ Treatment) 
summary(eclosionlm) 
#pvalue = 0.01984  

#mean pupal duration for each treatment
black.eclosion.mean <- mean(c.black.df$pupal.duration) 
#5.5
white.eclosion.mean <- mean(c.white.df$pupal.duration)
#4.33333


sex.e.ttest <- t.test(pupal.duration ~ sex) 
sex.e.ttest 
#pvalue = 0.6865

#for AIC 
elm1 <- lm(pupal.duration ~ X..black) 
summary(elm1) 
#pvalue = 0.02842

elm2 <- lm(pupal.duration ~ sex) 
summary(elm2) 
#pvalue = 0.6701

elm3 <- lm(pupal.duration ~ X..black*sex) 
summary(elm3) 
#pvalue = 0.1302 

elm4 <- lm(pupal.duration ~ larval.color) 
summary(elm4) 
#pvalue = 0.5476

elm5 <- lm(pupal.duration ~ X..black*sex*larval.color)  
summary(elm5)
#pvalue = 0.2939

elm6 <- lm(pupal.duration ~ sex*larval.color) 
summary(elm6) 
#pvalue = 0.5562 

elm7 <- lm(pupal.duration ~ X..black*larval.color)
summary(elm7) 
#pvalue = 0.1135

#define list of models
e.models <- list(elm1,elm2,elm3,elm4,elm5, elm6, elm7)

#specify model names
e.mod.names <- c('1','2','3','4','5','6','7') 

#calculate AIC of each model
eAIC <- aictab(cand.set = e.models, modnames = e.mod.names) 
eAIC
```

