R Markdown

library(ggplot2)
library(lme4)
## Loading required package: Matrix
library(MuMIn)
library(vegan)
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.5-7
library(emmeans)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(lmerTest)
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
setwd("C:/Users/jc130714/OneDrive - James Cook University/Documents/Yogi Yasutake/Honours/Post-honours/Statistics and data/Final data")
getwd()
## [1] "C:/Users/jc130714/OneDrive - James Cook University/Documents/Yogi Yasutake/Honours/Post-honours/Statistics and data/Final data"

#Adult physical differences

Adults = read.csv('F1 adult fish_FINAL.csv', strip.white=TRUE)
head(Adults)
##   Tank.. Treatment Dev.treatment Repro.treatment Sex   ID Family Family.number
## 1     17      CC-C             C               C   F CCC8      C             3
## 2     17      CC-C             C               C   M  DCC      D             4
## 3     38      CC-C             C               C   F BCC8      B             2
## 4     38      CC-C             C               C   M  FCC      F             6
## 5     84      CC-C             C               C   F BCC4      B             2
## 6     84      CC-C             C               C   M DCC8      D             4
##   Reproduced SL.cm Weight.g Condition  X X.1 X.2 X.3 X.4 X.5 X.6 X.7 X.8 X.9
## 1          N 9.534    33.90  3.911779 NA  NA  NA  NA  NA  NA  NA  NA  NA  NA
## 2          N 8.710    28.35  4.290408 NA  NA  NA  NA  NA  NA  NA  NA  NA  NA
## 3          N 9.874    43.38  4.506197 NA  NA  NA  NA  NA  NA  NA  NA  NA  NA
## 4          N 8.448    27.63  4.582676 NA  NA  NA  NA  NA  NA  NA  NA  NA  NA
## 5          N 9.504    42.30  4.927437 NA  NA  NA  NA  NA  NA  NA  NA  NA  NA
## 6          N 9.804    40.20  4.265954 NA  NA  NA  NA  NA  NA  NA  NA  NA  NA
##   X.10 X.11 X.12 X.13 X.14 X.15 X.16 X.17 X.18 X.19 X.20 X.21 X.22 X.23 X.24
## 1   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA
## 2   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA
## 3   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA
## 4   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA
## 5   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA
## 6   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA
##   X.25 X.26 X.27 X.28 X.29 X.30 X.31 X.32 Notes
## 1   NA   NA   NA   NA   NA   NA   NA   NA      
## 2   NA   NA   NA   NA   NA   NA   NA   NA      
## 3   NA   NA   NA   NA   NA   NA   NA   NA      
## 4   NA   NA   NA   NA   NA   NA   NA   NA      
## 5   NA   NA   NA   NA   NA   NA   NA   NA      
## 6   NA   NA   NA   NA   NA   NA   NA   NA

#PCA building

AdultsPCA = read.csv('F1 adult fish_dataPCA.csv', strip.white=TRUE)
head(AdultsPCA)
##   SL.cm Weight.g Condition
## 1 9.534    33.90  3.911779
## 2 8.710    28.35  4.290408
## 3 9.874    43.38  4.506197
## 4 8.448    27.63  4.582676
## 5 9.504    42.30  4.927437
## 6 9.804    40.20  4.265954
#standardising
Adults.std = stdize(AdultsPCA, centre=TRUE, scale = TRUE)
head(Adults.std)
##      z.SL.cm z.Weight.g z.Condition
## 1  0.9702193  0.5146636  -0.7918202
## 2 -0.4623866 -0.4018001   0.1814054
## 3  1.5613431  2.0800825   0.7360674
## 4 -0.9178996 -0.5206927   0.9326485
## 5  0.9180613  1.9017437   1.8188185
## 6  1.4396412  1.5549736   0.1185484
Adults.rda = rda(Adults.std)
summary(Adults.rda, display=NULL)
## 
## Call:
## rda(X = Adults.std) 
## 
## Partitioning of variance:
##               Inertia Proportion
## Total               3          1
## Unconstrained       3          1
## 
## Eigenvalues, and their contribution to the variance 
## 
## Importance of components:
##                          PC1    PC2      PC3
## Eigenvalue            1.8877 1.1065 0.005842
## Proportion Explained  0.6292 0.3688 0.001947
## Cumulative Proportion 0.6292 0.9981 1.000000
## 
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores:
screeplot(Adults.rda)

