Introduction: You definitely should go into these scalings here, both in terms of solving the original system and the automatic differentiation. The quadratic scaling of the the autodiff (really, the quadratic growth of the elements in the Jacobian) is what makes expanding the system so costly.
Yes, yes, I now see what you mean! I definitely want to discuss the algorithm more thouroughly. Per Daniel's recommendation, I've written some bare bones code with measurements of the CPU times required for each solver. I'm comparing the "dd" case, in which the init and the parameters are double, to the "vv" case, where init and parms are var and we compute the Jacobian. This should give us an "empirical" sense of where the gain in efficiency comes from.
Figure 2: Can you identify graphically the effect of the drug? Perhaps even a before/after with two figures, where you show the nominal feedback mechanism and then the effect of the drug.
Certainly. I can show how the drug causes a drop in neutrophil count and how the feedback mechanism brings it back up, once the body clears the drug. You even see a "spring" effect, where the feedback mechanism overcompensates before moving back to baseline.
Performance Analysis: Why not simulate the data directly in Stan? That would avoid any potential conflicts between the R package implementations and the model coded up in Stan.
To echo Bob, this avoids replicating mistakes. For testing purposes, it is easier to use a third party. If the two solvers disagree, but one agrees with a third program, this gives a clue as to where an error may lie (and it did actually help me debug a mistake).
Figure 3: The figure is fine, but readers will be tempted...
Good point. I also think having a plot with some summary statistics would be helpful. Bill G suggested plotting ratios. For the parameter names, I could probably use greek letters (instead of sigma, etc.)
Figure 4: The relatively large variation in performance is slightly worrying. Any idea why the runs are so variable?
I'm not sure, but I'd guess it is because I've only used 200 iterations. Even though the initial estimates are good and the priors informative, I think the chain can still wander off and land in pathological regions of the posterior. Because there are so few iterations, one such "expedition" would be enough to make a run relatively slow. To know for sure, I'd need to check the output of slower and faster runs. It may make sense to run the test with more iterations.
I haven't even looked, but until I get the word that Comic Sans is removed, I'm not going to!
Ok, I'll change to luminari.