Error: number of rows of matrices must match (see arg 2)

#1

I’m a relatively new greta user, and I’m trying to calculate deviance so that I can calculate DIC for a logistic regression model (see this related post).

I get the following error when trying to run my code (appended at the end of this post):

Error in (function (..., deparse.level = 1)  : 
  number of rows of matrices must match (see arg 2)

Any help would be much appreciated. I’m sure it’s something simple…

Thanks,
Vincent

Here is my code:

### Simulate
n <- 1000
beta_gen <- c(-1, 0.25)
covariates <- cbind(1, rnorm(n, 0, 1))

mu_gen <- covariates %*% beta_gen

y_r <- rbinom(n, 1, prob = boot::inv.logit(mu_gen))

y <- as_data(y_r)

# priors
beta <- normal(0, 1.5, dim = ncol(covariates))

# likelihood
prob <- covariates %*% beta
distribution(y) <- bernoulli(prob = ilogit(prob))

# Deviance
deviance <- -2 * sum(dbinom(y, size = 1, prob = boot::inv.logit(prob), log = TRUE))

# model
m <- model(beta, deviance)
n_chains <- 3

# inits
mean_int <- boot::logit(sum(y_r)/length(y_r))
beta_inits <- replicate(
  n_chains,
  initials(
    beta = c(rnorm(1, mean_int, 0.05), rnorm(ncol(covariates) - 1, -0.1, 0.001))
  ),
  simplify = FALSE
)

# Fit model
draws <- mcmc(m, chains = n_chains, initial_values = beta_inits, n_samples = 100, warmup = 100, pb_update = 10)
#2

You passed a greta array to some R functions (boot::inv.logit, dbinom) for which there aren’t greta equivalent functions. It should have errored at that point and I’m a bit surprised it didn’t!

I think the best way to achieve this is to remove the deviance bit from the greta model, do inference, then get the posterior samples for prob as a matrix with: as.matrix(calculate(prob, draws)) and do your DIC calculations on the R side.

#3

Ah, I see! Thanks very much, I’ll give it a shot and post back here with code.

#4

Just following up here with the solution:

## Calculate DIC
# Deviance
probs <- as.data.frame(t(do.call(rbind, calculate(prob, draws))))

deviance <- function(prob, y){
  -2 * sum(dbinom(y, size = 1, prob = boot::inv.logit(prob), log = TRUE))
}

cl <- makeCluster(detectCores() - 1)
dev <- unlist(parLapply(cl = cl, X = probs, fun = deviance, y = y_r))
stopCluster(cl)
rm(cl)

p_d <- mean(dev) - deviance(covariates %*% colMeans(do.call(rbind,draws)), y_r)

DIC <- p_d + mean(dev)
#5

Nice, thanks for posting!

BTW, you can do as.matrix(draws) to turn the coda mcmc.list object into a matrix (each row a different posterior sample). And plogis() from the stats package (so in base R) should give you the same result as boot::inv.logit().

1 Like