ulam fits giving odd errors

[jebyrnes]

ulam fits giving odd errors

[rmcelreath]

Some bug specific to binomial outcomes perhaps? I may not have time to figure it out today, but will make an issue of it.

[rmcelreath]

Quick assessment, two things going on:

(1) not enough samples in model, so subscript error. this is a bug, but goes away if you sample more chains.

(2) the logit(prob) <- p line generates only a single value, so postcheck generates only a single prediction. can hack it to work by doing something like:

logit(prob) <- p + 0*density

but will need to add special logic to somehow deal with this kind of thing in general.

[rmcelreath]

I have been thinking of replacing the old postcheck with Aki's PIT plots as default. But I need to make more time in my life to work through my dev notes.

[jebyrnes]

Huh - had not run across LOO-PIT yet. Damn this field keeps evolving. I had however, re-implemented Quantile Residual checks, much like in the DHARMa package.

quantile_residuals <- function(fit, n=1000){
  # get y - from rethinking::postcheck
  lik <- (fit@formula)[[1]]
  outcome <- as.character(lik[[2]])
  if (class(fit) == "ulam") 
    outcome <- undot(outcome)
  y <- fit@data[[outcome]]
  
  # get the ecdfs
  s <- sim(fit, n = n)
  dists <- apply(s, 1, ecdf)
  
  #get the quantile residuals
  quant_res <- numeric(ncol(s))
  for(i in 1:length(y)){
    quant_res[i] <- dists[[i]](y[i])
  }
  
  return(quant_res)
}

This is nice because then you can make QQ plots like so:

data(cars)
flist <- alist(
  dist ~ dnorm( mu , sigma ) ,
  mu <- a+b*speed ,
  c(a,b) ~ dnorm(0,1) , 
  sigma ~ dexp(1)
)
fit <- quap( flist , start=list(a=40,b=0.1,sigma=20) , data=cars )

quantile_residuals(fit) |>
  gap::qqunif(logscale=FALSE)

But - I'm guessing this will have similar issues for binomial models?
I'd imagine calculating LOO-PITs must be somewhat similar? Need to read, though.....

And, I suppose if you want to make a hot plot like the ones I now see in bayesplot - well, you could do it with quantil residuals as well

# plot against unif dists
get_rep_unifs <- function(samples, n = 100){
  data.frame(
    id = sort(rep(1:samples, n)),
    unif = runif(samples*n)
  )
}

library(ggplot2)
bayes_qr_plot <- function(fit, n = 1000, n_unif = 100, ...){
  q <- quantile_residuals(fit)
  get_rep_unifs(length(q), n_unif) |>
    ggplot(aes(x = unif, group = id)) +
    geom_density(color = alpha("black", 0.3), ...) +
    geom_density(data = data.frame(unif = q, id = 1),
                 linewidth = 3, ...)
}

bayes_qr_plot(fit)

Which I suppose could serve the basis of a plotting function. Yes, yes, I know, I prefer ggplot2. I'm very tidy. Or try to be. I'm sure this could be done in base, but, wow, my base brain is.... stale.