rbchan / unmarked

R package for hierarchical models in ecological research
https://rbchan.github.io/unmarked/
37 stars 25 forks source link

Adds predict method for unmarkedRanef class #170

Closed kenkellner closed 4 years ago

kenkellner commented 4 years ago

User supplies the ranef object and a custom function, and predict returns an estimate, SE, and 95% CI for the function's output parameters. Also adds a posteriorSamples function that just grabs the samples, so you can work on them separately. See example below. Fixes #166.

@rbchan, if you could quickly look over the implementation and the docs that would be great since I'm not sure I explained everything correctly.

library(unmarked)
set.seed(4564)
R <- 10
J <- 5
N <- rpois(R, 3)
y <- matrix(NA, R, J)
y[] <- rbinom(R*J, N, 0.5)
K <- 15
umf <- unmarkedFramePCount(y=y)
fm <- pcount(~1 ~1, umf, K=K)
re <- ranef(fm)

#Get posterior samples
(ps <- posteriorSamples(re, nsim=10))

## Posterior samples from unmarked model
## 10 sites x 1 primary periods x 10 sims
## Showing first 5 sites and first 3 simulations
## To see all samples, use print()
##      [,1] [,2] [,3]
## [1,]    1    2    1
## [2,]    6    5    6
## [3,]    2    4    3
## [4,]    5    5    4
## [5,]    5    5    5

#Custom function that calculates average abundance of sites 1-5 and sites 6-10
 myfunc <- function(x){ #x is 10 x 1
    c(mean(x[1:5,1]), mean(x[6:10,1]))
}

(pr <- predict(re, fun=myfunc, nsim=1000))

##   Predicted        SE lower upper
## 1    4.2168 0.4153677   3.4   5.0
## 2    4.1012 0.3682039   3.4   4.8

Source package: unmarked_0.13-2.9003.tar.gz Binary: unmarked_0.13-2.9003.zip

rbchan commented 4 years ago

Looks great. Just a couple of minor requests. Could you let predict return the entire set of samples from the posterior predictive distribution, instead of just the mean, SE, and CI? Also, for the SE, I think you want to use sd instead instead var, right?

I'll look at the docs later.

After you wrap this up, I'd say we're ready for 1.0. What do you think?

kenkellner commented 4 years ago

I did take sqrt() of var(), so I think the values are correct, but I don't know why I didn't just use sd().

Do you think the function should return the summary stats AND the samples, or just the samples?

kenkellner commented 4 years ago

Regarding 1.0 I think we're almost ready. In addition to this I want to fix #164 and #168 first but should get them done this weekend.

If possible could you go through the old issues quick and close any that aren't relevant anymore? There are a few where I can't tell if they're fixed or not, eg #62, #63, #81. I don't think all outstanding issues necessarily need to be addressed before 1.0 but just want to make sure all the important stuff is done.

rbchan commented 4 years ago

My bad. I didn't notice that sqrt.

I suppose it might be easier for users if predict returned the summary stats and samples.

Sounds like a good plan for 1.0. Thanks again.

kenkellner commented 4 years ago

predict now returns all the samples. I put an example of calculating summary stats on the output in the help file, rather than a dedicated function.