Skip to contents

Create an extended SEM with OpenMx (Boker et al., 2011) using a lavaan-style (Rosseel, 2012) syntax.

Usage

mxsem(
  model,
  data,
  scale_loadings = TRUE,
  scale_latent_variances = FALSE,
  add_intercepts = TRUE,
  add_variances = TRUE,
  add_exogenous_latent_covariances = TRUE,
  add_exogenous_manifest_covariances = TRUE,
  lbound_variances = TRUE,
  directed = unicode_directed(),
  undirected = unicode_undirected(),
  return_parameter_table = FALSE
)

Arguments

model

model syntax similar to lavaan's syntax

data

raw data used to fit the model. Alternatively, an object created with OpenMx::mxData can be used (e.g., OpenMx::mxData(observed = cov(OpenMx::Bollen), means = colMeans(OpenMx::Bollen), numObs = nrow(OpenMx::Bollen), type = "cov")).

scale_loadings

should the first loading of each latent variable be used for scaling?

scale_latent_variances

should the latent variances be used for scaling?

add_intercepts

should intercepts for manifest variables be added automatically? If set to false, intercepts must be added manually. If no intercepts are added, mxsem will automatically use just the observed covariances and not the observed means.

add_variances

should variances for manifest and latent variables be added automatically?

add_exogenous_latent_covariances

should covariances between exogenous latent variables be added automatically?

add_exogenous_manifest_covariances

should covariances between exogenous manifest variables be added automatically?

lbound_variances

should the lower bound for variances be set to 0.000001?

directed

symbol used to indicate directed effects (regressions and loadings)

undirected

symbol used to indicate undirected effects (variances and covariances)

return_parameter_table

if set to TRUE, the internal parameter table is returend together with the mxModel

Value

mxModel object that can be fitted with mxRun or mxTryHard. If return_parameter_table is TRUE, a list with the mxModel and the parameter table is returned.

Details

Setting up SEM can be tedious. The lavaan (Rosseel, 2012) package provides a great syntax to make the process easier. The objective of mxsem is to provide a similar syntax for OpenMx. OpenMx is a flexible R package for extended SEM. However, note that mxsem only covers a small part of the OpenMx framework by focusing on "standard" SEM. Similar to lavaan's sem()-function, mxsem tries to set up parts of the model automatically (e.g., adding variances automatically or scaling the latent variables automatically). If you want to unlock the full potential of OpenMx, mxsem may not be the best option.

Warning: The syntax and settings of mxsem may differ from lavaan in some cases. See vignette("Syntax", package = "mxsem") for more details on the syntax and the default arguments.

Alternatives

You will find similar functions in the following packages:

  • metaSEM (Cheung, 2015) provides a lavaan2RAM function that can be combined with the create.mxModel function. This combination offers more features than mxsem. For instance, constraints of the form a < b are supported. In mxsem such constraints require algebras (e.g., !diff; a := b - exp(diff)).

  • umx (Bates et al., 2019) provides the umxRAM and umxLav2RAM functions that can parse single lavaan-style statements (e.g., eta =~ y1 + y2 + y3) or an entire lavaan models to OpenMx models.

  • tidySEM (van Lissa, 2023) provides the as_ram function to translate lavaan syntax to OpenMx and also implements a unified syntax to specify both, lavaan and OpenMx models. Additionally, it works well with the tidyverse.

  • ezMx (Bates, et al. 2014) simplifies fitting SEM with OpenMx and also provides a translation of lavaan models to OpenMx with the lavaan.to.OpenMx function.

Because mxsem implements the syntax parser from scratch, it can extend the lavaan syntax to account for specific OpenMx features. This enables implicit transformations with curly braces.

Citation

Cite OpenMx (Boker et al., 2011) for the modeling and lavaan for the syntax (Rosseel, 2012). mxsem itself is just a very small package and lets OpenMx do all the heavy lifting.

Defaults

By default, mxsem scales latent variables by setting the loadings on the first item to 1. This can be changed by setting scale_loadings = FALSE in the function call. Setting scale_latent_variances = TRUE sets latent variances to 1 for scaling.

mxsem will add intercepts for all manifest variables as well as variances for all manifest and latent variables. A lower bound of 1e-6 will be added to all variances. Finally, covariances for all exogenous variables will be added. All of these options can be changed when calling mxsem.

