Load library

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.1     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.2     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors

library(lme4)

## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## 
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack

library(lmerTest)

## 
## Attaching package: 'lmerTest'
## 
## The following object is masked from 'package:lme4':
## 
##     lmer
## 
## The following object is masked from 'package:stats':
## 
##     step

library(emmeans)
library(plyr)

## ------------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## ------------------------------------------------------------------------------
## 
## Attaching package: 'plyr'
## 
## The following objects are masked from 'package:dplyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize
## 
## The following object is masked from 'package:purrr':
## 
##     compact

library(MuMIn)
library(performance)
library(ggstatsplot)

## You can cite this package as:
##      Patil, I. (2021). Visualizations with statistical details: The 'ggstatsplot' approach.
##      Journal of Open Source Software, 6(61), 3167, doi:10.21105/joss.03167

library(effectsize)

Set seed to 1234

set.seed(1234)

Data<- read.csv("Rotating_diet.csv")
head(Data)

##   year  Cat      host site.lab Mat..line Mother Father     eggs J_Eggs
## 1 2022 1431 B. cherry      Lab      143A    30b      6 6/9/2022    160
## 2 2022 1428 B. cherry      Lab      143A    30b      6 6/9/2022    160
## 3 2022 1438 B. cherry      Lab      143A    30b      6 6/9/2022    160
## 4 2022  893 B. cherry      Lab      136A      4      5 6/4/2022    155
## 5 2022 1316 B. cherry      Lab      138A      7    30b 6/7/2022    158
## 6 2022 1141 B. cherry      Lab      139A      5     15 6/5/2022    156
##   hatch.col J_Hatch Sex Color Type Pupa.mass Pupa.mass2 pupa.date J_Pupa
## 1 6/23/2022     174   f     B   FB    0.0864       85.4 7/30/2022    211
## 2 6/23/2022     174   f     B   FB    0.0875       87.0 7/30/2022    211
## 3 6/23/2022     174   f     B   FB    0.1012      100.9 7/30/2022    211
## 4 6/15/2022     166   f     B   FB    0.1101      110.0 7/21/2022    202
## 5 6/18/2022     169   f     B   FB    0.1166      116.3 7/27/2022    208
## 6 6/17/2022     168   f     B   FB    0.1166      116.3 7/30/2022    211
##   Pupation Mortality Survive Count lost death.date J_Death Notes
## 1        1         0     yes   266                      NA      
## 2        1         0     yes   266                      NA      
## 3        1         0     yes   266                      NA      
## 4        1         0     yes   266                      NA      
## 5        1         0     yes   266                      NA      
## 6        1         0     yes   266                      NA      
##   Genetic.analysis Death.host Dev_time    Log_T Mass_per_day
## 1               NA                  37 1.568202  0.000409479
## 2               NA                  37 1.568202  0.000414692
## 3               NA                  37 1.568202  0.000479621
## 4               NA                  36 1.556303  0.000545050
## 5               NA                  39 1.591065  0.000560577
## 6               NA                  43 1.633468  0.000552607

clean data by removing empty or missing data points, also removing none complete lineages

Data$Pupa.mass2<-as.numeric(Data$Pupa.mass2)

data.sex<- Data[Data$Sex != "",] #remove all rows missing sex

data.color<-data.sex[data.sex$Color !="",] #remove all rows missing color

data.pupa<-data.color[data.color$Pupation !="0",] #remove all rows of non-pupating individuals

data.mort<-data.pupa[data.pupa$Mortality !="1",] # remove all rows of individuals that died

data.host<-data.mort[data.mort$host !="",] # removing all rows missing host data

data.1<- data.host[data.host$Pupa.mass2 !="missing",]
data.1<-data.1[data.1$Pupa.mass != "",]


data <- data.1[data.1$pupa.date !="missing",]

data<-data[data$Mat..line != "138A",] #remove matrile lines that are not complete
data<-data[data$Mat..line != "143A",]#remove matrile lines that are not complete

head(data)

##    year  Cat      host site.lab Mat..line Mother Father     eggs J_Eggs
## 4  2022  893 B. cherry      Lab      136A      4      5 6/4/2022    155
## 6  2022 1141 B. cherry      Lab      139A      5     15 6/5/2022    156
## 10 2022  962 B. cherry      Lab      122B      5      4 6/3/2022    154
## 11 2022 1239 B. cherry      Lab      145A      6    18c 6/7/2022    158
## 12 2022  964 B. cherry      Lab      122B      5      4 6/3/2022    154
## 14 2022  961 B. cherry      Lab      122B      5      4 6/3/2022    154
##    hatch.col J_Hatch Sex Color Type Pupa.mass Pupa.mass2 pupa.date J_Pupa
## 4  6/15/2022     166   f     B   FB    0.1101      110.0 7/21/2022    202
## 6  6/17/2022     168   f     B   FB    0.1166      116.3 7/30/2022    211
## 10 6/16/2022     167   f     B   FB    0.1271      127.0 7/22/2022    203
## 11 6/18/2022     169   f     B   FB    0.1293      128.8 7/25/2022    206
## 12 6/16/2022     167   f     B   FB    0.1316      131.2 7/26/2022    207
## 14 6/16/2022     167   f     B   FB    0.1366      136.0 7/20/2022    201
##    Pupation Mortality Survive Count lost death.date J_Death Notes
## 4         1         0     yes   266                      NA      
## 6         1         0     yes   266                      NA      
## 10        1         0     yes   266                      NA      
## 11        1         0     yes   266                      NA      
## 12        1         0     yes   266                      NA      
## 14        1         0     yes   266                      NA      
##    Genetic.analysis Death.host Dev_time    Log_T Mass_per_day
## 4                NA                  36 1.556303  0.000545050
## 6                NA                  43 1.633468  0.000552607
## 10               NA                  36 1.556303  0.000626108
## 11               NA                  37 1.568202  0.000627670
## 12               NA                  40 1.602060  0.000635749
## 14               NA                  34 1.531479  0.000679602

