R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

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)

Import data

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

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

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.mass !="missing",]

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

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

development time analysis

data.pupa comes from data tab

dev<- data

table by sex and color

table(dev$Sex, dev$Color)

##    
##       B   R
##   f 484 169
##   m 549 140

Table 1: test of Biotypes and sex’s effect on larval development

"Test if sex and color interaction effect development time"

## [1] "Test if sex and color interaction effect development time"

mod.1<- lmer(Dev_time~Sex*Color+(1|Mat..line),data=dev)
summary(mod.1)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Dev_time ~ Sex * Color + (1 | Mat..line)
##    Data: dev
## 
## REML criterion at convergence: 5962.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8504 -0.5882 -0.0666  0.4567  6.3476 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Mat..line (Intercept) 10.076   3.174   
##  Residual               6.467   2.543   
## Number of obs: 1250, groups:  Mat..line, 19
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   37.3512     0.8281   17.3498  45.104  < 2e-16 ***
## Sexm          -1.2942     0.1628 1229.2348  -7.950 4.19e-15 ***
## ColorR        12.8859     1.8031   17.2837   7.147 1.49e-06 ***
## Sexm:ColorR   -0.1266     0.3592 1229.4014  -0.352    0.725    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Sexm   ColorR
## Sexm        -0.105              
## ColorR      -0.459  0.048       
## Sexm:ColorR  0.047 -0.453 -0.093

anova(mod.1)

## Type III Analysis of Variance Table with Satterthwaite's method
##           Sum Sq Mean Sq NumDF   DenDF F value    Pr(>F)    
## Sex       369.43  369.43     1 1229.40 57.1281 7.950e-14 ***
## Color     329.88  329.88     1   16.99 51.0126 1.656e-06 ***
## Sex:Color   0.80    0.80     1 1229.40  0.1241    0.7247    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r.squaredGLMM(mod.1)

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

##            R2m       R2c
## [1,] 0.6279265 0.8545507

random effect contribution

print("(Random effect size = conditional R squared) - (marginal R squared)" )

## [1] "(Random effect size = conditional R squared) - (marginal R squared)"

(0.853-0.626)*100

## [1] 22.7

print(" eta-squared for sex")

## [1] " eta-squared for sex"

F_to_eta2(58.428, 1, 1246)

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

print(" eta-squared for biotype")

## [1] " eta-squared for biotype"

F_to_eta2(50.664, 1, 17)

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

print(" eta-squared for interaction")

## [1] " eta-squared for interaction"

F_to_eta2(0.2306, 1, 1246)

## Eta2 (partial) |       95% CI
## -----------------------------
## 1.85e-04       | [0.00, 1.00]
## 
## - One-sided CIs: upper bound fixed at [1.00].

Tukeys test

emmeans(mod.1, list(pairwise~Sex*Color), adjust="tukey") # sex and color must be run as independent separate models

## $`emmeans of Sex, Color`
##  Sex Color emmean    SE   df lower.CL upper.CL
##  f   B       37.4 0.828 17.4     35.6     39.1
##  m   B       36.1 0.827 17.3     34.3     37.8
##  f   R       50.2 1.602 17.3     46.9     53.6
##  m   R       48.8 1.605 17.4     45.4     52.2
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $`pairwise differences of Sex, Color`
##  1         estimate    SE     df t.ratio p.value
##  f B - m B     1.29 0.163 1229.3   7.950  <.0001
##  f B - f R   -12.89 1.803   17.3  -7.147  <.0001
##  f B - m R   -11.47 1.806   17.4  -6.349  <.0001
##  m B - f R   -14.18 1.803   17.3  -7.866  <.0001
##  m B - m R   -12.76 1.805   17.4  -7.067  <.0001
##  f R - m R     1.42 0.320 1229.5   4.437  0.0001
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 4 estimates

Table 2, SD for development time by sex and color

