## set l_min
l_min <- 0

## get l_max
l_max <- test.l.max(X = X, n = 5)

## reduce l_max
l_max <- 0.6 * l_max

## make l vector
l <- seq(from = l_min, to = l_max, length.out = 20)

## get minimum q
tf <- test.factors(X = X, l = l, plot = TRUE)

q_min <- tf$q.min

## get maximum q
tp <- test.parameters(X = X, q = 2:20, l = l, plot = "mRt")

q_max <- tp$q.max

## build P
q_min <- ifelse(test = q_min < 2, yes = 2, no = q_min)
q_max <- ifelse(test = q_max < q_min, yes = NA, no = q_max)
P <- cbind(q_min, q_max, l)
P <- na.exclude(P)

## run all possible models
EM.pot <- test.robustness(X = X, P = P, plot = TRUE)

## plot all loadings
plot(NA, xlim = c(1, length(s)), ylim = range(EM.pot$loadings))
for(i in 1:nrow(EM.pot$loadings)) {lines(EM.pot$loadings[i,], col = adjustcolor(col = 1, alpha.f = 0.1))}

## plot histogram
hist(x = EM.pot$modes, breaks = 90)

## plot stem-and-leaf plot
stem(x = EM.pot$modes)

## define limits
limits <- cbind(c(63, 75, 98, 102),
                c(64, 77, 99, 103))

## extract robust loadings
loadings.robust <- robust.loadings(em = EM.pot, limits = limits, plot = TRUE)

## get optimum l
l_opt <- get.l.opt(X = X, l = l, quality = "mRt", Vqn = loadings.robust$Vqn$mean, plot = TRUE)

## model scores
scores_robust <- robust.scores(loadings = loadings.robust, l = l_opt, mc_n = 20, plot = TRUE)
