Comments on assess.sim.basic.pdf

[mpinsky]

Based on 22fb59b973a1145646e333705d822f156218570b @rBatt

Questions/thoughts:

1) Write out MSOM equations? Would help me differentiate Rhat_true from Rhat_obs.

2) Looks like MSOMs can mistakenly attribute low detectability to low probability of presence. Improved with more data to separate observation process from true process?

3) What kinds of simulation testing has already been published in the literature? What are the key novel messages in this simulation?

4) The mix of different kinds of replication becomes confusing. When we write this up, we'll have to be strategic about how to talk about it all.

5) It is worth focusing on aggregate measures of richness, not just the individual values of \psi.

6) Depending on your message a framing, it may be key to compare MSOM to Chao and other estimators.

7) What do we know about the distribution of p in real systems? Do we tend to have a few species with high p and many with low, or is the distribution closer to Gaussian?

8) Possible framework: Paper framing in terms of interest in biodiversity change at a global level seems reasonable (http://dx.doi.org/10.1126/science.1248484, http://dx.doi.org/10.1038/nclimate2769). But methodological challenges imposed by observation process. Simulations to test ability to account for common observation biases. Then apply to a dataset.

[mtingley]

2) Yes. All Occ-models can mix up psi/p if there's not enough information. Such models are generally unstable, but the typical situation is where the TRUTH is low occupancy and moderate detectability, but you have few sites / replicates, and model interprets this is high occupancy and VERY low detectability (e.g., all sites are occupied but you can never find the species). In these situations, this result is because you don't have enough info on detectability.

3) Very little, as far as I'm aware! Maybe Elise Zipkin has done some? I can't think of any good MSOM simulation papers.

4) Not quite sure what kind of reps you're referring to, but yes. There are some good ways to do this, so I'm happy to help standardize terminology (and give some examples).

5) The purpose of MSOMs is generally to predict richness, not to predict individual occurrence (psi). So, I agree with Malin. If you cared about individual psi, then you would just do a regular occupancy model (1-spp). Similarly, inference on covariates comes from the community- or group-level hyperparameters, not just from a single species. Whether a simulation or an actual study, inference is strongest for MSOMs on community level indices (hyper-params, richness, turnover, etc.)

6) I disagree. We discuss this a bit in our TREE paper (Iknayan et al. 2014). They're really different things, and even Colwell agrees with me that MSOM's should always be better than Chao, given adequate data to fit them. (and I don't see the point of simulations to decide how much data is adequate)

7) Well... it's going to depend on your system, your sampling, and your taxa! In one of the favorite figures I've ever created:

The idea is that each species has an intrinsic species-level detectability that is somewhere on the spectrum from 0 to 1. We can call this a species p-trait. However, the observation and sampling process can then distort this intrinsic p, which is why we add detection covariates (e.g., effort) into our models, because we may have very poor detectability of a highly detectable species if, e.g., we are using the wrong method (e.g., surveying for owls during the daytime).

Just to show the spread, from more of my work (with birds in CA), the following shows the 'p-trait' (as differing by grinnell vs modern) for 200+ species of birds: As you can see, there's a fair amount of clumping (but this could also just be a clumping around the hyper-parameter mean), but also a fair amount of spread, from 0 to 1. (and note that we've gotten systematically better at detecting over time)

[rBatt]

I appreciate these comments, but I'll wait for others to catch up before responding in full.

On Tuesday, September 1, 2015, Morgan Tingley notifications@github.com wrote:

2) Yes. All Occ-models can mix up psi/p if there's not enough information. Such models are generally unstable, but the typical situation is where the TRUTH is low occupancy and moderate detectability, but you have few sites / replicates, and model interprets this is high occupancy and VERY low detectability (e.g., all sites are occupied but you can never find the species). In these situations, this result is because you don't have enough info on detectability.

3) Very little, as far as I'm aware! Maybe Elise Zipkin has done some? I can't think of any good MSOM simulation papers.

4) Not quite sure what kind of reps you're referring to, but yes. There are some good ways to do this, so I'm happy to help standardize terminology (and give some examples).

5) The purpose of MSOMs is generally to predict richness, not to predict individual occurrence (psi). So, I agree with Malin. If you cared about individual psi, then you would just do a regular occupancy model (1-spp). Similarly, inference on covariates comes from the community- or group-level hyperparameters, not just from a single species. Whether a simulation or an actual study, inference is strongest for MSOMs on community level indices (hyper-params, richness, turnover, etc.)

