I really can’t alter the structure of the model itself as you suggest, for theoretical reasons it should be nested with at least a two-level model. This is actually a simplified model before I work with higher-order models. I’m also assessing these stan models with the ‘standard’ gamm4 models I’ve already built.

I’m not re-creating the exact error, but I suspect they may be related. Here is a worked example:

N <- 300 #Number of students per school

sigma <- 200

x1 <- runif(N, 10, 40)

x2 <- runif(N, 25, 55)

x3 <- runif(N, 40, 70)

x4 <- runif(N, 55, 85)

x5 <- runif(N, 70, 100)

y1 <- 20 + 0 * x1 + rnorm(N, 0, sigma)

y2 <- 40 + 10 * x2 + rnorm(N, 0, sigma)

y3 <- 60 + 20 * x3 + rnorm(N, 0, sigma)

y4 <- 80 + 30 * x4 + rnorm(N, 0, sigma)

y5 <- 100 + 40 * x5 + rnorm(N, 0, sigma)

ID <- rep(LETTERS[1:20], each = N)

test <- data.frame(SES = c(x1, x2, x3, x4, x5),

Score = c(y1, y2, y3, y4, y5), ID = ID)

HLM0 <- stan_gamm4(Score ~ s(SES), family = gaussian(), random = ~(1 | ID), data = test)

strat.test <- kfold_split_stratified(K=10, x= test$ID)

kfold(HLM0, K=10, folds = strat.test)

I’ve been unable to find a worked example of a multilevel model which uses kfold(), which would help with my doing it correctly. I can always write my own code to do cross-validation, I just assumed the built in loo functions would be more efficient.