table(data$Sex, data$Color)

##    
##       B   R
##   f 483 168
##   m 544 136

development time analysis

data.pupa comes from data tab

mass<- data


mass<- mass[mass$Pupa.mass2>40.0,] #removing dead pupa that are outliers

head(mass)

##    year  Cat      host site.lab Mat..line Mother Father     eggs J_Eggs
## 4  2022  893 B. cherry      Lab      136A      4      5 6/4/2022    155
## 6  2022 1141 B. cherry      Lab      139A      5     15 6/5/2022    156
## 10 2022  962 B. cherry      Lab      122B      5      4 6/3/2022    154
## 11 2022 1239 B. cherry      Lab      145A      6    18c 6/7/2022    158
## 12 2022  964 B. cherry      Lab      122B      5      4 6/3/2022    154
## 14 2022  961 B. cherry      Lab      122B      5      4 6/3/2022    154
##    hatch.col J_Hatch Sex Color Type Pupa.mass Pupa.mass2 pupa.date J_Pupa
## 4  6/15/2022     166   f     B   FB    0.1101      110.0 7/21/2022    202
## 6  6/17/2022     168   f     B   FB    0.1166      116.3 7/30/2022    211
## 10 6/16/2022     167   f     B   FB    0.1271      127.0 7/22/2022    203
## 11 6/18/2022     169   f     B   FB    0.1293      128.8 7/25/2022    206
## 12 6/16/2022     167   f     B   FB    0.1316      131.2 7/26/2022    207
## 14 6/16/2022     167   f     B   FB    0.1366      136.0 7/20/2022    201
##    Pupation Mortality Survive Count lost death.date J_Death Notes
## 4         1         0     yes   266                      NA      
## 6         1         0     yes   266                      NA      
## 10        1         0     yes   266                      NA      
## 11        1         0     yes   266                      NA      
## 12        1         0     yes   266                      NA      
## 14        1         0     yes   266                      NA      
##    Genetic.analysis Death.host Dev_time    Log_T Mass_per_day
## 4                NA                  36 1.556303  0.000545050
## 6                NA                  43 1.633468  0.000552607
## 10               NA                  36 1.556303  0.000626108
## 11               NA                  37 1.568202  0.000627670
## 12               NA                  40 1.602060  0.000635749
## 14               NA                  34 1.531479  0.000679602

Subset females by color

Female<-subset(mass, Sex == "f")
Fem.b<- subset(Female, Color == "B")
Fem.r<- subset(Female, Color == "R")

subset males by color

Male<-subset(mass, Sex == "m")
mal.b<- subset(Male, Color == "B")
mal.r<- subset(Male, Color == "R")

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

Black females

subset black females by host

print("       black female pupal mass  (mg)        ")

## [1] "       black female pupal mass  (mg)        "

print(".... black cherry [1] mean, [2]SD, [3] count, [4] SE   ......")

## [1] ".... black cherry [1] mean, [2]SD, [3] count, [4] SE   ......"

f.b.pru<-subset(Fem.b, host == "B. cherry")
          mean(f.b.pru$Pupa.mass2, na.rm = TRUE)

## [1] 175.1129

          sd(f.b.pru$Pupa.mass2, na.rm = TRUE)

## [1] 28.21661

          length(f.b.pru$Pupa.mass2)

## [1] 93

          sd(f.b.pru$Pupa.mass2, na.rm = TRUE)/sqrt(length(f.b.pru$Pupa.mass2))

## [1] 2.925926

 print(".... elm [1] mean, [2]SD, [3] count, [4] SE   ......")

## [1] ".... elm [1] mean, [2]SD, [3] count, [4] SE   ......"

f.b.ulm<-subset(Fem.b, host == "Elm")
          mean(f.b.ulm$Pupa.mass2, na.rm = TRUE)

## [1] 159.7584

           sd(f.b.ulm$Pupa.mass2)

## [1] 26.53074

           length(f.b.ulm$Pupa.mass2)

## [1] 101

           sd(f.b.ulm$Pupa.mass2, na.rm = TRUE)/sqrt(length(f.b.ulm$Pupa.mass2))

## [1] 2.639907

print(".... hickory [1] mean, [2]SD, [3] count, [4] SE   ......")

## [1] ".... hickory [1] mean, [2]SD, [3] count, [4] SE   ......"

f.b.car<-subset(Fem.b, host == "Hickory")
            mean(f.b.car$Pupa.mass2, na.rm = TRUE)

## [1] 152.997

            sd(f.b.car$Pupa.mass2)