dev.black<-subset(dev, Color == "B")
dev.red<-subset(dev, Color=="R")

print(".........Black females ........")

## [1] ".........Black females ........"

dev.fb<-subset(dev.black, Sex== "f")
print("balck female [1] mean, [2] SD, [3] count, [4] SE")

## [1] "balck female [1] mean, [2] SD, [3] count, [4] SE"

mean(na.omit(dev.fb$Dev_time))

## [1] 37.35498

a<-sd(na.omit(dev.fb$Dev_time))
a

## [1] 2.702848

b<-length(na.omit(dev.fb$Dev_time))
b

## [1] 462

a/sqrt(b)

## [1] 0.1257479

print(".........Black males ........")

## [1] ".........Black males ........"

dev.mb<-subset(dev.black, Sex== "m")
print("black male [1] mean, [2] SD, [3] count")

## [1] "black male [1] mean, [2] SD, [3] count"

mean(na.omit(dev.mb$Dev_time))

## [1] 35.99811

sd(na.omit(dev.mb$Dev_time))

## [1] 2.39831

length(na.omit(dev.mb$Dev_time))

## [1] 529

(sd(na.omit(dev.mb$Dev_time)))/sqrt(length(na.omit(dev.mb$Dev_time)))

## [1] 0.1042743

print(".........Red females ........")

## [1] ".........Red females ........"

dev.fr<-subset(dev.red, Sex== "f")
print("red female [1] mean, [2] SD, [3] count")

## [1] "red female [1] mean, [2] SD, [3] count"

mean(na.omit(dev.fr$Dev_time))

## [1] 50.29787

sd(na.omit(dev.fr$Dev_time))

## [1] 7.109093

length(na.omit(dev.fr$Dev_time))

## [1] 141

sd(na.omit(dev.fr$Dev_time))/sqrt(length(na.omit(dev.fr$Dev_time)))

## [1] 0.5986937

print(".........Red males ........")

## [1] ".........Red males ........"

dev.mr<-subset(dev.red, Sex== "m")
print("red male [1] mean, [2] SD, [3] count")

## [1] "red male [1] mean, [2] SD, [3] count"

mean(na.omit(dev.mr$Dev_time))

## [1] 47.64407

sd(na.omit(dev.mr$Dev_time))

## [1] 6.691652

length(na.omit(dev.mr$Dev_time))

## [1] 118

sd(na.omit(dev.mr$Dev_time))/sqrt(length(na.omit(dev.mr$Dev_time)))

## [1] 0.6160165

Subsetting by demographic

subset by sex

dev.F<- subset(dev, Sex=="f") # females only


dev.M<- subset(dev, Sex =="m") # males only

subset females by color

dev.F.R<- subset(dev.F, Color=="R") # female Reds only

dev.F.B<- subset(dev.F, Color=="B") # female blacks only

subset males by color

dev.M.R<- subset(dev.M, Color=="R") # male Reds only

dev.M.B<- subset(dev.M, Color=="B") # male blacks only

black female

subset black females by host

Means, SD, Count, and SE by host

print("       Black female develoment time  (days)        ")

## [1] "       Black female develoment time  (days)        "

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

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

dev.fb.pru<-subset(dev.fb, host == "B. cherry")
       mean(dev.fb.pru$Dev_time, na.rm = TRUE)

## [1] 36.04819

       sd(dev.fb.pru$Dev_time, na.rm = TRUE)

## [1] 2.112123

       length(dev.fb.pru$Dev_time)

## [1] 94

       sd(dev.fb.pru$Dev_time, na.rm = TRUE)/sqrt(length(dev.fb.pru$Dev_time))

## [1] 0.2178489

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

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

dev.fb.ulm<-subset(dev.fb, host == "Elm")
      mean(dev.fb.ulm$Dev_time, na.rm = TRUE)

## [1] 37.34021

      sd(dev.fb.ulm$Dev_time, na.rm = TRUE)

