Simulated Annealing with simulatedAnnealing()

How it works

Simulated annealing mimics the physical process of annealing metals. Kirkpatrick et al. (1983) introduced the analogy; simulatedAnnealing() follows that demonstration closely, adapted for psychometric models.

At each step:

A “neighbor” candidate model is generated from the current model (by default, swapping one or more items for a short form, or freeing/fixing a parameter for a full-model search).
The neighbor’s criterion value is compared to the current model’s.
If the neighbor is better, it’s accepted. If it’s worse, it’s accepted anyway with some probability that depends on how much worse it is and the current “temperature” – which shrinks toward zero as the search progresses, so the algorithm accepts worse moves more readily early on (to escape local optima) and becomes greedier as it converges.
The best model found so far is tracked separately from the current model, since the current model can still wander to worse solutions.

A basic example

As with antColony(), every candidate item must already appear on its factor’s line in initialModel – factors and each factor’s candidate item pool are derived directly from that syntax.

set.seed(58310)

result <- suppressWarnings(simulatedAnnealing(
  initialModel = " visual  =~ x1 + x2 + x3
                   textual =~ x4 + x5 + x6
                   speed   =~ x7 + x8 + x9 ",
  originalData = lavaan::HolzingerSwineford1939,
  maxIterations = 3,
  criterion = "cfi",
  negateCriterion = TRUE,
  itemsPerFactor = c(2, 2, 2),
  items = paste0("x", 1:9)
))
#> Initializing short form creation.
#> The initial short form is:
#> visual =~ x2 + x1
#> textual =~ x6 + x4
#> speed =~ x8 + x7
#> 
#> Using the short form randomNeighbor function.
#> Finished initializing short form options.
#>  Current Progress: 
#> Old Fit: 0.97 New Fit: 0.985                                                                     Current Step = 2 of a maximum 3.   Current Step = 3 of a maximum 3.

result
#> Algorithm: Simulated Annealing
#> Total Run Time: 0.058 secs using 1 chains. 
#> 
#> Function call:
#> simulatedAnnealing(initialModel = " visual =~ x1 + x2 + x3\n textual =~ x4 + x5
#>   + x6\n speed =~ x7 + x8 + x9 ", originalData = lavaan::HolzingerSwineford1939,
#>   maxIterations = 3, criterion = "cfi", negateCriterion = TRUE, itemsPerFactor
#>   = c(2, 2, 2), items = paste0("x", 1:9), temperature = "linear", Kirkpatrick
#>   = TRUE, randomNeighbor = TRUE, lavaan.model.specs = list(model.type = "cfa",
#>   auto.var = TRUE, estimator = "default", ordered = NULL, int.ov.free = TRUE,
#>   int.lv.free = FALSE, std.lv = TRUE, auto.fix.first = FALSE, auto.fix.single
#>   = TRUE, auto.cov.lv.x = TRUE, auto.th = TRUE, auto.delta = TRUE, auto.cov.y =
#>   TRUE), maxChanges = 5, restartCriteria = "consecutive", maximumConsecutive =
#>   25, bifactor = NULL, setChains = 1, shortForm = T)
#> 
#> Final Model Syntax:
#> visual =~ x3 + x1
#> textual =~ x5 + x4
#> speed =~ x8 + x9
#> 
#> 
#> 
#> Criterion: "cfi" (maximized)
#> Final Model Value: 0.985

itemsPerFactor sets the target number of items to keep per factor (in the order factors appear in initialModel); items is the flat pool of candidate item names to draw from (defaulting to all column names in originalData if omitted).

Inspecting the result

summary(result)
#> Algorithm: Simulated Annealing
#> Total Run Time: 0.058 secs
#> 
#> lavaan 0.7-2 ended normally after 28 iterations
#> 
#>   Estimator                                         ML
#>   Optimization method                           NLMINB
#>   Number of model parameters                        21
#> 
#>   Number of observations                           301
#> 
#> Model Test User Model:
#>                                                       
#>   Test statistic                                13.109
#>   Degrees of freedom                                 6
#>   P-value (Chi-square)                           0.041
#> 
#> 
#> Final Model Syntax:
#> visual =~ x3 + x1
#> textual =~ x5 + x4
#> speed =~ x8 + x9
#> 
#> 
#> 
#> Criterion: "cfi" (maximized)
#> Final Model Value: 0.985

plot() shows the criterion value across chain steps. Since the algorithm can wander to worse solutions before recovering, this trace is not necessarily monotonic – the best model found is what’s returned in result@best_model, not necessarily the last one visited.

plot(result)

An early “burn-in” period (common practice for Monte Carlo-style methods) can be excluded from the plot:

plot(result, burn_in = 1)

The criterion

criterion accepts either a character fit-measure name recognized by lavaan::fitmeasures() (as above), or an arbitrary function that takes a fitted lavaan object and returns a single numeric value:

simulatedAnnealing(
  initialModel = "...",
  originalData = myData,
  maxIterations = 20,
  criterion = function(fit) AIC(fit),
  negateCriterion = FALSE, # smaller AIC is better
  itemsPerFactor = c(6, 6, 6)
)

negateCriterion controls the search direction: TRUE looks for the largest value of criterion (e.g. CFI, where larger is better), FALSE looks for the smallest (e.g. RMSEA or AIC, where smaller is better). criterion/negateCriterion is supplied on its natural scale either way – you never need to pre-negate your own criterion function.

Full-model (non-short-form) searches

Omitting itemsPerFactor switches from item-swap short-form search to a full-model search, where each step frees or fixes one parameter rather than swapping items:

fittedModel <- lavaan::cfa(
  model = " visual  =~ x1 + x2 + x3
            textual =~ x4 + x5 + x6
            speed   =~ x7 + x8 + x9",
  data = lavaan::HolzingerSwineford1939
)

simulatedAnnealing(
  initialModel = fittedModel,
  originalData = lavaan::HolzingerSwineford1939,
  maxIterations = 20,
  criterion = "cfi",
  negateCriterion = TRUE
)

Bifactor models

Pass the name of the general factor as bifactor to have all of the retained items across the other factors also load on it – as with antColony(), this only applies when creating a short form (itemsPerFactor supplied):

bifactorModel <- "
visual  =~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9
textual =~ x4 + x5 + x6
speed   =~ x7 + x8 + x9"

simulatedAnnealing(
  initialModel = bifactorModel,
  originalData = lavaan::HolzingerSwineford1939,
  maxIterations = 20,
  criterion = "cfi", negateCriterion = TRUE,
  itemsPerFactor = c(6, 3, 3),
  items = paste0("x", 1:9),
  bifactor = "visual"
)

Parallel chains

setChains runs multiple independent searches in parallel (each starting from the same initialModel but exploring different random neighbors), which is a useful check against any one chain getting stuck:

simulatedAnnealing(
  initialModel = "...",
  originalData = myData,
  maxIterations = 100,
  criterion = "cfi", negateCriterion = TRUE,
  itemsPerFactor = c(6, 6, 6),
  setChains = 4
)

With setChains > 1, plot() overlays each chain’s trace (up to 8 chains), and the best model/fit reported is the best across all chains.