diff options
author | Ken Kellner <ken@kenkellner.com> | 2023-03-21 16:23:27 -0400 |
---|---|---|
committer | Ken Kellner <ken@kenkellner.com> | 2023-03-21 16:23:27 -0400 |
commit | 81c0bb9bd490980daa49fe787eac788f1f9bc8cb (patch) | |
tree | 87b3f67fae03473e20225caacf12c0fdec904127 | |
parent | 8382a647fe6d1a57c2d5e2f2b9f2dc70bc26f1ad (diff) |
-rw-r--r-- | Makefile | 2 | ||||
-rw-r--r-- | renv.lock | 12 | ||||
-rw-r--r-- | unmarked_Paper_Analysis.Rmd | 77 |
3 files changed, 67 insertions, 24 deletions
@@ -4,5 +4,5 @@ docker: rm -rf docker-output mkdir -p docker-output docker cp $(ID):/unmarked-paper/unmarked_Paper_Analysis.html docker-output/ - docker cp $(ID):/unmarked-paper/Figure_4.tiff docker-output/ + docker cp $(ID):/unmarked-paper/Figure_2.tiff docker-output/ docker rm -v $(ID) @@ -116,10 +116,10 @@ }, "TMB": { "Package": "TMB", - "Version": "1.9.1", + "Version": "1.9.2", "Source": "Repository", "Repository": "CRAN", - "Hash": "99cfaf7c7447b6848ff204d22783897c", + "Hash": "73d30a9f18ac4e021aa5ad231939d932", "Requirements": [ "Matrix", "RcppEigen" @@ -1084,15 +1084,15 @@ }, "unmarked": { "Package": "unmarked", - "Version": "1.2.5.9008", + "Version": "1.2.5.9012", "Source": "GitHub", "RemoteType": "github", "RemoteHost": "api.github.com", "RemoteRepo": "unmarked", "RemoteUsername": "rbchan", - "RemoteRef": "2acf96a", - "RemoteSha": "2acf96a7b3feb55e8cc014ad13133fd610e89ff9", - "Hash": "b8191712b6941d47028b3a02f9a603c3", + "RemoteRef": "845a8e7", + "RemoteSha": "845a8e74b5e5eca078bd9e1f5581c8ae94ee266e", + "Hash": "33dba75aded983dbcd75045dcba996b7", "Requirements": [ "MASS", "Matrix", diff --git a/unmarked_Paper_Analysis.Rmd b/unmarked_Paper_Analysis.Rmd index a92a6ba..72e6cf3 100644 --- a/unmarked_Paper_Analysis.Rmd +++ b/unmarked_Paper_Analysis.Rmd @@ -1,7 +1,7 @@ --- title: "The unmarked R package: Twelve years of advances in occurrence and abundance modeling in ecology" author: Kenneth F. Kellner -date: 1 December 2022 +date: 21 March 2023 output: html_document: toc: true @@ -13,7 +13,7 @@ To replicate this analysis, you currently need the dev version of `unmarked` fro You can install it with the `remotes` package: ```{r, eval=FALSE} -remotes::install_github("rbchan/unmarked", ref="2acf96a") +remotes::install_github("rbchan/unmarked", ref="845a8e7") ``` ```{r} @@ -234,7 +234,9 @@ This can be 0, so no need to fill in the blanks. distdata <- acfl %>% group_by(TransectName, Point, Year, DOY, Habitat) %>% summarize(dist25=sum(Count[DistanceBin=="L25"], na.rm=TRUE), - dist50=sum(Count[DistanceBin=="G25"], na.rm=TRUE), .groups='drop') + dist50=sum(Count[DistanceBin=="G25"], na.rm=TRUE), .groups='drop') %>% + mutate(TransectName=as.factor(TransectName), Point=as.factor(Point), + Habitat=as.factor(Habitat)) head(distdata) ``` @@ -340,11 +342,11 @@ stopifnot(all(mods_table$AIC == c(4322.14, 4329.85, 4330.65, 4366.48, 4367.29))) mods_table ``` -Table 1 in paper: +### Build Table 4 ```{r} tab1 <- mods_table[,c("nPars", "AIC", "delta")] -colnames(tab1) <- c("Parameters", "AIC", "ΔAIC") +colnames(tab1) <- c("Parameters", "AIC", "deltaAIC") tab1 <- cbind(Model=rownames(tab1), tab1) rownames(tab1) <- NULL knitr::kable(tab1) @@ -356,7 +358,7 @@ Top model is HABxYEAR. Plot residuals, separately by submodel: -```{r} +```{r, fig.height=7} plot(mod_habxyear) ``` @@ -372,7 +374,9 @@ Model does not appear to fit the data well, but this is a very crude test. ## Inference from top model -Summary table for the top model (Table 2 in paper): +### Build Figure 1 + +This table (with annotations) is Figure 1 in the paper. ```{r} stopifnot(all(round(coef(mod_habxyear),2) == @@ -384,17 +388,26 @@ Statistics reported in the results follow. Hardwood plantations initially had higher abundance: ```{r} +# Get model parameter estimates est <- coef(mod_habxyear) -pct_higher <- round((exp(est[2])-1) * 100) +# Effect of hardwood relative to river levee/bottomland [e^beta] +hardwood_effect <- exp(est[2]) +# Convert the effect size to a percent increase +pct_higher <- round((hardwood_effect-1) * 100) +# Print result names(pct_higher) <- "Percent hardwood higher than bottomland" stopifnot(pct_higher == 55) pct_higher ``` -Percent yearly decline in bottomland hardwoods: +Percent yearly decline in river levee/bottomland hardwoods: ```{r} -pct_decline_bottom <- round((1 - exp(est[3]))*100, 1) +# Since bottomland is the reference level for habitat, this is just +# The effect size of the parameter for "year" +year_effect <- exp(est[3]) +# Convert to a percent decline +pct_decline_bottom <- round((1 - year_effect)*100, 1) names(pct_decline_bottom) <- "Percent yearly decline bottomland" stopifnot(pct_decline_bottom == 2.7) pct_decline_bottom @@ -403,15 +416,44 @@ pct_decline_bottom Percent yearly decline in hardwood plantations: ```{r} -pct_decline_plant <- round((1 - exp(est[3]+est[4]))*100, 1) +# For yearly decline in hardwoods, the effect size is +# year + interaction of year and hardwood +hardwood_year_effect <- exp(est[3] + est[4]) +# Convert to percent decline +pct_decline_plant <- round((1 - hardwood_year_effect)*100, 1) names(pct_decline_plant) <- "Percent yearly decline hardwood plantation" stopifnot(pct_decline_plant == 7.0) pct_decline_plant ``` -## Build Figure 4 +## Effects plots -Predicted abundance over time by habitat: +We can get simple covariate effects plots (i.e., marginal effects plots) using the `plotEffects` function. +This function only allows one covariate at a time to be varying, while the other covariates are held at their median values or reference levels. +For example, comparing predicted habitat: + +```{r} +plotEffects(mod_habxyear, 'lambda', covariate='Habitat') +``` + +Or changes by year: + +```{r} +plotEffects(mod_habxyear, 'lambda', covariate='Year') +``` + +Note that since the habitat covariate is held at its reference level (river levee/bottomland), the plot above represents changes in abundance over time for river levee/bottomland habitat only. + +You can get the data used to make this plot (e.g. to customize the plot yourself) using `plotEffectsData`: + +```{r} +head(plotEffectsData(mod_habxyear, 'lambda', covariate='Year')) +``` + +## Build Figure 2 + +For this model we have an interaction between habitat type and year, which `plotEffects` cannot currently handle. +To show this interaction, we'll need to build the plot manually using `predict` and `ggplot`. ```{r} nd <- data.frame(Habitat=factor(levels(umf@siteCovs$Habitat), @@ -439,7 +481,7 @@ ggplot(data=pr, aes(x=Year+2005, y=Predicted)) + Save plot to file: ```{r} -tiff("Figure_4.tiff", height=5, width=7, units='in', res=300, compression='lzw') +tiff("Figure_2.tiff", height=5, width=7, units='in', res=300, compression='lzw') ggplot(data=pr, aes(x=Year+2005, y=Predicted)) + geom_ribbon(aes(ymin=lower, ymax=upper, fill=Habitat), alpha=0.2) + geom_line(aes(col=Habitat)) + @@ -452,15 +494,16 @@ dev.off() ## Latent Abundance estimates for each site-year -Based on empirical Bayes methods. +The `ranef` function generates posterior distributions of latent abundance or occurrence (abundance in this example) for each site-year. +The `bup` function converts these posterior distributions into a point estimate of abundance for each site-year. See `?ranef` and `?bup` for more information. ```{r} r <- ranef(mod_habxyear) stopifnot(round(bup(r)[1],3) == 8.138) head(bup(r)) ``` - -Compare latent abundance estimates among habitat types: +You can apply the `predict` function to the output from `ranef` in order to generate derived posterior distributions - for example, posteriors for the mean latent abundance in each habitat type. +See `?predict` for more information on this functionality. ```{r} # Function to calculate the mean across sites of each habitat type |