lavaan provides a range of optimization methods with the optim.method argument (nlminb, BFGS, L-BFGS-B, GN, and nlminb.constr). `lav_export_estimation` allows exporting objects and functions necessary to pass a lavaan model into any optimizer that takes a combination of (1) starting values, (2) fit-function, (3) gradient-function, and (4) upper and lower bounds. This allows testing new optimization frameworks.
lav_export_estimation(lavaan_model)List with:
get_coef - When working with equality constraints, lavaan internally uses some transformations. get_coef is a functions that recreates the coef function for the parameters.
starting_values - starting_values to be used in the optimization
objective_function - objective function, expecting the current parameter values and the lavaan model
gradient_function - gradient function, expecting the current parameter values and the lavaan model
lower - lower bounds for parameters
upper - upper bound for parameters
library(lavaan)
model <- '
# latent variable definitions
ind60 =~ x1 + x2 + x3
dem60 =~ y1 + y2 + y3 + y4
dem65 =~ y5 + a*y6 + y7 + y8
# regressions
dem60 ~ ind60
dem65 ~ ind60 + dem60
'
fit <- sem(model,
data = PoliticalDemocracy,
do.fit = FALSE)
est <- lav_export_estimation(lavaan_model = fit)
# The starting values are:
est$starting_values
#> [1] 2.1933907 1.8236688 1.2960357 1.0553606 1.2937500 1.3028203 1.4031904
#> [8] 1.4011790 0.0000000 0.0000000 0.0000000 0.2649934 1.1258394 0.9748386
#> [15] 3.3934258 7.6860427 5.3103620 5.5346730 3.3673404 5.6118270 5.3276901
#> [22] 5.1967100 0.0500000 0.0500000 0.0500000
# Note that these do not have labels (and may also differ from coef(fit)
# in case of equality constraints):
coef(fit)
#> ind60=~x2 ind60=~x3 dem60=~y2 dem60=~y3 dem60=~y4 a
#> 2.193 1.824 1.296 1.055 1.294 1.303
#> dem65=~y7 dem65=~y8 dem60~ind60 dem65~ind60 dem65~dem60 x1~~x1
#> 1.403 1.401 0.000 0.000 0.000 0.265
#> x2~~x2 x3~~x3 y1~~y1 y2~~y2 y3~~y3 y4~~y4
#> 1.126 0.975 3.393 7.686 5.310 5.535
#> y5~~y5 y6~~y6 y7~~y7 y8~~y8 ind60~~ind60 dem60~~dem60
#> 3.367 5.612 5.328 5.197 0.050 0.050
#> dem65~~dem65
#> 0.050
# To get the same parameters, use:
est$get_coef(parameter_values = est$starting_values,
lavaan_model = fit)
#> ind60=~x2 ind60=~x3 dem60=~y2 dem60=~y3 dem60=~y4 a
#> 2.1933907 1.8236688 1.2960357 1.0553606 1.2937500 1.3028203
#> dem65=~y7 dem65=~y8 dem60~ind60 dem65~ind60 dem65~dem60 x1~~x1
#> 1.4031904 1.4011790 0.0000000 0.0000000 0.0000000 0.2649934
#> x2~~x2 x3~~x3 y1~~y1 y2~~y2 y3~~y3 y4~~y4
#> 1.1258394 0.9748386 3.3934258 7.6860427 5.3103620 5.5346730
#> y5~~y5 y6~~y6 y7~~y7 y8~~y8 ind60~~ind60 dem60~~dem60
#> 3.3673404 5.6118270 5.3276901 5.1967100 0.0500000 0.0500000
#> dem65~~dem65
#> 0.0500000
# The objective function can be used to compute the fit at the current estimates:
est$objective_function(parameter_values = est$starting_values,
lavaan_model = fit)
#> [1] 5.533355
#> attr(,"fx.group")
#> [1] 5.533355
# The gradient function can be used to compute the gradients at the current estimates:
est$gradient_function(parameter_values = est$starting_values,
lavaan_model = fit)
#> [1] -0.138617036 -0.155930082 -0.037670199 -0.044517915 -0.048072325
#> [6] -0.051920303 -0.054058104 -0.056808694 -0.230316536 -0.336105367
#> [11] -0.264571954 -0.102276584 0.004678712 -0.056401584 -0.126903723
#> [16] -0.057834891 -0.084105190 -0.077438158 -0.126235245 -0.075307846
#> [21] -0.077377232 -0.078555760 -8.836008864 -2.183835199 -2.884529355
# Together, these elements provide the means to estimate the parameters with a large
# range of optimizers. For simplicity, here is an example using optim:
est_fit <- optim(par = est$starting_values,
fn = est$objective_function,
gr = est$gradient_function,
lavaan_model = fit,
method = "BFGS")
est$get_coef(parameter_values = est_fit$par,
lavaan_model = fit)
#> ind60=~x2 ind60=~x3 dem60=~y2 dem60=~y3 dem60=~y4 a
#> 2.18171736 1.81874551 1.35416563 1.04409332 1.29964381 1.25860619
#> dem65=~y7 dem65=~y8 dem60~ind60 dem65~ind60 dem65~dem60 x1~~x1
#> 1.28255808 1.30986845 1.47374687 0.45323406 0.86442152 0.08181681
#> x2~~x2 x3~~x3 y1~~y1 y2~~y2 y3~~y3 y4~~y4
#> 0.11842731 0.46720776 1.94198161 6.49048382 5.34036763 2.88740523
#> y5~~y5 y6~~y6 y7~~y7 y8~~y8 ind60~~ind60 dem60~~dem60
#> 2.39016366 4.34001561 3.50896168 2.93958456 0.44816321 3.87064378
#> dem65~~dem65
#> 0.11490729
# This is identical to
coef(sem(model,
data = PoliticalDemocracy))
#> Warning: lavaan->lav_model_vcov():
#> The variance-covariance matrix of the estimated parameters (vcov) does not
#> appear to be positive definite! The smallest eigenvalue (= -1.487778e+00)
#> is smaller than zero. This may be a symptom that the model is not
#> identified.
#> ind60=~x2 ind60=~x3 dem60=~y2 dem60=~y3 dem60=~y4 a
#> 2.182 1.819 1.354 1.044 1.300 1.258
#> dem65=~y7 dem65=~y8 dem60~ind60 dem65~ind60 dem65~dem60 x1~~x1
#> 1.282 1.310 1.474 0.453 0.864 0.082
#> x2~~x2 x3~~x3 y1~~y1 y2~~y2 y3~~y3 y4~~y4
#> 0.118 0.467 1.942 6.490 5.340 2.887
#> y5~~y5 y6~~y6 y7~~y7 y8~~y8 ind60~~ind60 dem60~~dem60
#> 2.390 4.343 3.510 2.940 0.448 3.872
#> dem65~~dem65
#> 0.115
# Example using ridge regularization for parameter a
fn_ridge <- function(parameter_values, lavaan_model, est, lambda){
return(est$objective_function(parameter_values = parameter_values,
lavaan_model = lavaan_model) + lambda * parameter_values[6]^2)
}
ridge_fit <- optim(par = est$get_coef(est$starting_values,
lavaan_model = fit),
fn = fn_ridge,
lavaan_model = fit,
est = est,
lambda = 10)
est$get_coef(parameter_values = ridge_fit$par,
lavaan_model = fit)
#> ind60=~x2 ind60=~x3 dem60=~y2 dem60=~y3 dem60=~y4 a
#> 2.15620403 1.86034843 1.24391694 1.11021140 1.51437493 0.06675702
#> dem65=~y7 dem65=~y8 dem60~ind60 dem65~ind60 dem65~dem60 x1~~x1
#> 1.01131213 1.18600528 0.76634803 0.33283591 1.19359335 0.10041833
#> x2~~x2 x3~~x3 y1~~y1 y2~~y2 y3~~y3 y4~~y4
#> 0.13547008 0.38138650 3.05148185 9.39563065 4.16928052 6.22630508
#> y5~~y5 y6~~y6 y7~~y7 y8~~y8 ind60~~ind60 dem60~~dem60
#> 2.87932244 5.86926036 5.16553248 4.14051590 0.71796152 2.28951984
#> dem65~~dem65
#> -0.15685831