## [1] 22.61446

            length(f.b.car$Pupa.mass2)

## [1] 99

            sd(f.b.car$Pupa.mass2, na.rm = TRUE)/sqrt(length(f.b.car$Pupa.mass2))

## [1] 2.272839

print(".... apple [1] mean, [2]SD, [3] count, [4] SE   ......")

## [1] ".... apple [1] mean, [2]SD, [3] count, [4] SE   ......"

f.b.mal<-subset(Fem.b, host =="Malus")
           mean(f.b.mal$Pupa.mass2, na.rm = TRUE)

## [1] 182.6623

           sd(f.b.mal$Pupa.mass2)

## [1] 26.51229

           length(f.b.mal$Pupa.mass2)

## [1] 106

           sd(f.b.mal$Pupa.mass2, na.rm = TRUE)/sqrt(length(f.b.mal$Pupa.mass2))

## [1] 2.575101

print(".... rotating [1] mean, [2]SD, [3] count, [4] SE   ......")

## [1] ".... rotating [1] mean, [2]SD, [3] count, [4] SE   ......"

f.b.mix<-subset(Fem.b, host == "Mix")
           mean(f.b.mix$Pupa.mass2, na.rm = TRUE)

## [1] 161.9214

           sd(f.b.mix$Pupa.mass2)

## [1] 28.70121

           length(f.b.mix$Pupa.mass2)

## [1] 84

           sd(f.b.mix$Pupa.mass2, na.rm = TRUE)/sqrt(length(f.b.mix$Pupa.mass2))

## [1] 3.131559

Linear mixed model– pupa mass (black females)

"female black Pupal mass"

## [1] "female black Pupal mass"

mass.female.b<- lmer(Pupa.mass2~host+(1|Mother)+(1|Father), data=Fem.b)
summary(mass.female.b)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Pupa.mass2 ~ host + (1 | Mother) + (1 | Father)
##    Data: Fem.b
## 
## REML criterion at convergence: 4387.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.4422 -0.5858  0.0251  0.6370  2.6205 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Father   (Intercept) 290.41   17.041  
##  Mother   (Intercept)  59.99    7.745  
##  Residual             504.76   22.467  
## Number of obs: 483, groups:  Father, 9; Mother, 6
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)  175.599      7.020  12.945  25.015 2.41e-12 ***
## hostElm      -17.507      3.267 467.888  -5.358 1.32e-07 ***
## hostHickory  -23.427      3.253 466.991  -7.202 2.38e-12 ***
## hostMalus      4.610      3.219 467.610   1.432    0.153    
## hostMix      -14.324      3.434 468.993  -4.171 3.61e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) hstElm hstHck hstMls
## hostElm     -0.241                     
## hostHickory -0.243  0.515              
## hostMalus   -0.248  0.524  0.524       
## hostMix     -0.223  0.504  0.489  0.499

anova(mass.female.b)

## Type III Analysis of Variance Table with Satterthwaite's method
##      Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## host  56223   14056     4 467.95  27.846 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

random effect (black females)

print("BLACK FEMALES")

## [1] "BLACK FEMALES"

print("parental effect")

## [1] "parental effect"

r.squaredGLMM(mass.female.b)

## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.

##            R2m       R2c
## [1,] 0.1215493 0.4814904

(0.481-0.119)*100

## [1] 36.2

host effect size (black females)

F_to_eta2(27.327, 4, 469)

## Eta2 (partial) |       95% CI
## -----------------------------
## 0.19           | [0.13, 1.00]
## 
## - One-sided CIs: upper bound fixed at [1.00].

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

Black males

subset black males by host

print("       black male pupal mass  (mg)        ")

## [1] "       black male pupal mass  (mg)        "

print(".... black cherry [1] mean, [2]SD, [3] count, [4] SE  ......")

## [1] ".... black cherry [1] mean, [2]SD, [3] count, [4] SE  ......"

m.b.pru<-subset(mal.b, host == "B. cherry")
       mean(m.b.pru$Pupa.mass2)

## [1] 144.4331

        sd(m.b.pru$Pupa.mass2)

## [1] 21.34244

        length(m.b.pru$Pupa.mass2)

## [1] 124

        sd(m.b.pru$Pupa.mass2)/sqrt(length(m.b.pru$Pupa.mass2))

## [1] 1.916607

print(".... elm [1] mean, [2]SD, [3] count, [4] SE  ......")

## [1] ".... elm [1] mean, [2]SD, [3] count, [4] SE  ......"

m.b.ulm<-subset(mal.b, host == "Elm")
      mean(m.b.ulm$Pupa.mass2)

## [1] 130.6204

       sd(m.b.ulm$Pupa.mass2)

## [1] 17.10214

       length(m.b.ulm$Pupa.mass2)

## [1] 108

       sd(m.b.ulm$Pupa.mass2)/sqrt(length(m.b.ulm$Pupa.mass2))

## [1] 1.645655

print(".... hickory [1] mean, [2]SD, [3] count, [4] SE  ......")

## [1] ".... hickory [1] mean, [2]SD, [3] count, [4] SE  ......"

m.b.car<-subset(mal.b, host == "Hickory")
         mean(m.b.car$Pupa.mass2)

## [1] 131.6221

         sd(m.b.car$Pupa.mass2)