## [1] 2.014978

      length(dev.fb.ulm$Dev_time)

## [1] 101

      sd(dev.fb.ulm$Dev_time, na.rm = TRUE)/sqrt(length(dev.fb.ulm$Dev_time))

## [1] 0.2004978

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

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

dev.fb.car<-subset(dev.fb, host == "Hickory")
      mean(dev.fb.car$Dev_time, na.rm = TRUE)

## [1] 38.88421

      sd(dev.fb.car$Dev_time, na.rm = TRUE)

## [1] 2.76322

      length(dev.fb.car$Dev_time)

## [1] 99

      sd(dev.fb.car$Dev_time, na.rm = TRUE)/sqrt(length(dev.fb.car$Dev_time))

## [1] 0.2777141

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

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

dev.fb.mal<-subset(dev.fb, host =="Malus")
      mean(dev.fb.mal$Dev_time, na.rm = TRUE)

## [1] 36.64151

      sd(dev.fb.mal$Dev_time, na.rm = TRUE)

## [1] 3.117351

      length(dev.fb.mal$Dev_time)

## [1] 106

      sd(dev.fb.mal$Dev_time, na.rm = TRUE)/sqrt(length(dev.fb.mal$Dev_time))

## [1] 0.3027839

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

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

dev.fb.mix<-subset(dev.fb, host == "Mix")
       mean(dev.fb.mix$Dev_time, na.rm = TRUE )

## [1] 37.85185

       sd(dev.fb.mix$Dev_time, na.rm = TRUE )

## [1] 2.335118

       length(dev.fb.mix$Dev_time)

## [1] 84

       sd(dev.fb.mix$Dev_time, na.rm = TRUE)/sqrt(length(dev.fb.mix$Dev_time))

## [1] 0.2547823

Female black model – development time

"female black development time"

## [1] "female black development time"

mod.female.b<- lmer(Dev_time~host+(1|Mother)+(1|Father), data=dev.F.B)
summary(mod.female.b)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Dev_time ~ host + (1 | Mother) + (1 | Father)
##    Data: dev.F.B
## 
## REML criterion at convergence: 2128.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4970 -0.5668 -0.0442  0.5359  6.8319 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Father   (Intercept) 0.2177   0.4666  
##  Mother   (Intercept) 0.8457   0.9196  
##  Residual             5.6560   2.3782  
## Number of obs: 462, groups:  Father, 9; Mother, 6
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  36.0975     0.4985   9.3492  72.406 3.57e-14 ***
## hostElm       1.0718     0.3595 451.4804   2.982  0.00302 ** 
## hostHickory   2.8055     0.3581 446.7111   7.834 3.47e-14 ***
## hostMalus     0.4493     0.3523 449.4292   1.275  0.20290    
## hostMix       1.5854     0.3763 452.2771   4.213 3.04e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) hstElm hstHck hstMls
## hostElm     -0.391                     
## hostHickory -0.393  0.535              
## hostMalus   -0.406  0.548  0.549       
## hostMix     -0.372  0.523  0.511  0.526

anova(mod.female.b)

## Type III Analysis of Variance Table with Satterthwaite's method
##      Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## host  436.9  109.22     4 450.18  19.312 1.123e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

variation contributed by lineage to black female development time

r.squaredGLMM(mod.female.b)

##            R2m       R2c
## [1,] 0.1237883 0.2624592

0.265-0.121

## [1] 0.144

host effect size on black female development time

F_to_eta2(19.202, 4, 454)

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

tukey’s test of pairwise diet on female development time

"female black development time Tukey"

## [1] "female black development time Tukey"

emmeans(mod.female.b, list(pairwise~host), adjust="tukey")

