Date rendered: Tue Sep 24 08:45:17 2019

Objectives

The goal of this project is to investigate sex assortment. (1) Does sex assortment differ between the two resource dispersion treatments? (2) Does sex assortment impact population fitness? (3) How does sex assortment influence selection on male body size?

NOTE (1) This data is not fully error-checked. (2) The social network metrics are calcualted from social networks that include beetles that died and moved (3MZ) during the experiment. But the analyses do not include beetles that died or 3MZ. (3) The beetle size measurements are elytra measurements done by Ruth and Patrick (work study students). In the future, we will have horn measurements as well and can calculate PCA. (4) These analyses ignore the lack of indepdence in network metrics. In the future, I need to incorporate permutations.

Set-Up and Data

#clear memory
rm(list=ls())

#libraries
library(lme4)
## Loading required package: Matrix
library(car)
## Loading required package: carData
## Registered S3 methods overwritten by 'car':
##   method                          from
##   influence.merMod                lme4
##   cooks.distance.influence.merMod lme4
##   dfbeta.influence.merMod         lme4
##   dfbetas.influence.merMod        lme4
library(glmmTMB)
library(DHARMa)
## Registered S3 method overwritten by 'DHARMa':
##   method        from   
##   refit.glmmTMB glmmTMB
library(effects)
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
library(ggeffects)
library(tidyverse)
## ── Attaching packages ──────────────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──
## ✔ ggplot2 3.2.0     ✔ purrr   0.3.2
## ✔ tibble  2.1.3     ✔ dplyr   0.8.3
## ✔ tidyr   0.8.3     ✔ stringr 1.4.0
## ✔ readr   1.3.1     ✔ forcats 0.4.0
## ── Conflicts ─────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::expand() masks Matrix::expand()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ✖ dplyr::recode() masks car::recode()
## ✖ purrr::some()   masks car::some()
#function that excludes
`%sans%` <- Negate(`%in%`)

#Change R's defult to contrasts 
options(contrasts=c("contr.sum", "contr.poly")) 

#FUNCTION: check.assumptions
check.assumptions <- function(mod){
  # Multiple panels
  par(mfrow=c(2,2))
  
  # Normality
  hist(resid(mod))
  qqnorm(resid(mod))
  qqline(resid(mod))
  
  # Linearity
  # --> Plot residual versus fitted values with loess line
  plot(fitted(mod), resid(mod), main = "Residual vs fitted values")
  fr.fit <- loess(resid(mod)~fitted(mod))
  fr.vec <- seq(min(fitted(mod)), max(fitted(mod)), length.out=100)
  lines(fr.vec,predict(fr.fit,fr.vec), col="red")
  
  # Homoscedasticity
  # --> Plot the absolute value of the residuals against the fitted values
  plot(fitted(mod), sqrt(abs(resid(mod))), main="Scale-location")
  sl.fit <- loess(sqrt(abs(resid(mod)))~fitted(mod))
  sl.vec <- seq(min(fitted(mod)), max(fitted(mod)), length.out=100)
  lines(sl.vec,predict(sl.fit,sl.vec), col="red")
}

#upload individual-level social network data
IndividSN<-read.csv("socialnetworks_20190513_.csv")

#upload population-level social network data
PopSN<-read.csv("PopSN_20190513_.csv")

#upload surveys dataset
surveys_2018<-read.csv("20190422_OffMountainBeetleBase2.0.csv")
#NOTE: This is not a fully cleaned dataset.

#upload beetle measurements
elytra<-read.csv("BeetleMeasurements_20190324.csv")
#NOTE: This elytra data is not completely cleaned.
#(1) Need to check with Hannah Donald about how to measure elytra 3FM and 2MX. These are dirty scans. Ruth and Patrick measured these scans very differently. I already emailed Hannah about this (20190221) and just need to follow-up by looking at these measurements. Right now, these are wrong measurements!
#(2) Missing dorsal scan for 3NH. For now, I am using Rachel's measurement of the elytra. I do not have a pronotum measurement for 3NH. Luckily, this beetle died during Period 1. Rachel's measurement: 7.281

#upload bracket attributes dataset
#use this datset for total eggs laid
bracketattributes<-read.csv("bracketattributes_20190602_.csv")

#create objects for ease of coding
DeadBeetles <- c("2BO", "2BX", "2H2", "2M6", "2MS", "2N6", "2N7", "2OH", "2OW", "2PK", "2PW", "2T9", "2TS", "2VB", "2VO", "2WV", "3A6", "3AB", "3AW", "3HK", "3HX", "3KA", "3LV", "3LX", "3MN", "3NH", "3ON", "3VO", "3W9", "3XM", "3XN", "2FF", "2HF", "3F3", "3KL", "3YF")

Dispersed <- c("1A_Period1", "2A_Period1", "3A_Period1", "4A_Period1", "5A_Period1", "6A_Period1", "1B_Period2", "2B_Period2", "3B_Period2", "4B_Period2", "5B_Period2", "6B_Period2")
Clumped <- c("1B_Period1", "2B_Period1", "3B_Period1", "4B_Period1", "5B_Period1", "6B_Period1", "1A_Period2", "2A_Period2", "3A_Period2", "4A_Period2", "5A_Period2", "6A_Period2")

Clean Up Population Social Network Dataset

Note: These Pop SN Metrics are calculated from SRI matrices that include all males and females.

PopSN %>%
  mutate(Condo = substr(Condo_Period, 1, 2)) %>%
  mutate(Period = substr(Condo_Period, 4, 10)) %>%
  mutate(Treatment = ifelse(Condo_Period %in% Dispersed, "Dispersed", "Clumped")) %>%
  mutate(Condo = as.factor(Condo)) %>%
  mutate(Period = as.factor(Period)) %>%
  mutate(Treatment = as.factor(Treatment)) -> PopSN

Does sex assortment differ between the two resource dispersion treatments?

This is the only resource distribution effect. I tested the effect of resource dispersion on 6 other network-level metrics (average strength, average shortest path, global clustering coefficient, elytra assortment, coefficient of variation) Sexes are more negatively assorted when resources are dispersed.

hist(PopSN$sex.assortment)

sex.assortment.comparison <- lm(sex.assortment~Treatment+Period+Condo, data=PopSN)
summary(sex.assortment.comparison)
## 
## Call:
## lm(formula = sex.assortment ~ Treatment + Period + Condo, data = PopSN)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.04571 -0.01374  0.00000  0.01374  0.04571 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.1537828  0.0069027 -22.279 7.46e-10 ***
## Treatment1   0.0159241  0.0069027   2.307   0.0437 *  
## Period1     -0.0053747  0.0069027  -0.779   0.4542    
## Condo1      -0.0318395  0.0228936  -1.391   0.1945    
## Condo2       0.0516261  0.0228936   2.255   0.0478 *  
## Condo3      -0.0002418  0.0228936  -0.011   0.9918    
## Condo4      -0.0215107  0.0228936  -0.940   0.3696    
## Condo5       0.0078469  0.0228936   0.343   0.7389    
## Condo6      -0.0238871  0.0228936  -1.043   0.3213    
## Condo7       0.0310010  0.0228936   1.354   0.2055    
## Condo8       0.0403012  0.0228936   1.760   0.1088    
## Condo9      -0.0109768  0.0228936  -0.479   0.6419    
## Condo10     -0.0232233  0.0228936  -1.014   0.3343    
## Condo11     -0.0079482  0.0228936  -0.347   0.7357    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03382 on 10 degrees of freedom
## Multiple R-squared:  0.6698, Adjusted R-squared:  0.2405 
## F-statistic:  1.56 on 13 and 10 DF,  p-value: 0.2433
Anova(sex.assortment.comparison, type=3)
## Anova Table (Type III tests)
## 
## Response: sex.assortment
##              Sum Sq Df  F value    Pr(>F)    
## (Intercept) 0.56758  1 496.3423 7.456e-10 ***
## Treatment   0.00609  1   5.3220   0.04373 *  
## Period      0.00069  1   0.6063   0.45422    
## Condo       0.01641 11   1.3048   0.34125    
## Residuals   0.01144 10                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(sex.assortment.comparison))

check.assumptions(sex.assortment.comparison)

Analyze within periods

Subset by Period

PopSN %>%
  filter(Period == "Period1") %>%
  droplevels() -> PopSN1