Adults.sites.score = (scores(Adults.rda, display='sites'))
head(Adults.sites.score)
##             PC1         PC2
## sit1 -0.3259455  0.35980439
## sit2  0.1915765 -0.08145263
## sit3 -0.8170804 -0.31441361
## sit4  0.3148028 -0.40937491
## sit5 -0.6408565 -0.78517034
## sit6 -0.6687076 -0.03975941
Adults.species.score = as.data.frame(scores(Adults.rda, display='species'))
Adults.species.score$Species = rownames(Adults.species.score)
Adults.sites.score = data.frame(Adults.sites.score, AdultsPCA)

#getting and adding groupings

AdultGroup = read.csv('F1 adults groupingsPCA.csv')
head(AdultGroup)
##   Treatment Dev.treatment Repro.treatment Sex Family
## 1      CC-C             C               C   F      C
## 2      CC-C             C               C   M      D
## 3      CC-C             C               C   F      B
## 4      CC-C             C               C   M      F
## 5      CC-C             C               C   F      B
## 6      CC-C             C               C   M      D
Adults.sites.score2 = cbind(Adults.sites.score, AdultGroup)
head(Adults.sites.score2)
##             PC1         PC2 SL.cm Weight.g Condition Treatment Dev.treatment
## sit1 -0.3259455  0.35980439 9.534    33.90  3.911779      CC-C             C
## sit2  0.1915765 -0.08145263 8.710    28.35  4.290408      CC-C             C
## sit3 -0.8170804 -0.31441361 9.874    43.38  4.506197      CC-C             C
## sit4  0.3148028 -0.40937491 8.448    27.63  4.582676      CC-C             C
## sit5 -0.6408565 -0.78517034 9.504    42.30  4.927437      CC-C             C
## sit6 -0.6687076 -0.03975941 9.804    40.20  4.265954      CC-C             C
##      Repro.treatment Sex Family
## sit1               C   F      C
## sit2               C   M      D
## sit3               C   F      B
## sit4               C   M      F
## sit5               C   F      B
## sit6               C   M      D
#making a PCA figure 
#putting the sites on the graph
g = ggplot() + geom_point(data=Adults.sites.score2, aes(y=PC2, x=PC1, shape=Treatment, colour= Treatment),  show.legend = TRUE)   + theme_classic()

#adding the arrows
g = g + geom_segment(data=Adults.species.score, aes(y=0, x=0, yend=PC2, xend=PC1), arrow = arrow(length=unit(0.3, 'lines')), show.legend = FALSE)
hjust = ifelse(Adults.species.score$PC1>0,0,1)
vjust = ifelse(Adults.species.score$PC2>0,0,1)
g = g + geom_text(data=Adults.species.score, aes(y=PC2, x=PC1, label = Species))

#now adding cross hairs and axis lines 
g = g + geom_segment(data=NULL, aes(y=-Inf, x=-0, yend=Inf, xend=0))
g = g + geom_segment(data=NULL, aes(y=0, x=-Inf, yend=0, xend=Inf))
g

#lmer testing of PCA

