library(knitr)
library(rmarkdown)
library(ggplot2)
library(lme4)
## Loading required package: Matrix
library(MuMIn)
library(vegan)
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.6-4
library(emmeans)
library(lmerTest)
##
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
##
## lmer
## The following object is masked from 'package:stats':
##
## step
library(tibble)
library(plyr)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:plyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(readr)
library(janitor)
##
## Attaching package: 'janitor'
## The following objects are masked from 'package:stats':
##
## chisq.test, fisher.test
library(pbkrtest)
library(tidyr)
##
## Attaching package: 'tidyr'
## The following objects are masked from 'package:Matrix':
##
## expand, pack, unpack
library(ggpubr)
##
## Attaching package: 'ggpubr'
## The following object is masked from 'package:plyr':
##
## mutate
library(ggthemes)
library(stringr)
library(car)
## Loading required package: carData
##
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
##
## recode
library(ordinal)
##
## Attaching package: 'ordinal'
## The following object is masked from 'package:dplyr':
##
## slice
library(RVAideMemoire)
## *** Package RVAideMemoire v 0.9-81-2 ***
##
## Attaching package: 'RVAideMemoire'
## The following object is masked from 'package:lme4':
##
## dummy
library(scales)
##
## Attaching package: 'scales'
## The following object is masked from 'package:readr':
##
## col_factor
library(glmmTMB)
## Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
## TMB was built with Matrix version 1.5.3
## Current Matrix version is 1.5.4
## Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' package
library(nlme)
##
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
##
## collapse
## The following object is masked from 'package:lme4':
##
## lmList
library(moments)
library(extrafont)
## Registering fonts with R
extrafont::loadfonts(quiet=TRUE)
Input and sort data
library(readr)
Behaviour_and_resp <- read_csv("C:/Users/jc819096/OneDrive - James Cook University/DATA/Lab 3/Behaviour data/Behaviour_and_resp_JMD2.csv")
## Rows: 544 Columns: 37
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (17): Maternal, Maternal_GF, Maternal_GM, Paternal, Paternal_GF, Paterna...
## dbl (20): Tank, Density, Female, Male, Age, Wet_weight, Weight, RMR, Mo2_res...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
View(Behaviour_and_resp)
Behaviour_and_resp <- Behaviour_and_resp %>%
unite(Temp:Control_CO2, col="Juv_treat", sep="_", remove=FALSE)
head(Behaviour_and_resp)
## # A tibble: 6 × 38
## Tank Density Female Male Maternal Maternal_GF Maternal_GM Paternal
## <dbl> <dbl> <dbl> <dbl> <chr> <chr> <chr> <chr>
## 1 100 18 132 9 DC D C CA
## 2 30 17 132 9 DC D C CA
## 3 113 14 75 9 AE A E CA
## 4 107 15 190 139 AE A E FA
## 5 32 17 190 139 AE A E FA
## 6 117 19 42 138 DA D A CA
## # ℹ 30 more variables: Paternal_GF <chr>, Paternal_GM <chr>,
## # Parental_treatment <chr>, DOM <chr>, DOT <chr>, DOD <chr>, Age <dbl>,
## # am_pm <chr>, Wet_weight <dbl>, Weight <dbl>, Juv_treat <chr>,
## # Control_CO2 <chr>, Temp <chr>, Chamber_ID <chr>, Parental_number <chr>,
## # RMR <dbl>, Mo2_rest <dbl>, MMR <dbl>, AS <dbl>, FS <dbl>, Mo2_max <dbl>,
## # Mo2_AS <dbl>, `Behaviour day` <chr>, Cam <dbl>, `Video ID` <chr>,
## # Time_emerge_s <dbl>, Time_emerge_min <dbl>, Bold <dbl>, Sum_score <dbl>, …
str(Behaviour_and_resp)
## tibble [544 × 38] (S3: tbl_df/tbl/data.frame)
## $ Tank : num [1:544] 100 30 113 107 32 117 39 167 113 75 ...
## $ Density : num [1:544] 18 17 14 15 17 19 16 3 14 17 ...
## $ Female : num [1:544] 132 132 75 190 190 42 190 185 75 190 ...
## $ Male : num [1:544] 9 9 9 139 139 138 197 190 9 139 ...
## $ Maternal : chr [1:544] "DC" "DC" "AE" "AE" ...
## $ Maternal_GF : chr [1:544] "D" "D" "A" "A" ...
## $ Maternal_GM : chr [1:544] "C" "C" "E" "E" ...
## $ Paternal : chr [1:544] "CA" "CA" "CA" "FA" ...
## $ Paternal_GF : chr [1:544] "C" "C" "C" "F" ...
## $ Paternal_GM : chr [1:544] "A" "A" "A" "A" ...
## $ Parental_treatment: chr [1:544] "HHHC" "HHHC" "HHHC" "CCCC" ...
## $ DOM : chr [1:544] "8/01/2022" "8/01/2022" "6/03/2022" "16/02/2022" ...
## $ DOT : chr [1:544] "6/05/2022" "6/05/2022" "2/07/2022" "27/06/2022" ...
## $ DOD : chr [1:544] "7/05/2022" "7/05/2022" "3/07/2022" "28/06/2022" ...
## $ Age : num [1:544] 119 119 119 132 121 130 104 129 119 131 ...
## $ am_pm : chr [1:544] "pm" "am" "pm" "pm" ...
## $ Wet_weight : num [1:544] 1.55 1.45 2.04 1.59 1.42 1.26 1.24 2.61 1.9 1.86 ...
## $ Weight : num [1:544] 0.00155 0.00145 0.00204 0.00159 0.00142 0.00126 0.00124 0.00261 0.0019 0.00186 ...
## $ Juv_treat : chr [1:544] "zero_control" "elevated_control" "zero_control" "zero_control" ...
## $ Control_CO2 : chr [1:544] "control" "control" "control" "control" ...
## $ Temp : chr [1:544] "zero" "elevated" "zero" "zero" ...
## $ Chamber_ID : chr [1:544] "A" "D" "D" "C" ...
## $ Parental_number : chr [1:544] "65" "65" "9b" "92b" ...
## $ RMR : num [1:544] 216 205 219 244 205 ...
## $ Mo2_rest : num [1:544] 0.34 0.3 0.45 0.39 0.29 0.21 0.31 0.44 0.32 0.33 ...
## $ MMR : num [1:544] 592 512 528 719 585 ...
## $ AS : num [1:544] 376 307 309 475 380 ...
## $ FS : num [1:544] 2.74 2.5 2.41 2.95 2.85 ...
## $ Mo2_max : num [1:544] 0.918 0.742 1.078 1.144 0.831 ...
## $ Mo2_AS : num [1:544] 0.58 0.45 0.63 0.76 0.54 0.82 0.54 2.12 0.84 0.68 ...
## $ Behaviour day : chr [1:544] "7/05/2022" "7/05/2022" "3/07/2022" "28/06/2022" ...
## $ Cam : num [1:544] 3 4 4 4 4 3 4 4 4 2 ...
## $ Video ID : chr [1:544] "4349" "4352" "4734" "4674" ...
## $ Time_emerge_s : num [1:544] 3 23 0 0 0 0 0 0 0 0 ...
## $ Time_emerge_min : num [1:544] 0.05 0.383 0 0 0 ...
## $ Bold : num [1:544] 4 4 1 1 1 1 1 2 2 2 ...
## $ Sum_score : num [1:544] 8 8 2 2 2 2 2 4 4 4 ...
## $ Sum_score_v2 : num [1:544] 5 5 2 2 2 2 2 3 3 3 ...
Behaviour_and_resp$Tank= factor (Behaviour_and_resp$Tank)
Behaviour_and_resp$Parental_number = factor (Behaviour_and_resp$Parental_number)
Behaviour_and_resp$Parental_treatment = factor (Behaviour_and_resp$Parental_treatment, levels=c("CCCC", "CCCH", "HHCC", "CCHC","HHHC"))
Behaviour_and_resp$Temp = factor (Behaviour_and_resp$Temp)
Behaviour_and_resp$Control_CO2 = factor (Behaviour_and_resp$Control_CO2)
Behaviour_and_resp$am_pm = factor (Behaviour_and_resp$am_pm)
str(Behaviour_and_resp)
## tibble [544 × 38] (S3: tbl_df/tbl/data.frame)
## $ Tank : Factor w/ 162 levels "1","2","3","4",..: 89 23 102 96 25 106 31 147 102 64 ...
## $ Density : num [1:544] 18 17 14 15 17 19 16 3 14 17 ...
## $ Female : num [1:544] 132 132 75 190 190 42 190 185 75 190 ...
## $ Male : num [1:544] 9 9 9 139 139 138 197 190 9 139 ...
## $ Maternal : chr [1:544] "DC" "DC" "AE" "AE" ...
## $ Maternal_GF : chr [1:544] "D" "D" "A" "A" ...
## $ Maternal_GM : chr [1:544] "C" "C" "E" "E" ...
## $ Paternal : chr [1:544] "CA" "CA" "CA" "FA" ...
## $ Paternal_GF : chr [1:544] "C" "C" "C" "F" ...
## $ Paternal_GM : chr [1:544] "A" "A" "A" "A" ...
## $ Parental_treatment: Factor w/ 5 levels "CCCC","CCCH",..: 5 5 5 1 2 3 1 2 5 1 ...
## $ DOM : chr [1:544] "8/01/2022" "8/01/2022" "6/03/2022" "16/02/2022" ...
## $ DOT : chr [1:544] "6/05/2022" "6/05/2022" "2/07/2022" "27/06/2022" ...
## $ DOD : chr [1:544] "7/05/2022" "7/05/2022" "3/07/2022" "28/06/2022" ...
## $ Age : num [1:544] 119 119 119 132 121 130 104 129 119 131 ...
## $ am_pm : Factor w/ 2 levels "am","pm": 2 1 2 2 1 2 2 2 2 1 ...
## $ Wet_weight : num [1:544] 1.55 1.45 2.04 1.59 1.42 1.26 1.24 2.61 1.9 1.86 ...
## $ Weight : num [1:544] 0.00155 0.00145 0.00204 0.00159 0.00142 0.00126 0.00124 0.00261 0.0019 0.00186 ...
## $ Juv_treat : chr [1:544] "zero_control" "elevated_control" "zero_control" "zero_control" ...
## $ Control_CO2 : Factor w/ 2 levels "CO","control": 2 2 2 2 2 1 2 2 2 1 ...
## $ Temp : Factor w/ 2 levels "elevated","zero": 2 1 2 2 1 1 2 1 2 2 ...
## $ Chamber_ID : chr [1:544] "A" "D" "D" "C" ...
## $ Parental_number : Factor w/ 24 levels "13","15","4",..: 8 8 24 20 6 14 19 21 24 20 ...
## $ RMR : num [1:544] 216 205 219 244 205 ...
## $ Mo2_rest : num [1:544] 0.34 0.3 0.45 0.39 0.29 0.21 0.31 0.44 0.32 0.33 ...
## $ MMR : num [1:544] 592 512 528 719 585 ...
## $ AS : num [1:544] 376 307 309 475 380 ...
## $ FS : num [1:544] 2.74 2.5 2.41 2.95 2.85 ...
## $ Mo2_max : num [1:544] 0.918 0.742 1.078 1.144 0.831 ...
## $ Mo2_AS : num [1:544] 0.58 0.45 0.63 0.76 0.54 0.82 0.54 2.12 0.84 0.68 ...
## $ Behaviour day : chr [1:544] "7/05/2022" "7/05/2022" "3/07/2022" "28/06/2022" ...
## $ Cam : num [1:544] 3 4 4 4 4 3 4 4 4 2 ...
## $ Video ID : chr [1:544] "4349" "4352" "4734" "4674" ...
## $ Time_emerge_s : num [1:544] 3 23 0 0 0 0 0 0 0 0 ...
## $ Time_emerge_min : num [1:544] 0.05 0.383 0 0 0 ...
## $ Bold : num [1:544] 4 4 1 1 1 1 1 2 2 2 ...
## $ Sum_score : num [1:544] 8 8 2 2 2 2 2 4 4 4 ...
## $ Sum_score_v2 : num [1:544] 5 5 2 2 2 2 2 3 3 3 ...
Behaviour_and_resp$Juv_treat = factor (Behaviour_and_resp$Juv_treat, levels=c("zero_control", "elevated_control", "zero_CO", "elevated_CO"))
head(Behaviour_and_resp)
## # A tibble: 6 × 38
## Tank Density Female Male Maternal Maternal_GF Maternal_GM Paternal
## <fct> <dbl> <dbl> <dbl> <chr> <chr> <chr> <chr>
## 1 100 18 132 9 DC D C CA
## 2 30 17 132 9 DC D C CA
## 3 113 14 75 9 AE A E CA
## 4 107 15 190 139 AE A E FA
## 5 32 17 190 139 AE A E FA
## 6 117 19 42 138 DA D A CA
## # ℹ 30 more variables: Paternal_GF <chr>, Paternal_GM <chr>,
## # Parental_treatment <fct>, DOM <chr>, DOT <chr>, DOD <chr>, Age <dbl>,
## # am_pm <fct>, Wet_weight <dbl>, Weight <dbl>, Juv_treat <fct>,
## # Control_CO2 <fct>, Temp <fct>, Chamber_ID <chr>, Parental_number <fct>,
## # RMR <dbl>, Mo2_rest <dbl>, MMR <dbl>, AS <dbl>, FS <dbl>, Mo2_max <dbl>,
## # Mo2_AS <dbl>, `Behaviour day` <chr>, Cam <dbl>, `Video ID` <chr>,
## # Time_emerge_s <dbl>, Time_emerge_min <dbl>, Bold <dbl>, Sum_score <dbl>, …
Exploring raw data
ggplot(Behaviour_and_resp) + geom_bar(aes(x=Bold, y=(..count..)/sum(..count..)))
## Warning: The dot-dot notation (`..count..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(count)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
ggplot(Behaviour_and_resp, aes(x=Juv_treat, y=Bold, colour=Parental_treatment))+ geom_boxplot()
ggplot(Behaviour_and_resp) + geom_density(aes(x=Bold, y=((..count..)/sum(..count..)), fill = Parental_treatment, colour=Parental_treatment, alpha=0.3), position="dodge") + facet_grid(~Juv_treat)+ theme_classic()
## Warning: Width not defined
## ℹ Set with `position_dodge(width = ...)`
ggplot(data=Behaviour_and_resp, aes(x=Bold, color=Parental_treatment, fill= Parental_treatment)) + geom_density(aes(y=..density..), position="identity", alpha=0.3, adjust =3) + facet_grid(~Juv_treat) + theme_classic()
Raw data checks
hist(BehavCount$n)
mean(BehavCount$n)
## [1] 1.133333
var(BehavCount$n)
## [1] 1.547947
hist(BehavCount$Bold)
plot(n~Bold, BehavCount, log="y")
## Warning in xy.coords(x, y, xlabel, ylabel, log): 197 y values <= 0 omitted from
## logarithmic plot
with(BehavCount, lines(lowess(n~Bold)))
library(vcd)
## Loading required package: grid
Ord_plot(BehavCount$Bold)
Ord_plot(BehavCount$Bold, tol=0.2)
library(vcd)
fit <- goodfit(BehavCount$n, type='poisson')
summary(fit)
##
## Goodness-of-fit test for poisson distribution
##
## X^2 df P(> X^2)
## Likelihood Ratio 50.76845 5 9.646799e-10
rootogram(fit)
fit <- goodfit(BehavCount$n, type='nbinomial')
summary(fit)
##
## Goodness-of-fit test for nbinomial distribution
##
## X^2 df P(> X^2)
## Likelihood Ratio 21.69608 4 0.0002303572
rootogram(fit)
distplot(BehavCount$n, type='poisson')
distplot(BehavCount$n, type='nbinomial')
dat.nb.tab<-table(BehavCount$n==0)
dat.nb.tab/sum(dat.nb.tab)
##
## FALSE TRUE
## 0.5895833 0.4104167
mu <- mean(BehavCount$n)
cnts <- rpois(1000, mu)
dat.nb.tabE <- table(cnts == 0)
dat.nb.tabE/sum(dat.nb.tabE)
##
## FALSE TRUE
## 0.677 0.323
model
library(MASS)
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
negb.Bold = glm.nb(n ~ Bold * Parental_treatment * Control_CO2 * Temp, data = BehavCount)
post model checks
plot(negb.Bold)
library(DHARMa)
## This is DHARMa 0.4.6. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')
plot(simulateResiduals(negb.Bold))
dat.resid <- sum(resid(negb.Bold, type = "pearson")^2)
1 - pchisq(dat.resid, negb.Bold$df.resid)
## [1] 0.4618099
testUniformity(negb.Bold)
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.040057, p-value = 0.4244
## alternative hypothesis: two-sided
testDispersion(simulateResiduals(negb.Bold, refit=T))
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in glm.nb(formula = n ~ Bold * Parental_treatment * Control_CO2 * :
## alternation limit reached
##
## DHARMa nonparametric dispersion test via mean deviance residual fitted
## vs. simulated-refitted
##
## data: simulateResiduals(negb.Bold, refit = T)
## dispersion = 0.93276, p-value = 0.008
## alternative hypothesis: two.sided
m1 <- negb.Bold
dispfun <- function(m) {
r <- residuals(m,type="pearson")
n <- df.residual(m)
dsq <- sum(r^2)
c(dsq=dsq,n=n,disp=dsq/n)
}
options(digits=4)
dispfun(m1)
## dsq n disp
## 442.181 440.000 1.005
testZeroInflation(simulateResiduals(negb.Bold, refit=T))
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in glm.nb(formula = n ~ Bold * Parental_treatment * Control_CO2 * :
## alternation limit reached
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1, p-value = 0.6
## alternative hypothesis: two.sided
library(AER)
## Loading required package: lmtest
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## Loading required package: sandwich
## Loading required package: survival
check<- data.frame(obs=BehavCount$n, fitted=fitted(negb.Bold))
library(ggplot2)
ggplot(check, aes(y=fitted, x=obs)) + geom_point()
negb.Bold %>% performance::check_model()
## Model has interaction terms. VIFs might be inflated.
## You may check multicollinearity among predictors of a model without
## interaction terms.
negb.Bold %>% performance::check_overdispersion()
## # Overdispersion test
##
## dispersion ratio = 1.005
## Pearson's Chi-Squared = 442.181
## p-value = 0.462
## No overdispersion detected.
negb.Bold %>%performance::check_zeroinflation()
## # Check for zero-inflation
##
## Observed zeros: 197
## Predicted zeros: 189
## Ratio: 0.96
## Model seems ok, ratio of observed and predicted zeros is within the
## tolerance range.
negb.Bold %>%simulateResiduals(plot=TRUE)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.9276 0.8205 0.1794 0.1818 0.9608 0.6195 0.9416 0.8173 0.2709 0.9509 0.04781 0.8448 0.2309 0.6216 0.1631 0.8744 0.2694 0.6226 0.2293 0.9568 ...
model analysis
summary(negb.Bold)
##
## Call:
## glm.nb(formula = n ~ Bold * Parental_treatment * Control_CO2 *
## Temp, data = BehavCount, init.theta = 2.976457161, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.894 -1.282 -0.160 0.547 2.415
##
## Coefficients:
## Estimate Std. Error
## (Intercept) -0.73780 0.61455
## Bold 0.21852 0.17053
## Parental_treatmentCCCH -0.12468 0.87150
## Parental_treatmentCCHC 0.06242 0.86474
## Parental_treatmentHHCC 0.28263 0.83014
## Parental_treatmentHHHC 0.32483 0.87654
## Control_CO2control 0.50228 0.81707
## Tempzero 0.50073 0.82868
## Bold:Parental_treatmentCCCH 0.05856 0.23967
## Bold:Parental_treatmentCCHC -0.01954 0.24077
## Bold:Parental_treatmentHHCC -0.01835 0.23146
## Bold:Parental_treatmentHHHC -0.04323 0.24552
## Bold:Control_CO2control -0.08447 0.23009
## Parental_treatmentCCCH:Control_CO2control -0.07288 1.18206
## Parental_treatmentCCHC:Control_CO2control 0.17746 1.14186
## Parental_treatmentHHCC:Control_CO2control -1.22482 1.17350
## Parental_treatmentHHHC:Control_CO2control -0.41102 1.21311
## Bold:Tempzero -0.12873 0.23509
## Parental_treatmentCCCH:Tempzero 0.49332 1.16296
## Parental_treatmentCCHC:Tempzero -0.02695 1.16968
## Parental_treatmentHHCC:Tempzero -0.65194 1.15326
## Parental_treatmentHHHC:Tempzero 0.30998 1.18130
## Control_CO2control:Tempzero 0.20301 1.09990
## Bold:Parental_treatmentCCCH:Control_CO2control -0.08013 0.33197
## Bold:Parental_treatmentCCHC:Control_CO2control -0.04574 0.32374
## Bold:Parental_treatmentHHCC:Control_CO2control 0.29997 0.32390
## Bold:Parental_treatmentHHHC:Control_CO2control -0.00299 0.34521
## Bold:Parental_treatmentCCCH:Tempzero -0.19280 0.33228
## Bold:Parental_treatmentCCHC:Tempzero -0.00437 0.33319
## Bold:Parental_treatmentHHCC:Tempzero 0.16501 0.32435
## Bold:Parental_treatmentHHHC:Tempzero -0.15052 0.34347
## Bold:Control_CO2control:Tempzero -0.08067 0.31877
## Parental_treatmentCCCH:Control_CO2control:Tempzero 0.02315 1.56086
## Parental_treatmentCCHC:Control_CO2control:Tempzero -0.63386 1.55573
## Parental_treatmentHHCC:Control_CO2control:Tempzero 1.41109 1.57644
## Parental_treatmentHHHC:Control_CO2control:Tempzero -0.50968 1.62896
## Bold:Parental_treatmentCCCH:Control_CO2control:Tempzero 0.12155 0.45581
## Bold:Parental_treatmentCCHC:Control_CO2control:Tempzero 0.20006 0.45084
## Bold:Parental_treatmentHHCC:Control_CO2control:Tempzero -0.35425 0.44725
## Bold:Parental_treatmentHHHC:Control_CO2control:Tempzero 0.27210 0.47810
## z value Pr(>|z|)
## (Intercept) -1.20 0.23
## Bold 1.28 0.20
## Parental_treatmentCCCH -0.14 0.89
## Parental_treatmentCCHC 0.07 0.94
## Parental_treatmentHHCC 0.34 0.73
## Parental_treatmentHHHC 0.37 0.71
## Control_CO2control 0.61 0.54
## Tempzero 0.60 0.55
## Bold:Parental_treatmentCCCH 0.24 0.81
## Bold:Parental_treatmentCCHC -0.08 0.94
## Bold:Parental_treatmentHHCC -0.08 0.94
## Bold:Parental_treatmentHHHC -0.18 0.86
## Bold:Control_CO2control -0.37 0.71
## Parental_treatmentCCCH:Control_CO2control -0.06 0.95
## Parental_treatmentCCHC:Control_CO2control 0.16 0.88
## Parental_treatmentHHCC:Control_CO2control -1.04 0.30
## Parental_treatmentHHHC:Control_CO2control -0.34 0.73
## Bold:Tempzero -0.55 0.58
## Parental_treatmentCCCH:Tempzero 0.42 0.67
## Parental_treatmentCCHC:Tempzero -0.02 0.98
## Parental_treatmentHHCC:Tempzero -0.57 0.57
## Parental_treatmentHHHC:Tempzero 0.26 0.79
## Control_CO2control:Tempzero 0.18 0.85
## Bold:Parental_treatmentCCCH:Control_CO2control -0.24 0.81
## Bold:Parental_treatmentCCHC:Control_CO2control -0.14 0.89
## Bold:Parental_treatmentHHCC:Control_CO2control 0.93 0.35
## Bold:Parental_treatmentHHHC:Control_CO2control -0.01 0.99
## Bold:Parental_treatmentCCCH:Tempzero -0.58 0.56
## Bold:Parental_treatmentCCHC:Tempzero -0.01 0.99
## Bold:Parental_treatmentHHCC:Tempzero 0.51 0.61
## Bold:Parental_treatmentHHHC:Tempzero -0.44 0.66
## Bold:Control_CO2control:Tempzero -0.25 0.80
## Parental_treatmentCCCH:Control_CO2control:Tempzero 0.01 0.99
## Parental_treatmentCCHC:Control_CO2control:Tempzero -0.41 0.68
## Parental_treatmentHHCC:Control_CO2control:Tempzero 0.90 0.37
## Parental_treatmentHHHC:Control_CO2control:Tempzero -0.31 0.75
## Bold:Parental_treatmentCCCH:Control_CO2control:Tempzero 0.27 0.79
## Bold:Parental_treatmentCCHC:Control_CO2control:Tempzero 0.44 0.66
## Bold:Parental_treatmentHHCC:Control_CO2control:Tempzero -0.79 0.43
## Bold:Parental_treatmentHHHC:Control_CO2control:Tempzero 0.57 0.57
##
## (Dispersion parameter for Negative Binomial(2.9765) family taken to be 1)
##
## Null deviance: 557.11 on 479 degrees of freedom
## Residual deviance: 531.64 on 440 degrees of freedom
## AIC: 1451
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 2.976
## Std. Err.: 0.811
##
## 2 x log-likelihood: -1369.041
Anova(negb.Bold)
## Analysis of Deviance Table (Type II tests)
##
## Response: n
## LR Chisq Df Pr(>Chisq)
## Bold 6.36 1 0.012 *
## Parental_treatment 1.03 4 0.905
## Control_CO2 1.72 1 0.189
## Temp 0.78 1 0.377
## Bold:Parental_treatment 2.65 4 0.617
## Bold:Control_CO2 0.81 1 0.369
## Parental_treatment:Control_CO2 0.94 4 0.918
## Bold:Temp 5.58 1 0.018 *
## Parental_treatment:Temp 0.25 4 0.993
## Control_CO2:Temp 0.44 1 0.505
## Bold:Parental_treatment:Control_CO2 0.64 4 0.959
## Bold:Parental_treatment:Temp 0.97 4 0.914
## Bold:Control_CO2:Temp 0.07 1 0.789
## Parental_treatment:Control_CO2:Temp 0.67 4 0.955
## Bold:Parental_treatment:Control_CO2:Temp 2.05 4 0.726
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Bold = (emmeans(negb.Bold, ~ Parental_treatment * Temp * Control_CO2) %>% as.data.frame)
## NOTE: Results may be misleading due to involvement in interactions
Bold
## Parental_treatment Temp Control_CO2 emmean SE df asymp.LCL asymp.UCL
## CCCC elevated CO -0.0822 0.2427 Inf -0.5580 0.3935
## CCCH elevated CO -0.0312 0.2404 Inf -0.5025 0.4400
## CCHC elevated CO -0.0784 0.2417 Inf -0.5521 0.3952
## HHCC elevated CO 0.1453 0.2223 Inf -0.2903 0.5810
## HHHC elevated CO 0.1129 0.2506 Inf -0.3784 0.6042
## CCCC zero CO 0.0323 0.2291 Inf -0.4167 0.4813
## CCCH zero CO -0.0018 0.2315 Inf -0.4555 0.4519
## CCHC zero CO -0.0040 0.2319 Inf -0.4585 0.4505
## HHCC zero CO 0.1030 0.2270 Inf -0.3419 0.5478
## HHHC zero CO 0.0859 0.2514 Inf -0.4068 0.5785
## CCCC elevated control 0.1666 0.2189 Inf -0.2624 0.5956
## CCCH elevated control -0.0956 0.2409 Inf -0.5677 0.3764
## CCHC elevated control 0.2107 0.2145 Inf -0.2097 0.6310
## HHCC elevated control 0.0693 0.2387 Inf -0.3985 0.5370
## HHHC elevated control -0.0582 0.2650 Inf -0.5776 0.4611
## CCCC zero control 0.2422 0.2122 Inf -0.1737 0.6580
## CCCH zero control 0.2826 0.2108 Inf -0.1306 0.6958
## CCHC zero control 0.2124 0.2142 Inf -0.2074 0.6323
## HHCC zero control 0.3362 0.2050 Inf -0.0656 0.7380
## HHHC zero control 0.1823 0.2418 Inf -0.2916 0.6562
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
plot data
graph2 <- ggplot(data=Behaviour_and_resp, aes(x=Bold, color=Temp, fill= Temp)) + geom_density(aes(y=..density..), position="identity", alpha=0.2, adjust =3, linewidth = 0.8)+ theme_classic() + theme(text = element_text(size=20))
graph2
ggsave("Behaviour2_graph.eps", graph2, height = 6, width = 10, dpi = 320)
## Warning in grid.Call.graphics(C_polygon, x$x, x$y, index): semi-transparency is
## not supported on this device: reported only once per page
graph <- ggplot(data=Behaviour_and_resp, aes(x=Bold, color=Parental_treatment, fill= Parental_treatment)) + geom_density(aes(y=..density..), position="identity", alpha=0.2, adjust =3, linewidth = 0.8) +
facet_grid(~Juv_treat) +
scale_colour_manual(values=c("#56B4E9", "#E69F00", "#D55E00","#009E73", "#999999"), name = "Cross-generation exposure", labels = c("Control", "Parental development", "Grandparental development", "Grandparental post-maturation", "Continuous grandparent")) +
scale_fill_manual(values=c("#56B4E9", "#E69F00", "#D55E00","#009E73", "#999999"), name = "Cross-generation exposure", labels = c("Control", "Parental development", "Grandparental development", "Grandparental post-maturation", "Continuous grandparent")) + labs(x="Bold score", y="Density") +
theme(panel.spacing = unit(1, "cm", data = NULL)) + theme(text = element_text(size=20)) + theme_classic()
print(graph)
ggsave("Behaviour_graph.png", graph, height = 6, width = 14, dpi = 320)
max = which.max(density(Behaviour_and_resp$Bold)$y)
density(Behaviour_and_resp$Bold)$x[max]
## [1] 4.004
Control = Behaviour_and_resp %>% filter(Temp == "zero")
max = which.max(density(Control$Bold)$y)
density(Control$Bold)$x[max]
## [1] 3.006
Elevated = Behaviour_and_resp %>% filter(Temp == "elevated")
max = which.max(density(Elevated$Bold)$y)
density(Elevated$Bold)$x[max]
## [1] 4.005