Syntax

The syntax is, for the most part, identical to that of lavaan. The following specifies loadings of a latent variable eta on manifest variables y1-y4:

eta =~ y1 + y2 + y3

Regressions are specified with ~:

xi  =~ x1 + x2 + x3
eta =~ y1 + y2 + y3
# predict eta with xi:
eta ~  xi

Add covariances with ~~

xi  =~ x1 + x2 + x3
eta =~ y1 + y2 + y3
# predict eta with xi:
eta ~  xi
x1 ~~ x2

Intercepts are specified with ~1

xi  =~ x1 + x2 + x3
eta =~ y1 + y2 + y3
# predict eta with xi:
eta ~  xi
x1 ~~ x2

eta ~ 1

Parameter labels and constraints

Add labels to parameters as follows:

xi  =~ l1*x1 + l2*x2 + l3*x3
eta =~ l4*y1 + l5*y2 + l6*y3
# predict eta with xi:
eta ~  b*xi

Fix parameters by using numeric values instead of labels:

xi  =~ 1*x1 + l2*x2 + l3*x3
eta =~ 1*y1 + l5*y2 + l6*y3
# predict eta with xi:
eta ~  b*xi

Bounds

Lower and upper bounds allow for constraints on parameters. For instance, a lower bound can prevent negative variances.

xi  =~ 1*x1 + l2*x2 + l3*x3
eta =~ 1*y1 + l5*y2 + l6*y3
# predict eta with xi:
eta ~  b*xi
# residual variance for x1
x1 ~~ v*x1
# bound:
v > 0

Upper bounds are specified with v < 10. Note that the parameter label must always come first. The following is not allowed: 0 < v or 10 > v.

(Non-)linear constraints

Assume that latent construct eta was observed twice, where eta1 is the first observation and eta2 the second. We want to define the loadings of eta2 on its observations as l_1 + delta_l1. If delta_l1 is zero, we have measurement invariance.

eta1 =~ l1*y1 + l2*y2 + l3*y3
eta2 =~ l4*y4 + l5*y5 + l6*y6
# define new delta-parameter
!delta_1; !delta_2; !delta_3
# redefine l4-l6
l4 := l1 + delta_1
l5 := l2 + delta_2
l6 := l3 + delta_3

Alternatively, implicit transformations can be used as follows:

eta1 =~ l1*y1 + l2*y2 + l3*y3
eta2 =~ {l1 + delta_1} * y4 + {l2 + delta_2} * y5 + {l3 + delta_3} * y6

Specific labels for the transformation results can also be provided:

eta1 =~ l1*y1 + l2*y2 + l3*y3
eta2 =~ {l4 := l1 + delta_1} * y4 + {l5 := l2 + delta_2} * y5 + {l6 := l3 + delta_3} * y6

This is inspired by the approach in metaSEM (Cheung, 2015).

Definition variables

Definition variables allow for person-specific parameter constraints. Use the data.-prefix to specify definition variables.

I =~ 1*y1 + 1*y2 + 1*y3 + 1*y4 + 1*y5
S =~ data.t_1 * y1 + data.t_2 * y2 + data.t_3 * y3 + data.t_4 * y4 + data.t_5 * y5

I ~ int*1
S ~ slp*1

Starting Values

mxsem differs from lavaan in the specification of starting values. Instead of providing starting values in the model syntax, the set_starting_values function is used.

References

  • Bates, T. C., Maes, H., & Neale, M. C. (2019). umx: Twin and Path-Based Structural Equation Modeling in R. Twin Research and Human Genetics, 22(1), 27–41. https://doi.org/10.1017/thg.2019.2

  • Bates, T. C., Prindle, J. J. (2014). ezMx. https://github.com/OpenMx/ezMx

  • Boker, S. M., Neale, M., Maes, H., Wilde, M., Spiegel, M., Brick, T., Spies, J., Estabrook, R., Kenny, S., Bates, T., Mehta, P., & Fox, J. (2011). OpenMx: An Open Source Extended Structural Equation Modeling Framework. Psychometrika, 76(2), 306–317. https://doi.org/10.1007/s11336-010-9200-6

  • Cheung, M. W.-L. (2015). metaSEM: An R package for meta-analysis using structural equation modeling. Frontiers in Psychology, 5. https://doi.org/10.3389/fpsyg.2014.01521

  • Rosseel, Y. (2012). lavaan: An R package for structural equation modeling. Journal of Statistical Software, 48(2), 1–36. https://doi.org/10.18637/jss.v048.i02

  • van Lissa, C. J. (2023). tidySEM: Tidy Structural Equation Modeling. R package version 0.2.4, https://cjvanlissa.github.io/tidySEM/.