lmer.PC1a = lmer(PC1 ~ Dev.treatment * Repro.treatment + (1|Family) , data = Adults.sites.score2)
summary(lmer.PC1a)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: PC1 ~ Dev.treatment * Repro.treatment + (1 | Family)
##    Data: Adults.sites.score2
## 
## REML criterion at convergence: 100.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8369 -0.4298  0.0159  0.6142  1.9782 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Family   (Intercept) 0.03879  0.1970  
##  Residual             0.15702  0.3963  
## Number of obs: 86, groups:  Family, 6
## 
## Fixed effects:
##                                 Estimate Std. Error       df t value Pr(>|t|)  
## (Intercept)                     -0.03979    0.11258  9.40738  -0.353   0.7316  
## Dev.treatmentH                   0.23002    0.12796 78.46427   1.798   0.0761 .
## Repro.treatmentH                -0.04558    0.11084 77.31546  -0.411   0.6821  
## Dev.treatmentH:Repro.treatmentH -0.04218    0.17754 77.94350  -0.238   0.8129  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Dv.trH Rpr.tH
## Dev.trtmntH -0.422              
## Rpr.trtmntH -0.488  0.444       
## Dv.trtH:R.H  0.305 -0.722 -0.634
anova(lmer.PC1a)
## Type III Analysis of Variance Table with Satterthwaite's method
##                                Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
## Dev.treatment                 0.87177 0.87177     1 77.983  5.5520 0.02097 *
## Repro.treatment               0.09071 0.09071     1 76.547  0.5777 0.44956  
## Dev.treatment:Repro.treatment 0.00886 0.00886     1 77.943  0.0564 0.81285  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
PC1EMM = (emmeans(lmer.PC1a, ~ Dev.treatment * Repro.treatment) %>% as.data.frame) 
PC1EMM
##   Dev.treatment Repro.treatment      emmean        SE       df    lower.CL
## 1             C               C -0.03978547 0.1127554 10.78477 -0.28856412
## 2             H               C  0.19023873 0.1306568 17.54404 -0.08477335
## 3             C               H -0.08536295 0.1133612 10.83946 -0.33532095
## 4             H               H  0.10248618 0.1252105 15.50365 -0.16364046
##    upper.CL
## 1 0.2089932
## 2 0.4652508
## 3 0.1645950
## 4 0.3686128
#Plot
ggplot(PC1EMM, aes(x=Repro.treatment, y=emmean, colour=Dev.treatment))+ geom_point() + geom_linerange(aes (ymin = emmean-SE, ymax = emmean+SE))  +  theme_classic() 

lmer.PC2a = lmer(PC2 ~ Dev.treatment * Repro.treatment + (1|Family) , data = Adults.sites.score2)
summary(lmer.PC2a)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: PC2 ~ Dev.treatment * Repro.treatment + (1 | Family)
##    Data: Adults.sites.score2
## 
## REML criterion at convergence: 99.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3900 -0.5146  0.1451  0.6036  2.2319 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Family   (Intercept) 0.01344  0.1159  
##  Residual             0.16278  0.4035  
## Number of obs: 86, groups:  Family, 6
## 
## Fixed effects:
##                                 Estimate Std. Error       df t value Pr(>|t|)
## (Intercept)                     -0.13093    0.09289 17.58731  -1.410  0.17612
## Dev.treatmentH                   0.08752    0.12947 79.86097   0.676  0.50099
## Repro.treatmentH                 0.33757    0.11252 78.51255   3.000  0.00362
## Dev.treatmentH:Repro.treatmentH -0.36118    0.17987 79.49168  -2.008  0.04804
##                                   
## (Intercept)                       
## Dev.treatmentH                    
## Repro.treatmentH                **
## Dev.treatmentH:Repro.treatmentH * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Dv.trH Rpr.tH
## Dev.trtmntH -0.523              
## Rpr.trtmntH -0.603  0.441       
## Dv.trtH:R.H  0.377 -0.721 -0.632
anova(lmer.PC2a)
## Type III Analysis of Variance Table with Satterthwaite's method
##                                Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
## Dev.treatment                 0.17467 0.17467     1 79.366  1.0730 0.30341  
## Repro.treatment               0.50378 0.50378     1 77.416  3.0948 0.08249 .
## Dev.treatment:Repro.treatment 0.65634 0.65634     1 79.492  4.0320 0.04804 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
PC2EMM = (emmeans(lmer.PC2a, ~ Dev.treatment * Repro.treatment) %>% as.data.frame) 
PC2EMM
##   Dev.treatment Repro.treatment      emmean         SE       df     lower.CL
## 1             C               C -0.13092748 0.09340086 18.25290 -0.326960709
## 2             H               C -0.04340341 0.11441287 31.04538 -0.276736177
## 3             C               H  0.20663943 0.09398270 17.96771  0.009163668
## 4             H               H -0.06701742 0.10829319 27.98340 -0.288851884
##     upper.CL
## 1 0.06510575
## 2 0.18992935
## 3 0.40411520
## 4 0.15481704
#Plot
ggplot(PC2EMM, aes(x=Repro.treatment, y=emmean, colour=Dev.treatment))+ geom_point() + geom_linerange(aes (ymin = emmean-SE, ymax = emmean+SE))  +  theme_classic() 

