I'm estimating a mixed effects model using simulated data. The basis of this is a conjoint experiment: there are N number of countries in the study with P participants and each respondent is shown the experiment twice. This means that there are NxPx2 observations. Heterogeneity is introduced into the data at the country level and so I run a mixed effect model using lmer with random effects varying by country to account for this variance. However, because each respondent does the experiment twice, I also want to cluster my standard errors at the individual level. My data and model looks something like this:

library(lme4)
data(iris)
 # generating IDs for observations
iris <- iris %>% mutate(id = rep(1:(n()/2), each = 2))
#run model
mod <- lmer(Sepal.Length~Sepal.Width+Petal.Length+Petal.Width + (Sepal.Width+Petal.Length+Petal.Width || Species), data=iris, REML = F, control = lmerControl(optimizer = 'bobyqa'))

I then attempt to get clustered SEs using the parameters package:

library(parameters)
param <- model_parameters(
  mod,
  robust = TRUE,
  vcov_estimation = "CR",
  vcov_type = "CR1",
  vcov_args = list(cluster = iris$id)
)

This returns an error:

Error in vcovCR.lmerMod(obj = new("lmerModLmerTest", vcov_varpar = c(0.00740122363004, : Non-nested random effects detected. clubSandwich methods are not available for such models.

I'm not married to any one method or anything. I just want to return clustered SEs for this type of model specification. As of now I can't find any package that does this. Does anyone know how this can be done, or if such a model even makes sense? I'm new to MLMs but I was thinking if I were to run this as a simple linear model I would lm_robust and cluster by individual so it makes sense to me that I should do the same here as well.