Examples

# THE FOLLOWING EXAMPLE IS ADAPTED FROM LAVAAN
library(mxsem)

model <- '
  # latent variable definitions
     ind60 =~ x1 + x2 + x3
     dem60 =~ y1 + a1*y2 + b*y3 + c1*y4
     dem65 =~ y5 + a2*y6 + b*y7 + c2*y8

  # regressions
    dem60 ~ ind60
    dem65 ~ ind60 + dem60

  # residual correlations
    y1 ~~ y5
    y2 ~~ y4 + y6
    y3 ~~ y7
    y4 ~~ y8
    y6 ~~ y8
'

fit <- mxsem(model = model,
            data  = OpenMx::Bollen) |>
  mxTryHard()
#> Running untitled13 with 41 parameters
#> 
#> Beginning initial fit attempt
#> Running untitled13 with 41 parameters
#> 
#>  Lowest minimum so far:  3274.15170117083
#>  Not all eigenvalues of the Hessian are positive: 1002820.9693578, 4177.87917816189, 1100.79628533037, 715.798788568486, 466.997947118047, 315.164898674905, 191.210996173414, 124.723726935937, 113.088444201225, 104.256081240764, 84.1320091971086, 81.3076257177491, 69.4606363216076, 46.4275404178963, 46.0430712481736, 41.8253108289941, 41.6287312581773, 35.4717906795529, 32.9203894004398, 32.8766157598417, 27.9622813632953, 26.6406390383752, 25.2192983941879, 17.916483821106, 14.7320256689223, 13.9566137097917, 10.4036708003971, 9.66303168104416, 7.47548627862979, 5.35388822052549, 4.85071826227911, 3.97666777725942, 3.17327658704878, 2.7085959711609, 2.45460626023246, 1.86311038744241, 1.48654619095132, 0.867182184265156, 0.711663598845945, 0.00130784632992241, -0.317322692530756
#> 
#> Beginning fit attempt 1 of at maximum 10 extra tries
#> Running untitled13 with 41 parameters
#>  Not all eigenvalues of the Hessian are positive: 1002820.9693578, 4177.87917816189, 1100.79628533037, 715.798788568486, 466.997947118047, 315.164898674905, 191.210996173414, 124.723726935937, 113.088444201225, 104.256081240764, 84.1320091971086, 81.3076257177491, 69.4606363216076, 46.4275404178963, 46.0430712481736, 41.8253108289941, 41.6287312581773, 35.4717906795529, 32.9203894004398, 32.8766157598417, 27.9622813632953, 26.6406390383752, 25.2192983941879, 17.916483821106, 14.7320256689223, 13.9566137097917, 10.4036708003971, 9.66303168104416, 7.47548627862979, 5.35388822052549, 4.85071826227911, 3.97666777725942, 3.17327658704878, 2.7085959711609, 2.45460626023246, 1.86311038744241, 1.48654619095132, 0.867182184265156, 0.711663598845945, 0.00130784632992241, -0.317322692530756
#> 
#> Beginning fit attempt 2 of at maximum 10 extra tries
#> Running untitled13 with 41 parameters
#>  Not all eigenvalues of the Hessian are positive: 1095887.4576289, 4336.88859812291, 1101.60095584567, 715.741581089665, 466.998437149653, 315.194564297537, 191.214222532952, 124.73892882919, 113.088825054082, 104.256184668642, 84.1319306999197, 81.3074748088679, 69.4919938064738, 46.4274267129193, 46.042779484438, 41.8254181053709, 41.6255353605796, 35.4737582650051, 32.9198091671539, 32.8695838557629, 27.9513582497545, 26.6403953018149, 25.2229261942084, 17.9163618466597, 14.7230331580777, 13.9602445171162, 10.4029465882407, 9.66830000581792, 7.4764601631925, 5.35130724182903, 4.85510731933932, 3.97589558284252, 3.1704829839838, 2.70908859897041, 2.45409708182186, 1.86196187966812, 1.48615988724095, 0.863169838347198, 0.690872511605086, 0.00126393982350972, -0.301573566565726
#> 
#> Beginning fit attempt 3 of at maximum 10 extra tries
#> Running untitled13 with 41 parameters
#> 
#>  Lowest minimum so far:  3274.15163155368
#>  Not all eigenvalues of the Hessian are positive: 284798.769509955, 2549.7339611363, 1079.27039998206, 716.771194425571, 467.000720559717, 314.738340383315, 191.140661163423, 124.539638654959, 113.089880357975, 104.257004215137, 84.1321828908432, 81.3085153018582, 69.0018968616705, 46.4277217344138, 46.0429474323784, 42.0092498750708, 41.6714199906983, 35.4602515051097, 32.9212115246584, 32.888003738862, 28.0876766268602, 26.6411930020571, 25.2308687208095, 17.9169707979993, 14.9396095512264, 13.974346303053, 10.4331834050266, 9.68150043577825, 7.46689661248203, 5.41096702444768, 4.86554277268622, 4.0242657743474, 3.23299167546002, 2.70778327806949, 2.49858022640508, 1.88866518211077, 1.53414590989677, 1.23971867362244, 0.838322570143965, 0.00270159642946001, -0.475802037693299
#> 
#> Beginning fit attempt 4 of at maximum 10 extra tries
#> Running untitled13 with 41 parameters
#>  Not all eigenvalues of the Hessian are positive: 284798.769509955, 2549.7339611363, 1079.27039998206, 716.771194425571, 467.000720559717, 314.738340383315, 191.140661163423, 124.539638654959, 113.089880357975, 104.257004215137, 84.1321828908432, 81.3085153018582, 69.0018968616705, 46.4277217344138, 46.0429474323784, 42.0092498750708, 41.6714199906983, 35.4602515051097, 32.9212115246584, 32.888003738862, 28.0876766268602, 26.6411930020571, 25.2308687208095, 17.9169707979993, 14.9396095512264, 13.974346303053, 10.4331834050266, 9.68150043577825, 7.46689661248203, 5.41096702444768, 4.86554277268622, 4.0242657743474, 3.23299167546002, 2.70778327806949, 2.49858022640508, 1.88866518211077, 1.53414590989677, 1.23971867362244, 0.838322570143965, 0.00270159642946001, -0.475802037693299
#> 
#> Beginning fit attempt 5 of at maximum 10 extra tries
#> Running untitled13 with 41 parameters
#>  Not all eigenvalues of the Hessian are positive: 348677.73085087, 2745.5061438744, 1084.46232491371, 716.62417631028, 466.999851834883, 314.816467248177, 191.147275465012, 124.571426096468, 113.088825989555, 104.256430192424, 84.1320238704133, 81.3072623603106, 69.0885052137753, 46.4272997229159, 46.0427059622244, 41.9725327759057, 41.6681060032252, 35.4573644980286, 32.9198890396413, 32.8863043401259, 28.0616881551836, 26.6399994123142, 25.228529599107, 17.9163840232495, 14.898188647715, 13.9702946445312, 10.4259195306958, 9.6768959730061, 7.47025163865757, 5.40016478972039, 4.86217798680426, 4.01280021695399, 3.21813770251536, 2.70914615342392, 2.48851191520276, 1.88207249025303, 1.52384310326366, 1.1587837410884, 0.831955996214121, 0.00241703483254824, -0.447475094189825
#> 
#> Beginning fit attempt 6 of at maximum 10 extra tries
#> Running untitled13 with 41 parameters
#> 
#>  Lowest minimum so far:  3234.2876855728
#>  Not all eigenvalues of the Hessian are positive: 7598.57597817472, 1688.56359939695, 539.290535353664, 419.003575947642, 343.804573752009, 338.013425505306, 242.647822318344, 226.520961149572, 179.226984898953, 144.885970225671, 92.0293538950684, 81.0097659424557, 74.2518964424219, 73.5026364672453, 66.3328181163647, 42.7403096523063, 42.2691565112523, 36.0010704691878, 30.668714377437, 28.9584225606522, 22.9697436559588, 19.8982480528, 16.2268370934855, 14.1540607158276, 10.9222196334779, 9.84068882509725, 9.73988036580468, 9.04820257410117, 4.59963356501852, 4.3848690245792, 4.25184980918354, 3.10261191262817, 2.0561702557282, 1.96459038271992, 1.60223037518751, 1.29185272264826, 1.0085922850472, 0.668779485070543, 0.434913485764363, 0.327634180665386, -39.3029899422406
#> 
#> Beginning fit attempt 7 of at maximum 10 extra tries
#> Running untitled13 with 41 parameters
#> 
#>  Lowest minimum so far:  3234.28768556172
#>  Not all eigenvalues of the Hessian are positive: 7598.54465251359, 1688.56226687367, 539.299514018853, 419.003873887013, 343.816416161985, 338.020946238197, 242.654201343044, 226.541543024551, 179.229989013562, 144.893491147561, 92.0299023481529, 81.0105291174884, 74.2577001311838, 73.5029022821897, 66.3428232549056, 42.7408238931461, 42.2641494751532, 36.0038521667295, 30.6755726329355, 28.9587221973923, 22.9699043276952, 19.8988686999703, 16.2333610814914, 14.1543711603645, 10.9238654858736, 9.84132728108342, 9.74384396005334, 9.04747193372708, 4.6044876868251, 4.38396911488603, 4.25831726634215, 3.10420296260445, 2.05836341270424, 1.96572010852991, 1.60313797578899, 1.29368916054298, 1.00856658923608, 0.668641669932772, 0.43703238193563, 0.327851407063622, -39.2933247207047
#> 
#> Beginning fit attempt 8 of at maximum 10 extra tries
#> Running untitled13 with 41 parameters
#> 
#>  Lowest minimum so far:  3234.28768555776
#>  Not all eigenvalues of the Hessian are positive: 7598.5384191669, 1688.56219088878, 539.28559907923, 419.004541415532, 343.808437180686, 338.008215204671, 242.653642676101, 226.532900345163, 179.229359236613, 144.888570088225, 92.0294275720868, 81.010138845785, 74.2546620027886, 73.5025322298607, 66.33422191362, 42.7405414776084, 42.2707635255777, 36.002759645056, 30.6720357649246, 28.9584209198202, 22.9688169917519, 19.8985346257997, 16.2245011788703, 14.1537434356986, 10.9230374924946, 9.84084103968408, 9.74150559673712, 9.04625341131414, 4.60131725147045, 4.38460563123126, 4.25535754319832, 3.10269023446604, 2.05794148889311, 1.96486717013323, 1.60188936055615, 1.29298898564357, 1.00882158849929, 0.668918143856383, 0.436609405898277, 0.32783186823877, -39.3010815755573
#> 
#> Beginning fit attempt 9 of at maximum 10 extra tries
#> Running untitled13 with 41 parameters
#> 
#>  Lowest minimum so far:  3096.94452072231
#> 
#> Solution found
#> 


