Contextual Analysis 2 2018
Robin Costello
20190928
Date rendered: Mon Sep 30 08:46:26 2019
To Do
-In contextual analysis: should I include treatment? number of scans seen? interactions? -There is no relationship between female elytra and number of lays. This is strange. -What should I do with the observation that less negative sex assortment is associated with more eggs laid (but not more observed eggs laid). Is there a paper here?
Objectives
The goal of this project is to investigate how intersex interactions influence fitness.
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 independence 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: Matrixlibrary(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 lme4library(glmmTMB)
library(DHARMa)## Registered S3 method overwritten by 'DHARMa':
## method from
## refit.glmmTMB glmmTMBlibrary(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 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 Surveys Dataset
Rename X_fkfieldID_surveys
surveys_2018 %>%
dplyr::rename(ID = X_fkfieldID_surveys) -> surveys_2018Divide 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 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))) -> Period2FemalesCollapse Periods 1 and 2
surveydataM<-bind_rows(Period1Males, Period2Males)
surveydataF<-bind_rows(Period1Females, Period2Females)
surveydata<-bind_rows(surveydataM, surveydataF)Remove UKs
surveydata %>%
#remove UKM
filter(ID != "UKM") %>%
filter(ID != "UKF") %>%
filter(ID != "UK") -> surveydataRemove 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() -> surveydataAdd 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) -> surveydataBeetlesetRemove Extra Columns
surveydataBeetleset %>%
dplyr::select(ID, Elytra, Survey_Sex, ScansSeen, Lays, UniqueGuards, Period, Condo) -> surveydataBeetlesetStandardize 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 density
PopSN %>%
group_by(Period) %>%
mutate(StandardizedDensity = ((density-(mean(density)))/sd(density))) %>%
ungroup() -> PopSN
#standardize elytra assortment
PopSN %>%
group_by(Period) %>%
mutate(StandardizedElytraAssortment = ((elytra.assortment-(mean(elytra.assortment)))/sd(elytra.assortment))) %>%
ungroup() -> PopSN
#standardize global cc
PopSN %>%
group_by(Period) %>%
mutate(StandardizedGlobalCC = ((global.cc-(mean(global.cc)))/sd(global.cc))) %>%
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() -> surveydataBeetlesetMerge 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, StandardizedDensity, StandardizedElytraAssortment, StandardizedGlobalCC) -> PopSNSimple
#merge sex assortment to beetleset
surveydataBeetleset <- merge(surveydataBeetleset, PopSNSimple, by=("Condo_Period"), all.x=TRUE, all.y=FALSE)
#droplevels
surveydataBeetleset %>%
droplevels() -> surveydataBeetlesetSubset by Sex
#subset data for males
surveydataMBeetleset<-subset(surveydataBeetleset, Survey_Sex=="M")
#subset data for females
surveydataFBeetleset<-subset(surveydataBeetleset, Survey_Sex=="F")Do intersex interactions influence individual fitness?
Guards
guards.intersex<-glmer(UniqueGuards~StandardizedElytra+StandardizedSexAssortment+StandardizedScansSeen+Period+(1|Condo/ID), data=surveydataMBeetleset, family="poisson")## boundary (singular) fit: see ?isSingularsummary(guards.intersex)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: UniqueGuards ~ StandardizedElytra + StandardizedSexAssortment +
## StandardizedScansSeen + Period + (1 | Condo/ID)
## Data: surveydataMBeetleset
##
## AIC BIC logLik deviance df.resid
## 1631.3 1659.2 -808.6 1617.3 393
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.00202 -0.82284 -0.07646 0.61014 2.85851
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID:Condo (Intercept) 1.722e-01 4.150e-01
## Condo (Intercept) 7.067e-09 8.407e-05
## Number of obs: 400, groups: ID:Condo, 200; Condo, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.870634 0.001279 680.773 <2e-16 ***
## StandardizedElytra 0.202471 0.001276 158.621 <2e-16 ***
## StandardizedSexAssortment -0.034671 0.035841 -0.967 0.333
## StandardizedScansSeen 0.284315 0.001276 222.765 <2e-16 ***
## Period1 -0.029685 0.001276 -23.269 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) StndrE StndSA StndSS
## StndrdzdEly 0.000
## StndrdzdSxA 0.000 0.001
## StndrdzdScS 0.000 0.000 0.001
## Period1 0.000 0.000 0.000 0.000
## convergence code: 0
## boundary (singular) fit: see ?isSingularAnova(guards.intersex, type=3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: UniqueGuards
## Chisq Df Pr(>Chisq)
## (Intercept) 4.6345e+05 1 <2e-16 ***
## StandardizedElytra 2.5160e+04 1 <2e-16 ***
## StandardizedSexAssortment 9.3580e-01 1 0.3334
## StandardizedScansSeen 4.9624e+04 1 <2e-16 ***
## Period 5.4146e+02 1 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(allEffects(guards.intersex))
check.assumptions(guards.intersex)
#the model isSingular
#test for zero inflation and overdispersion with DHARMa
guards.intersex.resid<-simulateResiduals(guards.intersex, n=250)
plot(guards.intersex.resid)
testDispersion(guards.intersex.resid)##
## DHARMa nonparametric dispersion test via sd of residuals fitted
## vs. simulated
##
## data: simulationOutput
## ratioObsSim = 0.97147, p-value = 0.704
## alternative hypothesis: two.sidedtestZeroInflation(guards.intersex.resid)##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.4558, p-value < 2.2e-16
## alternative hypothesis: two.sided#the model is not overdispersed but it is zero inflated
#specify 0 inflation with glmmTMB
guards.intersex.0inflation<-glmmTMB(UniqueGuards~StandardizedElytra+StandardizedSexAssortment+StandardizedScansSeen+Period+(1|Condo/ID), data=surveydataMBeetleset, family="poisson", ziformula=~1)
summary(guards.intersex.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.9 1641.8 -796.9 1593.9 392
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## ID:Condo (Intercept) 6.653e-02 0.2579423
## Condo (Intercept) 6.350e-10 0.0000252
## Number of obs: 400, groups: ID:Condo, 200; Condo, 12
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.03194 0.04800 21.499 < 2e-16 ***
## StandardizedElytra 0.20558 0.04340 4.737 2.17e-06 ***
## StandardizedSexAssortment -0.02990 0.03502 -0.854 0.393
## StandardizedScansSeen 0.25743 0.04239 6.073 1.25e-09 ***
## Period1 -0.03330 0.03168 -1.051 0.293
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.1491 0.2497 -8.607 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Anova(guards.intersex.0inflation, type=3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: UniqueGuards
## Chisq Df Pr(>Chisq)
## (Intercept) 462.1858 1 < 2.2e-16 ***
## StandardizedElytra 22.4360 1 2.173e-06 ***
## StandardizedSexAssortment 0.7287 1 0.3933
## StandardizedScansSeen 36.8850 1 1.253e-09 ***
## Period 1.1051 1 0.2932
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(allEffects(guards.intersex.0inflation))
check.assumptions(guards.intersex.0inflation)
#test for zero inflation and overdispersion with DHARMa
guards.intersex.0inflation.resid<-simulateResiduals(guards.intersex.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(guards.intersex.0inflation.resid)
testDispersion(guards.intersex.0inflation.resid)##
## DHARMa nonparametric dispersion test via sd of residuals fitted
## vs. simulated
##
## data: simulationOutput
## ratioObsSim = 0.84456, p-value < 2.2e-16
## alternative hypothesis: two.sidedtestZeroInflation(guards.intersex.0inflation.resid)##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0502, p-value = 0.624
## alternative hypothesis: two.sided#the model is now underdispersed and is not zero inflated
Lays
lays.intersex<-glmer(Lays~StandardizedElytra+StandardizedSexAssortment+StandardizedScansSeen+Period+(1|Condo/ID), data=surveydataFBeetleset, family="poisson")
summary(lays.intersex)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula:
## Lays ~ StandardizedElytra + StandardizedSexAssortment + StandardizedScansSeen +
## Period + (1 | Condo/ID)
## Data: surveydataFBeetleset
##
## AIC BIC logLik deviance df.resid
## 1339.5 1367.4 -662.7 1325.5 393
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6289 -0.7115 -0.1356 0.5965 3.1678
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID:Condo (Intercept) 0.06903 0.2627
## Condo (Intercept) 0.01160 0.1077
## Number of obs: 400, groups: ID:Condo, 200; Condo, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.51735 0.05528 9.358 < 2e-16 ***
## StandardizedElytra -0.03587 0.04332 -0.828 0.408
## StandardizedSexAssortment 0.02635 0.04544 0.580 0.562
## StandardizedScansSeen 0.24928 0.04322 5.767 8.06e-09 ***
## Period1 -0.03593 0.03711 -0.968 0.333
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) StndrE StndSA StndSS
## StndrdzdEly 0.024
## StndrdzdSxA -0.005 0.003
## StndrdzdScS -0.163 -0.064 -0.003
## Period1 0.023 -0.004 0.021 -0.006Anova(lays.intersex, type=3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: Lays
## Chisq Df Pr(>Chisq)
## (Intercept) 87.5729 1 < 2.2e-16 ***
## StandardizedElytra 0.6857 1 0.4076
## StandardizedSexAssortment 0.3364 1 0.5619
## StandardizedScansSeen 33.2614 1 8.056e-09 ***
## Period 0.9372 1 0.3330
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(allEffects(lays.intersex))
check.assumptions(lays.intersex)
#test for zero inflation and overdispersion with DHARMa
lays.intersex.resid<-simulateResiduals(lays.intersex, n=250)
plot(lays.intersex.resid)
testDispersion(lays.intersex.resid)##
## DHARMa nonparametric dispersion test via sd of residuals fitted
## vs. simulated
##
## data: simulationOutput
## ratioObsSim = 0.99039, p-value = 0.864
## alternative hypothesis: two.sidedtestZeroInflation(lays.intersex.resid)##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0024, p-value = 0.976
## alternative hypothesis: two.sided#the model is not overdispersed nor zero inflated
Does population connectedness influence individual fitness?
Guards
guards.density<-glmer(UniqueGuards~StandardizedElytra+StandardizedDensity+StandardizedScansSeen+Period+(1|Condo/ID), data=surveydataMBeetleset, family="poisson")## boundary (singular) fit: see ?isSingularsummary(guards.density)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula:
## UniqueGuards ~ StandardizedElytra + StandardizedDensity + StandardizedScansSeen +
## Period + (1 | Condo/ID)
## Data: surveydataMBeetleset
##
## AIC BIC logLik deviance df.resid
## 1632.2 1660.1 -809.1 1618.2 393
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.00723 -0.81226 -0.09099 0.63159 2.91631
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID:Condo (Intercept) 1.732e-01 0.4162319
## Condo (Intercept) 6.603e-10 0.0000257
## Number of obs: 400, groups: ID:Condo, 200; Condo, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.871058 0.001275 683.04 <2e-16 ***
## StandardizedElytra 0.201769 0.001273 158.52 <2e-16 ***
## StandardizedDensity -0.002506 0.001272 -1.97 0.0489 *
## StandardizedScansSeen 0.284231 0.001273 223.33 <2e-16 ***
## Period1 -0.030295 0.001272 -23.82 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) StndrE StndrD StndSS
## StndrdzdEly 0.000
## StndrdzdDns 0.000 0.000
## StndrdzdScS 0.000 0.000 0.000
## Period1 0.000 0.000 0.000 0.000
## convergence code: 0
## boundary (singular) fit: see ?isSingularAnova(guards.density, type=3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: UniqueGuards
## Chisq Df Pr(>Chisq)
## (Intercept) 4.6654e+05 1 <2e-16 ***
## StandardizedElytra 2.5128e+04 1 <2e-16 ***
## StandardizedDensity 3.8789e+00 1 0.0489 *
## StandardizedScansSeen 4.9876e+04 1 <2e-16 ***
## Period 5.6715e+02 1 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(allEffects(guards.density))
check.assumptions(guards.density)
#the model isSingular
#test for zero inflation and overdispersion with DHARMa
guards.density.resid<-simulateResiduals(guards.density, n=250)
plot(guards.density.resid)
testDispersion(guards.density.resid)##
## DHARMa nonparametric dispersion test via sd of residuals fitted
## vs. simulated
##
## data: simulationOutput
## ratioObsSim = 0.97429, p-value = 0.792
## alternative hypothesis: two.sidedtestZeroInflation(guards.density.resid)##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.4688, p-value < 2.2e-16
## alternative hypothesis: two.sided#the model is not overdispersed but it is zero inflated
#specify 0 inflation with glmmTMB
guards.density.0inflation<-glmmTMB(UniqueGuards~StandardizedElytra+StandardizedDensity+StandardizedScansSeen+Period+(1|Condo/ID), data=surveydataMBeetleset, family="poisson", ziformula=~1)
summary(guards.density.0inflation)## Family: poisson ( log )
## Formula:
## UniqueGuards ~ StandardizedElytra + StandardizedDensity + StandardizedScansSeen +
## Period + (1 | Condo/ID)
## Zero inflation: ~1
## Data: surveydataMBeetleset
##
## AIC BIC logLik deviance df.resid
## 1610.5 1642.4 -797.2 1594.5 392
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## ID:Condo (Intercept) 6.717e-02 0.2591734
## Condo (Intercept) 7.077e-10 0.0000266
## Number of obs: 400, groups: ID:Condo, 200; Condo, 12
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.03363 0.04796 21.552 < 2e-16 ***
## StandardizedElytra 0.20384 0.04350 4.686 2.79e-06 ***
## StandardizedDensity -0.01259 0.03441 -0.366 0.714
## StandardizedScansSeen 0.25727 0.04245 6.060 1.36e-09 ***
## Period1 -0.03390 0.03170 -1.070 0.285
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.1401 0.2476 -8.644 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Anova(guards.density.0inflation, type=3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: UniqueGuards
## Chisq Df Pr(>Chisq)
## (Intercept) 464.4726 1 < 2.2e-16 ***
## StandardizedElytra 21.9553 1 2.791e-06 ***
## StandardizedDensity 0.1339 1 0.7144
## StandardizedScansSeen 36.7278 1 1.358e-09 ***
## Period 1.1440 1 0.2848
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(allEffects(guards.density.0inflation))
check.assumptions(guards.density.0inflation)
#test for zero inflation and overdispersion with DHARMa
guards.density.0inflation.resid<-simulateResiduals(guards.density.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(guards.density.0inflation.resid)
testDispersion(guards.density.0inflation.resid)##
## DHARMa nonparametric dispersion test via sd of residuals fitted
## vs. simulated
##
## data: simulationOutput
## ratioObsSim = 0.84222, p-value < 2.2e-16
## alternative hypothesis: two.sidedtestZeroInflation(guards.density.0inflation.resid)##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0404, p-value = 0.68
## alternative hypothesis: two.sided#the model is now underdispersed and is not zero inflated
Lays
lays.density<-glmer(Lays~StandardizedElytra+StandardizedDensity+StandardizedScansSeen+Period+(1|Condo/ID), data=surveydataFBeetleset, family="poisson")
summary(lays.density)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula:
## Lays ~ StandardizedElytra + StandardizedDensity + StandardizedScansSeen +
## Period + (1 | Condo/ID)
## Data: surveydataFBeetleset
##
## AIC BIC logLik deviance df.resid
## 1337.3 1365.2 -661.6 1323.3 393
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.5394 -0.7076 -0.1473 0.5687 3.0131
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID:Condo (Intercept) 0.06871 0.2621
## Condo (Intercept) 0.01078 0.1038
## Number of obs: 400, groups: ID:Condo, 200; Condo, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.51647 0.05467 9.447 < 2e-16 ***
## StandardizedElytra -0.03543 0.04330 -0.818 0.413
## StandardizedDensity -0.06583 0.04097 -1.607 0.108
## StandardizedScansSeen 0.24899 0.04323 5.760 8.4e-09 ***
## Period1 -0.03289 0.03721 -0.884 0.377
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) StndrE StndrD StndSS
## StndrdzdEly 0.024
## StndrdzdDns 0.032 -0.008
## StndrdzdScS -0.164 -0.063 0.004
## Period1 0.020 -0.003 -0.057 -0.007Anova(lays.density, type=3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: Lays
## Chisq Df Pr(>Chisq)
## (Intercept) 89.2544 1 < 2.2e-16 ***
## StandardizedElytra 0.6695 1 0.4132
## StandardizedDensity 2.5814 1 0.1081
## StandardizedScansSeen 33.1794 1 8.404e-09 ***
## Period 0.7812 1 0.3768
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(allEffects(lays.density))
check.assumptions(lays.density)
#test for zero inflation and overdispersion with DHARMa
lays.density.resid<-simulateResiduals(lays.density, n=250)
plot(lays.density.resid)
testDispersion(lays.density.resid)##
## DHARMa nonparametric dispersion test via sd of residuals fitted
## vs. simulated
##
## data: simulationOutput
## ratioObsSim = 0.98852, p-value = 0.84
## alternative hypothesis: two.sidedtestZeroInflation(lays.density.resid)##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0.99975, p-value = 1
## alternative hypothesis: two.sided#the model is not overdispersed nor zero inflated
Does elytra assortment influence fitness?
Guards
guards.elytra.assortment<-glmer(UniqueGuards~StandardizedElytra+StandardizedElytraAssortment+StandardizedScansSeen+Period+(1|Condo/ID), data=surveydataMBeetleset, family="poisson")## boundary (singular) fit: see ?isSingularsummary(guards.elytra.assortment)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula:
## UniqueGuards ~ StandardizedElytra + StandardizedElytraAssortment +
## StandardizedScansSeen + Period + (1 | Condo/ID)
## Data: surveydataMBeetleset
##
## AIC BIC logLik deviance df.resid
## 1632 1660 -809 1618 393
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.02700 -0.81199 -0.07933 0.62668 2.89997
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID:Condo (Intercept) 0.1722 0.4149
## Condo (Intercept) 0.0000 0.0000
## Number of obs: 400, groups: ID:Condo, 200; Condo, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.87258 0.04646 18.781 < 2e-16 ***
## StandardizedElytra 0.20134 0.04809 4.187 2.83e-05 ***
## StandardizedElytraAssortment -0.01599 0.03563 -0.449 0.654
## StandardizedScansSeen 0.28401 0.04370 6.499 8.09e-11 ***
## Period1 -0.03046 0.02976 -1.024 0.306
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) StndrE StndEA StndSS
## StndrdzdEly -0.110
## StndrdzdElA 0.008 -0.001
## StndrdzdScS -0.141 -0.286 -0.004
## Period1 0.031 0.002 0.014 -0.044
## convergence code: 0
## boundary (singular) fit: see ?isSingularAnova(guards.elytra.assortment, type=3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: UniqueGuards
## Chisq Df Pr(>Chisq)
## (Intercept) 352.7330 1 < 2.2e-16 ***
## StandardizedElytra 17.5275 1 2.832e-05 ***
## StandardizedElytraAssortment 0.2014 1 0.6536
## StandardizedScansSeen 42.2369 1 8.086e-11 ***
## Period 1.0479 1 0.3060
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(allEffects(guards.elytra.assortment))
check.assumptions(guards.elytra.assortment)
#the model isSingular
#test for zero inflation and overdispersion with DHARMa
guards.elytra.assortment.resid<-simulateResiduals(guards.elytra.assortment, n=250)
plot(guards.elytra.assortment.resid)
testDispersion(guards.elytra.assortment.resid)##
## DHARMa nonparametric dispersion test via sd of residuals fitted
## vs. simulated
##
## data: simulationOutput
## ratioObsSim = 0.97421, p-value = 0.752
## alternative hypothesis: two.sidedtestZeroInflation(guards.elytra.assortment.resid)##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.467, p-value = 0.008
## alternative hypothesis: two.sided#the model is not overdispersed but it is zero inflated
#specify 0 inflation with glmmTMB
guards.elytra.assortment.0inflation<-glmmTMB(UniqueGuards~StandardizedElytra+StandardizedElytraAssortment+StandardizedScansSeen+Period+(1|Condo/ID), data=surveydataMBeetleset, family="poisson", ziformula=~1)
summary(guards.elytra.assortment.0inflation)## Family: poisson ( log )
## Formula:
## UniqueGuards ~ StandardizedElytra + StandardizedElytraAssortment +
## StandardizedScansSeen + Period + (1 | Condo/ID)
## Zero inflation: ~1
## Data: surveydataMBeetleset
##
## AIC BIC logLik deviance df.resid
## 1609.8 1641.7 -796.9 1593.8 392
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## ID:Condo (Intercept) 6.677e-02 0.2584002
## Condo (Intercept) 1.069e-09 0.0000327
## Number of obs: 400, groups: ID:Condo, 200; Condo, 12
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.03303 0.04784 21.591 < 2e-16 ***
## StandardizedElytra 0.20417 0.04348 4.695 2.66e-06 ***
## StandardizedElytraAssortment -0.03360 0.03619 -0.928 0.353
## StandardizedScansSeen 0.25862 0.04246 6.090 1.13e-09 ***
## Period1 -0.03480 0.03170 -1.098 0.272
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.1437 0.2477 -8.654 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Anova(guards.elytra.assortment.0inflation, type=3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: UniqueGuards
## Chisq Df Pr(>Chisq)
## (Intercept) 466.1836 1 < 2.2e-16 ***
## StandardizedElytra 22.0464 1 2.661e-06 ***
## StandardizedElytraAssortment 0.8620 1 0.3532
## StandardizedScansSeen 37.0939 1 1.126e-09 ***
## Period 1.2054 1 0.2722
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(allEffects(guards.elytra.assortment.0inflation))
check.assumptions(guards.elytra.assortment.0inflation)
#test for zero inflation and overdispersion with DHARMa
guards.elytra.assortment.0inflation.resid<-simulateResiduals(guards.elytra.assortment.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(guards.elytra.assortment.0inflation.resid)
testDispersion(guards.elytra.assortment.0inflation.resid)##
## DHARMa nonparametric dispersion test via sd of residuals fitted
## vs. simulated
##
## data: simulationOutput
## ratioObsSim = 0.84443, p-value < 2.2e-16
## alternative hypothesis: two.sidedtestZeroInflation(guards.elytra.assortment.0inflation.resid)##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0489, p-value = 0.624
## alternative hypothesis: two.sided#the model is now underdispersed and is not zero inflated
Lays
lays.elytra.assortment<-glmer(Lays~StandardizedElytra+StandardizedElytraAssortment+StandardizedScansSeen+Period+(1|Condo/ID), data=surveydataFBeetleset, family="poisson")
summary(lays.elytra.assortment)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula:
## Lays ~ StandardizedElytra + StandardizedElytraAssortment + StandardizedScansSeen +
## Period + (1 | Condo/ID)
## Data: surveydataFBeetleset
##
## AIC BIC logLik deviance df.resid
## 1339.4 1367.3 -662.7 1325.4 393
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6173 -0.7161 -0.1328 0.5805 3.1199
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID:Condo (Intercept) 0.06854 0.2618
## Condo (Intercept) 0.01265 0.1125
## Number of obs: 400, groups: ID:Condo, 200; Condo, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.51744 0.05603 9.234 < 2e-16 ***
## StandardizedElytra -0.03517 0.04328 -0.812 0.417
## StandardizedElytraAssortment -0.02750 0.04395 -0.626 0.531
## StandardizedScansSeen 0.24895 0.04320 5.763 8.25e-09 ***
## Period1 -0.03478 0.03719 -0.935 0.350
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) StndrE StndEA StndSS
## StndrdzdEly 0.023
## StndrdzdElA 0.004 -0.029
## StndrdzdScS -0.160 -0.062 0.015
## Period1 0.022 -0.002 -0.070 -0.008Anova(lays.elytra.assortment, type=3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: Lays
## Chisq Df Pr(>Chisq)
## (Intercept) 85.2739 1 < 2.2e-16 ***
## StandardizedElytra 0.6600 1 0.4166
## StandardizedElytraAssortment 0.3916 1 0.5315
## StandardizedScansSeen 33.2156 1 8.248e-09 ***
## Period 0.8746 1 0.3497
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(allEffects(lays.elytra.assortment))
check.assumptions(lays.elytra.assortment)
#test for zero inflation and overdispersion with DHARMa
lays.elytra.assortment.resid<-simulateResiduals(lays.elytra.assortment, n=250)
plot(lays.elytra.assortment.resid)
testDispersion(lays.elytra.assortment.resid)##
## DHARMa nonparametric dispersion test via sd of residuals fitted
## vs. simulated
##
## data: simulationOutput
## ratioObsSim = 0.98973, p-value = 0.864
## alternative hypothesis: two.sidedtestZeroInflation(lays.elytra.assortment.resid)##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0029, p-value = 1
## alternative hypothesis: two.sided#the model is not overdispersed nor zero inflated
Does global cc influence fitness?
Guards
guards.global.cc<-glmer(UniqueGuards~StandardizedElytra+StandardizedGlobalCC+StandardizedScansSeen+Period+(1|Condo/ID), data=surveydataMBeetleset, family="poisson")## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.0553114
## (tol = 0.001, component 1)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?summary(guards.global.cc)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula:
## UniqueGuards ~ StandardizedElytra + StandardizedGlobalCC + StandardizedScansSeen +
## Period + (1 | Condo/ID)
## Data: surveydataMBeetleset
##
## AIC BIC logLik deviance df.resid
## 1631.6 1659.6 -808.8 1617.6 393
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.01675 -0.80925 -0.06796 0.61823 2.89221
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID:Condo (Intercept) 1.721e-01 0.4148562
## Condo (Intercept) 1.448e-07 0.0003806
## Number of obs: 400, groups: ID:Condo, 200; Condo, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.871172 0.001285 677.83 <2e-16 ***
## StandardizedElytra 0.202027 0.001283 157.49 <2e-16 ***
## StandardizedGlobalCC -0.027383 0.001282 -21.36 <2e-16 ***
## StandardizedScansSeen 0.283582 0.001283 221.09 <2e-16 ***
## Period1 -0.029289 0.001282 -22.84 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) StndrE StnGCC StndSS
## StndrdzdEly 0.000
## StndrdzdGCC 0.000 0.000
## StndrdzdScS 0.000 0.000 0.000
## Period1 0.000 0.000 0.000 0.000
## convergence code: 0
## Model failed to converge with max|grad| = 0.0553114 (tol = 0.001, component 1)
## Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?Anova(guards.global.cc, type=3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: UniqueGuards
## Chisq Df Pr(>Chisq)
## (Intercept) 459449.14 1 < 2.2e-16 ***
## StandardizedElytra 24802.60 1 < 2.2e-16 ***
## StandardizedGlobalCC 456.02 1 < 2.2e-16 ***
## StandardizedScansSeen 48882.31 1 < 2.2e-16 ***
## Period 521.87 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(allEffects(guards.global.cc))
check.assumptions(guards.global.cc)
#model failed to converge
#test for zero inflation and overdispersion with DHARMa
guards.global.cc.resid<-simulateResiduals(guards.global.cc, n=250)
plot(guards.global.cc.resid)
testDispersion(guards.global.cc.resid)##
## DHARMa nonparametric dispersion test via sd of residuals fitted
## vs. simulated
##
## data: simulationOutput
## ratioObsSim = 0.97107, p-value = 0.72
## alternative hypothesis: two.sidedtestZeroInflation(guards.global.cc.resid)##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.455, p-value < 2.2e-16
## alternative hypothesis: two.sided#the model is not overdispersed but it is zero inflated
#specify 0 inflation with glmmTMB
guards.global.cc.0inflation<-glmmTMB(UniqueGuards~StandardizedElytra+StandardizedGlobalCC+StandardizedScansSeen+Period+(1|Condo/ID), data=surveydataMBeetleset, family="poisson", ziformula=~1)
summary(guards.global.cc.0inflation)## Family: poisson ( log )
## Formula:
## UniqueGuards ~ StandardizedElytra + StandardizedGlobalCC + StandardizedScansSeen +
## Period + (1 | Condo/ID)
## Zero inflation: ~1
## Data: surveydataMBeetleset
##
## AIC BIC logLik deviance df.resid
## 1609.7 1641.7 -796.9 1593.7 392
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## ID:Condo (Intercept) 6.624e-02 2.574e-01
## Condo (Intercept) 9.830e-10 3.135e-05
## Number of obs: 400, groups: ID:Condo, 200; Condo, 12
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.03343 0.04792 21.568 < 2e-16 ***
## StandardizedElytra 0.20293 0.04345 4.670 3.01e-06 ***
## StandardizedGlobalCC -0.03277 0.03474 -0.943 0.345
## StandardizedScansSeen 0.25658 0.04239 6.053 1.42e-09 ***
## Period1 -0.03301 0.03170 -1.041 0.298
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.1382 0.2474 -8.641 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Anova(guards.global.cc.0inflation, type=3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: UniqueGuards
## Chisq Df Pr(>Chisq)
## (Intercept) 465.1598 1 < 2.2e-16 ***
## StandardizedElytra 21.8080 1 3.013e-06 ***
## StandardizedGlobalCC 0.8901 1 0.3455
## StandardizedScansSeen 36.6400 1 1.421e-09 ***
## Period 1.0843 1 0.2977
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(allEffects(guards.global.cc.0inflation))
check.assumptions(guards.global.cc.0inflation)
#test for zero inflation and overdispersion with DHARMa
guards.global.cc.0inflation.resid<-simulateResiduals(guards.global.cc.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(guards.global.cc.0inflation.resid)
testDispersion(guards.global.cc.0inflation.resid)##
## DHARMa nonparametric dispersion test via sd of residuals fitted
## vs. simulated
##
## data: simulationOutput
## ratioObsSim = 0.84229, p-value < 2.2e-16
## alternative hypothesis: two.sidedtestZeroInflation(guards.global.cc.0inflation.resid)##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.052, p-value = 0.648
## alternative hypothesis: two.sided#the model is now underdispersed and is not zero inflated
Lays
lays.global.cc<-glmer(Lays~StandardizedElytra+StandardizedGlobalCC+StandardizedScansSeen+Period+(1|Condo/ID), data=surveydataFBeetleset, family="poisson")
summary(lays.global.cc)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula:
## Lays ~ StandardizedElytra + StandardizedGlobalCC + StandardizedScansSeen +
## Period + (1 | Condo/ID)
## Data: surveydataFBeetleset
##
## AIC BIC logLik deviance df.resid
## 1339.5 1367.4 -662.7 1325.5 393
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6128 -0.7027 -0.1262 0.6087 2.9384
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID:Condo (Intercept) 0.06861 0.2619
## Condo (Intercept) 0.01308 0.1144
## Number of obs: 400, groups: ID:Condo, 200; Condo, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.51706 0.05637 9.173 < 2e-16 ***
## StandardizedElytra -0.03570 0.04328 -0.825 0.409
## StandardizedGlobalCC -0.02574 0.04595 -0.560 0.575
## StandardizedScansSeen 0.24914 0.04320 5.767 8.08e-09 ***
## Period1 -0.03615 0.03710 -0.974 0.330
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) StndrE StnGCC StndSS
## StndrdzdEly 0.023
## StndrdzdGCC 0.014 -0.010
## StndrdzdScS -0.160 -0.063 0.009
## Period1 0.024 -0.003 -0.012 -0.007Anova(lays.global.cc, type=3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: Lays
## Chisq Df Pr(>Chisq)
## (Intercept) 84.1426 1 < 2.2e-16 ***
## StandardizedElytra 0.6803 1 0.4095
## StandardizedGlobalCC 0.3137 1 0.5754
## StandardizedScansSeen 33.2556 1 8.081e-09 ***
## Period 0.9494 1 0.3299
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(allEffects(lays.global.cc))
check.assumptions(lays.global.cc)
#test for zero inflation and overdispersion with DHARMa
lays.global.cc.resid<-simulateResiduals(lays.global.cc, n=250)
plot(lays.global.cc.resid)
testDispersion(lays.global.cc.resid)##
## DHARMa nonparametric dispersion test via sd of residuals fitted
## vs. simulated
##
## data: simulationOutput
## ratioObsSim = 0.98959, p-value = 0.848
## alternative hypothesis: two.sidedtestZeroInflation(lays.global.cc.resid)##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0031, p-value = 1
## alternative hypothesis: two.sided#the model is not overdispersed nor zero inflated