## $`emmeans of host`
##  host      emmean    SE    df lower.CL upper.CL
##  B. cherry   36.1 0.504 10.47     35.0     37.2
##  Elm         37.2 0.492  9.66     36.1     38.3
##  Hickory     38.9 0.491  9.57     37.8     40.0
##  Malus       36.5 0.485  9.08     35.5     37.6
##  Mix         37.7 0.506 10.71     36.6     38.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       -1.072 0.360 452  -2.973  0.0257
##  B. cherry - Hickory   -2.806 0.358 447  -7.831  <.0001
##  B. cherry - Malus     -0.449 0.353 450  -1.273  0.7080
##  B. cherry - Mix       -1.585 0.378 453  -4.199  0.0003
##  Elm - Hickory         -1.734 0.347 451  -4.997  <.0001
##  Elm - Malus            0.622 0.339 452   1.835  0.3546
##  Elm - Mix             -0.514 0.360 450  -1.427  0.6108
##  Hickory - Malus        2.356 0.338 447   6.978  <.0001
##  Hickory - Mix          1.220 0.364 452   3.348  0.0078
##  Malus - Mix           -1.136 0.356 452  -3.189  0.0132
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 5 estimates

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

black Male

subset black males by host

Means, SD, Count, and SE by host

print("       Black male develoment time  (days)        ")

## [1] "       Black male develoment time  (days)        "

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

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

dev.mb.pru<-subset(dev.mb, host == "B. cherry")
       mean(dev.mb.pru$Dev_time, na.rm = TRUE)

## [1] 35.57895

       sd(dev.mb.pru$Dev_time, na.rm = TRUE)

## [1] 2.272945

       length(dev.mb.pru$Dev_time)

## [1] 125

       sd(dev.mb.pru$Dev_time, na.rm = TRUE)/sqrt(length(dev.mb.pru$Dev_time))

## [1] 0.2032984

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

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

dev.mb.ulm<-subset(dev.mb, host == "Elm")
      mean(dev.mb.ulm$Dev_time, na.rm = TRUE)

## [1] 35.4

      sd(dev.mb.ulm$Dev_time, na.rm = TRUE)

## [1] 2.293804

      length(dev.mb.ulm$Dev_time)

## [1] 109

      sd(dev.mb.ulm$Dev_time, na.rm = TRUE)/sqrt(length(dev.mb.ulm$Dev_time))

## [1] 0.2197066

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

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

dev.mb.car<-subset(dev.mb, host == "Hickory")
      mean(dev.mb.car$Dev_time, na.rm = TRUE)

## [1] 37.68142

      sd(dev.mb.car$Dev_time, na.rm = TRUE)

## [1] 2.56092

      length(dev.mb.car$Dev_time)

## [1] 115

      sd(dev.mb.car$Dev_time, na.rm = TRUE)/sqrt(length(dev.mb.car$Dev_time))

## [1] 0.238807

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

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

dev.mb.mal<-subset(dev.mb, host =="Malus")
      mean(dev.mb.mal$Dev_time, na.rm = TRUE)

## [1] 34.78261

      sd(dev.mb.mal$Dev_time, na.rm = TRUE)

## [1] 1.862446

      length(dev.mb.mal$Dev_time)

## [1] 95

      sd(dev.mb.mal$Dev_time, na.rm = TRUE)/sqrt(length(dev.mb.mal$Dev_time))

## [1] 0.1910829

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

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

dev.mb.mix<-subset(dev.mb, host == "Mix")
       mean(dev.mb.mix$Dev_time, na.rm = TRUE )

## [1] 36.30476

       sd(dev.mb.mix$Dev_time, na.rm = TRUE )

## [1] 1.787369

       length(dev.mb.mix$Dev_time)

## [1] 105

       sd(dev.mb.mix$Dev_time, na.rm = TRUE)/sqrt(length(dev.mb.mix$Dev_time))

## [1] 0.1744294

Male black model – development time

"male black development time"

## [1] "male black development time"