##Individual models ##Standard length

lmer.SL = lmer(SL.cm ~ Dev.treatment * Repro.treatment + (1|Family) , data = Adults)
summary(lmer.SL)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: SL.cm ~ Dev.treatment * Repro.treatment + (1 | Family)
##    Data: Adults
## 
## REML criterion at convergence: 144.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3610 -0.5419 -0.0211  0.4461  3.1791 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Family   (Intercept) 0.05792  0.2407  
##  Residual             0.27223  0.5218  
## Number of obs: 86, groups:  Family, 6
## 
## Fixed effects:
##                                 Estimate Std. Error       df t value Pr(>|t|)
## (Intercept)                      8.98408    0.14286  9.98394  62.886 2.62e-14
## Dev.treatmentH                  -0.26681    0.16836 78.59118  -1.585    0.117
## Repro.treatmentH                 0.17805    0.14589 77.35559   1.220    0.226
## Dev.treatmentH:Repro.treatmentH -0.08937    0.23363 78.05481  -0.383    0.703
##                                    
## (Intercept)                     ***
## Dev.treatmentH                     
## Repro.treatmentH                   
## Dev.treatmentH:Repro.treatmentH    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Dv.trH Rpr.tH
## Dev.trtmntH -0.438              
## Rpr.trtmntH -0.507  0.444       
## Dv.trtH:R.H  0.317 -0.721 -0.634
anova(lmer.SL)
## Type III Analysis of Variance Table with Satterthwaite's method
##                                Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)   
## Dev.treatment                 1.94002 1.94002     1 78.080  7.1263 0.009237 **
## Repro.treatment               0.36308 0.36308     1 76.513  1.3337 0.251743   
## Dev.treatment:Repro.treatment 0.03984 0.03984     1 78.055  0.1463 0.703106   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
hist(residuals(lmer.SL), col="darkgray")

plot(fitted(lmer.SL), residuals(lmer.SL))

SLEMM = (emmeans(lmer.SL, ~ Dev.treatment * Repro.treatment) %>% as.data.frame) 
SLEMM
##   Dev.treatment Repro.treatment   emmean        SE       df lower.CL upper.CL
## 1             C               C 8.984076 0.1431351 11.58062 8.670954 9.297198
## 2             H               C 8.717266 0.1673804 19.16347 8.367137 9.067395
## 3             C               H 9.162124 0.1439410 11.61816 8.847357 9.476892
## 4             H               H 8.805941 0.1600452 16.92288 8.468158 9.143725
#Plot
ggplot(SLEMM, aes(x=Repro.treatment, y=emmean, colour=Dev.treatment))+ geom_point() + geom_linerange(aes (ymin = emmean-SE, ymax = emmean+SE))  +  theme_classic() 

##Weight

lmer.W = lmer(Weight.g~ Dev.treatment * Repro.treatment + (1|Family) , data = Adults)

summary(lmer.W)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Weight.g ~ Dev.treatment * Repro.treatment + (1 | Family)
##    Data: Adults
## 
## REML criterion at convergence: 534.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8874 -0.6558 -0.0795  0.4815  4.2167 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Family   (Intercept)  8.028   2.833   
##  Residual             31.099   5.577   
## Number of obs: 86, groups:  Family, 6
## 
## Fixed effects:
##                                 Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)                      31.7260     1.6025  9.3061  19.798 6.37e-09
## Dev.treatmentH                   -3.3799     1.8013 78.4514  -1.876   0.0643
## Repro.treatmentH                 -0.4865     1.5600 77.3355  -0.312   0.7560
## Dev.treatmentH:Repro.treatmentH   1.9090     2.4990 77.9396   0.764   0.4472
##                                    
## (Intercept)                     ***
## Dev.treatmentH                  .  
## Repro.treatmentH                   
## Dev.treatmentH:Repro.treatmentH    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Dv.trH Rpr.tH
## Dev.trtmntH -0.417              
## Rpr.trtmntH -0.483  0.444       
## Dv.trtH:R.H  0.302 -0.722 -0.634
anova(lmer.W)
## Type III Analysis of Variance Table with Satterthwaite's method
##                                Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)  
## Dev.treatment                 117.431 117.431     1 77.983  3.7760 0.0556 .
## Repro.treatment                 4.470   4.470     1 76.593  0.1437 0.7056  
## Dev.treatment:Repro.treatment  18.148  18.148     1 77.940  0.5835 0.4472  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
hist(residuals(lmer.W), col="darkgray")