## [1] 19.42891

         length(m.b.car$Pupa.mass2)

## [1] 113

         sd(m.b.car$Pupa.mass2)/sqrt(length(m.b.car$Pupa.mass2))

## [1] 1.827718

print(".... apple [1] mean, [2]SD, [3] count, [4] SE  ......")

## [1] ".... apple [1] mean, [2]SD, [3] count, [4] SE  ......"

m.b.mal<-subset(mal.b, host =="Malus")
      mean(m.b.mal$Pupa.mass2)

## [1] 147.6798

       sd(m.b.mal$Pupa.mass2)

## [1] 20.07439

       length(m.b.mal$Pupa.mass2)

## [1] 94

       sd(m.b.mal$Pupa.mass2)/sqrt(length(m.b.mal$Pupa.mass2))

## [1] 2.070515

print(".... rotating [1] mean, [2]SD, [3] count, [4] SE  ......")

## [1] ".... rotating [1] mean, [2]SD, [3] count, [4] SE  ......"

m.b.mix<-subset(mal.b, host == "Mix")
      mean(m.b.mix$Pupa.mass2)

## [1] 135.7419

       sd(m.b.mix$Pupa.mass2)

## [1] 18.12751

       length(m.b.mix$Pupa.mass2)

## [1] 105

       sd(m.b.mix$Pupa.mass2)/sqrt(length(m.b.mix$Pupa.mass2))

## [1] 1.769064

Linear mixed model

"male black Pupal mass"

## [1] "male black Pupal mass"

mass.male.b<- lmer(Pupa.mass2~host+(1|Mother)+(1|Father),data=mal.b)
summary(mass.male.b)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Pupa.mass2 ~ host + (1 | Mother) + (1 | Father)
##    Data: mal.b
## 
## REML criterion at convergence: 4647
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.5912 -0.5961  0.0267  0.6374  3.5328 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Father   (Intercept) 142.86   11.952  
##  Mother   (Intercept)  61.35    7.833  
##  Residual             291.46   17.072  
## Number of obs: 544, groups:  Father, 9; Mother, 6
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)  143.971      5.413  10.408  26.597 6.57e-11 ***
## hostElm      -13.600      2.264 525.672  -6.007 3.53e-09 ***
## hostHickory  -12.502      2.227 525.155  -5.613 3.23e-08 ***
## hostMalus      3.452      2.352 525.647   1.468    0.143    
## hostMix       -9.841      2.282 525.531  -4.312 1.93e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) hstElm hstHck hstMls
## hostElm     -0.191                     
## hostHickory -0.193  0.471              
## hostMalus   -0.181  0.448  0.451       
## hostMix     -0.193  0.465  0.464  0.447

anova(mass.male.b)

## Type III Analysis of Variance Table with Satterthwaite's method
##      Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## host  24923  6230.8     4 525.43  21.378 2.363e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#random effect (black males)

print("black MALES")

## [1] "black MALES"

print("parental effect")

## [1] "parental effect"

r.squaredGLMM(mass.male.b)

##             R2m       R2c
## [1,] 0.08571004 0.4623855

(0.458-0.082)*100

## [1] 37.6

diet effect size (black males)

F_to_eta2(20.28, 4, 526)

## Eta2 (partial) |       95% CI
## -----------------------------
## 0.13           | [0.09, 1.00]
## 
## - One-sided CIs: upper bound fixed at [1.00].

Tukey’s test pairwise comparison (black males)

"male black pupal mass Tukey"

## [1] "male black pupal mass Tukey"

emmeans(mass.male.b, list(pairwise~host), adjust="tukey")

## $`emmeans of host`
##  host      emmean   SE   df lower.CL upper.CL
##  B. cherry    144 5.42 12.6      132      156
##  Elm          130 5.47 12.9      119      142
##  Hickory      131 5.45 12.8      120      143
##  Malus        147 5.51 13.3      136      159
##  Mix          134 5.46 12.9      122      146
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $`pairwise differences of host`
##  1                   estimate   SE  df t.ratio p.value
##  B. cherry - Elm        13.60 2.26 528   6.005  <.0001
##  B. cherry - Hickory    12.50 2.23 528   5.612  <.0001
##  B. cherry - Malus      -3.45 2.35 528  -1.467  0.5842
##  B. cherry - Mix         9.84 2.28 528   4.311  0.0002
##  Elm - Hickory          -1.10 2.31 528  -0.475  0.9896
##  Elm - Malus           -17.05 2.43 528  -7.029  <.0001
##  Elm - Mix              -3.76 2.35 528  -1.599  0.4988
##  Hickory - Malus       -15.95 2.40 528  -6.644  <.0001
##  Hickory - Mix          -2.66 2.34 528  -1.139  0.7856
##  Malus - Mix            13.29 2.44 528   5.454  <.0001
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 5 estimates

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

Red females

subset red females by host

print("       red female pupal mass  (mg)        ")

## [1] "       red female pupal mass  (mg)        "

print(".... black cherry [1] mean, [2]SD, [3] count, [4] SE   ......")

## [1] ".... black cherry [1] mean, [2]SD, [3] count, [4] SE   ......"

f.r.pru<-subset(Fem.r, host == "B. cherry")
            mean(f.r.pru$Pupa.mass2)