6) I disagree. We discuss this a bit in our TREE paper (Iknayan et al. 2014). They're really different things, and even Colwell agrees with me that MSOM's should always be better than Chao, given adequate data to fit them. (and I don't see the point of simulations to decide how much data is adequate)

7) Well... it's going to depend on your system, your sampling, and your taxa! In one of the favorite figures I've ever created:

[image: screen shot 2015-09-01 at 5 52 35 pm] https://cloud.githubusercontent.com/assets/6967321/9618115/462eab62-50d2-11e5-9de7-f5f5bd890649.png

The idea is that each species has an intrinsic species-level detectability that is somewhere on the spectrum from 0 to 1. We can call this a species p-trait. However, the observation and sampling process can then distort this intrinsic p, which is why we add detection covariates (e.g., effort) into our models, because we may have very poor detectability of a highly detectable species if, e.g., we are using the wrong method (e.g., surveying for owls during the daytime).

Just to show the spread, from more of my work (with birds in CA), the following shows the 'p-trait' (as differing by grinnell vs modern) for 200+ species of birds: [image: screen shot 2015-09-01 at 5 56 14 pm] https://cloud.githubusercontent.com/assets/6967321/9618193/c3317b9e-50d2-11e5-847b-7e5c1df9aeb9.png As you can see, there's a fair amount of clumping (but this could also just be a clumping around the hyper-parameter mean), but also a fair amount of spread, from 0 to 1. (and note that we've gotten systematically better at detecting over time)

— Reply to this email directly or view it on GitHub https://github.com/rBatt/trawl/issues/93#issuecomment-136876284.

[JWMorley]

Sorry for the delayed response, had to study up on the topic a bit more. Most of my thoughts relate to using an MSOM on trawl data in general, but for this first post I'll comment on the simulation modeling Ryan has done.

1. I have some difficulty with detectability being constant (within a year/taxa) if a species is occupying a site. Biomass greatly affects detectability and is spatially variable. Further, fish have much greater interannual variability in biomass than other verts. A solution would be If total biomass were annually drawn from a distribution, and then allocated to strata based on the environmental covariate. Detectability could them be refined per site based on biomass. Perhaps this is too complex for a simple sim. Also, you might end up with every species occupying all sites!
2. In Fig.1 and 2, it appears that the MSOM can give estimates of occupancy lower than what is observed...or am I interpreting this wrong?
3. The species identification issue could be more of a selling point for this simulation. In reading Morgan's paper, correct species ID is a basic assumption of MSOM. So this may pave new ground. Estimating diversity from trawl data may be increasingly used as an ecosystem indicator (e.g. Large, et al. (2015) MEPS 521: 1-17), but I'm guessing most folks that use survey data don't consider misidentification. Species richness estimates are probably driven by those rarer species, which are the ones more often misidentified. Plus, I think misidentification is probably more of an issue with trawl data than other surveys previously used with MSOMs. Perhaps run the model with and without that step for correct ID probability?
4. Could you better describe how you would use all the years together (e.g. line 511)? Does this mean pooling all the years to get one estimate of richness? If so, it seems questionable because abundance, especially of rare spp, is so variable year to year. This would make the assumption for a 'closed area for colonization/extinction' tough to meet.

[rBatt]

Thanks for the thoughts. @JWMorley -- I understand what you're saying for your second point, but now that you point this out, I'm confused by this too. If a species is observed, it should never be estimated as not-present. This might be due to the way I tallied things up. The Z matrix contains 1's and 0's; however, one averages (e.g.) each element of Z among all posterior draws, so that Z ends up having some fraction between 0 and 1, which basically represents the average probability of that species being present. But if a species is observed, I thought it should always be 1. Yeah, I'm confused now. I'll have to look at this later.

I'll need more time for the other points.

[JWMorley]

My other comments/questions regard using the MSOM on trawl survey data in general. These might be selling points for the paper, or potential criticisms.

1. Has an MSOM previously been used on survey data comparable to a bottom trawl? It seems like surveys for birds, amphibians, etc. are more common. So Ryan's application will really be looking at a lot more species than has previously been attempted, and a much greater variety of organisms including inverts. For example, the south atlantic survey has caught over 300 taxa during its history.
2. We should be explicit about what community we're analyzing. Are we estimating richness for this square of ocean in general, just the benthic community in this square, or probably more likely just the community that is vulnerable to this gear/method. If the latter, the survey will collect small numbers of species not vulnerable to a trawl (reef fish, pelagics), these will drive up the model estimates of how many rare and unsampled species are present.