PopSN %>%
  filter(Period == "Period2") %>%
  droplevels() -> PopSN2

Period 1

hist(PopSN1$sex.assortment)

sex.assortment.comparison.1 <- lm(sex.assortment~Treatment, data=PopSN1)
summary(sex.assortment.comparison.1)
## 
## Call:
## lm(formula = sex.assortment ~ Treatment, data = PopSN1)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.058617 -0.027016  0.007849  0.018012  0.062285 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.15916    0.01136  -14.01 6.72e-08 ***
## Treatment1   0.01795    0.01136    1.58    0.145    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03935 on 10 degrees of freedom
## Multiple R-squared:  0.1998, Adjusted R-squared:  0.1198 
## F-statistic: 2.497 on 1 and 10 DF,  p-value: 0.1451
Anova(sex.assortment.comparison.1, type=3)
## Anova Table (Type III tests)
## 
## Response: sex.assortment
##               Sum Sq Df  F value    Pr(>F)    
## (Intercept) 0.303973  1 196.3080 6.722e-08 ***
## Treatment   0.003867  1   2.4971    0.1451    
## Residuals   0.015484 10                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(sex.assortment.comparison.1))

check.assumptions(sex.assortment.comparison.1)
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.17729
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.03608
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 5.294e-19
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 0.0013018
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.17729
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.03608
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 5.294e-19
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 0.0013018

Period 2

hist(PopSN2$sex.assortment)

sex.assortment.comparison.2 <- lm(sex.assortment~Treatment, data=PopSN2)
summary(sex.assortment.comparison.2)
## 
## Call:
## lm(formula = sex.assortment ~ Treatment, data = PopSN2)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.071627 -0.012300 -0.001069  0.010347  0.059133 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.14841    0.01011 -14.680  4.3e-08 ***
## Treatment1   0.01390    0.01011   1.375    0.199    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03502 on 10 degrees of freedom
## Multiple R-squared:  0.1589, Adjusted R-squared:  0.07483 
## F-statistic:  1.89 on 1 and 10 DF,  p-value: 0.1993
Anova(sex.assortment.comparison.2, type=3)
## Anova Table (Type III tests)
## 
## Response: sex.assortment
##               Sum Sq Df  F value    Pr(>F)    
## (Intercept) 0.264300  1 215.4891 4.302e-08 ***
## Treatment   0.002318  1   1.8897    0.1993    
## Residuals   0.012265 10                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(sex.assortment.comparison.2))

check.assumptions(sex.assortment.comparison.2)
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.16244
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.027934
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 1.1927e-17
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 0.00078032
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.16244
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.027934
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 1.1927e-17
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 0.00078032

Clean Up Surveys Dataset

Rename X_fkfieldID_surveys

surveys_2018 %>%
  dplyr::rename(ID = X_fkfieldID_surveys) -> surveys_2018

Divide Dataset into two Periods

surveys_2018 %>%
  mutate(Period = ifelse(Date_Surveys < 20180719, "Period1", "Period2")) %>%
  mutate(Period = as.factor(Period)) -> surveys_2018

Period1<-subset(surveys_2018, Period=="Period1")
Period2<-subset(surveys_2018, Period=="Period2")

Subset by Sex

#subset data for males
Period1Males<-subset(Period1, Survey_Sex=="M")
Period2Males<-subset(Period2, Survey_Sex=="M")

#subset data for females
Period1Females<-subset(Period1, Survey_Sex=="F")
Period2Females<-subset(Period2, Survey_Sex=="F")

Collapse dataset so that each beetle has one row with the number of lays, number of guards, number of scans seen, condo_period, elytra

Number of Lays

Period1Females %>%
  #grab one row for each beetle ID during each Scan. need to grab one row so that don't overcount number of lays when beetles have multiple partners. 
  group_by(ID, ScanID) %>%
  filter(row_number() == 1) %>%
  #create a column that has number of lays
  ungroup() %>%
  group_by(ID) %>%
  dplyr::mutate(Lays = length(which(Behavior == "L"))) %>%
  #create dataset with one row per beetle with number of lays
  filter(row_number() == 1) %>%
  dplyr::select(ID, Lays) -> Period1FemalesLays

#merge number of lays into Period1Females
Period1Females<-merge(Period1Females, Period1FemalesLays, by="ID", all.x=TRUE, all.y=FALSE)

Period2Females %>%
  #grab one row for each beetle ID during each Scan. need to grab one row so that don't overcount number of lays when beetles have multiple partners. 
  group_by(ID, ScanID) %>%
  filter(row_number() == 1) %>%
  #create a column that has number of lays
  ungroup() %>%
  group_by(ID) %>%
  dplyr::mutate(Lays = length(which(Behavior == "L"))) %>% 
  #create dataset with one row per beetle with number of lays
  filter(row_number() == 1) %>%
  dplyr::select(ID, Lays) -> Period2FemalesLays

#merge number of lays into Period2Females
Period2Females<-merge(Period2Females, Period2FemalesLays, by="ID", all.x=TRUE, all.y=FALSE)

Number of Guards

Period1Males %>%
  #grab one row for each beetle ID during each Scan. need to grab one row so that don't overcount number of guards when beetles have multiple partners. 
  group_by(ID, ScanID) %>%
  filter(row_number() == 1) %>%
  #create a column that has number of guards
  ungroup() %>%
  group_by(ID) %>%
  dplyr::mutate(Guards = length(which(Behavior == "G"))) %>%
  #create dataset with one row per beetle with number of guards
  filter(row_number() == 1) %>%
  dplyr::select(ID, Guards) -> Period1MalesGuards

#merge number of guards into Period1Males
Period1Males<-merge(Period1Males, Period1MalesGuards, by="ID", all.x=TRUE, all.y=FALSE)

Period2Males %>%
  #grab one row for each beetle ID during each Scan. need to grab one row so that don't overcount number of guards when beetles have multiple partners. 
  group_by(ID, ScanID) %>%
  filter(row_number() == 1) %>%
  #create a column that has number of guards
  ungroup() %>%
  group_by(ID) %>%
  dplyr::mutate(Guards = length(which(Behavior == "G"))) %>%
  #create dataset with one row per beetle with number of guards
  filter(row_number() == 1) %>%
  dplyr::select(ID, Guards) -> Period2MalesGuards

#merge number of guards into Period2Males
Period2Males<-merge(Period2Males, Period2MalesGuards, by="ID", all.x=TRUE, all.y=FALSE)

Number of Unique Guards

#issue with this defintion of unqiue guards: it asks whether a beetle is guarding the same female that it was guarding during the last scan that the beetle was seen, not necessarily the last scan 
Period1Males %>%
  #grab one row for each beetle ID during each Scan. need to grab one row so that don't overcount number of guards when beetles have multiple partners. make sure to grab the row with the mating partner. need to grab the row with the mating partner to see if two guards on consecutive scans are acutally unique guards with different mating partners.
  group_by(ID, ScanID) %>%
  mutate(PartnerInteractionRank = ifelse(PartnerInteraction == "Mating Partners", 1, 2)) %>%
  dplyr::arrange(ID, ScanID, PartnerInteractionRank) %>%
  filter(row_number() == 1) %>%
  #create a column that has unique Gs (not counting guards with the same partner for two consecutive scans). to create the UniqueGuards column, compare the Behavior_Partner of one scan to the Behavior_Partner of the previous scan when that beetle was seen.
  ungroup() %>%
  group_by(ID) %>%
  mutate(Behavior_Partner = paste(Behavior, Partner, sep="_")) %>%
  mutate(UniqueBehaviors = ifelse(Behavior_Partner == lag(Behavior_Partner, default=first(Behavior)), "Same Behavior", as.character(Behavior))) %>%
  #count all unique Gs for each beetle ID
  dplyr::mutate(UniqueGuards = length(which(UniqueBehaviors == "G"))) %>%
  #create dataset with one row per beetle with number of unique guards
  filter(row_number() == 1) %>%
  dplyr::select(ID, UniqueGuards) -> Period1MalesUniqueGuards
#NOTE: this code doesn't account for a scenario where a beetle is seen guarding, then not seen for several scans, and then seen guarding the same beetle. VAF assigned a number to each scan and checked to see if the lag comparison increases by one or more than one.

