Updates based on reviewer commentsHEADmaster

3 files changed, 67 insertions, 24 deletions

diff --git a/Makefile b/Makefile
index 0dc5fbb..c223b65 100644
--- a/Makefile
+++ b/Makefile
@@ -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)
diff --git a/renv.lock b/renv.lock
index 3b9b7d9..ca0d558 100644
--- a/renv.lock
+++ b/renv.lock
@@ -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