## [1] 118.1889

            sd(f.r.pru$Pupa.mass2)

## [1] 17.47314

            length(f.r.pru$Pupa.mass2)

## [1] 36

            sd(f.r.pru$Pupa.mass2)/sqrt(length(f.r.pru$Pupa.mass2))

## [1] 2.912191

print(".... elm [1] mean, [2]SD, [3] count, [4] SE   ......")

## [1] ".... elm [1] mean, [2]SD, [3] count, [4] SE   ......"

f.r.ulm<-subset(Fem.r, host == "Elm")
           mean(f.r.ulm$Pupa.mass2)

## [1] 104.7667

           sd(f.r.ulm$Pupa.mass2)

## [1] 12.9736

           length(f.r.ulm$Pupa.mass2)

## [1] 36

           sd(f.r.ulm$Pupa.mass2)/sqrt(length(f.r.ulm$Pupa.mass2))

## [1] 2.162267

print(".... hickory [1] mean, [2]SD, [3] count, [4] SE   ......")

## [1] ".... hickory [1] mean, [2]SD, [3] count, [4] SE   ......"

f.r.car<-subset(Fem.r, host == "Hickory")
          mean(f.r.car$Pupa.mass2)

## [1] 108.5289

          sd(f.r.car$Pupa.mass2)

## [1] 11.30156

          length(f.r.car$Pupa.mass2)

## [1] 38

          sd(f.r.car$Pupa.mass2)/sqrt(length(f.r.car$Pupa.mass2))

## [1] 1.833356

print(".... apple [1] mean, [2]SD, [3] count, [4] SE   ......")

## [1] ".... apple [1] mean, [2]SD, [3] count, [4] SE   ......"

f.r.mal<-subset(Fem.r, host =="Malus")
          mean(f.r.mal$Pupa.mass2)

## [1] 99.4

          sd(f.r.mal$Pupa.mass2)

## [1] 18.18592

          length(f.r.mal$Pupa.mass2)

## [1] 29

          sd(f.r.mal$Pupa.mass2)/sqrt(length(f.r.mal$Pupa.mass2))

## [1] 3.377041

print(".... rotating [1] mean, [2]SD, [3] count, [4] SE   ......")

## [1] ".... rotating [1] mean, [2]SD, [3] count, [4] SE   ......"

f.r.mix<-subset(Fem.r, host == "Mix")
             mean(f.r.mix$Pupa.mass2)

## [1] 102.8621

               sd(f.r.mix$Pupa.mass2)

## [1] 14.85014

              length(f.r.mix$Pupa.mass2)

## [1] 29

              sd(f.r.mix$Pupa.mass2)/sqrt(length(f.r.mix$Pupa.mass2))

## [1] 2.757602

Linear mixed model– Pupa mass (red females)

"female red Pupal mass"

## [1] "female red Pupal mass"

mass.female.r<- lmer(Pupa.mass2~host+(1|Mother)+(1|Father),data=Fem.r)

## boundary (singular) fit: see help('isSingular')

summary(mass.female.r)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Pupa.mass2 ~ host + (1 | Mother) + (1 | Father)
##    Data: Fem.r
## 
## REML criterion at convergence: 1343.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3246 -0.5959 -0.0730  0.5448  3.4610 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Father   (Intercept)  40.22    6.342  
##  Mother   (Intercept)   0.00    0.000  
##  Residual             189.33   13.760  
## Number of obs: 168, groups:  Father, 6; Mother, 5
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)  119.607      3.598   9.553  33.240 3.38e-11 ***
## hostElm      -15.384      3.299 161.419  -4.663 6.50e-06 ***
## hostHickory  -10.829      3.226 159.126  -3.357 0.000987 ***
## hostMalus    -19.497      3.440 158.518  -5.668 6.66e-08 ***
## hostMix      -15.501      3.457 159.053  -4.484 1.40e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) hstElm hstHck hstMls
## hostElm     -0.477                     
## hostHickory -0.463  0.507              
## hostMalus   -0.427  0.466  0.478       
## hostMix     -0.419  0.461  0.464  0.437
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

anova(mass.female.r)

## Type III Analysis of Variance Table with Satterthwaite's method
##      Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## host 7585.9  1896.5     4 160.07  10.017 2.927e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

random effect (red females)

print("red FEMALES")

## [1] "red FEMALES"

print("parental effect")

## [1] "parental effect"

r.squaredGLMM(mass.female.r)

##            R2m       R2c
## [1,] 0.1659247 0.3120596

(0.312-0.166)*100

## [1] 14.6

diets effect size (red females)

F_to_eta2(10.017, 4, 160)

## Eta2 (partial) |       95% CI
## -----------------------------
## 0.20           | [0.10, 1.00]
## 
## - One-sided CIs: upper bound fixed at [1.00].

Tukey’s test pairwise comparison of diets (red female)

"female red pupal mass Tukey"

## [1] "female red pupal mass Tukey"

emmeans(mass.female.r, list(pairwise~host), adjust="tukey")