mod.male.b<- lmer(Dev_time~host+(1|Mother)+(1|Father), data=dev.M.B)
summary(mod.male.b)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Dev_time ~ host + (1 | Mother) + (1 | Father)
##    Data: dev.M.B
## 
## REML criterion at convergence: 2289.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4179 -0.6970 -0.1036  0.5410  5.3543 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Father   (Intercept) 0.5944   0.7710  
##  Mother   (Intercept) 0.1080   0.3286  
##  Residual             4.2580   2.0635  
## Number of obs: 529, groups:  Father, 9; Mother, 6
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  35.47958    0.35923  15.16874  98.765  < 2e-16 ***
## hostElm      -0.03843    0.28090 515.69460  -0.137  0.89122    
## hostHickory   2.03802    0.27454 514.20809   7.423 4.77e-13 ***
## hostMalus    -0.67396    0.29045 514.98945  -2.320  0.02071 *  
## hostMix       0.88074    0.28086 515.14275   3.136  0.00181 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) hstElm hstHck hstMls
## hostElm     -0.373                     
## hostHickory -0.378  0.487              
## hostMalus   -0.354  0.460  0.468       
## hostMix     -0.377  0.481  0.485  0.463

anova(mod.male.b)

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

variation contributed by lineage to black male development time

r.squaredGLMM(mod.male.b)

##            R2m       R2c
## [1,] 0.1516366 0.2717649

0.274-0.158

## [1] 0.116

host effect size on black male development time

F_to_eta2(29.007, 4, 525)

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

tukey’s test of pairwise diet on female development time

"male black development time Tukey"

## [1] "male black development time Tukey"

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

## $`emmeans of host`
##  host      emmean    SE   df lower.CL upper.CL
##  B. cherry   35.5 0.364 15.4     34.7     36.3
##  Elm         35.4 0.371 16.0     34.7     36.2
##  Hickory     37.5 0.366 15.5     36.7     38.3
##  Malus       34.8 0.381 17.6     34.0     35.6
##  Mix         36.4 0.369 16.0     35.6     37.1
## 
## 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       0.0384 0.281 516   0.137  0.9999
##  B. cherry - Hickory  -2.0380 0.275 514  -7.418  <.0001
##  B. cherry - Malus     0.6740 0.291 515   2.318  0.1409
##  B. cherry - Mix      -0.8807 0.281 515  -3.133  0.0156
##  Elm - Hickory        -2.0764 0.281 515  -7.377  <.0001
##  Elm - Malus           0.6355 0.297 516   2.137  0.2062
##  Elm - Mix            -0.9192 0.286 515  -3.211  0.0122
##  Hickory - Malus       2.7120 0.292 516   9.278  <.0001
##  Hickory - Mix         1.1573 0.282 515   4.104  0.0005
##  Malus - Mix          -1.5547 0.297 515  -5.242  <.0001
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 5 estimates

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

Red Females

subset red females by host

Means, SD, Count, and SE by host

print("       Red female develoment time  (days)        ")

## [1] "       Red female develoment time  (days)        "

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

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

dev.fr.pru<-subset(dev.fr, host == "B. cherry")
       mean(dev.fr.pru$Dev_time, na.rm = TRUE)

## [1] 51.16129

       sd(dev.fr.pru$Dev_time, na.rm = TRUE)

## [1] 7.95444

       length(dev.fr.pru$Dev_time)

## [1] 36

       sd(dev.fr.pru$Dev_time, na.rm = TRUE)/sqrt(length(dev.fr.pru$Dev_time))

## [1] 1.32574

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

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

dev.fr.ulm<-subset(dev.fr, host == "Elm")
      mean(dev.fr.ulm$Dev_time, na.rm = TRUE)

## [1] 49.03571

      sd(dev.fr.ulm$Dev_time, na.rm = TRUE)

## [1] 5.620611

      length(dev.fr.ulm$Dev_time)

## [1] 36

      sd(dev.fr.ulm$Dev_time, na.rm = TRUE)/sqrt(length(dev.fr.ulm$Dev_time))