plot(fitted(lmer.W), residuals(lmer.W))

WEMM = (emmeans(lmer.W, ~ Dev.treatment * Repro.treatment) %>% as.data.frame) 
WEMM
##   Dev.treatment Repro.treatment   emmean       SE       df lower.CL upper.CL
## 1             C               C 31.72604 1.604837 10.56602 28.17603 35.27604
## 2             H               C 28.34613 1.854547 17.09140 24.43497 32.25729
## 3             C               H 31.23953 1.613327 10.62448 27.67325 34.80581
## 4             H               H 29.76864 1.778452 15.11047 25.98038 33.55691
#Plot
ggplot(WEMM, aes(x=Repro.treatment, y=emmean, colour=Dev.treatment))+ geom_point() + geom_linerange(aes (ymin = emmean-SE, ymax = emmean+SE))  +  theme_classic() 

##Condition

lmer.Cond = lmer(Condition~ Dev.treatment * Repro.treatment + (1|Family), data= Adults)

summary(lmer.Cond)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Condition ~ Dev.treatment * Repro.treatment + (1 | Family)
##    Data: Adults
## 
## REML criterion at convergence: 82.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2628 -0.6172 -0.1655  0.5387  3.4198 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Family   (Intercept) 0.01099  0.1048  
##  Residual             0.13188  0.3632  
## Number of obs: 86, groups:  Family, 6
## 
## Fixed effects:
##                                 Estimate Std. Error       df t value Pr(>|t|)
## (Intercept)                      4.33816    0.08371 17.51595  51.826  < 2e-16
## Dev.treatmentH                  -0.08360    0.11654 79.85392  -0.717  0.47527
## Repro.treatmentH                -0.29878    0.10128 78.50611  -2.950  0.00418
## Dev.treatmentH:Repro.treatmentH  0.31607    0.16191 79.48247   1.952  0.05444
##                                    
## (Intercept)                     ***
## Dev.treatmentH                     
## Repro.treatmentH                ** 
## Dev.treatmentH:Repro.treatmentH .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Dv.trH Rpr.tH
## Dev.trtmntH -0.523              
## Rpr.trtmntH -0.602  0.441       
## Dv.trtH:R.H  0.377 -0.721 -0.632
anova(lmer.Cond)
## Type III Analysis of Variance Table with Satterthwaite's method
##                                Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
## Dev.treatment                 0.11173 0.11173     1 79.359  0.8472 0.36014  
## Repro.treatment               0.40499 0.40499     1 77.412  3.0709 0.08367 .
## Dev.treatment:Repro.treatment 0.50259 0.50259     1 79.482  3.8110 0.05444 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
hist(residuals(lmer.Cond), col="darkgray")

plot(fitted(lmer.Cond), residuals(lmer.Cond))

FKEMM = (emmeans(lmer.Cond, ~ Dev.treatment * Repro.treatment) %>% as.data.frame) 
FKEMM
##   Dev.treatment Repro.treatment   emmean         SE       df lower.CL upper.CL
## 1             C               C 4.338163 0.08416424 18.18172 4.161467 4.514859
## 2             H               C 4.254564 0.10306517 30.93334 4.044343 4.464785
## 3             C               H 4.039382 0.08468873 17.90143 3.861387 4.217377
## 4             H               H 4.271853 0.09755838 27.87370 4.071972 4.471733
#Plot
ggplot(FKEMM, aes(x=Repro.treatment, y=emmean, colour=Dev.treatment))+ geom_point() + geom_linerange(aes (ymin = emmean-SE, ymax = emmean+SE))  +  theme_classic() 

##Physical condition, weight for a given length

lmer.Cond = lmer(Weight.g ~ Dev.treatment * Repro.treatment + (1|Family) + SL.cm , data = Adults)