## $`emmeans of host`
##  host      emmean   SE   df lower.CL upper.CL
##  B. cherry    120 4.44 5.41    108.4      131
##  Elm          104 4.32 5.65     93.5      115
##  Hickory      109 4.33 5.46     97.9      120
##  Malus        100 4.57 6.56     89.2      111
##  Mix          104 4.66 6.43     92.9      115
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $`pairwise differences of host`
##  1                   estimate   SE  df t.ratio p.value
##  B. cherry - Elm       15.384 3.31 161   4.642  0.0001
##  B. cherry - Hickory   10.829 3.23 159   3.348  0.0089
##  B. cherry - Malus     19.497 3.44 159   5.665  <.0001
##  B. cherry - Mix       15.501 3.46 159   4.476  0.0001
##  Elm - Hickory         -4.555 3.25 161  -1.401  0.6277
##  Elm - Malus            4.113 3.50 161   1.177  0.7648
##  Elm - Mix              0.117 3.53 162   0.033  1.0000
##  Hickory - Malus        8.668 3.41 159   2.539  0.0872
##  Hickory - Mix          4.672 3.48 160   1.341  0.6661
##  Malus - Mix           -3.996 3.67 160  -1.089  0.8119
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 5 estimates

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

Red males

subset red males by host

print("       red male pupal mass  (mg)        ")

## [1] "       red male pupal mass  (mg)        "

print(".... black cherry [1] mean, [2]SD, [3] count, [4] SE   ......")

## [1] ".... black cherry [1] mean, [2]SD, [3] count, [4] SE   ......"

m.r.pru<-subset(mal.r, host == "B. cherry")
         mean(m.r.pru$Pupa.mass2)

## [1] 91.70385

        sd(m.r.pru$Pupa.mass2)

## [1] 12.11873

        length(m.r.pru$Pupa.mass2)

## [1] 26

        sd(m.r.pru$Pupa.mass2)/sqrt(length(m.r.pru$Pupa.mass2))

## [1] 2.376678

print(".... elm [1] mean, [2]SD, [3] count, [4] SE   ......")

## [1] ".... elm [1] mean, [2]SD, [3] count, [4] SE   ......"

m.r.ulm<-subset(mal.r, host == "Elm")
         mean(m.r.ulm$Pupa.mass2)

## [1] 83.93333

         sd(m.r.ulm$Pupa.mass2)

## [1] 13.23118

         length(m.r.ulm$Pupa.mass2)

## [1] 24

         sd(m.r.ulm$Pupa.mass2)/sqrt(length(m.r.ulm$Pupa.mass2))

## [1] 2.700803

print(".... hickory [1] mean, [2]SD, [3] count, [4] SE   ......")

## [1] ".... hickory [1] mean, [2]SD, [3] count, [4] SE   ......"

m.r.car<-subset(mal.r, host == "Hickory")
          mean(m.r.car$Pupa.mass2)

## [1] 91.61935

          sd(m.r.car$Pupa.mass2)

## [1] 20.41982

          length(m.r.car$Pupa.mass2)

## [1] 31

          sd(m.r.car$Pupa.mass2)/sqrt(length(m.r.car$Pupa.mass2))

## [1] 3.667507

print(".... apple [1] mean, [2]SD, [3] count, [4] SE   ......")

## [1] ".... apple [1] mean, [2]SD, [3] count, [4] SE   ......"

m.r.mal<-subset(mal.r, host =="Malus")
          mean(m.r.mal$Pupa.mass2)

## [1] 79.6697

          sd(m.r.mal$Pupa.mass2)

## [1] 10.07567

          length(m.r.mal$Pupa.mass2)

## [1] 33

          sd(m.r.mal$Pupa.mass2)/sqrt(length(m.r.mal$Pupa.mass2))

## [1] 1.753948

print(".... rotating [1] mean, [2]SD, [3] count, [4] SE   ......")

## [1] ".... rotating [1] mean, [2]SD, [3] count, [4] SE   ......"

m.r.mix<-subset(mal.r, host == "Mix")
          mean(m.r.mix$Pupa.mass2)

## [1] 84.73636

          sd(m.r.mix$Pupa.mass2)

## [1] 13.88413

          length(m.r.mix$Pupa.mass2)

## [1] 22

          sd(m.r.mix$Pupa.mass2)/sqrt(length(m.r.mix$Pupa.mass2))

## [1] 2.960107

Linear mixed model – pupa mass (red male)

"male red Pupal mass"

## [1] "male red Pupal mass"

mass.male.r<- lmer(Pupa.mass2~host+(1|Mother)+(1|Father),data=mal.r)

## boundary (singular) fit: see help('isSingular')

summary(mass.male.r)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Pupa.mass2 ~ host + (1 | Mother) + (1 | Father)
##    Data: mal.r
## 
## REML criterion at convergence: 1074.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1502 -0.6769 -0.1127  0.5850  5.4289 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Father   (Intercept)  36.5     6.042  
##  Mother   (Intercept)   0.0     0.000  
##  Residual             177.9    13.339  
## Number of obs: 136, groups:  Father, 5; Mother, 4
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)   90.937      3.768  10.732  24.134 1.07e-10 ***
## hostElm       -7.346      3.838 128.823  -1.914 0.057831 .  
## hostHickory    1.618      3.568 127.598   0.453 0.651065    
## hostMalus    -12.881      3.528 128.114  -3.651 0.000379 ***
## hostMix       -7.642      3.914 128.378  -1.952 0.053061 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) hstElm hstHck hstMls
## hostElm     -0.480                     
## hostHickory -0.514  0.510              
## hostMalus   -0.513  0.520  0.543       
## hostMix     -0.458  0.454  0.487  0.510
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