#> 
#>  Solution found!  Final fit=3096.9445 (started at 226053.29)  (10 attempt(s): 10 valid, 0 errors)
#>  Start values from best fit:
#> 2.17951987761656,1.81811323972771,1.44904166752005,0.604498514994528,1.29147147179014,1.17388124141934,1.30214898520524,0.898492748165293,1.13247286893006,1.2095781350104,1.91458653549528,7.40452642743505,4.99236397082348,1.32053675160167,3.15117526647134,2.17541393973224,5.01524171743309,0.08135248335329,0.120528462906145,0.466700510150328,0.590969841918841,2.30230216516881,0.731346168316912,3.52500932064948,0.353176392179849,1.41225146626662,3.32139966290673,0.448634263127453,3.71721712002252,0.164479917403726,5.46466666924916,4.25644288590711,6.56311025221454,4.45253304298167,2.97807408488281,5.05438384026312,4.7921946276716,3.55768978705078,5.13625192384593,6.1962638920508,4.04338968447061
omxGetParameters(fit)
#>    ind60→x2    ind60→x3 ind60→dem60 ind60→dem65          a1           b 
#>  2.17951988  1.81811324  1.44904167  0.60449851  1.29147147  1.17388124 
#>          c1 dem60→dem65          a2          c2       y1↔y1       y2↔y2 
#>  1.30214899  0.89849275  1.13247287  1.20957814  1.91458654  7.40452643 
#>       y3↔y3       y2↔y4       y4↔y4       y2↔y6       y6↔y6       x1↔x1 
#>  4.99236397  1.32053675  3.15117527  2.17541394  5.01524172  0.08135248 
#>       x2↔x2       x3↔x3       y1↔y5       y5↔y5       y3↔y7       y7↔y7 
#>  0.12052846  0.46670051  0.59096984  2.30230217  0.73134617  3.52500932 
#>       y4↔y8       y6↔y8       y8↔y8 ind60↔ind60 dem60↔dem60 dem65↔dem65 
#>  0.35317639  1.41225147  3.32139966  0.44863426  3.71721712  0.16447992 
#>      one→y1      one→y2      one→y3      one→y4      one→y6      one→x1 
#>  5.46466667  4.25644289  6.56311025  4.45253304  2.97807408  5.05438384 
#>      one→x2      one→x3      one→y5      one→y7      one→y8 
#>  4.79219463  3.55768979  5.13625192  6.19626389  4.04338968 