## [1] 0.9367685

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

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

dev.fr.car<-subset(dev.fr, host == "Hickory")
      mean(dev.fr.car$Dev_time, na.rm = TRUE)

## [1] 48.53333

      sd(dev.fr.car$Dev_time, na.rm = TRUE)

## [1] 6.463336

      length(dev.fr.car$Dev_time)

## [1] 39

      sd(dev.fr.car$Dev_time, na.rm = TRUE)/sqrt(length(dev.fr.car$Dev_time))

## [1] 1.034962

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

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

dev.fr.mal<-subset(dev.fr, host =="Malus")
      mean(dev.fr.mal$Dev_time, na.rm = TRUE)

## [1] 52

      sd(dev.fr.mal$Dev_time, na.rm = TRUE)

## [1] 8.151966

      length(dev.fr.mal$Dev_time)

## [1] 29

      sd(dev.fr.mal$Dev_time, na.rm = TRUE)/sqrt(length(dev.fr.mal$Dev_time))

## [1] 1.513782

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

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

dev.fr.mix<-subset(dev.fr, host == "Mix")
       mean(dev.fr.mix$Dev_time, na.rm = TRUE )

## [1] 51.06897

       sd(dev.fr.mix$Dev_time, na.rm = TRUE )

## [1] 7.085856

       length(dev.fr.mix$Dev_time)

## [1] 29

       sd(dev.fr.mix$Dev_time, na.rm = TRUE)/sqrt(length(dev.fr.mix$Dev_time))

## [1] 1.31581

Female red model – development time

"female red development time"

## [1] "female red development time"

mod.female.r<- lmer(Dev_time~host+(1|Mother)+(1|Father), data=dev.F.R)
summary(mod.female.r)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Dev_time ~ host + (1 | Mother) + (1 | Father)
##    Data: dev.F.R
## 
## REML criterion at convergence: 725.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8817 -0.6696 -0.0809  0.6123  3.9384 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. 
##  Father   (Intercept) 5.207e+01 7.2159214
##  Mother   (Intercept) 2.719e-07 0.0005215
##  Residual             9.541e+00 3.0888177
## Number of obs: 141, groups:  Father, 4; Mother, 3
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  51.1813     3.6509   3.1120  14.019 0.000647 ***
## hostElm      -1.5958     0.8190 133.0528  -1.949 0.053454 .  
## hostHickory  -2.7668     0.7976 133.0260  -3.469 0.000705 ***
## hostMalus     0.2737     0.8506 133.0025   0.322 0.748133    
## hostMix      -0.5742     0.7988 133.0035  -0.719 0.473500    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) hstElm hstHck hstMls
## hostElm     -0.105                     
## hostHickory -0.108  0.486              
## hostMalus   -0.099  0.442  0.453       
## hostMix     -0.106  0.477  0.487  0.454

anova(mod.female.r)

## Type III Analysis of Variance Table with Satterthwaite's method
##      Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)   
## host 171.68  42.921     4 133.03  4.4987 0.001931 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

variation contributed by lineage to red female development time

r.squaredGLMM(mod.female.r)

##            R2m       R2c
## [1,] 0.0201185 0.8482584

0.85-0.02

## [1] 0.83

host effect size on red female development time

r.squaredGLMM(mod.female.r)

##            R2m       R2c
## [1,] 0.0201185 0.8482584

F_to_eta2(4.6096, 4, 135)

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

tukey’s test of pairwise diet on red female development time

"female black development time Tukey"

## [1] "female black development time Tukey"

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

