index_recovery
Marina
2022-12-18
load packages and data
{code not included}
EDA
Table showing # of occurrences reported per 1,000km^2; table
## `summarise()` has grouped output by 'treatmentGroup', 'desert'. You can
## override using the `.groups` argument.## # A tibble: 12 × 6
## # Groups: treatmentGroup, desert [6]
## treatmentGroup desert ez_taxa area_km2 totalOccReported totalORPer1000…¹
## <chr> <chr> <chr> <dbl> <int> <dbl>
## 1 burned Mojave Aves 722 144 200
## 2 burned Mojave Other 722 16 23
## 3 burned San Joaquin Aves 399 116 291
## 4 burned San Joaquin Other 399 4 11
## 5 burned Sonoran Aves 219 16 74
## 6 burned Sonoran Other 219 51 233
## 7 never burned Mojave Aves 73161 8105 111
## 8 never burned Mojave Other 73161 586 9
## 9 never burned San Joaquin Aves 26537 8798 332
## 10 never burned San Joaquin Other 26537 226 9
## 11 never burned Sonoran Aves 27789 6086 220
## 12 never burned Sonoran Other 27789 788 29
## # … with abbreviated variable name ¹totalORPer1000km2## `summarise()` has grouped output by 'treatmentGroup'. You can override using
## the `.groups` argument.## # A tibble: 6 × 4
## # Groups: treatmentGroup [2]
## treatmentGroup desert mean_ndvi sd_ndvi
## <fct> <fct> <dbl> <dbl>
## 1 burned Mojave 0.23 0.02
## 2 burned San Joaquin 0.34 0.04
## 3 burned Sonoran 0.23 0.02
## 4 never burned Mojave 0.13 0.02
## 5 never burned San Joaquin 0.37 0.03
## 6 never burned Sonoran 0.14 0.01## `summarise()` has grouped output by 'treatmentGroup', 'desert', 'ez_taxa'. You
## can override using the `.groups` argument.
## `summarise()` has grouped output by 'treatmentGroup', 'desert'. You can
## override using the `.groups` argument.## # A tibble: 12 × 7
## # Groups: treatmentGroup, desert [6]
## treatmentGroup desert ez_taxa yrlyAvgOccReported yrlyAv…¹ sd_av…² SEM
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 burned Mojave Aves 6 9 6.91 2.82
## 2 burned Mojave Other 2 4 1.41 0.996
## 3 burned San Joaquin Aves 8 19 15.7 7.84
## 4 burned San Joaquin Other 2 4 1.73 2
## 5 burned Sonoran Aves 2 10 4.17 2.95
## 6 burned Sonoran Other 5 20 19.4 11.2
## 7 never burned Mojave Aves 312 5 3.36 1.32
## 8 never burned Mojave Other 28 2 0.359 0.156
## 9 never burned San Joaquin Aves 339 14 12.7 5.00
## 10 never burned San Joaquin Other 14 2 0.588 0.285
## 11 never burned Sonoran Aves 235 10 7.05 2.77
## 12 never burned Sonoran Other 53 3 2.44 1.26
## # … with abbreviated variable names ¹yrlyAvgORPer1000km2, ²sd_avgORper1000km2NDVI and occurrences
##
## burned never burned
## Mojave 21 21
## San Joaquin 21 21
## Sonoran 21 21## Df Sum Sq Mean Sq F value Pr(>F)
## desert 2 0.8119 0.4060 290.84 < 2e-16 ***
## treatmentGroup 1 0.0997 0.0997 71.45 7.16e-14 ***
## Residuals 122 0.1703 0.0014
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = ndvi ~ desert + treatmentGroup, data = ndvi)
##
## $desert
## diff lwr upr p adj
## San Joaquin-Mojave 0.174805861 0.155462394 0.19414933 0.0000000
## Sonoran-Mojave 0.009435481 -0.009907987 0.02877895 0.4809407
## Sonoran-San Joaquin -0.165370380 -0.184713848 -0.14602691 0.0000000##
## Shapiro-Wilk normality test
##
## data: ndvi$ndvi
## W = 0.9166, p-value = 9.167e-07## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = ndvi ~ desert + treatmentGroup, data = ndvi)
##
## $treatmentGroup
## diff lwr upr p adj
## never burned-burned -0.05626842 -0.06944584 -0.043091 0## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
##
## Shapiro-Wilk normality test
##
## data: test$ndvi
## W = 0.93761, p-value = 0.02359## Df Sum Sq Mean Sq F value Pr(>F)
## treatmentGroup 1 0.00621 0.006213 6.099 0.0179 *
## Residuals 40 0.04075 0.001019
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = ndvi ~ treatmentGroup, data = test)
##
## $treatmentGroup
## diff lwr upr p adj
## never burned-burned 0.02432483 0.004418128 0.04423154 0.01789##
## Shapiro-Wilk normality test
##
## data: test$ndvi
## W = 0.89743, p-value = 6.038e-06## Df Sum Sq Mean Sq F value Pr(>F)
## treatmentGroup 1 0.19582 0.19582 551.6 <2e-16 ***
## Residuals 82 0.02911 0.00035
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = ndvi ~ treatmentGroup, data = test)
##
## $treatmentGroup
## diff lwr upr p adj
## never burned-burned -0.09656505 -0.104744 -0.08838612 0##
## Shapiro-Wilk normality test
##
## data: test$ndvi
## W = 0.97545, p-value = 0.4937## Df Sum Sq Mean Sq F value Pr(>F)
## desert 1 0.000025 0.0000246 0.053 0.819
## Residuals 40 0.018554 0.0004638## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = ndvi ~ desert, data = test)
##
## $desert
## diff lwr upr p adj
## Sonoran-Mojave 0.001530101 -0.0119029 0.0149631 0.8190997##
## Shapiro-Wilk normality test
##
## data: test$ndvi
## W = 0.95208, p-value = 0.0767## Df Sum Sq Mean Sq F value Pr(>F)
## desert 1 0.003157 0.0031574 17.13 0.000175 ***
## Residuals 40 0.007373 0.0001843
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = ndvi ~ desert, data = test)
##
## $desert
## diff lwr upr p adj
## Sonoran-Mojave 0.01734086 0.008873093 0.02580863 0.0001746##
## Shapiro-Wilk normality test
##
## data: test$ndvi
## W = 0.89743, p-value = 6.038e-06## Df Sum Sq Mean Sq F value Pr(>F)
## treatmentGroup 1 0.19582 0.19582 551.6 <2e-16 ***
## Residuals 82 0.02911 0.00035
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = ndvi ~ treatmentGroup, data = test)
##
## $treatmentGroup
## diff lwr upr p adj
## never burned-burned -0.09656505 -0.104744 -0.08838612 0## `summarise()` has grouped output by 'treatmentGroup'. You can override using
## the `.groups` argument.## Warning: `fun.y` is deprecated. Use `fun` instead.## `summarise()` has grouped output by 'obsYear', 'treatmentGroup', 'ez_taxa'. You
## can override using the `.groups` argument.
## `summarise()` has grouped output by 'obsYear', 'treatmentGroup', 'ez_taxa'. You
## can override using the `.groups` argument.
## `summarise()` has grouped output by 'obsYear', 'treatmentGroup', 'ez_taxa'. You
## can override using the `.groups` argument.
## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 23 rows containing missing values (geom_smooth).## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 16 rows containing missing values (geom_smooth).
## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 36 rows containing missing values (geom_smooth).
Stats
glm average occurrences reported
## `summarise()` has grouped output by 'treatmentGroup', 'desert', 'ez_taxa'. You
## can override using the `.groups` argument.## # Overdispersion test
##
## dispersion ratio = 1.194
## Pearson's Chi-Squared = 233.974
## p-value = 0.033## Overdispersion detected.##
## Call:
## glm.nb(formula = totalORPer1000km2 ~ desert * treatmentGroup,
## data = test2, init.theta = 1.289509539, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5448 -0.9992 -0.4995 0.3743 2.1358
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.15380 0.15235 7.573 3.64e-14 ***
## desertSan Joaquin 0.99308 0.20967 4.736 2.18e-06 ***
## desertSonoran 0.74576 0.21395 3.486 0.000491 ***
## treatmentGroupburned 0.84853 0.22730 3.733 0.000189 ***
## desertSan Joaquin:treatmentGroupburned -0.19364 0.34121 -0.568 0.570364
## desertSonoran:treatmentGroupburned 0.01193 0.34083 0.035 0.972082
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.2895) family taken to be 1)
##
## Null deviance: 268.63 on 201 degrees of freedom
## Residual deviance: 207.15 on 196 degrees of freedom
## AIC: 1237.3
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.290
## Std. Err.: 0.142
##
## 2 x log-likelihood: -1223.296## Warning in anova.negbin(occ_glm, test = "Chisq"): tests made without re-
## estimating 'theta'## Analysis of Deviance Table
##
## Model: Negative Binomial(1.2895), link: log
##
## Response: totalORPer1000km2
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 201 268.63
## desert 2 28.039 199 240.59 8.154e-07 ***
## treatmentGroup 1 33.012 198 207.57 9.159e-09 ***
## desert:treatmentGroup 2 0.421 196 207.15 0.8102
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## NOTE: Results may be misleading due to involvement in interactions## $emmeans
## treatmentGroup emmean SE df z.ratio p.value
## never burned 1.73 0.086 Inf 20.162 <.0001
## burned 2.52 0.113 Inf 22.363 <.0001
##
## Results are averaged over the levels of: desert
## Results are given on the log (not the response) scale.
##
## $contrasts
## contrast estimate SE df z.ratio p.value
## never burned - burned -0.788 0.142 Inf -5.557 <.0001
##
## Results are averaged over the levels of: desert
## Results are given on the log (not the response) scale.## $emmeans
## desert = Mojave:
## treatmentGroup emmean SE df z.ratio p.value
## never burned 1.15 0.152 Inf 7.573 <.0001
## burned 2.00 0.169 Inf 11.870 <.0001
##
## desert = San Joaquin:
## treatmentGroup emmean SE df z.ratio p.value
## never burned 2.15 0.144 Inf 14.903 <.0001
## burned 2.80 0.210 Inf 13.355 <.0001
##
## desert = Sonoran:
## treatmentGroup emmean SE df z.ratio p.value
## never burned 1.90 0.150 Inf 12.646 <.0001
## burned 2.76 0.205 Inf 13.477 <.0001
##
## Results are given on the log (not the response) scale.
##
## $contrasts
## desert = Mojave:
## contrast estimate SE df z.ratio p.value
## never burned - burned -0.849 0.227 Inf -3.733 0.0002
##
## desert = San Joaquin:
## contrast estimate SE df z.ratio p.value
## never burned - burned -0.655 0.254 Inf -2.573 0.0101
##
## desert = Sonoran:
## contrast estimate SE df z.ratio p.value
## never burned - burned -0.860 0.254 Inf -3.388 0.0007
##
## Results are given on the log (not the response) scale.## # Overdispersion test
##
## dispersion ratio = 0.983
## Pearson's Chi-Squared = 118.017
## p-value = 0.534## No overdispersion detected.##
## Call:
## glm.nb(formula = totalORPer1000km2 ~ desert * treatmentGroup,
## data = filter(test2, ez_taxa == "Aves"), init.theta = 1.842130059,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.7966 -1.0985 -0.4326 0.5600 1.8207
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.5702 0.1699 9.240 < 2e-16 ***
## desertSan Joaquin 1.0181 0.2295 4.437 9.11e-06 ***
## desertSonoran 0.6313 0.2324 2.716 0.0066 **
## treatmentGroupburned 0.5988 0.2372 2.525 0.0116 *
## desertSan Joaquin:treatmentGroupburned -0.2527 0.3425 -0.738 0.4607
## desertSonoran:treatmentGroupburned -0.5360 0.4028 -1.331 0.1833
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.8421) family taken to be 1)
##
## Null deviance: 163.68 on 125 degrees of freedom
## Residual deviance: 129.87 on 120 degrees of freedom
## AIC: 827.71
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.842
## Std. Err.: 0.266
##
## 2 x log-likelihood: -813.708## Warning in anova.negbin(occA_glm, test = "Chisq"): tests made without re-
## estimating 'theta'## Analysis of Deviance Table
##
## Model: Negative Binomial(1.8421), link: log
##
## Response: totalORPer1000km2
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 125 163.68
## desert 2 25.3953 123 138.28 3.058e-06 ***
## treatmentGroup 1 6.6375 122 131.65 0.009985 **
## desert:treatmentGroup 2 1.7767 120 129.87 0.411329
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## NOTE: Results may be misleading due to involvement in interactions## $emmeans
## treatmentGroup emmean SE df z.ratio p.value
## never burned 2.12 0.093 Inf 22.805 <.0001
## burned 2.46 0.127 Inf 19.318 <.0001
##
## Results are averaged over the levels of: desert
## Results are given on the log (not the response) scale.
##
## $contrasts
## contrast estimate SE df z.ratio p.value
## never burned - burned -0.336 0.158 Inf -2.133 0.0329
##
## Results are averaged over the levels of: desert
## Results are given on the log (not the response) scale.## NOTE: Results may be misleading due to involvement in interactions## $emmeans
## desert emmean SE df z.ratio p.value
## Mojave 1.87 0.119 Inf 15.765 <.0001
## San Joaquin 2.76 0.124 Inf 22.358 <.0001
## Sonoran 2.23 0.163 Inf 13.718 <.0001
##
## Results are averaged over the levels of: treatmentGroup
## Results are given on the log (not the response) scale.
##
## $contrasts
## contrast estimate SE df z.ratio p.value
## Mojave - San Joaquin -0.892 0.171 Inf -5.208 <.0001
## Mojave - Sonoran -0.363 0.201 Inf -1.804 0.1683
## San Joaquin - Sonoran 0.529 0.204 Inf 2.587 0.0262
##
## Results are averaged over the levels of: treatmentGroup
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimates## $emmeans
## desert = Mojave:
## treatmentGroup emmean SE df z.ratio p.value
## never burned 1.57 0.170 Inf 9.240 <.0001
## burned 2.17 0.165 Inf 13.108 <.0001
##
## desert = San Joaquin:
## treatmentGroup emmean SE df z.ratio p.value
## never burned 2.59 0.154 Inf 16.789 <.0001
## burned 2.93 0.193 Inf 15.204 <.0001
##
## desert = Sonoran:
## treatmentGroup emmean SE df z.ratio p.value
## never burned 2.20 0.159 Inf 13.886 <.0001
## burned 2.26 0.284 Inf 7.964 <.0001
##
## Results are given on the log (not the response) scale.
##
## $contrasts
## desert = Mojave:
## contrast estimate SE df z.ratio p.value
## never burned - burned -0.5988 0.237 Inf -2.525 0.0116
##
## desert = San Joaquin:
## contrast estimate SE df z.ratio p.value
## never burned - burned -0.3462 0.247 Inf -1.401 0.1611
##
## desert = Sonoran:
## contrast estimate SE df z.ratio p.value
## never burned - burned -0.0629 0.326 Inf -0.193 0.8468
##
## Results are given on the log (not the response) scale.## $emmeans
## desert = Mojave:
## treatmentGroup response SE df asymp.LCL asymp.UCL
## never burned 4.81 0.817 Inf 3.45 6.71
## burned 8.75 1.448 Inf 6.33 12.10
##
## desert = San Joaquin:
## treatmentGroup response SE df asymp.LCL asymp.UCL
## never burned 13.31 2.052 Inf 9.84 18.00
## burned 18.81 3.631 Inf 12.89 27.46
##
## desert = Sonoran:
## treatmentGroup response SE df asymp.LCL asymp.UCL
## never burned 9.04 1.433 Inf 6.62 12.33
## burned 9.62 2.737 Inf 5.51 16.80
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## desert = Mojave:
## contrast ratio SE df asymp.LCL asymp.UCL
## never burned / burned 0.549 0.130 Inf 0.345 0.875
##
## desert = San Joaquin:
## contrast ratio SE df asymp.LCL asymp.UCL
## never burned / burned 0.707 0.175 Inf 0.436 1.148
##
## desert = Sonoran:
## contrast ratio SE df asymp.LCL asymp.UCL
## never burned / burned 0.939 0.306 Inf 0.496 1.777
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale## # Overdispersion test
##
## dispersion ratio = 0.767
## Pearson's Chi-Squared = 53.668
## p-value = 0.926## No overdispersion detected.## Warning in anova.negbin(occO_glm, test = "Chisq"): tests made without re-
## estimating 'theta'## Analysis of Deviance Table
##
## Model: Negative Binomial(3.1901), link: log
##
## Response: totalORPer1000km2
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 75 223.651
## desert 2 103.633 73 120.019 < 2.2e-16 ***
## treatmentGroup 1 65.048 72 54.971 7.309e-16 ***
## desert:treatmentGroup 2 5.363 70 49.607 0.06846 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## $emmeans
## desert = Mojave:
## treatmentGroup emmean SE df z.ratio p.value
## never burned 0.134 0.238 Inf 0.561 0.5746
## burned 1.216 0.276 Inf 4.406 <.0001
##
## desert = San Joaquin:
## treatmentGroup emmean SE df z.ratio p.value
## never burned 0.258 0.253 Inf 1.020 0.3077
## burned 1.386 0.433 Inf 3.199 0.0014
##
## desert = Sonoran:
## treatmentGroup emmean SE df z.ratio p.value
## never burned 0.956 0.216 Inf 4.429 <.0001
## burned 2.992 0.174 Inf 17.184 <.0001
##
## Results are given on the log (not the response) scale.
##
## $contrasts
## desert = Mojave:
## contrast estimate SE df z.ratio p.value
## never burned - burned -1.08 0.364 Inf -2.971 0.0030
##
## desert = San Joaquin:
## contrast estimate SE df z.ratio p.value
## never burned - burned -1.13 0.502 Inf -2.249 0.0245
##
## desert = Sonoran:
## contrast estimate SE df z.ratio p.value
## never burned - burned -2.04 0.277 Inf -7.345 <.0001
##
## Results are given on the log (not the response) scale.## $emmeans
## desert = Mojave:
## treatmentGroup response SE df asymp.LCL asymp.UCL
## never burned 1.14 0.272 Inf 0.717 1.82
## burned 3.38 0.932 Inf 1.965 5.80
##
## desert = San Joaquin:
## treatmentGroup response SE df asymp.LCL asymp.UCL
## never burned 1.29 0.327 Inf 0.789 2.12
## burned 4.00 1.734 Inf 1.711 9.35
##
## desert = Sonoran:
## treatmentGroup response SE df asymp.LCL asymp.UCL
## never burned 2.60 0.561 Inf 1.704 3.97
## burned 19.92 3.467 Inf 14.159 28.02
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## desert = Mojave:
## contrast ratio SE df asymp.LCL asymp.UCL
## never burned / burned 0.339 0.1234 Inf 0.1658 0.692
##
## desert = San Joaquin:
## contrast ratio SE df asymp.LCL asymp.UCL
## never burned / burned 0.324 0.1623 Inf 0.1210 0.865
##
## desert = Sonoran:
## contrast ratio SE df asymp.LCL asymp.UCL
## never burned / burned 0.131 0.0362 Inf 0.0758 0.225
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale## Warning in plot_grid(p2, p1a, cols = 2, labels = "auto", align = "h", axis =
## "tb", : Argument 'cols' is deprecated. Use 'ncol' instead.
occurrences annually time-series
## Warning: package 'car' was built under R version 4.2.3## Loading required package: carData## Warning: package 'carData' was built under R version 4.2.3##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recode## The following object is masked from 'package:purrr':
##
## some## `summarise()` has grouped output by 'treatmentGroup', 'obsYear', 'desert'. You
## can override using the `.groups` argument.## Warning in (function (start, objective, gradient = NULL, hessian = NULL, : NA/
## NaN function evaluation## Warning in (function (start, objective, gradient = NULL, hessian = NULL, : NA/
## NaN function evaluation## Family: nbinom1 ( log )
## Formula: stzd_Occ ~ obsYear + treatmentGroup + (1 | desert)
## Data: data
##
## AIC BIC logLik deviance df.resid
## 1240.5 1257.0 -615.2 1230.5 197
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## desert (Intercept) 0.02697 0.1642
## Number of obs: 202, groups: desert, 3
##
## Dispersion parameter for nbinom1 family (): 3.67
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.361e+02 1.595e+01 -8.529 < 2e-16 ***
## obsYear 6.869e-02 7.928e-03 8.664 < 2e-16 ***
## treatmentGroupburned 3.458e-01 9.499e-02 3.640 0.000273 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## # Overdispersion test
##
## dispersion ratio = 1.063
## Pearson's Chi-Squared = 209.411
## p-value = 0.259## No overdispersion detected.## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: stzd_Occ
## Chisq Df Pr(>Chisq)
## (Intercept) 72.751 1 < 2.2e-16 ***
## obsYear 75.063 1 < 2.2e-16 ***
## treatmentGroup 13.249 1 0.0002728 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## $emmeans
## treatmentGroup emmean SE df lower.CL upper.CL
## never burned 1.97 0.118 197 1.74 2.20
## burned 2.31 0.123 197 2.07 2.56
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## never burned - burned -0.346 0.095 197 -3.640 0.0003
##
## Results are given on the log (not the response) scale.## Warning in (function (start, objective, gradient = NULL, hessian = NULL, : NA/
## NaN function evaluation## Warning in (function (start, objective, gradient = NULL, hessian = NULL, : NA/
## NaN function evaluation## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')## Warning in sqrt(diag(vcovs)): NaNs produced## Family: nbinom1 ( log )
## Formula: stzd_Occ ~ obsYear + treatmentGroup + (1 | desert)
## Data: filter(data, ez_taxa == "Aves")
##
## AIC BIC logLik deviance df.resid
## NA NA NA NA 121
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## desert (Intercept) 2.12e-10 1.456e-05
## Number of obs: 126, groups: desert, 3
##
## Dispersion parameter for nbinom1 family (): 1.59
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.813e+02 3.500e-02 -5179 <2e-16 ***
## obsYear 9.135e-02 NaN NaN NaN
## treatmentGroupburned 1.514e-03 8.524e-02 0 0.986
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## # Overdispersion test
##
## dispersion ratio = 1.161
## Pearson's Chi-Squared = 140.465
## p-value = 0.109## No overdispersion detected.## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: stzd_Occ
## Chisq Df Pr(>Chisq)
## (Intercept) 2.6827e+07 1 <2e-16 ***
## obsYear -3.1735e+07 1 1.0000
## treatmentGroup 3.0000e-04 1 0.9858
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## $emmeans
## treatmentGroup emmean SE df lower.CL upper.CL
## never burned 2.2 0.0549 121 2.09 2.31
## burned 2.2 0.0633 121 2.08 2.33
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## never burned - burned -0.00151 0.0852 121 -0.018 0.9859
##
## Results are given on the log (not the response) scale.## Warning in (function (start, objective, gradient = NULL, hessian = NULL, : NA/
## NaN function evaluation## Warning in (function (start, objective, gradient = NULL, hessian = NULL, : NA/
## NaN function evaluation## Family: nbinom1 ( log )
## Formula: stzd_Occ ~ obsYear + treatmentGroup + (1 | desert)
## Data: filter(data, ez_taxa == "Other")
##
## AIC BIC logLik deviance df.resid
## 386.8 398.4 -188.4 376.8 71
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## desert (Intercept) 0.1967 0.4435
## Number of obs: 76, groups: desert, 3
##
## Dispersion parameter for nbinom1 family (): 0.848
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -238.8409 32.4675 -7.356 1.89e-13 ***
## obsYear 0.1192 0.0161 7.404 1.32e-13 ***
## treatmentGroupburned 0.7832 0.1346 5.819 5.91e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## # Overdispersion test
##
## dispersion ratio = 1.233
## Pearson's Chi-Squared = 87.523
## p-value = 0.089## No overdispersion detected.## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: stzd_Occ
## Chisq Df Pr(>Chisq)
## (Intercept) 54.115 1 1.891e-13 ***
## obsYear 54.820 1 1.321e-13 ***
## treatmentGroup 33.863 1 5.913e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## $emmeans
## treatmentGroup emmean SE df lower.CL upper.CL
## never burned 1.06 0.286 71 0.488 1.63
## burned 1.84 0.280 71 1.283 2.40
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## never burned - burned -0.783 0.135 71 -5.819 <.0001
##
## Results are given on the log (not the response) scale.## Family: nbinom2 ( log )
## Formula: stzd_Occ ~ obsYear + treatmentGroup + ez_taxa + (1 | desert)
## Data: data
##
## AIC BIC logLik deviance df.resid
## 1078.1 1098.0 -533.1 1066.1 196
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## desert (Intercept) 0.06196 0.2489
## Number of obs: 202, groups: desert, 3
##
## Dispersion parameter for nbinom2 family (): 5.62
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.122e+02 1.816e+01 -11.683 < 2e-16 ***
## obsYear 1.066e-01 9.037e-03 11.801 < 2e-16 ***
## treatmentGroupburned 4.275e-01 8.498e-02 5.031 4.87e-07 ***
## ez_taxaOther -1.044e+00 9.602e-02 -10.875 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## # Overdispersion test
##
## dispersion ratio = 0.993
## Pearson's Chi-Squared = 194.550
## p-value = 0.516## No overdispersion detected.## $emmeans
## ez_taxa = Aves:
## treatmentGroup emmean SE df lower.CL upper.CL
## never burned 2.12 0.156 196 1.817 2.43
## burned 2.55 0.161 196 2.233 2.87
##
## ez_taxa = Other:
## treatmentGroup emmean SE df lower.CL upper.CL
## never burned 1.08 0.169 196 0.747 1.41
## burned 1.51 0.174 196 1.165 1.85
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## ez_taxa = Aves:
## contrast estimate SE df t.ratio p.value
## never burned - burned -0.428 0.085 196 -5.031 <.0001
##
## ez_taxa = Other:
## contrast estimate SE df t.ratio p.value
## never burned - burned -0.428 0.085 196 -5.031 <.0001
##
## Results are given on the log (not the response) scale.Recovery
Occurrences before and after fire
GLMM model w/ emmeans - pre/post
all species
library(car)
data_pp_5yrAS <- data_pp %>%
group_by(desert, fireYear, pre_or_post, funcGroup, treatmentGroup, fireID) %>%
summarize(mean_occurrence = ceiling(mean(stzdObs))) %>%
mutate(ez_taxa = ifelse(grepl("Birds", funcGroup), "Aves", "Other"))## `summarise()` has grouped output by 'desert', 'fireYear', 'pre_or_post',
## 'funcGroup', 'treatmentGroup'. You can override using the `.groups` argument.pp_nbAS <- glmmTMB(mean_occurrence ~ pre_or_post*treatmentGroup + (1|ez_taxa) + (1|desert), family = nbinom1, data = data_pp_5yrAS)
check_overdispersion(pp_nbAS)## # Overdispersion test
##
## dispersion ratio = 0.904
## Pearson's Chi-Squared = 456.384
## p-value = 0.941## No overdispersion detected.#validate model using simulated residuals
simulationOutput <- simulateResiduals(fittedModel = pp_nbAS, plot = F)
#plot(simulationOutput)
# testDispersion -- A significant ratio > 1 indicates overdispersion, a significant ratio < 1 underdispersion.
testDispersion(simulationOutput, plot = FALSE)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.7186, p-value = 0.256
## alternative hypothesis: two.sided# testZeroInflation -- A value < 1 means that the observed data has less zeros than expected, a value > 1 means that it has more zeros than expected (aka zero-inflation)
testZeroInflation(simulationOutput, plot = FALSE) ##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0.89343, p-value = 0.312
## alternative hypothesis: two.sidedAnova(pp_nbAS, type = "III")## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: mean_occurrence
## Chisq Df Pr(>Chisq)
## (Intercept) 0.0003 1 0.986008
## pre_or_post 8.5289 1 0.003495 **
## treatmentGroup 159.2265 1 < 2.2e-16 ***
## pre_or_post:treatmentGroup 3.5436 1 0.059777 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#stepwise regression to select best predictors; both summary above and stepAIC show that funcGroup not significant
null <- glmmTMB(mean_occurrence ~ 1, family = nbinom1, data = data_pp_5yrAS)
stepAIC(pp_nbAS, scope=list(lower=null, upper=pp_nbAS), data=data_pp_5yrAS, direction='backward')## Start: AIC=1725.43
## mean_occurrence ~ pre_or_post * treatmentGroup
##
## Df AIC
## <none> 1725.4
## - pre_or_post:treatmentGroup 3 1783.1## Formula: mean_occurrence ~ pre_or_post * treatmentGroup
## Data: data_pp_5yrAS
## AIC BIC logLik df.resid
## 1725.4307 1755.0990 -855.7154 505
## Random-effects (co)variances:
##
## Conditional model:
## Groups Name Std.Dev.
## ez_taxa (Intercept) 0.6060
## desert (Intercept) 0.3236
##
## Number of obs: 512 / Conditional model: ez_taxa, 2; desert, 3
##
## Dispersion parameter for nbinom1 family (): 25.4
##
## Fixed Effects:
##
## Conditional model:
## (Intercept) pre_or_postpre
## -0.009194 -1.191218
## treatmentGroupcontrol pre_or_postpre:treatmentGroupcontrol
## 2.885037 0.807334emm_pp <- emmeans(pp_nbAS, pairwise ~ pre_or_post|treatmentGroup)
test(emm_pp)## $emmeans
## treatmentGroup = burned:
## pre_or_post emmean SE df t.ratio p.value
## post -0.00919 0.524 505 -0.018 0.9860
## pre -1.20041 0.600 505 -2.000 0.0460
##
## treatmentGroup = control:
## pre_or_post emmean SE df t.ratio p.value
## post 2.87584 0.480 505 5.996 <.0001
## pre 2.49196 0.486 505 5.125 <.0001
##
## Results are given on the log (not the response) scale.
##
## $contrasts
## treatmentGroup = burned:
## contrast estimate SE df t.ratio p.value
## post - pre 1.191 0.408 505 2.920 0.0037
##
## treatmentGroup = control:
## contrast estimate SE df t.ratio p.value
## post - pre 0.384 0.134 505 2.871 0.0043
##
## Results are given on the log (not the response) scale.confint(emm_pp, level = .95, type = "response")## $emmeans
## treatmentGroup = burned:
## pre_or_post response SE df lower.CL upper.CL
## post 0.991 0.519 505 0.3537 2.776
## pre 0.301 0.181 505 0.0926 0.979
##
## treatmentGroup = control:
## pre_or_post response SE df lower.CL upper.CL
## post 17.740 8.509 505 6.9138 45.521
## pre 12.085 5.876 505 4.6491 31.413
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## treatmentGroup = burned:
## contrast ratio SE df lower.CL upper.CL
## post / pre 3.29 1.342 505 1.48 7.33
##
## treatmentGroup = control:
## contrast ratio SE df lower.CL upper.CL
## post / pre 1.47 0.196 505 1.13 1.91
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scaleemm_pp <- emmeans(pp_nbAS, pairwise ~ treatmentGroup|pre_or_post)
test(emm_pp)## $emmeans
## pre_or_post = post:
## treatmentGroup emmean SE df t.ratio p.value
## burned -0.00919 0.524 505 -0.018 0.9860
## control 2.87584 0.480 505 5.996 <.0001
##
## pre_or_post = pre:
## treatmentGroup emmean SE df t.ratio p.value
## burned -1.20041 0.600 505 -2.000 0.0460
## control 2.49196 0.486 505 5.125 <.0001
##
## Results are given on the log (not the response) scale.
##
## $contrasts
## pre_or_post = post:
## contrast estimate SE df t.ratio p.value
## burned - control -2.89 0.229 505 -12.618 <.0001
##
## pre_or_post = pre:
## contrast estimate SE df t.ratio p.value
## burned - control -3.69 0.370 505 -9.980 <.0001
##
## Results are given on the log (not the response) scale.non-avian species
data_pp_5yrNA <- filter(data_pp, !grepl("Birds", funcGroup)) %>%
group_by(desert, fireYear, pre_or_post, funcGroup, treatmentGroup, fireID) %>%
summarize(mean_occurrence = ceiling(mean(stzdObs)))## `summarise()` has grouped output by 'desert', 'fireYear', 'pre_or_post',
## 'funcGroup', 'treatmentGroup'. You can override using the `.groups` argument.pp_nbNA <- glmmTMB(mean_occurrence ~ pre_or_post*treatmentGroup + (1|funcGroup) + (1|desert), family = nbinom2, data = data_pp_5yrNA)
check_overdispersion(pp_nbNA)## # Overdispersion test
##
## dispersion ratio = 2.820
## Pearson's Chi-Squared = 659.835
## p-value = < 0.001## Overdispersion detected.#validate model using simulated residuals
simulationOutput <- simulateResiduals(fittedModel = pp_nbNA, plot = F)
#plot(simulationOutput)
testDispersion(simulationOutput, plot = FALSE)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.34175, p-value = 0.744
## alternative hypothesis: two.sided# testZeroInflation -- A value < 1 means that the observed data has less zeros than expected, a value > 1 means that it has more zeros than expected (aka zero-inflation)
testZeroInflation(simulationOutput, plot = FALSE) ##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0.9791, p-value = 0.792
## alternative hypothesis: two.sidedAnova(pp_nbNA, type = "III")## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: mean_occurrence
## Chisq Df Pr(>Chisq)
## (Intercept) 2.3630 1 0.1242438
## pre_or_post 12.1845 1 0.0004819 ***
## treatmentGroup 28.1024 1 1.151e-07 ***
## pre_or_post:treatmentGroup 1.7508 1 0.1857812
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#stepwise regression to select best predictors; both summary above and stepAIC show that funcGroup not significant
null <- glmmTMB(mean_occurrence ~ 1, family = nbinom1, data = data_pp_5yrNA)
stepAIC(pp_nbNA, scope=list(lower=null, upper=pp_nbNA), data=data_pp_5yr_aves, direction='backward')## Start: AIC=535.04
## mean_occurrence ~ pre_or_post * treatmentGroup
##
## Df AIC
## <none> 535.04
## - pre_or_post:treatmentGroup 3 558.93## Formula: mean_occurrence ~ pre_or_post * treatmentGroup
## Data: data_pp_5yrNA
## AIC BIC logLik df.resid
## 535.0389 559.4325 -260.5195 234
## Random-effects (co)variances:
##
## Conditional model:
## Groups Name Std.Dev.
## funcGroup (Intercept) 0.6404
## desert (Intercept) 0.4465
##
## Number of obs: 241 / Conditional model: funcGroup, 4; desert, 3
##
## Dispersion parameter for nbinom2 family (): 0.542
##
## Fixed Effects:
##
## Conditional model:
## (Intercept) pre_or_postpre
## -0.7216 -1.2967
## treatmentGroupcontrol pre_or_postpre:treatmentGroupcontrol
## 1.8194 0.7194emm_pp <- emmeans(pp_nbNA, pairwise ~ pre_or_post|treatmentGroup)
test(emm_pp)## $emmeans
## treatmentGroup = burned:
## pre_or_post emmean SE df t.ratio p.value
## post -0.722 0.469 234 -1.537 0.1256
## pre -2.018 0.523 234 -3.861 0.0001
##
## treatmentGroup = control:
## pre_or_post emmean SE df t.ratio p.value
## post 1.098 0.477 234 2.299 0.0224
## pre 0.520 0.524 234 0.994 0.3213
##
## Results are given on the log (not the response) scale.
##
## $contrasts
## treatmentGroup = burned:
## contrast estimate SE df t.ratio p.value
## post - pre 1.297 0.371 234 3.491 0.0006
##
## treatmentGroup = control:
## contrast estimate SE df t.ratio p.value
## post - pre 0.577 0.394 234 1.465 0.1441
##
## Results are given on the log (not the response) scale.confint(emm_pp, level = .95, type = "response")## $emmeans
## treatmentGroup = burned:
## pre_or_post response SE df lower.CL upper.CL
## post 0.486 0.2281 234 0.1927 1.225
## pre 0.133 0.0695 234 0.0474 0.372
##
## treatmentGroup = control:
## pre_or_post response SE df lower.CL upper.CL
## post 2.998 1.4311 234 1.1702 7.678
## pre 1.683 0.8813 234 0.5997 4.722
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## treatmentGroup = burned:
## contrast ratio SE df lower.CL upper.CL
## post / pre 3.66 1.359 234 1.76 7.60
##
## treatmentGroup = control:
## contrast ratio SE df lower.CL upper.CL
## post / pre 1.78 0.702 234 0.82 3.87
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scaleemm_pp <- emmeans(pp_nbNA, pairwise ~ treatmentGroup|pre_or_post)
summary(emm_pp)## $emmeans
## pre_or_post = post:
## treatmentGroup emmean SE df lower.CL upper.CL
## burned -0.722 0.469 234 -1.646 0.203
## control 1.098 0.477 234 0.157 2.038
##
## pre_or_post = pre:
## treatmentGroup emmean SE df lower.CL upper.CL
## burned -2.018 0.523 234 -3.048 -0.988
## control 0.520 0.524 234 -0.511 1.552
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## pre_or_post = post:
## contrast estimate SE df t.ratio p.value
## burned - control -1.82 0.343 234 -5.301 <.0001
##
## pre_or_post = pre:
## contrast estimate SE df t.ratio p.value
## burned - control -2.54 0.476 234 -5.328 <.0001
##
## Results are given on the log (not the response) scale.all aves
data_pp_5yrA <- filter(data_pp, grepl("Birds", funcGroup)) %>%
group_by(desert, fireYear, pre_or_post, funcGroup, treatmentGroup, fireID) %>%
summarize(mean_occurrence = ceiling(mean(stzdObs)))## `summarise()` has grouped output by 'desert', 'fireYear', 'pre_or_post',
## 'funcGroup', 'treatmentGroup'. You can override using the `.groups` argument.pp_nbA <- glmmTMB(mean_occurrence ~ pre_or_post*treatmentGroup + (1|funcGroup) + (1|desert), family = nbinom1, data = data_pp_5yrA)
check_overdispersion(pp_nbA)## # Overdispersion test
##
## dispersion ratio = 0.763
## Pearson's Chi-Squared = 201.536
## p-value = 0.998## No overdispersion detected.#validate model using simulated residuals
simulationOutput <- simulateResiduals(fittedModel = pp_nbA, plot = F)
#plot(simulationOutput)
testDispersion(simulationOutput, plot = FALSE)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.3302, p-value = 0.344
## alternative hypothesis: two.sided# testZeroInflation -- A value < 1 means that the observed data has less zeros than expected, a value > 1 means that it has more zeros than expected (aka zero-inflation)
testZeroInflation(simulationOutput, plot = FALSE)##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0.91144, p-value = 0.328
## alternative hypothesis: two.sided#stepwise regression to select best predictors; both summary above and stepAIC show that funcGroup not significant
null <- glmmTMB(mean_occurrence ~ 1, family = nbinom1, data = data_pp_5yrA)
stepAIC(pp_nbA, scope=list(lower=null, upper=pp_nbA), data=data_pp_5yrA, direction='backward')## Start: AIC=1178.5
## mean_occurrence ~ pre_or_post * treatmentGroup
##
## Df AIC
## <none> 1178.5
## - pre_or_post:treatmentGroup 3 1221.6## Formula: mean_occurrence ~ pre_or_post * treatmentGroup
## Data: data_pp_5yrA
## AIC BIC logLik df.resid
## 1178.5029 1203.7177 -582.2515 264
## Random-effects (co)variances:
##
## Conditional model:
## Groups Name Std.Dev.
## funcGroup (Intercept) 0.6380
## desert (Intercept) 0.3575
##
## Number of obs: 271 / Conditional model: funcGroup, 4; desert, 3
##
## Dispersion parameter for nbinom1 family (): 26.2
##
## Fixed Effects:
##
## Conditional model:
## (Intercept) pre_or_postpre
## 0.4445 -1.0522
## treatmentGroupcontrol pre_or_postpre:treatmentGroupcontrol
## 2.8554 0.4898Anova(pp_nbA, type = "III")## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: mean_occurrence
## Chisq Df Pr(>Chisq)
## (Intercept) 0.9027 1 0.34206
## pre_or_post 5.5162 1 0.01884 *
## treatmentGroup 112.1569 1 < 2e-16 ***
## pre_or_post:treatmentGroup 1.0745 1 0.29993
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1emm_pp <- emmeans(pp_nbA, pairwise ~ pre_or_post|treatmentGroup)
test(emm_pp)## $emmeans
## treatmentGroup = burned:
## pre_or_post emmean SE df t.ratio p.value
## post 0.444 0.468 264 0.950 0.3429
## pre -0.608 0.555 264 -1.095 0.2744
##
## treatmentGroup = control:
## pre_or_post emmean SE df t.ratio p.value
## post 3.300 0.401 264 8.225 <.0001
## pre 2.737 0.412 264 6.648 <.0001
##
## Results are given on the log (not the response) scale.
##
## $contrasts
## treatmentGroup = burned:
## contrast estimate SE df t.ratio p.value
## post - pre 1.052 0.448 264 2.349 0.0196
##
## treatmentGroup = control:
## contrast estimate SE df t.ratio p.value
## post - pre 0.562 0.153 264 3.673 0.0003
##
## Results are given on the log (not the response) scale.confint(emm_pp, level = 0.95, type = "response")## $emmeans
## treatmentGroup = burned:
## pre_or_post response SE df lower.CL upper.CL
## post 1.560 0.730 264 0.621 3.92
## pre 0.545 0.302 264 0.183 1.62
##
## treatmentGroup = control:
## pre_or_post response SE df lower.CL upper.CL
## post 27.108 10.875 264 12.304 59.73
## pre 15.448 6.361 264 6.867 34.75
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## treatmentGroup = burned:
## contrast ratio SE df lower.CL upper.CL
## post / pre 2.86 1.283 264 1.19 6.92
##
## treatmentGroup = control:
## contrast ratio SE df lower.CL upper.CL
## post / pre 1.75 0.269 264 1.30 2.37
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scaleemm_pp <- emmeans(pp_nbA, pairwise ~ treatmentGroup|pre_or_post)
test(emm_pp)## $emmeans
## pre_or_post = post:
## treatmentGroup emmean SE df t.ratio p.value
## burned 0.444 0.468 264 0.950 0.3429
## control 3.300 0.401 264 8.225 <.0001
##
## pre_or_post = pre:
## treatmentGroup emmean SE df t.ratio p.value
## burned -0.608 0.555 264 -1.095 0.2744
## control 2.737 0.412 264 6.648 <.0001
##
## Results are given on the log (not the response) scale.
##
## $contrasts
## pre_or_post = post:
## contrast estimate SE df t.ratio p.value
## burned - control -2.86 0.27 264 -10.590 <.0001
##
## pre_or_post = pre:
## contrast estimate SE df t.ratio p.value
## burned - control -3.35 0.40 264 -8.365 <.0001
##
## Results are given on the log (not the response) scale.passerine and raptor only
data_pp_5yrPR <- filter(data_pp, funcGroup == "Birds:Passerines" | funcGroup == "Birds:Raptors") %>%
group_by(desert, fireYear, pre_or_post, funcGroup, treatmentGroup, fireID) %>%
summarize(mean_occurrence = ceiling(mean(stzdObs)))## `summarise()` has grouped output by 'desert', 'fireYear', 'pre_or_post',
## 'funcGroup', 'treatmentGroup'. You can override using the `.groups` argument.pp_nbPR <- glmmTMB(mean_occurrence ~ pre_or_post * treatmentGroup + (1|funcGroup) + (1|desert), family = nbinom1, data = data_pp_5yrPR)
check_overdispersion(pp_nbPR)## # Overdispersion test
##
## dispersion ratio = 0.882
## Pearson's Chi-Squared = 131.367
## p-value = 0.848## No overdispersion detected.#validate model using simulated residuals
simulationOutput <- simulateResiduals(fittedModel = pp_nbPR, plot = F)
#plot(simulationOutput)
testDispersion(simulationOutput, plot = FALSE)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.4606, p-value = 0.272
## alternative hypothesis: two.sided# testZeroInflation -- A value < 1 means that the observed data has less zeros than expected, a value > 1 means that it has more zeros than expected (aka zero-inflation)
testZeroInflation(simulationOutput, plot = FALSE)##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0.9572, p-value = 0.8
## alternative hypothesis: two.sidedsummary(pp_nbPR)## Family: nbinom1 ( log )
## Formula:
## mean_occurrence ~ pre_or_post * treatmentGroup + (1 | funcGroup) +
## (1 | desert)
## Data: data_pp_5yrPR
##
## AIC BIC logLik deviance df.resid
## 913.5 934.9 -449.8 899.5 149
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## funcGroup (Intercept) 0.0135 0.1162
## desert (Intercept) 0.1208 0.3475
## Number of obs: 156, groups: funcGroup, 2; desert, 3
##
## Dispersion parameter for nbinom1 family (): 26.2
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.3059 0.3497 3.735 0.000188 ***
## pre_or_postpre -1.2065 0.4764 -2.533 0.011324 *
## treatmentGroupcontrol 2.5211 0.2773 9.092 < 2e-16 ***
## pre_or_postpre:treatmentGroupcontrol 0.5589 0.5053 1.106 0.268621
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1emm_pp <- emmeans(pp_nbPR, pairwise ~ pre_or_post|treatmentGroup)
summary(emm_pp)## $emmeans
## treatmentGroup = burned:
## pre_or_post emmean SE df lower.CL upper.CL
## post 1.3059 0.350 149 0.615 2.00
## pre 0.0994 0.483 149 -0.855 1.05
##
## treatmentGroup = control:
## pre_or_post emmean SE df lower.CL upper.CL
## post 3.8270 0.248 149 3.337 4.32
## pre 3.1794 0.267 149 2.653 3.71
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## treatmentGroup = burned:
## contrast estimate SE df t.ratio p.value
## post - pre 1.206 0.476 149 2.533 0.0124
##
## treatmentGroup = control:
## contrast estimate SE df t.ratio p.value
## post - pre 0.648 0.172 149 3.772 0.0002
##
## Results are given on the log (not the response) scale.emm_pp <- emmeans(pp_nbPR, pairwise ~ treatmentGroup|pre_or_post)
summary(emm_pp)## $emmeans
## pre_or_post = post:
## treatmentGroup emmean SE df lower.CL upper.CL
## burned 1.3059 0.350 149 0.615 2.00
## control 3.8270 0.248 149 3.337 4.32
##
## pre_or_post = pre:
## treatmentGroup emmean SE df lower.CL upper.CL
## burned 0.0994 0.483 149 -0.855 1.05
## control 3.1794 0.267 149 2.653 3.71
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## pre_or_post = post:
## contrast estimate SE df t.ratio p.value
## burned - control -2.52 0.277 149 -9.092 <.0001
##
## pre_or_post = pre:
## contrast estimate SE df t.ratio p.value
## burned - control -3.08 0.433 149 -7.117 <.0001
##
## Results are given on the log (not the response) scale.graphs
Overall occurrences of raptors and passerines between 1995-2020
## `summarise()` has grouped output by 'desert', 'pre_or_post', 'treatmentGroup'.
## You can override using the `.groups` argument.
## `summarise()` has grouped output by 'pre_or_post', 'treatmentGroup'. You can
## override using the `.groups` argument.
exclude passerines and raptors
The average number of occurrences reported 5 years before fire and 5 years after fire. Zeros were added to years when no occurrences were reported.
Recovery over-time - SpadeR
Chao Two Community Similarity
shapiro.test(df_tsf$SI)##
## Shapiro-Wilk normality test
##
## data: df_tsf$SI
## W = 0.91556, p-value = 0.1897leveneTest(SI ~ desert, data = df_tsf)## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 0.4745 0.6344
## 11sim_mod <- lm(SI ~ yearsSinceFire + desert, data = df_tsf)
summary(sim_mod)##
## Call:
## lm(formula = SI ~ yearsSinceFire + desert, data = df_tsf)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.40340 -0.07529 0.01818 0.07235 0.29812
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.02284 0.16406 6.235 9.7e-05 ***
## yearsSinceFire -0.01069 0.01048 -1.020 0.33176
## desertSan Joaquin -0.14988 0.16216 -0.924 0.37712
## desertSonoran -0.57406 0.17301 -3.318 0.00777 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1994 on 10 degrees of freedom
## Multiple R-squared: 0.6153, Adjusted R-squared: 0.4999
## F-statistic: 5.332 on 3 and 10 DF, p-value: 0.01878anova(sim_mod)## Analysis of Variance Table
##
## Response: SI
## Df Sum Sq Mean Sq F value Pr(>F)
## yearsSinceFire 1 0.01730 0.017301 0.4351 0.524397
## desert 2 0.61877 0.309384 7.7805 0.009164 **
## Residuals 10 0.39764 0.039764
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(aov(SI~desert, data = df_tsf))## Df Sum Sq Mean Sq F value Pr(>F)
## desert 2 0.5947 0.29735 7.45 0.009 **
## Residuals 11 0.4390 0.03991
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1sim_emm <- emmeans(sim_mod, pairwise ~ desert)
test(sim_emm)## $emmeans
## desert emmean SE df t.ratio p.value
## Mojave 0.912 0.1432 10 6.372 0.0001
## San Joaquin 0.762 0.0717 10 10.624 <.0001
## Sonoran 0.338 0.1007 10 3.356 0.0073
##
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## Mojave - San Joaquin 0.150 0.162 10 0.924 0.6381
## Mojave - Sonoran 0.574 0.173 10 3.318 0.0194
## San Joaquin - Sonoran 0.424 0.125 10 3.389 0.0173
##
## P value adjustment: tukey method for comparing a family of 3 estimatesggplot(df_tsf, aes(x=yearsSinceFire, y=SI * 100, color = desert)) +
geom_point(aes(shape = desert), size = 3) +
#geom_errorbar(aes(ymin=SI_CIlower*100, ymax=ifelse(SI_CIupper*100 > 100, 100, SI_CIupper*100), color = desert), width=.25, show.legend = FALSE) +
geom_smooth(method = "lm", se = FALSE, show.legend = FALSE) +
ylab("Similarity (%)") +
xlab("Time since fire (years)") +
scale_y_continuous(limits = c(0, 105), breaks = c(0, 20, 40, 60, 80, 100)) +
scale_x_continuous(breaks = c(0, 5, 10, 15, 20, 25)) +
scale_color_manual(values = c("#d95f02", "#7570b3", "#1b9e77")) +
guides(linetype = "none") +
theme_bw() +
theme(axis.text.y = element_text(size = 16),
axis.text.x = element_text(size = 16),
axis.title.x = element_text(size = 16),
axis.title.y = element_text(size = 16),
legend.text=element_text(size=16),
axis.ticks.x = element_blank(),
panel.grid.major.y = element_blank(),
legend.position = "top",
legend.title=element_blank())## `geom_smooth()` using formula 'y ~ x'
print(df_tsf)## desert yearsSinceFire SpRich_burned SpRich_control SI SI_CIlower
## 1 Mojave 7 5 13 0.8746 0.622936
## 2 Mojave 9 5 17 1.0000 0.564292
## 3 San Joaquin 1 5 14 0.9039 0.387636
## 4 San Joaquin 4 2 15 0.7227 0.606080
## 5 San Joaquin 8 3 12 0.7115 0.287748
## 6 San Joaquin 12 4 16 0.9372 0.711996
## 7 San Joaquin 15 3 16 0.7818 0.536800
## 8 San Joaquin 16 4 16 1.0000 0.777932
## 9 San Joaquin 18 2 11 0.2771 0.181844
## 10 San Joaquin 19 2 14 0.6551 0.248204
## 11 Sonoran 6 3 18 0.1762 0.106816
## 12 Sonoran 7 4 20 0.5460 0.153608
## 13 Sonoran 8 3 17 0.3959 0.164228
## 14 Sonoran 15 3 19 0.2921 0.131380
## SI_CIupper
## 1 1.126264
## 2 1.435708
## 3 1.420164
## 4 0.839320
## 5 1.135252
## 6 1.162404
## 7 1.026800
## 8 1.222068
## 9 0.372356
## 10 1.061996
## 11 0.245584
## 12 0.938392
## 13 0.627572
## 14 0.452820species-specific
## `summarise()` has grouped output by 'species', 'treatmentGroup', 'ez_taxa'. You
## can override using the `.groups` argument.## # A tibble: 36 × 5
## # Groups: species, ez_taxa [36]
## species ez_taxa funcGroup burned control
## <chr> <chr> <chr> <dbl> <dbl>
## 1 Agelaius tricolor Aves Birds:Passerines 38 4432
## 2 Ambystoma californiense Other Salamanders 0 1
## 3 Ammospermophilus nelsoni Other Small Mammals 3 107
## 4 Anaxyrus californicus Other Frogs and Toads 2 0
## 5 Batrachoseps stebbinsi Other Salamanders 0 1
## 6 Bombus crotchii Other Bees 1 15
## 7 Branta hutchinsii Aves Birds:Rails and Geese 0 4
## 8 Buteo swainsoni Aves Birds:Raptors 43 6366
## 9 Colaptes chrysoides Aves Birds:Woodpeckers 74 113
## 10 Coleonyx switaki Other Lizards and Geckos 0 3
## # … with 26 more rows