Perhaps this isn't important, because the same gear is used throughout the time series. However, it is still worth noting the difference b/w a trawl survey and a bird/amphibian/ant survey, which have a bit more of an actual ceiling to total possible number of species. In a trawl, the number of undetected species is probably similar to the number of detected species....even if you were to restrict it to fish alone. Also, the sampling method is probably much less effective than terrestrial surveys, or even surveys of coral reef fish. As an example, I plotted the number of times each taxon in the southeast was collected over 25 years (1000s of hauls). Taxa collected only 1-8 times had the highest frequency.

[mtingley]

In response to @JWMorley (cc @rBatt ):

1. It's definitely odd that you'd get estimated occupancy (psi-hat) lower than naive occupancy (psi-naive). However, @rBatt is right that the Z-matrix can have 1 or 0 even if the observed state was 1. This is a model choice and depends on how you parameterize Z. Most people structure Z as: Z[i,j,k] ~ dbern(psi[i,j,k[) and then parameterize psi as: psi[i,j,k] ~ InvLogit(*params*) In doing this, the probability of occupancy at site Z is independent of the observed state! Psi[i,j,k] is purely a function of the covariates used, which means that dbern(psi) can be 1 or 0, and there will likely be zeros given all the MCMC runs. And this is OK! Remember, in a Bayesian analysis, everything is treated as a random variable. If your model doesn't reproduce your data well (i.e., via a posterior predictive check) then it's just not a very good model. If you really care about site-specific occurrence structure, then you can force the Z-matrix to include 1's where occurrence was known. This won't automatically be done for you.
2. There are a flavor of occupancy-type models that account for false-presences in addition to false-absences. David Miller has done a good job recently at showing how important this is. False-positive models have not been widely adopted yet however due to resistance within various communities (either because people are lazy or they feel offended at the suggestion that they would falsely identify something). I don't think any MSOM-MSAM type analysis has even considered false positives (it would, however, vastly slow down convergence).
3. I'm not sure about total numbers for other studies. My Sierra Nevada paper had 212 bird species in the MSOM. But, I agree with @JWMorley that I don't know of any other study that mixes broad taxonomic groups. I think the ONLY reason that this is justified in a MSOM is because 'detectability' here depends so much on sampling and density, and not on species-specific classical detectability characteristics (e.g., how loud does that salmon sing?).
4. The final point by @JWMorley is also good. The solution is appropriate filtering of your species pool. For example, with birds, EVEN when using a MSOM to predict richness, if you surveyed in forests, you would remove the rare and casual fly-over detections of aquatic birds (e.g., gulls, pelicans) because they do not represent the 'community' that you are trying to model. How much species filtering have you done, @rBatt?

[rBatt]

I'm lovin' the frank perspective on 1. @mtingley :)

I was confused about my confusion; when I initially read @JWMorley 's comment, I thought of this line from the Zipkin et al. 2010 paper (my mind jumps to Zipkin et al b/c it was referenced by Iknayan et al [the first MSOM-esque paper I ever found], and thus Zipkin et al is my first carefully-read MSOM paper):

The above text from Zipkin et al. corresponds to these lines in the JAGS model for the 1-covariate MSOM model, and to this line which sets w to 1 for all "real" (not 0-augmented) species.

Anyway, my fuzzy memory had caused me to incorrectly think that if a species was observed it was known to be present (it would be a 1 in Z), but that's not true; there's nothing about that which prevents a false positive. Really, it acts to ensure that if a species doesn't exist (it has a 0 in w), then it isn't present (has a 0 in Z) and therefore isn't observed (probability of detection, mu.p is 0, leading to a 0 in X, which contains the observed data). So it works to limit the count, not bolster it (speaking loosely here).

---

I don't include the augmented 0's in my tally of richness. Therefore, I'm not too concerned with the "richness of what community? what's the definition of 'community'?" issue; I view our objective as "figure out how many of these trawl critters were actually present in a given place in a given year". The term "richness" is a convenient term to describe the enumeration of these presences; if we're willing to consider macrofauna that either permanently or briefly inhabit coastal bottom habitat, then I think the sampling gear provides a meaningful (and convenient) definition of the community for us, and in that case "richness" refers to the diversity of the the "trawl community". We're using the widest scope for "community richness" that we can given our data.

Semantics aside, what I see as important for statistical reasons is that our "community" is defined in such a way that the distributional assumptions regarding how common (or, as is the more common case, rare) species are. I'm not sure how to validate this assumption. But so long as our definition of "community" doesn't result in all species being very common or very rare, my gut says we're fine.

I also think that it's a good thing that most species are extremely rare. This pattern (right-skewed species abundance curve) would be a standard expectation of macroecology (i.e., that "most species are rare") dating back to Darwin, if not earlier. Note that the skew-right form is not merely an artifact of sampling design (as I understand it). In fact, I'm tempted to use the "most of the species in the trawl data set are quite rare" observation made my @JWMorley as evidence that we have a valid "community".

Also, @mtingley , is this comment of yours:

The final point by @JWMorley is also good. The solution is appropriate filtering of your species pool. For example, with birds, EVEN when using a MSOM to predict richness, if you surveyed in forests, you would remove the rare and casual fly-over detections of aquatic birds (e.g., gulls, pelicans) because they do not represent the 'community' that you are trying to model. How much species filtering have you done, @rBatt?

address by my above comments? Also, was that in response to @JWMorley 's

We should be explicit about what community we're analyzing. Are we estimating richness for this square of ocean in general, just the benthic community in this square, or probably more likely just the community that is vulnerable to this gear/method. If the latter, the survey will collect small numbers of species not vulnerable to a trawl (reef fish, pelagics), these will drive up the model estimates of how many rare and unsampled species are present.

Thoughts on my thoughts (and on it goes, all the way down ...)?

---

The species identification issue could be more of a selling point for this simulation. In reading Morgan's paper, correct species ID is a basic assumption of MSOM.

I think that assumption refers to making sure that things are ID'd correctly when they are ID'd. Deviations from this aren't scaring me much, for some reason. I guess a lot of the improper identification might fall into the category of naming a species as "unid fish" type stuff, which I throw out of the data set, which would result in a lower value of p (probability of detection) for that species. If species are getting mixed up, then yeah, that sucks, and is related to needing to estimate false positives. However, if there isn't a consistent bias (e.g., Ryan always writes down "cod" when he catches a "flounder" or a "bluefish" or a "skate", but when he catches a "cod", he always writes down "cod"), then the consequences shouldn't be too drastic. I could be wrong here, but as of now, I'm not convinced this is of primary concern.

---

In his first comment in this thread, @mtingley said

5) The purpose of MSOMs is generally to predict richness, not to predict individual occurrence (psi). So, I agree with Malin. If you cared about individual psi, then you would just do a regular occupancy model (1-spp). Similarly, inference on covariates comes from the community- or group-level hyperparameters, not just from a single species. Whether a simulation or an actual study, inference is strongest for MSOMs on community level indices (hyper-params, richness, turnover, etc.)