## $`emmeans of host`
##  host      emmean   SE   df lower.CL upper.CL
##  B. cherry   51.2 4.96 1.45     19.9     82.5
##  Elm         49.6 4.96 1.46     18.5     80.7
##  Hickory     48.4 4.95 1.45     17.3     79.6
##  Malus       51.5 4.97 1.47     20.8     82.1
##  Mix         50.6 4.96 1.45     19.4     81.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        1.596 0.819 133   1.948  0.2974
##  B. cherry - Hickory    2.767 0.798 133   3.469  0.0062
##  B. cherry - Malus     -0.274 0.851 133  -0.322  0.9977
##  B. cherry - Mix        0.574 0.799 133   0.719  0.9518
##  Elm - Hickory          1.171 0.820 133   1.429  0.6103
##  Elm - Malus           -1.869 0.882 133  -2.119  0.2181
##  Elm - Mix             -1.022 0.828 133  -1.234  0.7316
##  Hickory - Malus       -3.040 0.864 133  -3.521  0.0053
##  Hickory - Mix         -2.193 0.809 133  -2.711  0.0577
##  Malus - Mix            0.848 0.863 133   0.983  0.8627
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 5 estimates

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

Male red

subset red males by host

Means, SD, Count, and SE by host

print("       Red male develoment time  (days)        ")

## [1] "       Red male develoment time  (days)        "

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

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

dev.mr.pru<-subset(dev.mr, host == "B. cherry")
       mean(dev.mr.pru$Dev_time, na.rm = TRUE)

## [1] 49.68182

       sd(dev.mr.pru$Dev_time, na.rm = TRUE)

## [1] 7.147336

       length(dev.mr.pru$Dev_time)

## [1] 26

       sd(dev.mr.pru$Dev_time, na.rm = TRUE)/sqrt(length(dev.mr.pru$Dev_time))

## [1] 1.401708

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

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

dev.mr.ulm<-subset(dev.mr, host == "Elm")
      mean(dev.mr.ulm$Dev_time, na.rm = TRUE)

## [1] 46.77778

      sd(dev.mr.ulm$Dev_time, na.rm = TRUE)

## [1] 6.34828

      length(dev.mr.ulm$Dev_time)

## [1] 25

      sd(dev.mr.ulm$Dev_time, na.rm = TRUE)/sqrt(length(dev.mr.ulm$Dev_time))

## [1] 1.269656

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

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

dev.mr.car<-subset(dev.mr, host == "Hickory")
      mean(dev.mr.car$Dev_time, na.rm = TRUE)

## [1] 46.40741

      sd(dev.mr.car$Dev_time, na.rm = TRUE)

## [1] 7.023364

      length(dev.mr.car$Dev_time)

## [1] 32

      sd(dev.mr.car$Dev_time, na.rm = TRUE)/sqrt(length(dev.mr.car$Dev_time))

## [1] 1.241567

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

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

dev.mr.mal<-subset(dev.mr, host =="Malus")
      mean(dev.mr.mal$Dev_time, na.rm = TRUE)

## [1] 47.43333

      sd(dev.mr.mal$Dev_time, na.rm = TRUE)

## [1] 7.280978

      length(dev.mr.mal$Dev_time)

## [1] 35

      sd(dev.mr.mal$Dev_time, na.rm = TRUE)/sqrt(length(dev.mr.mal$Dev_time))

## [1] 1.23071

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

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

dev.mr.mix<-subset(dev.mr, host == "Mix")
       mean(dev.mr.mix$Dev_time, na.rm = TRUE )

## [1] 48.14286

       sd(dev.mr.mix$Dev_time, na.rm = TRUE )

## [1] 5.042675

       length(dev.mr.mix$Dev_time)

## [1] 22

       sd(dev.mr.mix$Dev_time, na.rm = TRUE)/sqrt(length(dev.mr.mix$Dev_time))

## [1] 1.075102

Male Red model– development time

"male red development time"

## [1] "male red development time"

mod.male.r<- lmer(Dev_time~host+(1|Mother)+(1|Father), data=dev.M.R)

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

