Sanity check wanted on model testing

[infotroph]

Your mission

All of the scripts in the stan/ directory are fitting hierarchical Bayesian models to a ~15000-image minirhizotron dataset collected from four different bioenergy cropping systems. To test them, I want to run each model on synthetic data generated from known parameter values.

If you can help with this, please start by commenting on sim_mixdist.r and text_mix_tube_depth.r. All feedback from math corrections to nitpicky style suggestions are welcome, but especially:

• Can anyone besides me understand this? What parts should I document better?
• Does it looks like it produces data that in fact provide a reasonable test of the models I describe below?
• As currently implemented, this is a pretty monolithic testing approach. Are there ways it could be broken out into units?
• Are there obvious ways I'm shooting myself in the foot for future extensions to the model?

  Background

The main quantity of interest is total standing root mass (eventually through time, but time is not yet implemented). There are 96 minirhizotron tubes in total (24 per crop), sampled ~monthly during the growing season from 2010 through 2014. The experiment is laid out in 4 blocks, but I assume here that variance between tubes overwhelms any block effect.

Needed for an effective simulation

• [x] Root volume, in slightly mind-bending units of "mm^3 root per mm^2 image", declines log-linearly with depth, i.e. log(y) ~ a + b*log(depth) + e, e~N(0,sigma).
• [x] Between-tube variance is N(0,sig_tube) and sig_tube is large enough to justify ignoring the experimental blocking.
• [x] Many images contain zero visible roots, probably because each image samples a much smaller area (~ 200 mm^2) than the dominant scale of heterogeneity. To model this, I assume detection is a logistic function of mean root volume: p(y>0) ~ logistic(aa+bb*mu), p(y|y>0) ~ logN(mu, sigma). Credit to Sonderegger et al 2013 for this model formulation.
• [ ] Comparing to destructive harvest methods, minirhizotrons seem to underdetect roots very near the surface. The reasons that remain handwavy but probably involve poor contact between tube and soil. I handle this by assuming a purely empirical fudge factor, the "surface effect", assumed to be sigmoid and declining to zero effect by depth=30 cm or so. Implemented in mix_tube_depth.stan and mix_crop_tube_depth.stan, but not yet in the test function.
• [ ] Differences between crops are the main thing I'm interested in, but I haven't implemented simulating them yet! This should be as simple as adding another rnorm() term, yes?
• [ ] Time effects? Will require that I consider covariance more carefully. Not going to simulate them until I decide how the model should treat them.

[infotroph]

The above wall of text isn't very clear about which parts I'm soliciting help with and which are for me to fix. Looking for feedback on the existing versions of sim_mixdist.R and test_mix_tube_depth.R, anything else is gravy.

[bpbond]

sim_mixdist.R seems straightforward and clear. Two suggestions:

• Preface it with a comment describing what the function does, in general and succinct language. Even better (if you're willing to do a bit more work) would be to write up a full Roxygen header, documenting all parameters and return values.
• Defensive programming. Have a few guard statements at the beginning of the function, checking (via stopifnot, assertthat::assert_is, etc.) types and sizes of incoming parameters.

[bpbond]

I like the way you end test_mix_tube_depth.R with a sessionInfo() call - this is good documentation, particularly if the script's output is mirrored to a log file*. Consider also documenting, in the code, which version of rstan you're currently using.

* Shameless self-promotion: check out luzlogr.