model_transformations <- '
  # latent variable definitions
     ind60 =~ x1 + x2 + x3
     dem60 =~ y1 + a1*y2 + b1*y3 + c1*y4
     dem65 =~ y5 + {a2 := a1 + delta_a}*y6 + {b2 := b1 + delta_b}*y7 + c2*y8

  # regressions
    dem60 ~ ind60
    dem65 ~ ind60 + dem60

  # residual correlations
    y1 ~~ y5
    y2 ~~ y4 + y6
    y3 ~~ y7
    y4 ~~ y8
    y6 ~~ y8
'

fit <- mxsem(model = model_transformations,
            data  = OpenMx::Bollen) |>
  mxTryHard()
#> Running untitled16 with 42 parameters
#> 
#> Beginning initial fit attempt
#> Running untitled16 with 42 parameters
#> 
#>  Lowest minimum so far:  3095.58188680662
#> 
#> Solution found
#> 


#> 
#>  Solution found!  Final fit=3095.5819 (started at 225616.79)  (1 attempt(s): 1 valid, 0 errors)
#>  Start values from best fit:
#> 2.18036735281242,1.8185114387558,1.48299990554601,0.572335346018354,1.2567470003478,1.0577157901271,1.26478686906097,0.837345248817426,1.2659469908364,1.89139993018466,7.37285349889048,5.06747633124564,1.31310790496682,3.14789891490124,2.15285745862293,4.95396914989032,0.0815493442403247,0.119806458733075,0.46670268393101,0.623670135827627,2.35097075483806,0.79496133883091,3.43137460759556,0.348227254504968,1.356174285358,3.25409453796309,0.448437394617253,3.95602585309738,0.172476656580233,5.46466693034299,4.25644300148278,6.56311103615158,4.45253380690424,2.97807331745104,5.05438379082776,4.79219440394796,3.55768962874225,5.13625181790305,6.19626364513509,4.04339168196963,-0.0710501581032182,0.221796133564053
omxGetParameters(fit)
#>    ind60→x2    ind60→x3 ind60→dem60 ind60→dem65          a1          b1 
#>  2.18036735  1.81851144  1.48299991  0.57233535  1.25674700  1.05771579 
#>          c1 dem60→dem65          c2       y1↔y1       y2↔y2       y3↔y3 
#>  1.26478687  0.83734525  1.26594699  1.89139993  7.37285350  5.06747633 
#>       y2↔y4       y4↔y4       y2↔y6       y6↔y6       x1↔x1       x2↔x2 
#>  1.31310790  3.14789891  2.15285746  4.95396915  0.08154934  0.11980646 
#>       x3↔x3       y1↔y5       y5↔y5       y3↔y7       y7↔y7       y4↔y8 
#>  0.46670268  0.62367014  2.35097075  0.79496134  3.43137461  0.34822725 
#>       y6↔y8       y8↔y8 ind60↔ind60 dem60↔dem60 dem65↔dem65      one→y1 
#>  1.35617429  3.25409454  0.44843739  3.95602585  0.17247666  5.46466693 
#>      one→y2      one→y3      one→y4      one→y6      one→x1      one→x2 
#>  4.25644300  6.56311104  4.45253381  2.97807332  5.05438379  4.79219440 
#>      one→x3      one→y5      one→y7      one→y8     delta_a     delta_b 
#>  3.55768963  5.13625182  6.19626365  4.04339168 -0.07105016  0.22179613