Help merging the stancode_opt branch

[ScottClaessens]

This morning, I tried editing the stancode_opt branch to allow for measurement error. You can see the change that I made in 62a881e. However, I think that I haven't quite done the right thing (see below), so I'm opening this up for discussion.

There are a number of outstanding issues that need to be resolved before we can merge this branch.

1. Measurement error model does not run

When I run the measurement error model, I get the following error:

# fit model with measurement error
n <- 5
tree <- ape::rcoal(n)
d <- data.frame(
  id = tree$tip.label,
  x = rnorm(n),
  y = rnorm(n),
  x_se = rexp(n, 5),
  y_se = rexp(n, 5)
)
fit <-
  coev_fit(
    data = d,
    variables = list(
      x = "normal",
      y = "normal"
    ),
    id = "id",
    tree = tree,
    measurement_error = list(
      x = "x_se",
      y = "y_se"
    )
  )

Chain 1 Rejecting initial value:
Chain 1   Error evaluating the log probability at the initial value.
Chain 1 Exception: cholesky_decompose: Matrix m is not positive definite (in '/var/folders/1c/18978_5s35n79f0fw4_d2lyc0000gp/T/RtmpNwagPq/model-16a9ec9fe221.stan', line 246, column 8 to column 149)

The error points to the following line of the Stan code, which is the line I edited in 62a881e:

lp[t] = multi_normal_cholesky_lpdf(tdrift | rep_vector(0.0, J), cholesky_decompose(add_diag(VCV_tips[t, tip_to_seg[t, tip_id[i]]], se[i,])));

My reasoning for this change was that when we have measurement error, we want to do the Cholesky decomposition in the model block after adding the standard errors to the diagonal of VCV_tips. But my guess is that I've done the indexing wrong here?

2. Indexing in generated quantities block

Relatedly, I also have a question about the indexing of VCV_tips in the generated quantities block of this model. In the model block, we have the following Stan code:

for (i in 1:N_tips) {
  for (t in 1:N_tree) {
    for (j in 1:J) terminal_drift_rep[t,i][j] = normal_rng(0, 1);
    terminal_drift_rep[t,i] = cholesky_decompose(add_diag(VCV_tips[t, i], se[i,])) * terminal_drift_rep[t,i];
  }
}

But does VCV_tips[t, i] make sense here, since i refers to tips rather than tree segments?

Apologies @ErikRingen if these are silly questions - I think I'm just not quite following the indexing that was added by the caching.

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.

📢 Thoughts on this report? Let us know!

[ScottClaessens]

My hunch is that because VCV_tips and L_VCV_tips end up being ordered by node, and because the tips are always first, I should be doing:

VCV_tips[t, tip_id[i]]

in the model block rather than:

VCV_tips[t, tip_to_seg[t, tip_id[i]]]

I did the latter because that's what was implemented in f799c8f.

@ErikRingen could you confirm this?

[ScottClaessens]

I've given this a go in 99f0b92. There were two places in the code where I changed the indexing for VCV_tips and L_VCV_tips. With these changes, the measurement error model above now runs without errors. All the test fixture models run successfully too.

@ErikRingen please let me know if you agree with these changes - I don't want to perpetuate indexing issues in the code!

[ErikRingen]

I've given this a go in 99f0b92. There were two places in the code where I changed the indexing for VCV_tips and L_VCV_tips. With these changes, the measurement error model above now runs without errors. All the test fixture models run successfully too.

@ErikRingen please let me know if you agree with these changes - I don't want to perpetuate indexing issues in the code!

Ah yes this makes sense--my original indexing pattern was incorrect. Disturbed that this wasn't caught by the test suite so I add new tests in 93fb665 to ensure that tip vs internal indexing is being used in the right parts of the Stan program. Hopefully this will prevent regressions in the future.

[ScottClaessens]

Thanks @ErikRingen for checking that and for the new tests! I've updated the test fixtures and will merge this into the main branch now.