I'm not sure I share this perspective. I think that an MSOM is really good for getting at psi and generating plausible, data-baed response curves for entire communities, because most species in a community are rare, and therefore it would be hard to estimate their thermal niche. I realize that, for well-observed species, you might take a more efficient approach.

It makes sense that MSOMs are most valuable at the community level. That's a good point. But I don't see the general framework of an MSOM as being inherently limiting --- the model structure is flexible. I know it could get unwieldy, possibly, but anything you could do for a 1-spp model, couldn't you incorporate that into an MSOM? Also, even though there are hyperparameters, with lots of data for the well-observed species, it seems like the posterior should arrive at an answer that would be similar to the single-species answer (i.e., the hierarchy doesn't necessarily taint the species-specific estimates).

---

Next Steps

There's a lot more here that I could comment on. I know I'll be referring back to this discussion, and I'd be glad to see it continue to grow. However, this has me thinking that we might be due for a group meeting or Skype call. I'll be writing an abstract for the ASLO/ Ocean Sciences meeting tomorrow. I'll be revisiting possible angles to take with this work, and that thinking might provide some good fodder for group discussion.

Tomorrow I'll send out an email that will survey your dispositions and availabilities to a call or meeting.

[rBatt]

Just read through this again. There's a lot of good material here. I think it might be appropriate to split some of this into new issues. But I'll leave it open for now.

Thanks again for all these comments, everyone. I think this on-line discussion, our phone call, and subsequent discussion over the past several months have been very helpful. I really think everyone in our little group has done a great job of stepping up and making strong intellectual contributions.

Now that the Stan model is together and verified, I need to tweak it to properly handle the many years (namely, make it dynamic #83 ), and see if we can get meaningful results for real data.