summary(lmer.Cond)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Weight.g ~ Dev.treatment * Repro.treatment + (1 | Family) + SL.cm
##    Data: Adults
## 
## REML criterion at convergence: 403.3
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.21058 -0.62911 -0.09576  0.45871  2.92659 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Family   (Intercept) 0.8586   0.9266  
##  Residual             6.6231   2.5735  
## Number of obs: 86, groups:  Family, 6
## 
## Fixed effects:
##                                 Estimate Std. Error       df t value Pr(>|t|)
## (Intercept)                     -54.5204     4.8786  80.0104 -11.175  < 2e-16
## Dev.treatmentH                   -0.8054     0.8396  79.3540  -0.959  0.34034
## Repro.treatmentH                 -2.1754     0.7255  76.7488  -2.998  0.00366
## SL.cm                             9.6000     0.5384  80.5186  17.832  < 2e-16
## Dev.treatmentH:Repro.treatmentH   2.7506     1.1512  77.7467   2.389  0.01931
##                                    
## (Intercept)                     ***
## Dev.treatmentH                     
## Repro.treatmentH                ** 
## SL.cm                           ***
## Dev.treatmentH:Repro.treatmentH *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Dv.trH Rpr.tH SL.cm 
## Dev.trtmntH -0.226                     
## Rpr.trtmntH  0.064  0.409              
## SL.cm       -0.991  0.165 -0.137       
## Dv.trtH:R.H -0.003 -0.702 -0.633  0.049
anova(lmer.Cond)
## Type III Analysis of Variance Table with Satterthwaite's method
##                                Sum Sq Mean Sq NumDF  DenDF  F value  Pr(>F)    
## Dev.treatment                    6.01    6.01     1 78.787   0.9080 0.34355    
## Repro.treatment                 12.88   12.88     1 76.609   1.9441 0.16726    
## SL.cm                         2105.95 2105.95     1 80.519 317.9713 < 2e-16 ***
## Dev.treatment:Repro.treatment   37.81   37.81     1 77.747   5.7085 0.01931 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
hist(residuals(lmer.Cond), col="darkgray")

plot(fitted(lmer.Cond), residuals(lmer.Cond))

CondEMM = (emmeans(lmer.Cond, ~ Dev.treatment * Repro.treatment) %>% as.data.frame) 
CondEMM
##   Dev.treatment Repro.treatment   emmean        SE       df lower.CL upper.CL
## 1             C               C 31.64888 0.6378392 14.39579 30.28437 33.01339
## 2             H               C 30.84350 0.7805864 24.84596 29.23535 32.45165
## 3             C               H 29.47351 0.6490274 15.25794 28.09217 30.85485
## 4             H               H 31.41871 0.7342391 22.32097 29.89726 32.94016
#Plot
ggplot(CondEMM, aes(x=Repro.treatment, y=emmean, colour=Dev.treatment))+ geom_point() + geom_linerange(aes (ymin = emmean-SE, ymax = emmean+SE))  +  theme_classic() 

##Proportion of pairs breeding

Breeding = read.csv('Breeding proportions FINAL.csv', strip.white=TRUE)
head(Breeding)
##   Development Post.maturity Treatment.combo Reproduced
## 1     Control       Control              CC          N
## 2     Control       Control              CC          N
## 3     Control       Control              CC          N
## 4     Control       Control              CC          N
## 5     Control       Control              CC          N
## 6     Control       Control              CC          N
table(Breeding$Treatment.combo, Breeding$Reproduced)
##     
##      N Y
##   CC 7 6
##   CH 7 6
##   HC 3 5
##   HH 6 3
chisq.test(table(Breeding$Treatment.combo, Breeding$Reproduced))
## Warning in chisq.test(table(Breeding$Treatment.combo, Breeding$Reproduced)):
## Chi-squared approximation may be incorrect
## 
##  Pearson's Chi-squared test
## 
## data:  table(Breeding$Treatment.combo, Breeding$Reproduced)
## X-squared = 1.4516, df = 3, p-value = 0.6935
chisq.test(table(Breeding$Development, Breeding$Reproduced))
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  table(Breeding$Development, Breeding$Reproduced)
## X-squared = 0, df = 1, p-value = 1
chisq.test(table(Breeding$Post.maturity, Breeding$Reproduced))
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  table(Breeding$Post.maturity, Breeding$Reproduced)
## X-squared = 0.20077, df = 1, p-value = 0.6541