Right now, if the model includes other smooth terms, the indexing for the random effects won't work. This can probably be solved by creating a smooths list object that are just the model$smooth that are constructed with bs = 're' from the outset.