anova(mass.male.r)

## Type III Analysis of Variance Table with Satterthwaite's method
##      Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
## host 4167.1  1041.8     4 128.4  5.8548 0.0002324 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

random effect (red male)

print("red MALES")

## [1] "red MALES"

print("parental effect")

## [1] "parental effect"

r.squaredGLMM(mass.male.r)

##            R2m       R2c
## [1,] 0.1301509 0.2782195

(0.278-0.13)*100

## [1] 14.8

diet effect size (red male)

F_to_eta2(5.8548, 4, 128.4)

## Eta2 (partial) |       95% CI
## -----------------------------
## 0.15           | [0.05, 1.00]
## 
## - One-sided CIs: upper bound fixed at [1.00].

Tukey’s test pairwise diet (red males)

"male red pupal mass Tukey"

## [1] "male red pupal mass Tukey"

emmeans(mass.male.r, list(pairwise~host), adjust="tukey")

## $`emmeans of host`
##  host      emmean   SE   df lower.CL upper.CL
##  B. cherry   90.9 4.51 5.92     79.9    102.0
##  Elm         83.6 4.66 6.31     72.3     94.9
##  Hickory     92.6 4.44 4.86     81.0    104.1
##  Malus       78.1 4.43 4.74     66.5     89.6
##  Mix         83.3 4.82 6.77     71.8     94.8
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $`pairwise differences of host`
##  1                   estimate   SE  df t.ratio p.value
##  B. cherry - Elm        7.346 3.86 129   1.904  0.3207
##  B. cherry - Hickory   -1.618 3.58 127  -0.452  0.9913
##  B. cherry - Malus     12.881 3.54 128   3.641  0.0035
##  B. cherry - Mix        7.642 3.93 128   1.945  0.2993
##  Elm - Hickory         -8.964 3.69 129  -2.432  0.1134
##  Elm - Malus            5.534 3.63 128   1.524  0.5487
##  Elm - Mix              0.296 4.08 129   0.072  1.0000
##  Hickory - Malus       14.499 3.40 129   4.259  0.0004
##  Hickory - Mix          9.260 3.82 129   2.423  0.1157
##  Malus - Mix           -5.239 3.71 128  -1.413  0.6205
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 5 estimates

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

x Figures x

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

# using package ggstatsplot 

bf<-ggstatsplot::ggbetweenstats(
  data = Fem.b,
    x = host,
  y = Pupa.mass2,
  color=host,
   messages=FALSE,
  centrality.label.args = list(size  = 2),
 mean.plotting = FALSE,
  mean.ci=FALSE,
  mean.label.size=1,
 bf.message = FALSE, 
  pairwise.comparisons = FALSE,
    results.subtitle=FALSE,)+
  ylim(50, 300)+
 
  
  labs(
    x = "Caterpillar diet",
    y = "Pupal mass (mg)",
    #title = "                    Distribution of black-headed males by diet"
  )+
  scale_color_manual(values= c("pink", "green", "lightblue", "yellow", "orange"))+
  theme(text = element_text( size = 1, color = "black"),
    plot.title = element_text(
      size = 20,
      face = "bold",
      color = "#2a475e",
     
    ),
    plot.subtitle = element_text(
       
      size = 20, 
      face = "bold",
      color="red"
    ),
    plot.title.position = "plot", # slightly different from default
    axis.text = element_text(size = 10, color = "black"),
    axis.title = element_text(size = 15)
  )+
  theme(
    axis.ticks = element_blank(),
    axis.line = element_line(colour = "grey50"),
    panel.grid = element_line(color = "#b4aea9"),
    panel.grid.minor = element_blank(),
    panel.grid.major.x = element_blank(),
    panel.grid.major.y = element_line(linetype = "dashed"),
    panel.background = element_rect(fill = "white", color = "white"),
    plot.background = element_rect(fill = "white", color = "white"),
      # panel.background = element_rect(fill = "#fbf9f4", color = "#fbf9f4"),
    #plot.background = element_rect(fill = "#fbf9f4", color = "#fbf9f4"),
     )

## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

# using package ggstatsplot 

bm<-ggstatsplot::ggbetweenstats(
  data = mal.b,
    x = host,
  y = Pupa.mass2,
  color=host,
   messages=FALSE,
  centrality.label.args = list(size  = 2),
 mean.plotting = FALSE,
  mean.ci=FALSE,
  mean.label.size=1,
 bf.message = FALSE, 
  pairwise.comparisons = FALSE,
    results.subtitle=FALSE,)+
  ylim(50, 300)+
 
  
  labs(
    x = "Caterpillar diet",
    y = "Pupal mass (mg)",
    #title = "                    Distribution of black-headed males by diet"
  )+
  scale_color_manual(values= c("pink", "green", "lightblue", "yellow", "orange"))+
  theme(text = element_text( size = 1, color = "black"),
    plot.title = element_text(
      size = 20,
      face = "bold",
      color = "#2a475e",
     
    ),
    plot.subtitle = element_text(
       
      size = 20, 
      face = "bold",
      color="red"
    ),
    plot.title.position = "plot", # slightly different from default
    axis.text = element_text(size = 10, color = "black"),
    axis.title = element_text(size = 15)
  )+
  theme(
    axis.ticks = element_blank(),
    axis.line = element_line(colour = "grey50"),
    panel.grid = element_line(color = "#b4aea9"),
    panel.grid.minor = element_blank(),
    panel.grid.major.x = element_blank(),
    panel.grid.major.y = element_line(linetype = "dashed"),
    panel.background = element_rect(fill = "white", color = "white"),
    plot.background = element_rect(fill = "white", color = "white"),
      # panel.background = element_rect(fill = "#fbf9f4", color = "#fbf9f4"),
    #plot.background = element_rect(fill = "#fbf9f4", color = "#fbf9f4"),
     )

## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

# using package ggstatsplot 

rf<-ggstatsplot::ggbetweenstats(
  data = Fem.r,
    x = host,
  y = Pupa.mass2,
  color=host,
   messages=FALSE,
  centrality.label.args = list(size  = 2),
 mean.plotting = FALSE,
  mean.ci=FALSE,
  mean.label.size=1,
 bf.message = FALSE, 
  pairwise.comparisons = FALSE,
    results.subtitle=FALSE,)+
  ylim(50, 300)+
 
  
  labs(
    x = "Caterpillar diet",
    y = "Pupal mass (mg)",
    #title = "                    Distribution of black-headed males by diet"
  )+
  scale_color_manual(values= c("pink", "green", "lightblue", "yellow", "orange"))+
  theme(text = element_text( size = 1, color = "black"),
    plot.title = element_text(
      size = 20,
      face = "bold",
      color = "#2a475e",
     
    ),
    plot.subtitle = element_text(
       
      size = 20, 
      face = "bold",
      color="red"
    ),
    plot.title.position = "plot", # slightly different from default
    axis.text = element_text(size = 10, color = "black"),
    axis.title = element_text(size = 15)
  )+
  theme(
    axis.ticks = element_blank(),
    axis.line = element_line(colour = "grey50"),
    panel.grid = element_line(color = "#b4aea9"),
    panel.grid.minor = element_blank(),
    panel.grid.major.x = element_blank(),
    panel.grid.major.y = element_line(linetype = "dashed"),
    panel.background = element_rect(fill = "white", color = "white"),
    plot.background = element_rect(fill = "white", color = "white"),
      # panel.background = element_rect(fill = "#fbf9f4", color = "#fbf9f4"),
    #plot.background = element_rect(fill = "#fbf9f4", color = "#fbf9f4"),
     )

## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

# using package ggstatsplot 

rm<-ggstatsplot::ggbetweenstats(
  data = mal.r,
    x = host,
  y = Pupa.mass2,
  color=host,
   messages=FALSE,
  centrality.label.args = list(size  = 2),
 mean.plotting = FALSE,
  mean.ci=FALSE,
  mean.label.size=1,
 bf.message = FALSE, 
  pairwise.comparisons = FALSE,
    results.subtitle=FALSE,)+
  ylim(50, 300)+
 
  
  labs(
    x = "Caterpillar diet",
    y = "Pupal mass (mg)",
    #title = "                    Distribution of black-headed males by diet"
  )+
  scale_color_manual(values= c("pink", "green", "lightblue", "yellow", "orange"))+
  theme(text = element_text( size = 1, color = "black"),
    plot.title = element_text(
      size = 20,
      face = "bold",
      color = "#2a475e",
     
    ),
    plot.subtitle = element_text(
       
      size = 20, 
      face = "bold",
      color="red"
    ),
    plot.title.position = "plot", # slightly different from default
    axis.text = element_text(size = 10, color = "black"),
    axis.title = element_text(size = 15)
  )+
  theme(
    axis.ticks = element_blank(),
    axis.line = element_line(colour = "grey50"),
    panel.grid = element_line(color = "#b4aea9"),
    panel.grid.minor = element_blank(),
    panel.grid.major.x = element_blank(),
    panel.grid.major.y = element_line(linetype = "dashed"),
    panel.background = element_rect(fill = "white", color = "white"),
    plot.background = element_rect(fill = "white", color = "white"),
      # panel.background = element_rect(fill = "#fbf9f4", color = "#fbf9f4"),
    #plot.background = element_rect(fill = "#fbf9f4", color = "#fbf9f4"),
     )

## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

black female mass

bf

black male mass

bm

red female mass

rf

red male mass

rm

mass knit

Mattheau Comerford

2023-10-13

Load library

Set seed to 1234

clean data by removing empty or missing data points, also removing none complete lineages

development time analysis

data.pupa comes from data tab

Subset females by color

subset males by color

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

Black females

subset black females by host

Linear mixed model– pupa mass (black females)

random effect (black females)

host effect size (black females)

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

Black males

subset black males by host

Linear mixed model

Tukey’s test pairwise comparison (black males)

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

Red females

subset red females by host

Linear mixed model– Pupa mass (red females)

random effect (red females)

diets effect size (red females)

Tukey’s test pairwise comparison of diets (red female)

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

Red males

subset red males by host

Linear mixed model – pupa mass (red male)

random effect (red male)

diet effect size (red male)

Tukey’s test pairwise diet (red males)

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

x Figures x

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

black female mass

black male mass

red female mass

red male mass