Lutra_lutra_INPA
Ron Chen
2024-03-18
Data validation and exploration
## [1] "Lutra lutra"## [1] "Golan Heights" "Harod and Beit Shean valley"
## [3] "Hula valley and Higher Yarden river" "Kineret Basin and Yarden valley"## [1] "Golan Heights" "Hula valley and Higher Yarden river"
## [3] "Yarden valley"| Site | Area | Year | Year_rescaled | |
|---|---|---|---|---|
| Golan Heights :21 | Golan Heights :21 | Min. :2002 | Min. :-10 | |
| Harod and Beit Shean valley :21 | Hula valley and Higher Yarden river:21 | 1st Qu.:2007 | 1st Qu.: -5 | |
| Hula valley and Higher Yarden river:21 | Yarden valley :42 | Median :2012 | Median : 0 | |
| Kineret Basin and Yarden valley :21 | NA | Mean :2012 | Mean : 0 | |
| NA | NA | 3rd Qu.:2017 | 3rd Qu.: 5 | |
| NA | NA | Max. :2022 | Max. : 10 |
| Site | Area | Year | Year_rescaled | |
|---|---|---|---|---|
| Site | 1.0000000 | 0.6741999 | 0 | 0 |
| Area | 0.6741999 | 1.0000000 | 0 | 0 |
| Year | 0.0000000 | 0.0000000 | 1 | 1 |
| Year_rescaled | 0.0000000 | 0.0000000 | 1 | 1 |
Surveys occured once a year between 2002 - 2022
Harod and Beit Shean valley - no NAs, zero counts between 2012 - 2019 and in 2022
Golan Heights - no NAs, zero counts between 2016 - 2022
Kineret Basin and Yarden valley - no NAs, no zero counts
Hula valley and Higer Yarden river - no NAs, no zero counts
Explore data
Try fixed model with interaction between site and year
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!| Observations | 84 |
| Dependent variable | Count_proportion |
| Type | Generalized linear model |
| Family | binomial |
| Link | logit |
| χ²(7) | 26.17 |
| Pseudo-R² (Cragg-Uhler) | 0.57 |
| Pseudo-R² (McFadden) | 0.41 |
| AIC | 81.09 |
| BIC | 100.54 |
| Est. | S.E. | z val. | p | |
|---|---|---|---|---|
| (Intercept) | -2.25 | 0.91 | -2.49 | 0.01 |
| Year_rescaled | -0.20 | 0.14 | -1.43 | 0.15 |
| SiteHarod and Beit Shean valley | 0.97 | 1.16 | 0.84 | 0.40 |
| SiteHula valley and Higher Yarden river | 2.87 | 1.02 | 2.81 | 0.00 |
| SiteKineret Basin and Yarden valley | 1.99 | 1.01 | 1.97 | 0.05 |
| Year_rescaled:SiteHarod and Beit Shean valley | -0.12 | 0.20 | -0.60 | 0.55 |
| Year_rescaled:SiteHula valley and Higher Yarden river | 0.12 | 0.16 | 0.74 | 0.46 |
| Year_rescaled:SiteKineret Basin and Yarden valley | 0.14 | 0.16 | 0.86 | 0.39 |
| Standard errors: MLE |
All interaction terms are not significant. Test if model w/o interaction term is preferable:
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!## Single term deletions
##
## Model:
## Count_proportion ~ Year_rescaled * Site
## Df Deviance AIC
## <none> 7.4436 81.089
## Year_rescaled:Site 3 11.4718 79.117Lutra - Model w/o interaction term is preferable.
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!| Observations | 84 |
| Dependent variable | Count_proportion |
| Type | Generalized linear model |
| Family | binomial |
| Link | logit |
| χ²(4) | 22.15 |
| Pseudo-R² (Cragg-Uhler) | 0.54 |
| Pseudo-R² (McFadden) | 0.38 |
| AIC | 78.30 |
| BIC | 90.45 |
| Est. | S.E. | z val. | p | |
|---|---|---|---|---|
| (Intercept) | -2.00 | 0.66 | -3.03 | 0.00 |
| Year_rescaled | -0.14 | 0.05 | -2.95 | 0.00 |
| SiteHarod and Beit Shean valley | 1.16 | 0.82 | 1.42 | 0.15 |
| SiteHula valley and Higher Yarden river | 2.68 | 0.83 | 3.21 | 0.00 |
| SiteKineret Basin and Yarden valley | 1.71 | 0.81 | 2.12 | 0.03 |
| Standard errors: MLE |
Apply mixed model, using site as random variable, and compare with fixed model
## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm!| Observations | 84 |
| Dependent variable | Count_proportion |
| Type | Mixed effects generalized linear model |
| Family | binomial |
| Link | logit |
| AIC | 80.73 |
| BIC | 88.02 |
| Pseudo-R² (fixed effects) | 0.32 |
| Pseudo-R² (total) | 0.62 |
| Est. | S.E. | z val. | p | |
|---|---|---|---|---|
| (Intercept) | -1.00 | 0.87 | -1.15 | 0.25 |
| Year_rescaled | -0.27 | 0.07 | -3.98 | 0.00 |
| Group | Parameter | Std. Dev. |
|---|---|---|
| Site | (Intercept) | 1.60 |
| Group | # groups | ICC |
|---|---|---|
| Site | 4 | 0.44 |
Compare the two models
AIC(Lutra_fixed, Lutra_mixed)## df AIC
## Lutra_fixed 5 78.29554
## Lutra_mixed 3 80.73104Fixed model is preferable by AIC; however values are close and we have a rule to prefer mixed model whenever possible. Since we ruled out interaction with site, we Keep mixed model. Final model includes year as fixed factor and site as random factor.
Diagnostics of mixed model:
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Count_proportion ~ Year_rescaled + (1 | Site)
## Data: D
##
## AIC BIC logLik deviance df.resid
## 80.7 88.0 -37.4 74.7 81
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.32634 -0.32206 -0.07228 0.32916 1.59342
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 2.549 1.597
## Number of obs: 84, groups: Site, 4
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.00001 0.87050 -1.149 0.251
## Year_rescaled -0.27276 0.06857 -3.978 6.96e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Year_rescld 0.157## $Site
## (Intercept)
## Golan Heights -1.3481307
## Harod and Beit Shean valley -0.1239352
## Hula valley and Higher Yarden river 1.7112316
## Kineret Basin and Yarden valley 0.5371529
##
## with conditional variances for "Site"## $Site
Residuals vs. fitted plot is quite BAD. Revert to fixed model. Diagnostics of fixed model:
Diagnostics of fixed model are a little better. Attempted probit link, same diagnostics. Keep fixed model with logit link:
##
## Call:
## glm(formula = Count_proportion ~ Year_rescaled + Site, family = binomial(link = "logit"),
## data = D)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.01793 -0.28255 -0.02017 0.24063 0.72013
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.99939 0.65940 -3.032 0.00243 **
## Year_rescaled -0.13601 0.04604 -2.954 0.00314 **
## SiteHarod and Beit Shean valley 1.16146 0.81597 1.423 0.15462
## SiteHula valley and Higher Yarden river 2.67551 0.83349 3.210 0.00133 **
## SiteKineret Basin and Yarden valley 1.71032 0.80764 2.118 0.03420 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 33.617 on 83 degrees of freedom
## Residual deviance: 11.472 on 79 degrees of freedom
## AIC: 78.296
##
## Number of Fisher Scoring iterations: 5Lutra - Average response - response averaged across all four sites.
For a given year, if I were to pick a site at random, what would be the response I would get on average?
## Data points may overlap. Use the `jitter` argument to add some amount of
## random variation to the location of data points and avoid overplotting.
Lutra - Predicted response for all four sites:
We do not calculate mean annual decrease, as this decrease is not exponential but rather logistic. Calculate total decrease over entire period from 2002 to 2022 (percent decline):
tot_decline <- diff(range(Lutra_response_mean$predicted))/Lutra_response_mean$predicted[1]*100
print(tot_decline)## [1] 82.00716Average decrease per generation time (7.6 years) based on total decrease, and assuming exponential growth rate (NOT true for binomial model):
## [1] 47.88799Calculate total decrease over entire period from 2002 to 2022, for each of the four sites (percent decline):
## group percent_decline
## 1: Golan Heights 90.27600
## 2: Harod and Beit Shean valley 84.07868
## 3: Hula valley and Higher Yarden river 62.08450
## 4: Kineret Basin and Yarden valley 78.35292Plot average response - with proper years:
Plot average response - for publication of state of nature report
## png
## 2Plot average response without points, with site specific predictions:
Plot average response without points, with site specific predictions - for publication of state of nature report
## png
## 2