log_lik fails with matrices in models fit by ulam

[jebyrnes]

log_lik fails with matrices in models fit by ulam

[rmcelreath]

This is tricky, because WAIC will need the multi_normal likelihood decomposed, not summed up as Stan does it. So the code in the gq will need to be completely different than the code in the model block. I don't believe there is a way to get multi_normal_lpdf to return a vector, which is what we need.

[jebyrnes]

I'm always here to find fun bugs for you!

Alternately, if lppd could calculate the Log Likelihood matrix directly from a model (extract data and likelihood function) that could solve it.... but, I can see from playing with that, it's a PITA. But, a thought.

[jebyrnes]

Ah, yes, I see even using sim() on the above with log_lik = FALSE fails. Interesting.

Error in Imat * sigma_sq : non-conformable arrays

[jebyrnes]

BTW: you get a much more workable model if you move where you insert the covariance matrix and make it more similar to the GP models earlier in the chapter. Then again, I don't love using variance again - but, the coefs are the same otherwise.

AND now you can get simulated values and compare models. Haven't gone much further, but, perhaps worth looking at?

dat_list <- list(
  N_spp = nrow(dstan),
  sp_id = 1:nrow(dstan),
  M = standardize(log(dstan$body)),
  B = standardize(log(dstan$brain)),
  G = standardize(log(dstan$group_size)),
  Imat = diag(nrow(dstan)) )

no_phylo_withlik <- ulam(
  alist(
    #likelihood
    B ~ dnorm( mu , sigma ),
    
    # DGP
    mu <- a[sp_id] + bM*M + bG*G,
    
    # Autocorr
    vector[N_spp]:a ~ multi_normal(0, K),
    matrix[N_spp,N_spp]: K <- Imat * sigma_sq,
    
    #priors
    a ~ normal( 0 , 1 ),
    c(bM,bG) ~ normal( 0 , 0.5 ),
    sigma_sq ~ exponential( 1 ),
    sigma ~ exponential (1)
  ), data=dat_list , chains=4 , cores=4, log_lik = TRUE)

[rmcelreath]

Yeah that form is great for computation and extension to GLMs. But less transparent for teaching, because it has an extra distribution. Tradeoffs.

[jebyrnes]

I mean, yes and no. It's the same form as the GP stuff above (i.e., like a random intercept), so, if anything, I think there is some good scaffolding there. Just have to explain the extra distribution. Potato Potahto. I'll try it with my students today and see if it works!