#merge number of unique guards into Period1Males
Period1Males<-merge(Period1Males, Period1MalesUniqueGuards, by="ID", all.x=TRUE, all.y=FALSE)

Period2Males %>%
  #grab one row for each beetle ID during each Scan. need to grab one row so that don't overcount number of guards when beetles have multiple partners. make sure to grab the row with the mating partner. need to grab the row with the mating partner to see if two guards on consecutive scans are acutally unique guards with different mating partners.
  group_by(ID, ScanID) %>%
  mutate(PartnerInteractionRank = ifelse(PartnerInteraction == "Mating Partners", 1, 2)) %>%
  dplyr::arrange(ID, ScanID, PartnerInteractionRank) %>%
  filter(row_number() == 1) %>%
  #create a column that has unique Gs (not counting guards with the same partner for two consecutive scans). to create the UniqueGuards column, compare the Behavior_Partner of one scan to the Behavior_Partner of the previous scan when that beetle was seen.
  ungroup() %>%
  group_by(ID) %>%
  mutate(Behavior_Partner = paste(Behavior, Partner, sep="_")) %>%
  mutate(UniqueBehaviors = ifelse(Behavior_Partner == lag(Behavior_Partner, default=first(Behavior)), "Same Behavior", as.character(Behavior))) %>%
  #count all unique Gs for each beetle ID
  dplyr::mutate(UniqueGuards = length(which(UniqueBehaviors == "G"))) %>%
  #create dataset with one row per beetle with number of unique guards
  filter(row_number() == 1) %>%
  dplyr::select(ID, UniqueGuards) -> Period2MalesUniqueGuards

#merge number of unique guards into Period1Males
Period2Males<-merge(Period2Males, Period2MalesUniqueGuards, by="ID", all.x=TRUE, all.y=FALSE)

Number of scans seen

Period1Males %>%
  group_by(ID) %>%
  mutate(ScansSeen = length(unique(ScanID)))-> Period1Males

Period2Males %>%
  group_by(ID) %>%
  mutate(ScansSeen = length(unique(ScanID))) -> Period2Males

Period1Females %>%
  group_by(ID) %>%
  mutate(ScansSeen = length(unique(ScanID)))-> Period1Females

Period2Females %>%
  group_by(ID) %>%
  mutate(ScansSeen = length(unique(ScanID))) -> Period2Females

Collapse Periods 1 and 2

surveydataM<-bind_rows(Period1Males, Period2Males)
surveydataF<-bind_rows(Period1Females, Period2Females)
surveydata<-bind_rows(surveydataM, surveydataF)

Remove UKM

surveydata %>%
  #remove UKM
  filter(ID != "UKM") %>%
  filter(ID != "UKF") %>%
  filter(ID != "UK") -> surveydata

Remove beetles that died from the analyses

These beetles are currently included in the social networks, just removed from the analyses

surveydata %>%
  #remove dead beetles
  filter(ID %sans% DeadBeetles) %>%
  #remove 3MZ
  filter(ID != "3MZ") %>%
  #remove 3XY from 1B
  mutate(ID_Condo = paste(ID, Condo, sep="_")) %>%
  filter(ID_Condo != "3XY_1B") %>%
  #droplevels
  droplevels() -> surveydata

Add Elytra to Dataset

In the future, can use a PCA for body size. For now, use elytra measurements done by Ruth and Patrick (work-study students).

#remove pronotum and horn measurements
elytra %>%
  filter(Item.Measured == "Elytra") -> elytra

#add Rachel's measurement for 3NH. 3NH is missing a dorsal scan.
MissingDorsal <- c("3NH-D.jpg")
elytra %>%
  mutate(Average = replace(Average, Image.Name %in% MissingDorsal, 7.281)) -> elytra

#create simplified elytra dataset: just ID and Elytra Length
elytra %>%
  #create ID column
  mutate(ID = as.factor(substr(Image.Name, 1, 3))) %>%
  #rename elytra column
  dplyr::rename(Elytra = Average) %>%
  #simplify the dataset to just ID and Elytra Length
  dplyr::select(ID, Elytra) -> elytra

#merge body size measures to raw dataset
surveydata <- merge(surveydata, elytra, by=("ID"), all.x=TRUE, all.y=FALSE)

Collapse data by beetle

#select only one row for each beetle_period
surveydata %>%
  #make a ID_Period column
  mutate(ID_Period = as.factor(paste(ID, Period, sep="_"))) %>%
  #select only one row for each beetle_period
  distinct(ID_Period, .keep_all = TRUE) -> surveydataBeetleset

Remove Extra Columns

surveydataBeetleset %>%
  dplyr::select(ID, Elytra, Survey_Sex, ScansSeen, Lays, UniqueGuards, Period, Condo) -> surveydataBeetleset

Standardize variables for selection analyses

I’m standardizing body size and number of scans seen within condo and period I’m standardizing sex assortment within period Because these are poisson models, I have not standardized the number of unique guards. But I will need to standardize this fitness metric to calculate selection gradients.

#standardize sex assortment
PopSN %>%
  group_by(Period) %>%
  mutate(StandardizedSexAssortment = ((sex.assortment-(mean(sex.assortment)))/sd(sex.assortment))) %>%
  ungroup() -> PopSN

#standardize elytra
surveydataBeetleset %>%
  group_by(Period, Condo, Survey_Sex) %>%
  mutate(StandardizedElytra = ((Elytra-(mean(Elytra)))/sd(Elytra))) %>%
  ungroup() -> surveydataBeetleset

#standardize scans seen
surveydataBeetleset %>%
  group_by(Period, Condo, Survey_Sex) %>%
  mutate(StandardizedScansSeen = ((ScansSeen-(mean(ScansSeen)))/sd(ScansSeen))) %>%
  ungroup() -> surveydataBeetleset

Merge PopSN to surveydataBeetleset

#create Condo_Period column in surveydataMBeetleset
surveydataBeetleset %>%
  mutate(Condo_Period = as.factor(paste(Condo, Period, sep="_"))) -> surveydataBeetleset

#select for only Condo_Period and sex assortment
PopSN %>%
  dplyr::select(Condo_Period, Treatment, StandardizedSexAssortment) -> PopSNSimple

#merge sex assortment to beetleset
surveydataBeetleset <- merge(surveydataBeetleset, PopSNSimple, by=("Condo_Period"), all.x=TRUE, all.y=FALSE)

#droplevels
surveydataBeetleset %>%
  droplevels() -> surveydataBeetleset

Subset by Sex

#subset data for males
surveydataMBeetleset<-subset(surveydataBeetleset, Survey_Sex=="M")

#subset data for females
surveydataFBeetleset<-subset(surveydataBeetleset, Survey_Sex=="F")

Is there group level selection?

Calculate Total Unique Guards and Total Lays for each Condo

surveydataMBeetleset %>%
  #group by Condo and Period
  group_by(Condo, Period) %>%
  #sum up all Unique Guards
  mutate(TotalUniqueGuards = sum(UniqueGuards)) -> surveydataMBeetleset
 
surveydataFBeetleset %>%
  #group by Condo and Period
  group_by(Condo, Period) %>%
  #sum up all Lays
  mutate(TotalLays = sum(Lays)) -> surveydataFBeetleset

Select only one row from datasets to get a population-level dataset

surveydataMBeetleset %>%
  mutate(Condo_Period = paste(Condo, Period, sep="_")) %>%
  group_by(Condo_Period) %>%
  distinct(Condo_Period, .keep_all=TRUE) -> PopMBeetleset

surveydataFBeetleset %>%
  mutate(Condo_Period = paste(Condo, Period, sep="_")) %>%
  group_by(Condo_Period) %>%
  distinct(Condo_Period, .keep_all=TRUE) -> PopFBeetleset

Merge in Egg Counts

#create dataset with number of eggs for each Condo_Period
bracketattributes %>%
  group_by(Condo, Period) %>%
  mutate(TotalEggs = sum(Number.of.Eggs)) %>%
  #select only one row for each Condo_Period
  filter(row_number() == 1) %>%
  ungroup() -> Eggs

#merge TotalEggs to dataset
Eggs %>%
  mutate(Period = ifelse(Period==1, "Period1", "Period2")) %>%
  mutate(Condo_Period = paste(Condo, Period, sep="_"))  %>%   
  dplyr::select(Condo_Period, TotalEggs) -> Eggs

