Cij
These files are made from very very old script that Robin and I wrote for 2-15 survey analysis. The files used were called: FFB2016SurveyAnalysis.R and FFB2017SurveyAnalysis_fixed.R
## ── Attaching packages ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──## ✔ ggplot2 3.1.0 ✔ purrr 0.3.3
## ✔ tibble 2.1.3 ✔ dplyr 0.8.5
## ✔ tidyr 1.0.0 ✔ stringr 1.3.1
## ✔ readr 1.3.1 ✔ forcats 0.3.0## Warning: package 'tibble' was built under R version 3.5.2## Warning: package 'tidyr' was built under R version 3.5.2## Warning: package 'purrr' was built under R version 3.5.2## Warning: package 'dplyr' was built under R version 3.5.2## ── Conflicts ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()## -------------------------------------------------------------------------## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)## -------------------------------------------------------------------------##
## Attaching package: 'plyr'## The following objects are masked from 'package:dplyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize## The following object is masked from 'package:purrr':
##
## compact##2016 pooled
##
## Pearson's product-moment correlation
##
## data: MaleBeetleset_2016$beetle_elytra and MaleBeetleset_2016$weighted.avg.male.partner.elytra
## t = 1.378, df = 204, p-value = 0.1697
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.04120595 0.22972177
## sample estimates:
## cor
## 0.09603632
##
## Pearson's product-moment correlation
##
## data: MaleBeetleset_2016$beetle_thoracichorn and MaleBeetleset_2016$weighted.avg.partner.thoracichorn
## t = 0.61281, df = 204, p-value = 0.5407
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.09438887 0.17852075
## sample estimates:
## cor
## 0.04286554
## Warning: `...` is not empty.
##
## We detected these problematic arguments:
## * `needs_dots`
##
## These dots only exist to allow future extensions and should be empty.
## Did you misspecify an argument?## # A tibble: 13 x 2
## Pop n
## <fct> <int>
## 1 PDR-331 33
## 2 PDR-607 40
## 3 PDR-613 15
## 4 PDR-614 10
## 5 PDR-619 18
## 6 PDR-633 4
## 7 PDR-670 11
## 8 PDR-674 2
## 9 PDR-680 22
## 10 PDR-681 6
## 11 PDR-769 8
## 12 PDR-771 2
## 13 PDR-795 10##2017 Pooled 
##
## Pearson's product-moment correlation
##
## data: MaleBeetleset_2017$beetle_thoracichorn and MaleBeetleset_2017$weighted.avg.partner.thoracichorn
## t = -0.66526, df = 171, p-value = 0.5068
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.19850324 0.09914406
## sample estimates:
## cor
## -0.05080775
2016
## Pop corr up low df pval
## 1 PDR-331 0.02 -0.21 0.25 70 0.86
## 2 PDR-607 -0.01 -0.41 0.39 22 0.96
## 3 PDR-611 0.24 -0.81 0.93 3 0.69
## 4 PDR-613 0.06 -0.41 0.50 17 0.81
## 5 PDR-614 0.12 -0.42 0.60 13 0.66
## 6 PDR-617 -0.16 -0.86 0.75 4 0.76
## 7 PDR-618 0.69 0.10 0.92 8 0.03
## 8 PDR-619 -0.34 -0.70 0.15 16 0.16
## 9 PDR-633 0.08 -0.78 0.84 4 0.88
## 10 PDR-680 -0.40 -0.82 0.31 8 0.25
## 11 PDR-681 -0.54 -0.82 -0.03 13 0.04
## 12 PDR-748 -0.46 -0.93 0.56 4 0.36#2017
## Pop corr up low df pval
## 1 PDR-331 0.07 -0.28 0.41 31 0.68
## 2 PDR-607 -0.07 -0.38 0.24 38 0.65
## 3 PDR-613 0.43 -0.11 0.77 13 0.11
## 4 PDR-614 -0.36 -0.81 0.35 8 0.31
## 5 PDR-619 -0.49 -0.78 -0.03 16 0.04
## 6 PDR-670 -0.55 -0.87 0.07 9 0.08
## 7 PDR-680 -0.36 -0.68 0.08 20 0.10
## 8 PDR-681 0.13 -0.76 0.85 4 0.81
## 9 PDR-769 -0.54 -0.90 0.27 6 0.17
## 10 PDR-795 0.17 -0.51 0.72 8 0.64#### Social Selection
library(glmmTMB)
library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recode## The following object is masked from 'package:purrr':
##
## somelibrary(visreg)
library(DHARMa)## Warning: package 'DHARMa' was built under R version 3.5.2library(plyr)
mod2016<-glmmTMB(num.scans.guards~ beetle_thoracichorn+weighted.avg.partner.thoracichorn +num.scans.seen + (1|Pop), data=MaleBeetleset_2016, family="poisson", dispformula = ~num.scans.seen)
res.mod2016 <- simulateResiduals(mod2016, refit = F, n= 500)## 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.newdata = MaleBeetleset_2016
newdata$Pop = NA
newdata <- droplevels(newdata)
pred = rank(predict(mod2016, newdata = newdata))
plotResiduals(pred, res.mod2016$scaledResiduals)
testUniformity(res.mod2016)
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.10437, p-value = 0.02248
## alternative hypothesis: two-sidedtestZeroInflation(res.mod2016)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.1693, p-value = 0.068
## alternative hypothesis: two.sidedtestDispersion(res.mod2016)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted
## vs. simulated
##
## data: simulationOutput
## ratioObsSim = 0.92741, p-value = 0.688
## alternative hypothesis: two.sidedAnova(mod2016, type=3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: num.scans.guards
## Chisq Df Pr(>Chisq)
## (Intercept) 14.2269 1 0.000162 ***
## beetle_thoracichorn 11.9312 1 0.000552 ***
## weighted.avg.partner.thoracichorn 5.3439 1 0.020794 *
## num.scans.seen 215.3997 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## 2017
mod2017<-glmmTMB(num.scans.guards~ beetle_thoracichorn+weighted.avg.partner.thoracichorn +num.scans.seen + (1|Pop), data=MaleBeetleset_2017, family="nbinom2", dispformula = ~num.scans.seen)
res.mod2017 <- simulateResiduals(mod2017, refit = F, n= 500)## 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.newdata = MaleBeetleset_2017
newdata$Pop = NA
newdata <- droplevels(newdata)
pred = rank(predict(mod2017, newdata = newdata))
plotResiduals(pred, res.mod2017$scaledResiduals)
testUniformity(res.mod2017)
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.071746, p-value = 0.3353
## alternative hypothesis: two-sidedtestZeroInflation(res.mod2017)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0884, p-value = 0.34
## alternative hypothesis: two.sidedtestDispersion(res.mod2017)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted
## vs. simulated
##
## data: simulationOutput
## ratioObsSim = 0.76105, p-value = 0.092
## alternative hypothesis: two.sidedAnova(mod2017, type=3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: num.scans.guards
## Chisq Df Pr(>Chisq)
## (Intercept) 9.4660 1 0.0020931 **
## beetle_thoracichorn 13.2767 1 0.0002687 ***
## weighted.avg.partner.thoracichorn 1.3546 1 0.2444830
## num.scans.seen 88.7211 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(mod2017)## Family: nbinom2 ( log )
## Formula:
## num.scans.guards ~ beetle_thoracichorn + weighted.avg.partner.thoracichorn +
## num.scans.seen + (1 | Pop)
## Dispersion: ~num.scans.seen
## Data: MaleBeetleset_2017
##
## AIC BIC logLik deviance df.resid
## 517.0 539.1 -251.5 503.0 166
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Pop (Intercept) 0.1078 0.3284
## Number of obs: 173, groups: Pop, 10
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.460019 0.474542 -3.077 0.002093
## beetle_thoracichorn 0.291572 0.080020 3.644 0.000269
## weighted.avg.partner.thoracichorn -0.145254 0.124805 -1.164 0.244483
## num.scans.seen 0.023860 0.002533 9.419 < 2e-16
##
## (Intercept) **
## beetle_thoracichorn ***
## weighted.avg.partner.thoracichorn
## num.scans.seen ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Dispersion model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.91708 0.71308 -1.286 0.1984
## num.scans.seen 0.02464 0.01043 2.362 0.0182 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# #PC1 assortment ---> Vince doesn't have this now.
# cor.test(MaleBeetleset_2016_331$PC1.mean, MaleBeetleset_2016_331$WeightedPC1Touch5CM, method=c("pearson"))
# plot(MaleBeetleset_2016_331$WeightedPC1Touch5CM~MaleBeetleset_2016_331$Focal.PC1, xlab="Focal PC1", ylab="Partner PC1")
# #r=-0.13, P=0.011