summary(mod.male.r)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Dev_time ~ host + (1 | Mother) + (1 | Father)
##    Data: dev.M.R
## 
## REML criterion at convergence: 568.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.35249 -0.52109  0.03821  0.49558  2.26504 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Father   (Intercept) 56.552   7.520   
##  Mother   (Intercept)  0.000   0.000   
##  Residual              6.754   2.599   
## Number of obs: 118, groups:  Father, 4; Mother, 3
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  49.56375    3.80094   3.10184  13.040 0.000822 ***
## hostElm       0.52059    0.84241 110.03999   0.618 0.537871    
## hostHickory  -2.65724    0.75202 110.00949  -3.533 0.000601 ***
## hostMalus    -0.51785    0.73328 110.00963  -0.706 0.481552    
## hostMix       0.05447    0.79853 110.01508   0.068 0.945738    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) hstElm hstHck hstMls
## hostElm     -0.097                     
## hostHickory -0.109  0.501              
## hostMalus   -0.110  0.512  0.557       
## hostMix     -0.100  0.469  0.507  0.535
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

anova(mod.male.r)

## Type III Analysis of Variance Table with Satterthwaite's method
##      Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## host 149.92  37.481     4 110.02  5.5497 0.0004161 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

variation contributed by lineage to red male development time

r.squaredGLMM(mod.male.r)

##             R2m       R2c
## [1,] 0.02032802 0.8954843

0.8954-0.203

## [1] 0.6924

host effect size on red males

F_to_eta2(5.5497, 4, 110)

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

tukey’s test pairwise comparison of diet on red males

"male red development time Tukey"

## [1] "male red development time Tukey"

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

## $`emmeans of host`
##  host      emmean   SE   df lower.CL upper.CL
##  B. cherry   49.6 5.16 1.45     17.0     82.2
##  Elm         50.1 5.17 1.46     17.8     82.4
##  Hickory     46.9 5.16 1.44     13.9     80.0
##  Malus       49.0 5.15 1.43     15.9     82.2
##  Mix         49.6 5.16 1.45     17.1     82.1
## 
## 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      -0.5206 0.843 110  -0.618  0.9720
##  B. cherry - Hickory   2.6572 0.752 110   3.533  0.0053
##  B. cherry - Malus     0.5178 0.733 110   0.706  0.9547
##  B. cherry - Mix      -0.0545 0.799 110  -0.068  1.0000
##  Elm - Hickory         3.1778 0.801 110   3.970  0.0012
##  Elm - Malus           1.0384 0.785 110   1.324  0.6771
##  Elm - Mix             0.4661 0.847 110   0.550  0.9817
##  Hickory - Malus      -2.1394 0.700 110  -3.058  0.0229
##  Hickory - Mix        -2.7117 0.771 110  -3.519  0.0056
##  Malus - Mix          -0.5723 0.741 110  -0.773  0.9379
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 5 estimates

knit file

Mattheau Comerford

2023-10-11

R Markdown

Load library

Set seed to 1234

Import data

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

development time analysis

data.pupa comes from data tab

table by sex and color

Table 1: test of Biotypes and sex’s effect on larval development

random effect contribution

Tukeys test

Table 2, SD for development time by sex and color

Subsetting by demographic

subset by sex

subset females by color

subset males by color

black female

subset black females by host

Means, SD, Count, and SE by host

Female black model – development time

variation contributed by lineage to black female development time

host effect size on black female development time

tukey’s test of pairwise diet on female development time

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

black Male

subset black males by host

Means, SD, Count, and SE by host

variation contributed by lineage to black male development time

host effect size on black male development time

tukey’s test of pairwise diet on female development time

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

Red Females

subset red females by host

subset red females by host

Means, SD, Count, and SE by host

Female red model – development time

variation contributed by lineage to red female development time

host effect size on red female development time

tukey’s test of pairwise diet on red female development time

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

Male red

subset red males by host

Means, SD, Count, and SE by host

Male Red model– development time

variation contributed by lineage to red male development time

host effect size on red males

tukey’s test pairwise comparison of diet on red males