PopFBeetleset<-merge(PopFBeetleset, Eggs, by="Condo_Period", all.x=TRUE, all.y=FALSE)

Selection on Sex Assortment

Guards

guards.sex.assortment.selection<-glmer(TotalUniqueGuards~StandardizedSexAssortment+(1|Condo), data=PopMBeetleset, family="poisson")
summary(guards.sex.assortment.selection)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: TotalUniqueGuards ~ StandardizedSexAssortment + (1 | Condo)
##    Data: PopMBeetleset
## 
##      AIC      BIC   logLik deviance df.resid 
##    164.9    168.4    -79.4    158.9       21 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.4629 -0.4309  0.0537  0.4003  1.3652 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Condo  (Intercept) 0.006871 0.08289 
## Number of obs: 24, groups:  Condo, 12
## 
## Fixed effects:
##                           Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                3.83735    0.03845   99.81   <2e-16 ***
## StandardizedSexAssortment -0.03055    0.03434   -0.89    0.374    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## StndrdzdSxA 0.018
Anova(guards.sex.assortment.selection, type=3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: TotalUniqueGuards
##                               Chisq Df Pr(>Chisq)    
## (Intercept)               9961.1457  1     <2e-16 ***
## StandardizedSexAssortment    0.7913  1     0.3737    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(guards.sex.assortment.selection))

check.assumptions(guards.sex.assortment.selection)

#test for zero inflation and overdispersion with DHARMa
guards.sex.assortment.selection.resid<-simulateResiduals(guards.sex.assortment.selection, n=250)
plot(guards.sex.assortment.selection.resid)

testDispersion(guards.sex.assortment.selection.resid)
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted
##  vs. simulated
## 
## data:  simulationOutput
## ratioObsSim = 0.91302, p-value = 0.608
## alternative hypothesis: two.sided
testZeroInflation(guards.sex.assortment.selection.resid)
## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sided
#the model is not overdispersed nor zero-inflated.

Lays

lays.sex.assortment.selection<-glmer(TotalLays~StandardizedSexAssortment+(1|Condo), data=PopFBeetleset, family="poisson")
summary(lays.sex.assortment.selection)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: TotalLays ~ StandardizedSexAssortment + (1 | Condo)
##    Data: PopFBeetleset
## 
##      AIC      BIC   logLik deviance df.resid 
##    164.7    168.3    -79.4    158.7       21 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.28825 -0.52249 -0.09117  0.28264  2.72172 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Condo  (Intercept) 0.0253   0.1591  
## Number of obs: 24, groups:  Condo, 12
## 
## Fixed effects:
##                           Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                3.38767    0.05950  56.934   <2e-16 ***
## StandardizedSexAssortment  0.02653    0.04745   0.559    0.576    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## StndrdzdSxA -0.006
Anova(lays.sex.assortment.selection, type=3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: TotalLays
##                               Chisq Df Pr(>Chisq)    
## (Intercept)               3241.5176  1     <2e-16 ***
## StandardizedSexAssortment    0.3126  1     0.5761    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(lays.sex.assortment.selection))

check.assumptions(lays.sex.assortment.selection)

#test for zero inflation and overdispersion with DHARMa
lays.sex.assortment.selection.resid<-simulateResiduals(lays.sex.assortment.selection, n=250)
plot(lays.sex.assortment.selection.resid)

testDispersion(lays.sex.assortment.selection.resid)
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted
##  vs. simulated
## 
## data:  simulationOutput
## ratioObsSim = 1.0824, p-value = 0.624
## alternative hypothesis: two.sided
testZeroInflation(lays.sex.assortment.selection.resid)
## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sided
#the model is not overdispersed nor zero-inflated.

Eggs

eggs.sex.assortment.selection<-glmer(TotalEggs~StandardizedSexAssortment+(1|Condo), data=PopFBeetleset, family="poisson")
summary(eggs.sex.assortment.selection)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: TotalEggs ~ StandardizedSexAssortment + (1 | Condo)
##    Data: PopFBeetleset
## 
##      AIC      BIC   logLik deviance df.resid 
##    262.4    266.0   -128.2    256.4       21 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.88731 -1.26284 -0.08057  1.26477  2.16905 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Condo  (Intercept) 0.009455 0.09724 
## Number of obs: 24, groups:  Condo, 12
## 
## Fixed effects:
##                           Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                5.64102    0.03061 184.316  < 2e-16 ***
## StandardizedSexAssortment  0.05817    0.01674   3.475 0.000511 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## StndrdzdSxA -0.018
Anova(eggs.sex.assortment.selection, type=3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: TotalEggs
##                               Chisq Df Pr(>Chisq)    
## (Intercept)               33972.499  1  < 2.2e-16 ***
## StandardizedSexAssortment    12.074  1  0.0005113 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(eggs.sex.assortment.selection))

check.assumptions(eggs.sex.assortment.selection)

#test for zero inflation and overdispersion with DHARMa
eggs.sex.assortment.selection.resid<-simulateResiduals(eggs.sex.assortment.selection, n=250)
plot(eggs.sex.assortment.selection.resid)

testDispersion(eggs.sex.assortment.selection.resid)
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted
##  vs. simulated
## 
## data:  simulationOutput
## ratioObsSim = 1.2394, p-value = 0.176
## alternative hypothesis: two.sided
testZeroInflation(eggs.sex.assortment.selection.resid)
## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sided
#the model is not overdispersed nor zero-inflated.

Is there selection on male body size?

#male.elytra.mod<-glmer(UniqueGuards~StandardizedElytra+StandardizedScansSeen+Period+(1|Condo/ID), data=surveydataMBeetleset, family="poisson")
#summary(male.elytra.mod)
#Anova(male.elytra.mod, type=3)
#plot(allEffects(male.elytra.mod))
#check.assumptions(male.elytra.mod)
#the model isSingular

#test for zero inflation and overdispersion with DHARMa
#male.elytra.mod.resid<-simulateResiduals(male.elytra.mod, n=250)
#plot(male.elytra.mod.resid)
#testDispersion(male.elytra.mod.resid)
#testZeroInflation(male.elytra.mod.resid)
#the model is not overdispersed but it is zero inflated

#specify 0 inflation with glmmTMB
male.elytra.0inflation<-glmmTMB(UniqueGuards~StandardizedElytra+StandardizedScansSeen+Period+(1|Condo/ID), data=surveydataMBeetleset, family="poisson", ziformula=~1)
summary(male.elytra.0inflation)
##  Family: poisson  ( log )
## Formula:          
## UniqueGuards ~ StandardizedElytra + StandardizedScansSeen + Period +  
##     (1 | Condo/ID)
## Zero inflation:                ~1
## Data: surveydataMBeetleset
## 
##      AIC      BIC   logLik deviance df.resid 
##   1608.6   1636.6   -797.3   1594.6      393 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance  Std.Dev. 
##  ID:Condo (Intercept) 6.690e-02 2.586e-01
##  Condo    (Intercept) 8.188e-10 2.861e-05
## Number of obs: 400, groups:  ID:Condo, 200; Condo, 12
## 
## Conditional model:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            1.03335    0.04799  21.531  < 2e-16 ***
## StandardizedElytra     0.20419    0.04348   4.696 2.65e-06 ***
## StandardizedScansSeen  0.25706    0.04247   6.053 1.42e-09 ***
## Period1               -0.03362    0.03168  -1.061    0.289    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.1420     0.2483  -8.626   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(male.elytra.0inflation, type=3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: UniqueGuards
##                          Chisq Df Pr(>Chisq)    
## (Intercept)           463.5931  1  < 2.2e-16 ***
## StandardizedElytra     22.0570  1  2.647e-06 ***
## StandardizedScansSeen  36.6381  1  1.422e-09 ***
## Period                  1.1257  1     0.2887    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(male.elytra.0inflation))

check.assumptions(male.elytra.0inflation)

#test for zero inflation and overdispersion with DHARMa
male.elytra.0inflation.resid<-simulateResiduals(male.elytra.0inflation, n=250)
## It seems you are diagnosing a glmmTBM model. There are still a few minor limitations associatd with this package. The most important is that glmmTMB doesn't implement an option to create unconditional predictions from the model, which means that predicted values (in res ~ pred) plots include the random effects. With strong random effects, this can sometimes create diagonal patterns from bottom left to top right in the res ~ pred plot. All other tests and plots should work as desired. Please see https://github.com/florianhartig/DHARMa/issues/16 for further details.
plot(male.elytra.0inflation.resid)

testDispersion(male.elytra.0inflation.resid)
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted
##  vs. simulated
## 
## data:  simulationOutput
## ratioObsSim = 0.84407, p-value < 2.2e-16
## alternative hypothesis: two.sided
testZeroInflation(male.elytra.0inflation.resid)
## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.0305, p-value = 0.768
## alternative hypothesis: two.sided
#the model is now underdispersed but it is not zero inflated

How does treatment influence selection on male body size?

#male.elytra.treatment.mod<-glmer(UniqueGuards~StandardizedElytra*Treatment+StandardizedScansSeen+Period+(1|Condo/ID), data=surveydataMBeetleset, family="poisson")
#summary(male.elytra.treatment.mod)
#Anova(male.elytra.treatment.mod, type=3)
#plot(allEffects(male.elytra.treatment.mod))
#check.assumptions(male.elytra.treatment.mod)
#model failed to converge. Model is nearly unidentifiable: very large eigenvalue - Rescale variables?

#test for zero inflation and overdispersion with DHARMa
#male.elytra.treatment.mod.resid<-simulateResiduals(male.elytra.treatment.mod, n=250)
#plot(male.elytra.treatment.mod.resid)
#testDispersion(male.elytra.treatment.mod.resid)
#testZeroInflation(male.elytra.treatment.mod.resid)
#the model is not overdispersed but it is zero inflated

#specify 0 inflation with glmmTMB
male.elytra.treatment.0inflation<-glmmTMB(UniqueGuards~StandardizedElytra*Treatment+StandardizedScansSeen+Period+(1|Condo/ID), data=surveydataMBeetleset, family="poisson", ziformula=~1)
summary(male.elytra.treatment.0inflation)
##  Family: poisson  ( log )
## Formula:          
## UniqueGuards ~ StandardizedElytra * Treatment + StandardizedScansSeen +  
##     Period + (1 | Condo/ID)
## Zero inflation:                ~1
## Data: surveydataMBeetleset
## 
##      AIC      BIC   logLik deviance df.resid 
##   1609.2   1645.1   -795.6   1591.2      391 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance  Std.Dev. 
##  ID:Condo (Intercept) 6.623e-02 2.573e-01
##  Condo    (Intercept) 7.030e-10 2.651e-05
## Number of obs: 400, groups:  ID:Condo, 200; Condo, 12
## 
## Conditional model:
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    1.02704    0.04820  21.306  < 2e-16 ***
## StandardizedElytra             0.21130    0.04341   4.867 1.13e-06 ***
## Treatment1                    -0.03155    0.03315  -0.952   0.3412    
## StandardizedScansSeen          0.25586    0.04235   6.042 1.52e-09 ***
## Period1                       -0.03152    0.03167  -0.995   0.3197    
## StandardizedElytra:Treatment1  0.06307    0.03507   1.799   0.0721 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.1751     0.2555  -8.512   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(male.elytra.treatment.0inflation, type=3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: UniqueGuards
##                                 Chisq Df Pr(>Chisq)    
## (Intercept)                  453.9650  1  < 2.2e-16 ***
## StandardizedElytra            23.6897  1  1.132e-06 ***
## Treatment                      0.9058  1     0.3412    
## StandardizedScansSeen         36.5046  1  1.523e-09 ***
## Period                         0.9903  1     0.3197    
## StandardizedElytra:Treatment   3.2346  1     0.0721 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(male.elytra.treatment.0inflation))

check.assumptions(male.elytra.treatment.0inflation)

#test for zero inflation and overdispersion with DHARMa
male.elytra.treatment.0inflation.resid<-simulateResiduals(male.elytra.treatment.0inflation, n=250)
## It seems you are diagnosing a glmmTBM model. There are still a few minor limitations associatd with this package. The most important is that glmmTMB doesn't implement an option to create unconditional predictions from the model, which means that predicted values (in res ~ pred) plots include the random effects. With strong random effects, this can sometimes create diagonal patterns from bottom left to top right in the res ~ pred plot. All other tests and plots should work as desired. Please see https://github.com/florianhartig/DHARMa/issues/16 for further details.
plot(male.elytra.treatment.0inflation.resid)

testDispersion(male.elytra.treatment.0inflation.resid)
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted
##  vs. simulated
## 
## data:  simulationOutput
## ratioObsSim = 0.83948, p-value < 2.2e-16
## alternative hypothesis: two.sided
testZeroInflation(male.elytra.treatment.0inflation.resid)
## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.0526, p-value = 0.624
## alternative hypothesis: two.sided
#the model is now underdispersed but it is not zero inflated

How does sex assortment influence selection on male body size?

#male.elytra.sexassort.mod<-glmer(UniqueGuards~StandardizedElytra*StandardizedSexAssortment+StandardizedScansSeen+Period+(1|Condo/ID), data=surveydataMBeetleset, family="poisson")
#summary(male.elytra.sexassort.mod)
#Anova(male.elytra.sexassort.mod, type=3)
#plot(allEffects(male.elytra.sexassort.mod))
#check.assumptions(male.elytra.sexassort.mod)
#model failed to converge. Model is nearly unidentifiable: very large eigenvalue - Rescale variables?

#test for zero inflation and overdispersion with DHARMa
#male.elytra.sexassort.mod.resid<-simulateResiduals(male.elytra.sexassort.mod, n=250)
#plot(male.elytra.sexassort.mod.resid)
#testDispersion(male.elytra.sexassort.mod.resid)
#testZeroInflation(male.elytra.sexassort.mod.resid)
#the model is not overdispersed but it is zero inflated

#specify 0 inflation with glmmTMB
male.elytra.sexassort.0inflation<-glmmTMB(UniqueGuards~StandardizedElytra*StandardizedSexAssortment+StandardizedScansSeen+Period+(1|Condo/ID), data=surveydataMBeetleset, family="poisson", ziformula=~1)
summary(male.elytra.sexassort.0inflation)
##  Family: poisson  ( log )
## Formula:          
## UniqueGuards ~ StandardizedElytra * StandardizedSexAssortment +  
##     StandardizedScansSeen + Period + (1 | Condo/ID)
## Zero inflation:                ~1
## Data: surveydataMBeetleset
## 
##      AIC      BIC   logLik deviance df.resid 
##   1609.2   1645.2   -795.6   1591.2      391 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance  Std.Dev. 
##  ID:Condo (Intercept) 6.052e-02 2.460e-01
##  Condo    (Intercept) 9.684e-10 3.112e-05
## Number of obs: 400, groups:  ID:Condo, 200; Condo, 12
## 
## Conditional model:
##                                              Estimate Std. Error z value
## (Intercept)                                   1.03140    0.04747  21.729
## StandardizedElytra                            0.21356    0.04307   4.958
## StandardizedSexAssortment                    -0.04451    0.03600  -1.236
## StandardizedScansSeen                         0.25358    0.04206   6.030
## Period1                                      -0.03279    0.03164  -1.036
## StandardizedElytra:StandardizedSexAssortment  0.06607    0.04047   1.633
##                                              Pr(>|z|)    
## (Intercept)                                   < 2e-16 ***
## StandardizedElytra                           7.13e-07 ***
## StandardizedSexAssortment                       0.216    
## StandardizedScansSeen                        1.64e-09 ***
## Period1                                         0.300    
## StandardizedElytra:StandardizedSexAssortment    0.103    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.1534     0.2488  -8.654   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(male.elytra.sexassort.0inflation, type=3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: UniqueGuards
##                                                 Chisq Df Pr(>Chisq)    
## (Intercept)                                  472.1383  1  < 2.2e-16 ***
## StandardizedElytra                            24.5805  1  7.127e-07 ***
## StandardizedSexAssortment                      1.5282  1     0.2164    
## StandardizedScansSeen                         36.3565  1  1.643e-09 ***
## Period                                         1.0741  1     0.3000    
## StandardizedElytra:StandardizedSexAssortment   2.6660  1     0.1025    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(male.elytra.sexassort.0inflation))

check.assumptions(male.elytra.sexassort.0inflation)

#test for zero inflation and overdispersion with DHARMa
male.elytra.sexassort.0inflation.resid<-simulateResiduals(male.elytra.sexassort.0inflation, n=250)
## It seems you are diagnosing a glmmTBM model. There are still a few minor limitations associatd with this package. The most important is that glmmTMB doesn't implement an option to create unconditional predictions from the model, which means that predicted values (in res ~ pred) plots include the random effects. With strong random effects, this can sometimes create diagonal patterns from bottom left to top right in the res ~ pred plot. All other tests and plots should work as desired. Please see https://github.com/florianhartig/DHARMa/issues/16 for further details.
plot(male.elytra.sexassort.0inflation.resid)

testDispersion(male.elytra.sexassort.0inflation.resid)
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted
##  vs. simulated
## 
## data:  simulationOutput
## ratioObsSim = 0.85293, p-value = 0.008
## alternative hypothesis: two.sided
testZeroInflation(male.elytra.sexassort.0inflation.resid)
## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.0412, p-value = 0.664
## alternative hypothesis: two.sided
#the model is now underdispersed but it is not zero inflated

Period 1 Analyses

Subset by Period

surveydataMBeetleset %>%
  filter(Period == "Period1") -> surveydataMBeetleset1

Is there selection on male body size?

#male.elytra.mod.1<-glmer(UniqueGuards~StandardizedElytra+StandardizedScansSeen+(1|Condo), data=surveydataMBeetleset1, family="poisson")
#summary(male.elytra.mod.1)
#Anova(male.elytra.mod.1, type=3)
#plot(allEffects(male.elytra.mod.1))
#check.assumptions(male.elytra.mod.1)
#the model isSingular

#test for zero inflation and overdispersion with DHARMa
#male.elytra.mod.1.resid<-simulateResiduals(male.elytra.mod.1, n=250)
#plot(male.elytra.mod.1.resid)
#testDispersion(male.elytra.mod.1.resid)
#testZeroInflation(male.elytra.mod.1.resid)
#the model is both overdispersed and zero inflated

#specify 0 inflation with glmmTMB
male.elytra.0inflation.1<-glmmTMB(UniqueGuards~StandardizedElytra+StandardizedScansSeen+(1|Condo), data=surveydataMBeetleset1, family="poisson", ziformula=~1)
summary(male.elytra.0inflation.1)
##  Family: poisson  ( log )
## Formula:          
## UniqueGuards ~ StandardizedElytra + StandardizedScansSeen + (1 |  
##     Condo)
## Zero inflation:                ~1
## Data: surveydataMBeetleset1
## 
##      AIC      BIC   logLik deviance df.resid 
##    782.4    798.9   -386.2    772.4      195 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Condo  (Intercept) 5.2e-10  2.28e-05
## Number of obs: 200, groups:  Condo, 12
## 
## Conditional model:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            1.01279    0.05290  19.146  < 2e-16 ***
## StandardizedElytra     0.20167    0.05327   3.786 0.000153 ***
## StandardizedScansSeen  0.32605    0.05568   5.856 4.73e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.2288     0.3439  -6.481 9.14e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(male.elytra.0inflation.1, type=3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: UniqueGuards
##                         Chisq Df Pr(>Chisq)    
## (Intercept)           366.558  1  < 2.2e-16 ***
## StandardizedElytra     14.332  1  0.0001532 ***
## StandardizedScansSeen  34.296  1  4.733e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(male.elytra.0inflation.1))

check.assumptions(male.elytra.0inflation.1)

#test for zero inflation and overdispersion with DHARMa
male.elytra.0inflation.1.resid<-simulateResiduals(male.elytra.0inflation.1, n=250)
## It seems you are diagnosing a glmmTBM model. There are still a few minor limitations associatd with this package. The most important is that glmmTMB doesn't implement an option to create unconditional predictions from the model, which means that predicted values (in res ~ pred) plots include the random effects. With strong random effects, this can sometimes create diagonal patterns from bottom left to top right in the res ~ pred plot. All other tests and plots should work as desired. Please see https://github.com/florianhartig/DHARMa/issues/16 for further details.
plot(male.elytra.0inflation.1.resid)

testDispersion(male.elytra.0inflation.1.resid)
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted
##  vs. simulated
## 
## data:  simulationOutput
## ratioObsSim = 1.022, p-value = 0.712
## alternative hypothesis: two.sided
testZeroInflation(male.elytra.0inflation.1.resid)
## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.0775, p-value = 0.608
## alternative hypothesis: two.sided
#the model is now neither overdispersed nor zero inflated

How does treatment influence selection on male body size?

#male.elytra.treatment.mod.1<-glmer(UniqueGuards~StandardizedElytra*Treatment+StandardizedScansSeen+(1|Condo), data=surveydataMBeetleset1, family="poisson")
#summary(male.elytra.treatment.mod.1)
#Anova(male.elytra.treatment.mod.1, type=3)
#plot(allEffects(male.elytra.treatment.mod.1))
#check.assumptions(male.elytra.treatment.mod.1)
#the model ?isSingular

#test for zero inflation and overdispersion with DHARMa
#male.elytra.treatment.mod.1.resid<-simulateResiduals(male.elytra.treatment.mod.1, n=250)
#plot(male.elytra.treatment.mod.1.resid)
#testDispersion(male.elytra.treatment.mod.1.resid)
#testZeroInflation(male.elytra.treatment.mod.1.resid)
#the model is both overdispersed and zero inflated

#specify 0 inflation with glmmTMB
male.elytra.treatment.0inflation.1<-glmmTMB(UniqueGuards~StandardizedElytra*Treatment+StandardizedScansSeen+(1|Condo), data=surveydataMBeetleset1, family="poisson", ziformula=~1)
summary(male.elytra.treatment.0inflation.1)
##  Family: poisson  ( log )
## Formula:          
## UniqueGuards ~ StandardizedElytra * Treatment + StandardizedScansSeen +  
##     (1 | Condo)
## Zero inflation:                ~1
## Data: surveydataMBeetleset1
## 
##      AIC      BIC   logLik deviance df.resid 
##    785.4    808.5   -385.7    771.4      193 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance  Std.Dev. 
##  Condo  (Intercept) 5.787e-10 2.406e-05
## Number of obs: 200, groups:  Condo, 12
## 
## Conditional model:
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    1.01375    0.05298  19.135  < 2e-16 ***
## StandardizedElytra             0.20274    0.05316   3.814 0.000137 ***
## Treatment1                     0.01375    0.04672   0.294 0.768597    
## StandardizedScansSeen          0.32250    0.05585   5.775  7.7e-09 ***
## StandardizedElytra:Treatment1  0.03991    0.04887   0.817 0.414151    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.2243     0.3435  -6.476 9.41e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(male.elytra.treatment.0inflation.1, type=3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: UniqueGuards
##                                 Chisq Df Pr(>Chisq)    
## (Intercept)                  366.1660  1  < 2.2e-16 ***
## StandardizedElytra            14.5447  1  0.0001369 ***
## Treatment                      0.0866  1  0.7685967    
## StandardizedScansSeen         33.3497  1  7.699e-09 ***
## StandardizedElytra:Treatment   0.6669  1  0.4141507    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(male.elytra.treatment.0inflation.1))

check.assumptions(male.elytra.treatment.0inflation.1)

#test for zero inflation and overdispersion with DHARMa
male.elytra.treatment.0inflation.1.resid<-simulateResiduals(male.elytra.treatment.0inflation.1, n=250)
## It seems you are diagnosing a glmmTBM model. There are still a few minor limitations associatd with this package. The most important is that glmmTMB doesn't implement an option to create unconditional predictions from the model, which means that predicted values (in res ~ pred) plots include the random effects. With strong random effects, this can sometimes create diagonal patterns from bottom left to top right in the res ~ pred plot. All other tests and plots should work as desired. Please see https://github.com/florianhartig/DHARMa/issues/16 for further details.
plot(male.elytra.treatment.0inflation.1.resid)

testDispersion(male.elytra.treatment.0inflation.1.resid)
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted
##  vs. simulated
## 
## data:  simulationOutput
## ratioObsSim = 1.0228, p-value = 0.664
## alternative hypothesis: two.sided
testZeroInflation(male.elytra.treatment.0inflation.1.resid)
## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.08, p-value = 0.632
## alternative hypothesis: two.sided
#the model is now neither overdispersed nor zero inflated

How does sex assortment influence selection on male body size?

#male.elytra.sexassort.mod.1<-glmer(UniqueGuards~StandardizedElytra*StandardizedSexAssortment+StandardizedScansSeen+(1|Condo), data=surveydataMBeetleset1, family="poisson")
#summary(male.elytra.sexassort.mod.1)
#Anova(male.elytra.sexassort.mod.1, type=3)
#plot(allEffects(male.elytra.sexassort.mod.1))
#check.assumptions(male.elytra.sexassort.mod.1)
#the model ?isSingular

#test for zero inflation and overdispersion with DHARMa
#male.elytra.sexassort.mod.1.resid<-simulateResiduals(male.elytra.sexassort.mod.1, n=250)
#plot(male.elytra.sexassort.mod.1.resid)
#testDispersion(male.elytra.sexassort.mod.1.resid)
#testZeroInflation(male.elytra.sexassort.mod.1.resid)
#the model is both overdispersed and zero inflated

#specify 0 inflation with glmmTMB
male.elytra.sexassort.0inflation.1<-glmmTMB(UniqueGuards~StandardizedElytra*StandardizedSexAssortment+StandardizedScansSeen+(1|Condo), data=surveydataMBeetleset1, family="poisson", ziformula=~1)
summary(male.elytra.sexassort.0inflation.1)
##  Family: poisson  ( log )
## Formula:          
## UniqueGuards ~ StandardizedElytra * StandardizedSexAssortment +  
##     StandardizedScansSeen + (1 | Condo)
## Zero inflation:                ~1
## Data: surveydataMBeetleset1
## 
##      AIC      BIC   logLik deviance df.resid 
##    782.7    805.8   -384.3    768.7      193 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance  Std.Dev. 
##  Condo  (Intercept) 5.528e-10 2.351e-05
## Number of obs: 200, groups:  Condo, 12
## 
## Conditional model:
##                                              Estimate Std. Error z value
## (Intercept)                                   1.00201    0.05313  18.859
## StandardizedElytra                            0.21326    0.05339   3.994
## StandardizedSexAssortment                    -0.04659    0.04789  -0.973
## StandardizedScansSeen                         0.32412    0.05571   5.818
## StandardizedElytra:StandardizedSexAssortment  0.09706    0.05122   1.895
##                                              Pr(>|z|)    
## (Intercept)                                   < 2e-16 ***
## StandardizedElytra                           6.49e-05 ***
## StandardizedSexAssortment                      0.3307    
## StandardizedScansSeen                        5.97e-09 ***
## StandardizedElytra:StandardizedSexAssortment   0.0581 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.2822     0.3554  -6.422 1.34e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(male.elytra.sexassort.0inflation.1, type=3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: UniqueGuards
##                                                 Chisq Df Pr(>Chisq)    
## (Intercept)                                  355.6593  1  < 2.2e-16 ***
## StandardizedElytra                            15.9542  1  6.490e-05 ***
## StandardizedSexAssortment                      0.9462  1    0.33070    
## StandardizedScansSeen                         33.8437  1  5.972e-09 ***
## StandardizedElytra:StandardizedSexAssortment   3.5915  1    0.05807 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(male.elytra.sexassort.0inflation.1))

check.assumptions(male.elytra.sexassort.0inflation.1)

#test for zero inflation and overdispersion with DHARMa
male.elytra.sexassort.0inflation.1.resid<-simulateResiduals(male.elytra.sexassort.0inflation.1, n=250)
## It seems you are diagnosing a glmmTBM model. There are still a few minor limitations associatd with this package. The most important is that glmmTMB doesn't implement an option to create unconditional predictions from the model, which means that predicted values (in res ~ pred) plots include the random effects. With strong random effects, this can sometimes create diagonal patterns from bottom left to top right in the res ~ pred plot. All other tests and plots should work as desired. Please see https://github.com/florianhartig/DHARMa/issues/16 for further details.
plot(male.elytra.sexassort.0inflation.1.resid)

testDispersion(male.elytra.sexassort.0inflation.1.resid)
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted
##  vs. simulated
## 
## data:  simulationOutput
## ratioObsSim = 1.0174, p-value = 0.744
## alternative hypothesis: two.sided
testZeroInflation(male.elytra.sexassort.0inflation.1.resid)
## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.0848, p-value = 0.608
## alternative hypothesis: two.sided
#the model is now neither overdispersed nor zero inflated

Visualize the Interaction

#set ggplot theme
theme_set(theme_classic(base_size=18))

#sex assortment by elytra
geff.out1 <- ggeffect(male.elytra.sexassort.0inflation.1, c("StandardizedElytra","StandardizedSexAssortment[quart]"))
plot(geff.out1, ci=TRUE, ci.style=c("ribbon"), rawdata=TRUE, alpha=0.5, use.theme=FALSE, dot.alpha=1, jitter=0.01, show.title=FALSE, show.x.title=TRUE, show.y.title=TRUE, dot.size=0.5) + 
  labs(x ="Elytra", y="Mating Success")

Period 2 Analyses

Subset by Period

surveydataMBeetleset %>%
  filter(Period == "Period2") -> surveydataMBeetleset2

Is there selection on male body size?

#male.elytra.mod.2<-glmer(UniqueGuards~StandardizedElytra+StandardizedScansSeen+(1|Condo), data=surveydataMBeetleset2, family="poisson")
#summary(male.elytra.mod.2)
#Anova(male.elytra.mod.2, type=3)
#plot(allEffects(male.elytra.mod.2))
#check.assumptions(male.elytra.mod.2)
#the model isSingular

#test for zero inflation and overdispersion with DHARMa
#male.elytra.mod.2.resid<-simulateResiduals(male.elytra.mod.2, n=250)
#plot(male.elytra.mod.2.resid)
#testDispersion(male.elytra.mod.2.resid)
#testZeroInflation(male.elytra.mod.2.resid)
#the model is both overdispersed and zero inflated

#specify 0 inflation with glmmTMB
male.elytra.0inflation.2<-glmmTMB(UniqueGuards~StandardizedElytra+StandardizedScansSeen+(1|Condo), data=surveydataMBeetleset2, family="poisson", ziformula=~1)
summary(male.elytra.0inflation.2)
##  Family: poisson  ( log )
## Formula:          
## UniqueGuards ~ StandardizedElytra + StandardizedScansSeen + (1 |  
##     Condo)
## Zero inflation:                ~1
## Data: surveydataMBeetleset2
## 
##      AIC      BIC   logLik deviance df.resid 
##    835.1    851.6   -412.5    825.1      195 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance  Std.Dev. 
##  Condo  (Intercept) 1.712e-09 4.138e-05
## Number of obs: 200, groups:  Condo, 12
## 
## Conditional model:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            1.15592    0.04848  23.843  < 2e-16 ***
## StandardizedElytra     0.17725    0.05070   3.496 0.000472 ***
## StandardizedScansSeen  0.18481    0.04998   3.697 0.000218 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.7853     0.2457  -7.266  3.7e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(male.elytra.0inflation.2, type=3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: UniqueGuards
##                         Chisq Df Pr(>Chisq)    
## (Intercept)           568.484  1  < 2.2e-16 ***
## StandardizedElytra     12.223  1  0.0004721 ***
## StandardizedScansSeen  13.671  1  0.0002177 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(male.elytra.0inflation.2))

check.assumptions(male.elytra.0inflation.2)

#test for zero inflation and overdispersion with DHARMa
male.elytra.0inflation.2.resid<-simulateResiduals(male.elytra.0inflation.2, n=250)
## It seems you are diagnosing a glmmTBM model. There are still a few minor limitations associatd with this package. The most important is that glmmTMB doesn't implement an option to create unconditional predictions from the model, which means that predicted values (in res ~ pred) plots include the random effects. With strong random effects, this can sometimes create diagonal patterns from bottom left to top right in the res ~ pred plot. All other tests and plots should work as desired. Please see https://github.com/florianhartig/DHARMa/issues/16 for further details.
plot(male.elytra.0inflation.2.resid)

testDispersion(male.elytra.0inflation.2.resid)
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted
##  vs. simulated
## 
## data:  simulationOutput
## ratioObsSim = 1.0581, p-value = 0.28
## alternative hypothesis: two.sided
testZeroInflation(male.elytra.0inflation.2.resid)
## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.0402, p-value = 0.824
## alternative hypothesis: two.sided
#the model is now neither overdispersed nor zero inflated

How does treatment influence selection on male body size?

#male.elytra.treatment.mod.2<-glmer(UniqueGuards~StandardizedElytra*Treatment+StandardizedScansSeen+(1|Condo), data=surveydataMBeetleset2, family="poisson")
#summary(male.elytra.treatment.mod.2)
#Anova(male.elytra.treatment.mod.2, type=3)
#plot(allEffects(male.elytra.treatment.mod.2))
#check.assumptions(male.elytra.treatment.mod.2)
#the model ?isSingular

#test for zero inflation and overdispersion with DHARMa
#male.elytra.treatment.mod.2.resid<-simulateResiduals(male.elytra.treatment.mod.2, n=250)
#plot(male.elytra.treatment.mod.2.resid)
#testDispersion(male.elytra.treatment.mod.2.resid)
#testZeroInflation(male.elytra.treatment.mod.2.resid)
#the model is both overdispersed and zero inflated

#specify 0 inflation with glmmTMB
male.elytra.treatment.0inflation.2<-glmmTMB(UniqueGuards~StandardizedElytra*Treatment+StandardizedScansSeen+(1|Condo), data=surveydataMBeetleset2, family="poisson", ziformula=~1)
summary(male.elytra.treatment.0inflation.2)
##  Family: poisson  ( log )
## Formula:          
## UniqueGuards ~ StandardizedElytra * Treatment + StandardizedScansSeen +  
##     (1 | Condo)
## Zero inflation:                ~1
## Data: surveydataMBeetleset2
## 
##      AIC      BIC   logLik deviance df.resid 
##    834.8    857.9   -410.4    820.8      193 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance  Std.Dev. 
##  Condo  (Intercept) 5.214e-10 2.283e-05
## Number of obs: 200, groups:  Condo, 12
## 
## Conditional model:
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    1.14136    0.04916  23.217  < 2e-16 ***
## StandardizedElytra             0.19582    0.05110   3.832 0.000127 ***
## Treatment1                    -0.07888    0.04666  -1.690 0.090946 .  
## StandardizedScansSeen          0.18192    0.05008   3.633 0.000280 ***
## StandardizedElytra:Treatment1  0.08038    0.04881   1.647 0.099622 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.8329     0.2544  -7.206 5.78e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(male.elytra.treatment.0inflation.2, type=3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: UniqueGuards
##                                 Chisq Df Pr(>Chisq)    
## (Intercept)                  539.0310  1  < 2.2e-16 ***
## StandardizedElytra            14.6877  1  0.0001269 ***
## Treatment                      2.8575  1  0.0909462 .  
## StandardizedScansSeen         13.1985  1  0.0002802 ***
## StandardizedElytra:Treatment   2.7116  1  0.0996217 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(male.elytra.treatment.0inflation.2))

check.assumptions(male.elytra.treatment.0inflation.2)

#test for zero inflation and overdispersion with DHARMa
male.elytra.treatment.0inflation.2.resid<-simulateResiduals(male.elytra.treatment.0inflation.2, n=250)
## It seems you are diagnosing a glmmTBM model. There are still a few minor limitations associatd with this package. The most important is that glmmTMB doesn't implement an option to create unconditional predictions from the model, which means that predicted values (in res ~ pred) plots include the random effects. With strong random effects, this can sometimes create diagonal patterns from bottom left to top right in the res ~ pred plot. All other tests and plots should work as desired. Please see https://github.com/florianhartig/DHARMa/issues/16 for further details.
plot(male.elytra.treatment.0inflation.2.resid)

testDispersion(male.elytra.treatment.0inflation.2.resid)
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted
##  vs. simulated
## 
## data:  simulationOutput
## ratioObsSim = 1.0548, p-value = 0.272
## alternative hypothesis: two.sided
testZeroInflation(male.elytra.treatment.0inflation.2.resid)
## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.0433, p-value = 0.792
## alternative hypothesis: two.sided
#the model is now neither overdispersed nor zero inflated

How does sex assortment influence selection on male body size?

#male.elytra.sexassort.mod.2<-glmer(UniqueGuards~StandardizedElytra*StandardizedSexAssortment+StandardizedScansSeen+(1|Condo), data=surveydataMBeetleset2, family="poisson")
#summary(male.elytra.sexassort.mod.2)
#Anova(male.elytra.sexassort.mod.2, type=3)
#plot(allEffects(male.elytra.sexassort.mod.2))
#check.assumptions(male.elytra.sexassort.mod.2)
#the model ?isSingular

#test for zero inflation and overdispersion with DHARMa
#male.elytra.sexassort.mod.2.resid<-simulateResiduals(male.elytra.sexassort.mod.2, n=250)
#plot(male.elytra.sexassort.mod.2.resid)
#testDispersion(male.elytra.sexassort.mod.2.resid)
#testZeroInflation(male.elytra.sexassort.mod.2.resid)
#the model is both overdispersed and zero inflated

#specify 0 inflation with glmmTMB
male.elytra.sexassort.0inflation.2<-glmmTMB(UniqueGuards~StandardizedElytra*StandardizedSexAssortment+StandardizedScansSeen+(1|Condo), data=surveydataMBeetleset2, family="poisson", ziformula=~1)
summary(male.elytra.sexassort.0inflation.2)
##  Family: poisson  ( log )
## Formula:          
## UniqueGuards ~ StandardizedElytra * StandardizedSexAssortment +  
##     StandardizedScansSeen + (1 | Condo)
## Zero inflation:                ~1
## Data: surveydataMBeetleset2
## 
##      AIC      BIC   logLik deviance df.resid 
##    837.2    860.3   -411.6    823.2      193 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance  Std.Dev.
##  Condo  (Intercept) 7.955e-10 2.82e-05
## Number of obs: 200, groups:  Condo, 12
## 
## Conditional model:
##                                              Estimate Std. Error z value
## (Intercept)                                   1.14988    0.04872  23.604
## StandardizedElytra                            0.18798    0.05131   3.664
## StandardizedSexAssortment                    -0.05755    0.04758  -1.210
## StandardizedScansSeen                         0.18078    0.04982   3.628
## StandardizedElytra:StandardizedSexAssortment  0.05177    0.05185   0.998
##                                              Pr(>|z|)    
## (Intercept)                                   < 2e-16 ***
## StandardizedElytra                           0.000249 ***
## StandardizedSexAssortment                    0.226447    
## StandardizedScansSeen                        0.000285 ***
## StandardizedElytra:StandardizedSexAssortment 0.318040    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -1.802      0.248  -7.263  3.8e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(male.elytra.sexassort.0inflation.2, type=3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: UniqueGuards
##                                                Chisq Df Pr(>Chisq)    
## (Intercept)                                  557.132  1  < 2.2e-16 ***
## StandardizedElytra                            13.423  1  0.0002486 ***
## StandardizedSexAssortment                      1.463  1  0.2264474    
## StandardizedScansSeen                         13.165  1  0.0002853 ***
## StandardizedElytra:StandardizedSexAssortment   0.997  1  0.3180405    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(male.elytra.sexassort.0inflation.2))

check.assumptions(male.elytra.sexassort.0inflation.2)

#test for zero inflation and overdispersion with DHARMa
male.elytra.sexassort.0inflation.2.resid<-simulateResiduals(male.elytra.sexassort.0inflation.2, n=250)
## It seems you are diagnosing a glmmTBM model. There are still a few minor limitations associatd with this package. The most important is that glmmTMB doesn't implement an option to create unconditional predictions from the model, which means that predicted values (in res ~ pred) plots include the random effects. With strong random effects, this can sometimes create diagonal patterns from bottom left to top right in the res ~ pred plot. All other tests and plots should work as desired. Please see https://github.com/florianhartig/DHARMa/issues/16 for further details.
plot(male.elytra.sexassort.0inflation.2.resid)

testDispersion(male.elytra.sexassort.0inflation.2.resid)
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted
##  vs. simulated
## 
## data:  simulationOutput
## ratioObsSim = 1.0549, p-value = 0.28
## alternative hypothesis: two.sided
testZeroInflation(male.elytra.sexassort.0inflation.2.resid)
## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.0383, p-value = 0.872
## alternative hypothesis: two.sided
#the model is now neither overdispersed nor zero inflated