# Requires lme4; no automatic installation. Synthetic crossed design.
if (!requireNamespace('lme4', quietly=TRUE)) stop('Install lme4 before running this example.')
library(lme4)
set.seed(6006)
d <- expand.grid(subject=factor(1:30),item=factor(1:20))
# Between-subject predictor: no within-subject variation. A subject-specific
# slope for it is confounded with that subject's intercept.
d$x <- rep(rep(c(-.5,.5),each=15), times=20)
d$y <- .7*d$x + rnorm(30)[as.integer(d$subject)] + rnorm(20,sd=.5)[as.integer(d$item)] + rnorm(nrow(d))
within <- tapply(d$x,d$subject,function(z) length(unique(z)))
cat(sprintf('Rows=%d; subjects=%d; items=%d; max within-subject x values=%d\n',nrow(d),nlevels(d$subject),nlevels(d$item),max(within)))
bad <- lmer(y ~ x + (1+x|subject) + (1|item),d,REML=FALSE)
fit <- lmer(y ~ x + (1|subject) + (1|item),d,REML=FALSE)
other <- update(fit,control=lmerControl(optimizer='bobyqa',optCtrl=list(maxfun=100000)))
cat(sprintf('Unsupported subject slope singular=%s; supported model singular=%s\n',isSingular(bad),isSingular(fit)))
cat(sprintf('Supported x=%.4f; alternative optimizer x=%.4f; absolute difference=%.8f\n',fixef(fit)['x'],fixef(other)['x'],abs(fixef(fit)['x']-fixef(other)['x'])))
cat('Supported model convergence messages:\n'); print(fit@optinfo$conv$lme4$messages)
cat('lme4 version:',as.character(packageVersion('lme4')),'\n')
stopifnot(all(within==1),abs(fixef(fit)['x']-fixef(other)['x'])<1e-4)

# Maintainer-only export; run from repository root. Synthetic data only.
if (Sys.getenv("CLINIC_EXPORT_DATA") == "1") {
write.csv(d,'analysis-clinic/examples/case-006-shared-data.csv',row.names=FALSE,quote=FALSE)
}