Let's say we are perusing the Journal of Statistical Software (who isn't?) and we find an article on a new package for structural equation modeling – lavaan. To obtain the package we go to the description file and see:
Package:lavaan
Title: Latent Variable Analysis
Version: 0.4-14
Author@R: person(“Yves”, “Rosseel”, email=“Yves.Rosseel@UGent.be”)
Author: Yves Rosseel Yves.Rosseel@UGent.be with contributions from
many people. See the lavaan website for a list of contributors.
Maintainer: Yves Rosseel Yves.Rosseel@UGent.be
Description: Fit a variety of latent variable models, including
confirmatory factor analysis, structural equation modeling and
latent growth curve models.
Depends: methods, R(>= 2.14.0)
Imports: stats4, stats, graphics
Suggests: psych, qgraph, quadprog, boot
License: GPL (>= 2)
LazyData: yes
URL: http://lavaan.org
Packaged: 2012-05-11 09:26:32 UTC; yves
Repository: CRAN
Date/Publication: 2012-05-11 12:46:59
Built: R 2.14.0; ; 2012-06-27 15:50:01 UTC; unix
The website is likely to provide a more readable description of the package and if it is available, it would be the best place to start. The R Reference Manual is technical and we will eventually need to bite the bullet and look at it, but there are some other things we can do first. The best way to evaluate a package is to go ahead and install it and experiment.
Look for a Demo, Vignette, Data or some samples for the lavaan package.
library(lavaan)
## This is lavaan 0.4-14
## lavaan is BETA software! Please report any bugs.
demo(package = "lavaan")
## no demos found
vignette(package = "lavaan")
## no vignettes found
data(package = "lavaan")
(Note: See the 'zoo' package for example of demo vs vignette)
If data is included with the package, the 'data' function will open a tab in the upper left pane with a very brief description of the data.
Although there are no demos or vignettes for lavaan, there are three sample datasets. Now it's time to RTFM.
Where to start??
Looking at the sample data will help to determine whether the package will handle your own data. The first sample dataset is Demo.growth. If your goal is to use lavaan to fit a growth curve model, you are on the right track. The documentation will tell you how to access the data and explain the variables.
Sample datasets are 'preloaded', but you might want to make a copy of the data in your workspace so that you can examine it with RStudio's data.
myDemo.growth<-Demo.growth
The description's See Also points you to another part of the documention on 'growth'
'growth' is a function to fit a growth curve model. There are many arguments to the function to examine and it is easier to do this if you skip down to locate the 'Examples' section.
str(growth)
## function (model = NULL, fixed.x = "default", orthogonal = FALSE, std.lv = FALSE,
## data = NULL, std.ov = FALSE, missing = "default", sample.cov = NULL,
## sample.mean = NULL, sample.nobs = NULL, group = NULL, group.label = NULL,
## group.equal = "", group.partial = "", constraints = "", estimator = "default",
## likelihood = "default", information = "default", se = "default",
## test = "default", bootstrap = 1000L, mimic = "default", representation = "default",
## do.fit = TRUE, control = list(), start = "default", verbose = FALSE,
## warn = TRUE, debug = FALSE)
## linear growth model with a time-varying covariate
model.syntax <- "
# intercept and slope with fixed coefficients
i =~ 1*t1 + 1*t2 + 1*t3 + 1*t4
s =~ 0*t1 + 1*t2 + 2*t3 + 3*t4
# regressions
i ~ x1 + x2
s ~ x1 + x2
# time-varying covariates
t1 ~ c1
t2 ~ c2
t3 ~ c3
t4 ~ c4
"
fit <- growth(model.syntax, data=Demo.growth)
summary(fit)
## lavaan (0.4-14) converged normally after 41 iterations
##
## Number of observations 400
##
## Estimator ML
## Minimum Function Chi-square 26.059
## Degrees of freedom 21
## P-value 0.204
##
## Parameter estimates:
##
## Information Expected
## Standard Errors Standard
##
## Estimate Std.err Z-value P(>|z|)
## Latent variables:
## i =~
## t1 1.000
## t2 1.000
## t3 1.000
## t4 1.000
## s =~
## t1 0.000
## t2 1.000
## t3 2.000
## t4 3.000
##
## Regressions:
## i ~
## x1 0.608 0.060 10.134 0.000
## x2 0.604 0.064 9.412 0.000
## s ~
## x1 0.262 0.029 9.198 0.000
## x2 0.522 0.031 17.083 0.000
## t1 ~
## c1 0.143 0.050 2.883 0.004
## t2 ~
## c2 0.289 0.046 6.295 0.000
## t3 ~
## c3 0.328 0.044 7.361 0.000
## t4 ~
## c4 0.330 0.058 5.655 0.000
##
## Covariances:
## i ~~
## s 0.075 0.040 1.855 0.064
##
## Intercepts:
## t1 0.000
## t2 0.000
## t3 0.000
## t4 0.000
## i 0.580 0.062 9.368 0.000
## s 0.958 0.029 32.552 0.000
##
## Variances:
## t1 0.580 0.080
## t2 0.596 0.054
## t3 0.481 0.055
## t4 0.535 0.098
## i 1.079 0.112
## s 0.224 0.027
##
Now you can begin to play around with the various arguments, for example mimic=“Mplus”
fit <- growth(model.syntax, data = Demo.growth, mimic = "Mplus")
summary(fit)
## lavaan (0.4-14) converged normally after 38 iterations
##
## Number of observations 400
##
## Number of missing patterns 1
##
## Estimator ML
## Minimum Function Chi-square 26.059
## Degrees of freedom 21
## P-value 0.204
##
## Parameter estimates:
##
## Information Observed
## Standard Errors Standard
##
## Estimate Std.err Z-value P(>|z|)
## Latent variables:
## i =~
## t1 1.000
## t2 1.000
## t3 1.000
## t4 1.000
## s =~
## t1 0.000
## t2 1.000
## t3 2.000
## t4 3.000
##
## Regressions:
## i ~
## x1 0.608 0.060 10.133 0.000
## x2 0.604 0.064 9.408 0.000
## s ~
## x1 0.262 0.029 9.194 0.000
## x2 0.522 0.031 17.077 0.000
## t1 ~
## c1 0.143 0.050 2.870 0.004
## t2 ~
## c2 0.289 0.046 6.294 0.000
## t3 ~
## c3 0.328 0.045 7.354 0.000
## t4 ~
## c4 0.330 0.059 5.646 0.000
##
## Covariances:
## i ~~
## s 0.075 0.040 1.879 0.060
##
## Intercepts:
## t1 0.000
## t2 0.000
## t3 0.000
## t4 0.000
## i 0.580 0.062 9.361 0.000
## s 0.958 0.029 32.492 0.000
##
## Variances:
## t1 0.580 0.079 7.309 0.000
## t2 0.596 0.054 10.964 0.000
## t3 0.481 0.055 8.786 0.000
## t4 0.535 0.096 5.573 0.000
## i 1.079 0.112 9.628 0.000
## s 0.224 0.026 8.542 0.000
##
Be sure to look at the Values section of the documentation to see that this function returns:
Value
An object of class lavaan, for which several methods are available including a summary method.
You will notice that after running the growth example two objects were created in your workspace:
| Values | (class) |
|---|---|
| fit | lavaan[1] |
| model.syntax | character[1] |
We can see that the model.syntax is simply the character string containing our model definition.
But what about fit? If we click on it in RStudio we get..
S4 object of class structure(“lavaan”, package =
“lavaan”)
So a 'lavaan' is a class created by the package lavaan.
Back to the documentation, we see a section for 'lavaan' and another for 'lavaan-class'. We want the class.
class(fit)
## [1] "lavaan"
## attr(,"package")
## [1] "lavaan"
We know about int, character, dataframes, but what is a 'lavaan'? We can explore it's structure.
str(fit)
## Formal class 'lavaan' [package "lavaan"] with 8 slots
## ..@ call : language lavaan(model = model.syntax, model.type = "growth", int.ov.free = FALSE, int.lv.free = TRUE, auto.fix.first = TRUE, auto.fix.single = TRUE, ...
## ..@ timing :List of 10
## .. ..$ InitOptions: Named num 0.005
## .. .. ..- attr(*, "names")= chr "elapsed"
## .. ..$ InitData : Named num 0.009
## .. .. ..- attr(*, "names")= chr "elapsed"
## .. ..$ ParTable : Named num 0.005
## .. .. ..- attr(*, "names")= chr "elapsed"
## .. ..$ Sample : Named num 0.005
## .. .. ..- attr(*, "names")= chr "elapsed"
## .. ..$ Start : Named num 0.001
## .. .. ..- attr(*, "names")= chr "elapsed"
## .. ..$ Model : Named num 0.006
## .. .. ..- attr(*, "names")= chr "elapsed"
## .. ..$ Estimate : Named num 0.085
## .. .. ..- attr(*, "names")= chr "elapsed"
## .. ..$ VCOV : Named num 0.075
## .. .. ..- attr(*, "names")= chr "elapsed"
## .. ..$ TEST : Named num 0.001
## .. .. ..- attr(*, "names")= chr "elapsed"
## .. ..$ total : Named num 0.196
## .. .. ..- attr(*, "names")= chr "elapsed"
## ..@ Options :List of 29
## .. ..$ model : chr "\n # intercept and slope with fixed coefficients\n i =~ 1*t1 + 1*t2 + 1*t3 + 1*t4\n s =~ 0*t1 + 1*t2 + 2*t3 + 3*t4\n\n "| __truncated__
## .. ..$ model.type : chr "growth"
## .. ..$ meanstructure : logi TRUE
## .. ..$ int.ov.free : logi FALSE
## .. ..$ int.lv.free : logi TRUE
## .. ..$ fixed.x : logi TRUE
## .. ..$ orthogonal : logi FALSE
## .. ..$ std.lv : logi FALSE
## .. ..$ auto.fix.first : logi TRUE
## .. ..$ auto.fix.single: logi TRUE
## .. ..$ auto.var : logi TRUE
## .. ..$ auto.cov.lv.x : logi TRUE
## .. ..$ auto.cov.y : logi TRUE
## .. ..$ missing : chr "ml"
## .. ..$ group.equal : chr(0)
## .. ..$ group.partial : chr(0)
## .. ..$ constraints : chr ""
## .. ..$ estimator : chr "ML"
## .. ..$ likelihood : chr "normal"
## .. ..$ information : chr "observed"
## .. ..$ se : chr "standard"
## .. ..$ test : chr "standard"
## .. ..$ bootstrap : int 1000
## .. ..$ mimic : chr "Mplus"
## .. ..$ representation : chr "LISREL"
## .. ..$ do.fit : logi TRUE
## .. ..$ verbose : logi FALSE
## .. ..$ warn : logi TRUE
## .. ..$ debug : logi FALSE
## ..@ ParTable :List of 12
## .. ..$ id : int [1:56] 1 2 3 4 5 6 7 8 9 10 ...
## .. ..$ lhs : chr [1:56] "i" "i" "i" "i" ...
## .. ..$ op : chr [1:56] "=~" "=~" "=~" "=~" ...
## .. ..$ rhs : chr [1:56] "t1" "t2" "t3" "t4" ...
## .. ..$ user : int [1:56] 1 1 1 1 1 1 1 1 1 1 ...
## .. ..$ group : int [1:56] 1 1 1 1 1 1 1 1 1 1 ...
## .. ..$ free : int [1:56] 0 0 0 0 0 0 0 0 1 2 ...
## .. ..$ ustart: num [1:56] 1 1 1 1 0 1 2 3 NA NA ...
## .. ..$ exo : int [1:56] 0 0 0 0 0 0 0 0 0 0 ...
## .. ..$ label : chr [1:56] "" "" "" "" ...
## .. ..$ eq.id : int [1:56] 0 0 0 0 0 0 0 0 0 0 ...
## .. ..$ unco : int [1:56] 0 0 0 0 0 0 0 0 1 2 ...
## ..@ Data :Formal class 'lavData' [package "lavaan"] with 13 slots
## .. .. ..@ data.type : chr "full"
## .. .. ..@ ngroups : int 1
## .. .. ..@ group.label: chr(0)
## .. .. ..@ std.ov : logi FALSE
## .. .. ..@ nobs :List of 1
## .. .. .. ..$ : int 400
## .. .. ..@ norig :List of 1
## .. .. .. ..$ : int 400
## .. .. ..@ ov.names :List of 1
## .. .. .. ..$ : chr [1:10] "t1" "t2" "t3" "t4" ...
## .. .. ..@ ov.types :List of 1
## .. .. .. ..$ : chr "numeric"
## .. .. ..@ ov.idx :List of 1
## .. .. .. ..$ : int [1:10] 1 2 3 4 5 6 7 8 9 10
## .. .. ..@ case.idx :List of 1
## .. .. .. ..$ : int [1:400] 1 2 3 4 5 6 7 8 9 10 ...
## .. .. ..@ missing : chr "ml"
## .. .. ..@ Mp :List of 1
## .. .. .. ..$ :List of 8
## .. .. .. .. ..$ nobs : int 400
## .. .. .. .. ..$ nvar : int 10
## .. .. .. .. ..$ coverage : num [1:10, 1:10] 1 1 1 1 1 1 1 1 1 1 ...
## .. .. .. .. ..$ id : chr [1:400] "0000000000" "0000000000" "0000000000" "0000000000" ...
## .. .. .. .. ..$ npatterns: int 1
## .. .. .. .. ..$ order : chr "0000000000"
## .. .. .. .. ..$ pat : logi [1, 1:10] TRUE TRUE TRUE TRUE TRUE TRUE ...
## .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. .. .. .. ..$ : chr "400"
## .. .. .. .. .. .. ..$ : NULL
## .. .. .. .. ..$ empty.idx: int(0)
## .. .. ..@ X :List of 1
## .. .. .. ..$ : num [1:400, 1:10] 1.726 -1.984 0.32 0.777 0.449 ...
## ..@ SampleStats:Formal class 'lavSampleStats' [package "lavaan"] with 12 slots
## .. .. ..@ cov :List of 1
## .. .. .. ..$ : num [1:10, 1:10] 2.494 2.58 3.081 3.771 0.732 ...
## .. .. ..@ mean :List of 1
## .. .. .. ..$ : num [1:10] 0.5947 1.6733 2.5932 3.6387 -0.0921 ...
## .. .. ..@ nobs :List of 1
## .. .. .. ..$ : int 400
## .. .. ..@ ntotal : int 400
## .. .. ..@ ngroups : int 1
## .. .. ..@ icov :List of 1
## .. .. .. ..$ : num [1:10, 1:10] 1.0569 -0.4821 -0.1198 -0.02 -0.0638 ...
## .. .. ..@ cov.log.det :List of 1
## .. .. .. ..$ : num 0.432
## .. .. ..@ cov.vecs :List of 1
## .. .. .. ..$ : num [1:55] 2.494 2.58 3.081 3.771 0.732 ...
## .. .. ..@ WLS.V :List of 1
## .. .. .. ..$ : NULL
## .. .. ..@ missing.flag: logi TRUE
## .. .. ..@ missing :List of 1
## .. .. .. ..$ :List of 1
## .. .. .. .. ..$ :List of 5
## .. .. .. .. .. ..$ X : num [1:400, 1:10] 1.726 -1.984 0.32 0.777 0.449 ...
## .. .. .. .. .. ..$ SX : num [1:10, 1:10] 2.494 2.58 3.081 3.771 0.732 ...
## .. .. .. .. .. ..$ MX : num [1:10] 0.5947 1.6733 2.5932 3.6387 -0.0921 ...
## .. .. .. .. .. ..$ nobs : int 400
## .. .. .. .. .. ..$ var.idx: logi [1:10] TRUE TRUE TRUE TRUE TRUE TRUE ...
## .. .. ..@ missing.h1 :List of 1
## .. .. .. ..$ :List of 3
## .. .. .. .. ..$ sigma: num [1:10, 1:10] 2.494 2.58 3.081 3.771 0.732 ...
## .. .. .. .. ..$ mu : num [1:10] 0.5947 1.6733 2.5932 3.6387 -0.0921 ...
## .. .. .. .. ..$ h1 : num 10.4
## ..@ Model :Formal class 'Model' [package "lavaan"] with 30 slots
## .. .. ..@ GLIST :List of 6
## .. .. .. ..$ lambda: num [1:10, 1:12] 0 0 0 0 0 0 0 0 0 0 ...
## .. .. .. ..$ theta : num [1:10, 1:10] 0 0 0 0 0 0 0 0 0 0 ...
## .. .. .. ..$ psi : num [1:12, 1:12] 1.0795 0.0748 0 0 0 ...
## .. .. .. ..$ beta : num [1:12, 1:12] 0 0 1 1 1 1 0 0 0 0 ...
## .. .. .. ..$ nu : num [1:10, 1] 0 0 0 0 0 0 0 0 0 0
## .. .. .. ..$ alpha : num [1:12, 1] 0.58 0.958 0 0 0 ...
## .. .. ..@ dimNames :List of 6
## .. .. .. ..$ :List of 2
## .. .. .. .. ..$ : chr [1:10] "t1" "t2" "t3" "t4" ...
## .. .. .. .. ..$ : chr [1:12] "i" "s" "t1" "t2" ...
## .. .. .. ..$ :List of 2
## .. .. .. .. ..$ : chr [1:10] "t1" "t2" "t3" "t4" ...
## .. .. .. .. ..$ : chr [1:10] "t1" "t2" "t3" "t4" ...
## .. .. .. ..$ :List of 2
## .. .. .. .. ..$ : chr [1:12] "i" "s" "t1" "t2" ...
## .. .. .. .. ..$ : chr [1:12] "i" "s" "t1" "t2" ...
## .. .. .. ..$ :List of 2
## .. .. .. .. ..$ : chr [1:12] "i" "s" "t1" "t2" ...
## .. .. .. .. ..$ : chr [1:12] "i" "s" "t1" "t2" ...
## .. .. .. ..$ :List of 2
## .. .. .. .. ..$ : chr [1:10] "t1" "t2" "t3" "t4" ...
## .. .. .. .. ..$ : chr "intercept"
## .. .. .. ..$ :List of 2
## .. .. .. .. ..$ : chr [1:12] "i" "s" "t1" "t2" ...
## .. .. .. .. ..$ : chr "intercept"
## .. .. ..@ isSymmetric : logi [1:6] FALSE TRUE TRUE FALSE FALSE FALSE
## .. .. ..@ mmSize : int [1:6] 120 55 78 144 10 12
## .. .. ..@ representation : chr "LISREL"
## .. .. ..@ meanstructure : logi TRUE
## .. .. ..@ ngroups : int 1
## .. .. ..@ nmat : int 6
## .. .. ..@ nvar : int 10
## .. .. ..@ nx.free : int 17
## .. .. ..@ nx.unco : int 17
## .. .. ..@ nx.user : int 56
## .. .. ..@ m.free.idx :List of 6
## .. .. .. ..$ : int(0)
## .. .. .. ..$ : int(0)
## .. .. .. ..$ : int [1:8] 1 2 13 14 27 40 53 66
## .. .. .. ..$ : int [1:8] 73 74 85 86 99 112 125 138
## .. .. .. ..$ : int(0)
## .. .. .. ..$ : int [1:2] 1 2
## .. .. ..@ x.free.idx :List of 6
## .. .. .. ..$ : int(0)
## .. .. .. ..$ : int(0)
## .. .. .. ..$ : int [1:8] 13 15 15 14 9 10 11 12
## .. .. .. ..$ : int [1:8] 1 3 2 4 5 6 7 8
## .. .. .. ..$ : int(0)
## .. .. .. ..$ : int [1:2] 16 17
## .. .. ..@ m.unco.idx :List of 6
## .. .. .. ..$ : int(0)
## .. .. .. ..$ : int(0)
## .. .. .. ..$ : int [1:8] 1 2 13 14 27 40 53 66
## .. .. .. ..$ : int [1:8] 73 74 85 86 99 112 125 138
## .. .. .. ..$ : int(0)
## .. .. .. ..$ : int [1:2] 1 2
## .. .. ..@ x.unco.idx :List of 6
## .. .. .. ..$ : int(0)
## .. .. .. ..$ : int(0)
## .. .. .. ..$ : int [1:8] 13 15 15 14 9 10 11 12
## .. .. .. ..$ : int [1:8] 1 3 2 4 5 6 7 8
## .. .. .. ..$ : int(0)
## .. .. .. ..$ : int [1:2] 16 17
## .. .. ..@ m.user.idx :List of 6
## .. .. .. ..$ : int(0)
## .. .. .. ..$ : int(0)
## .. .. .. ..$ : int [1:44] 1 2 13 14 27 40 53 66 79 80 ...
## .. .. .. ..$ : int [1:16] 3 4 5 6 15 16 17 18 73 74 ...
## .. .. .. ..$ : int(0)
## .. .. .. ..$ : int [1:12] 1 2 3 4 5 6 7 8 9 10 ...
## .. .. ..@ x.user.idx :List of 6
## .. .. .. ..$ : int(0)
## .. .. .. ..$ : int(0)
## .. .. .. ..$ : int [1:44] 21 23 23 22 17 18 19 20 24 25 ...
## .. .. .. ..$ : int [1:16] 1 2 3 4 5 6 7 8 9 11 ...
## .. .. .. ..$ : int(0)
## .. .. .. ..$ : int [1:12] 55 56 45 46 47 48 49 50 51 52 ...
## .. .. ..@ x.def.idx : int(0)
## .. .. ..@ x.ceq.idx : int(0)
## .. .. ..@ x.cin.idx : int(0)
## .. .. ..@ eq.constraints : logi FALSE
## .. .. ..@ eq.constraints.K: num[0 , 0 ]
## .. .. ..@ def.function :function ()
## .. .. ..@ ceq.function :function ()
## .. .. ..@ ceq.jacobian :function ()
## .. .. ..@ cin.function :function ()
## .. .. ..@ cin.jacobian :function ()
## .. .. ..@ con.jac : num[0 , 0 ]
## .. .. ..@ fixed.x : logi TRUE
## ..@ Fit :Formal class 'Fit' [package "lavaan"] with 13 slots
## .. .. ..@ npar : int 17
## .. .. ..@ x : num [1:17] 0.608 0.604 0.262 0.522 0.143 ...
## .. .. ..@ start : num [1:56] 1 1 1 1 0 1 2 3 0 0 ...
## .. .. ..@ est : num [1:56] 1 1 1 1 0 ...
## .. .. ..@ se : num [1:56] 0 0 0 0 0 ...
## .. .. ..@ fx : num 0.0326
## .. .. ..@ fx.group : num 0.0326
## .. .. ..@ iterations: int 38
## .. .. ..@ converged : logi TRUE
## .. .. ..@ control :List of 7
## .. .. .. ..$ eval.max: int 20000
## .. .. .. ..$ iter.max: int 10000
## .. .. .. ..$ trace : int 0
## .. .. .. ..$ abs.tol : num 2.22e-15
## .. .. .. ..$ rel.tol : num 1e-10
## .. .. .. ..$ x.tol : num 1.5e-08
## .. .. .. ..$ step.min: num 2.2e-14
## .. .. ..@ Sigma.hat :List of 1
## .. .. .. ..$ : num [1:10, 1:10] 2.51 2.526 3.158 3.709 0.729 ...
## .. .. ..@ Mu.hat :List of 1
## .. .. .. ..$ : num [1:10, 1] 0.6089 1.6217 2.6411 3.6186 -0.0921 ...
## .. .. ..@ test :List of 1
## .. .. .. ..$ :List of 5
## .. .. .. .. ..$ test : chr "standard"
## .. .. .. .. ..$ stat : num 26.1
## .. .. .. .. ..$ stat.group: num 26.1
## .. .. .. .. ..$ df : int 21
## .. .. .. .. ..$ pvalue : num 0.204
Hmmn….
The documentation will explain what some of this is.
Some standard methods will work with custom classes, like the summary method we used above. Look for accessor-type methods that will extract values from the class.
logLik(fit)
## 'log Lik.' -5775 (df=17)
fitted.values(fit)
## $cov
## t1 t2 t3 t4 x1 x2 c1 c2 c3 c4
## t1 2.510
## t2 2.526 4.403
## t3 3.158 4.990 7.465
## t4 3.709 6.083 8.553 11.331
## x1 0.729 1.099 1.459 1.797 1.056
## x2 0.649 1.171 1.743 2.184 0.153 0.924
## c1 0.105 -0.031 -0.084 -0.052 -0.038 -0.019 0.972
## c2 0.023 0.274 0.018 0.031 0.026 -0.007 0.080 0.899
## c3 0.104 0.193 0.563 0.373 0.033 0.145 -0.030 0.001 0.872
## c4 -0.058 -0.127 -0.184 0.026 -0.025 -0.102 0.128 0.032 0.035 0.851
##
## $mean
## t1 t2 t3 t4 x1 x2 c1 c2 c3 c4
## 0.609 1.622 2.641 3.619 -0.092 0.138 0.008 0.029 0.068 -0.018
##
# sometimes there is more than one way to get the same values
partable1 <- parameterTable(fit)
partable2 <- inspect(fit, "list")
Putting the results in a workspace object will allow you to access all the values using standard naming.
partable1[, "id"]
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
## [24] 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46
## [47] 47 48 49 50 51 52 53 54 55 56
Methods may do things like plot or perform calculations on the values stored in the class.
Visiting the [lavaan.org website][http://lavaan.ugent.be/?q=node/16] we are provided with code to replicate several sample programs included with the package Mplus. Let's try to run an Mplus sample analysis using R
# Mplus example 5.3
Data <- read.table("ex5.3.dat")
names(Data) <- c("u1", "u2", "u3", "y4", "y5", "y6")
Data$u1 <- ordered(Data$u1)
Data$u2 <- ordered(Data$u2)
Data$u3 <- ordered(Data$u3)
model <- " f1 =~ u1 + u2 + u3\nf2 =~ y4 + y5 + y6 "
fit <- cfa(model, data = Data, mimic = "Mplus")
summary(fit, standardized = TRUE, fit.measures = TRUE, rsquare = TRUE)
## lavaan (0.4-14) converged normally after 4 iterations
##
## Number of observations 500
##
## Number of missing patterns 1
##
## Estimator ML
## Minimum Function Chi-square 0.000
## Degrees of freedom 8
## P-value 1.000
##
## Chi-square test baseline model:
##
## Minimum Function Chi-square 939.551
## Degrees of freedom 15
## P-value 0.000
##
## Full model versus baseline model:
##
## Comparative Fit Index (CFI) 1.000
## Tucker-Lewis Index (TLI) 1.016
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -2756.816
## Loglikelihood unrestricted model (H1) -2756.816
##
## Number of free parameters 19
## Akaike (AIC) 5551.631
## Bayesian (BIC) 5631.709
## Sample-size adjusted Bayesian (BIC) 5571.402
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.000
## 90 Percent Confidence Interval 0.000 0.000
## P-value RMSEA <= 0.05 1.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.271
##
## Parameter estimates:
##
## Information Observed
## Standard Errors Standard
##
## Estimate Std.err Z-value P(>|z|) Std.lv Std.all
## Latent variables:
## f1 =~
## u1 1.000 0.284 0.638
## u2 1.104 0.078 14.077 0.000 0.314 0.686
## u3 1.003 0.077 13.082 0.000 0.285 0.639
## f2 =~
## y4 1.000 0.374 0.462
## y5 1.056 0.084 12.558 0.000 0.395 0.482
## y6 1.002 0.087 11.481 0.000 0.375 0.469
##
## Covariances:
## f1 ~~
## f2 0.000 0.005 0.102 0.919 0.004 0.004
##
## Intercepts:
## u1 1.496 0.020 75.027 0.000 1.496 3.355
## u2 1.494 0.020 73.067 0.000 1.494 3.268
## u3 1.504 0.020 75.352 0.000 1.504 3.370
## y4 -0.016 0.036 -0.438 0.661 -0.016 -0.020
## y5 -0.016 0.037 -0.444 0.657 -0.016 -0.020
## y6 0.014 0.036 0.388 0.698 0.014 0.017
## f1 0.000 0.000 0.000
## f2 0.000 0.000 0.000
##
## Variances:
## u1 0.118 0.011 10.456 0.000 0.118 0.593
## u2 0.111 0.016 6.980 0.000 0.111 0.529
## u3 0.118 0.013 9.117 0.000 0.118 0.592
## y4 0.516 0.036 14.523 0.000 0.516 0.787
## y5 0.515 0.035 14.688 0.000 0.515 0.768
## y6 0.497 0.038 13.100 0.000 0.497 0.780
## f1 0.081 0.007 12.008 0.000 1.000 1.000
## f2 0.140 0.007 20.093 0.000 1.000 1.000
##
## R-Square:
##
## u1 0.407
## u2 0.471
## u3 0.408
## y4 0.213
## y5 0.232
## y6 0.220
Use the various methods to explore your results.
inspect(fit, "rsquare")
## u1 u2 u3 y4 y5 y6
## 0.4067 0.4711 0.4082 0.2133 0.2324 0.2203
# Inspecting fit gives
inspect(fit, "fit")
## chisq df pvalue baseline.chisq
## 0.000 8.000 1.000 939.551
## baseline.df baseline.pvalue cfi tli
## 15.000 0.000 1.000 1.016
## logl unrestricted.logl npar aic
## -2756.816 -2756.816 19.000 5551.631
## bic ntotal bic2 rmsea
## 5631.709 500.000 5571.402 0.000
## rmsea.ci.lower rmsea.ci.upper rmsea.pvalue srmr
## 0.000 0.000 1.000 0.271
# We can pull just the BIC from the fit measures using
fitMeasures(fit, "bic")
## bic
## 5632
fitted.values(fit)
## $cov
## u1 u2 u3 y4 y5 y6
## u1 0.199
## u2 0.089 0.209
## u3 0.081 0.089 0.199
## y4 0.000 0.001 0.000 0.655
## y5 0.001 0.001 0.001 0.148 0.671
## y6 0.000 0.001 0.000 0.140 0.148 0.637
##
## $mean
## u1 u2 u3 y4 y5 y6
## 1.496 1.494 1.504 -0.016 -0.016 0.014
##
coef(fit)
## f1=~u2 f1=~u3 f2=~y5 f2=~y6 u1~~u1 u2~~u2 u3~~u3 y4~~y4 y5~~y5 y6~~y6
## 1.104 1.003 1.056 1.002 0.118 0.111 0.118 0.516 0.515 0.497
## f1~~f1 f2~~f2 f1~~f2 u1~1 u2~1 u3~1 y4~1 y5~1 y6~1
## 0.081 0.140 0.000 1.496 1.494 1.504 -0.016 -0.016 0.014
resid(fit, type = "normalized")
## $cov
## u1 u2 u3 y4 y5 y6
## u1 3.238
## u2 5.526 2.589
## u3 5.151 5.509 3.212
## y4 0.220 -0.028 0.198 5.701
## y5 -0.081 -0.658 -0.529 11.380 5.468
## y6 0.324 -0.312 0.092 11.144 11.473 5.613
##
## $mean
## u1 u2 u3 y4 y5 y6
## 0 0 0 0 0 0
##
On the first page of the documentation, there was a line “Suggests: psych, qgraph”. This is telling you about other packages that may be aware of lavaan and are able to utilize it in some way.
The psych package has lavaan.diagram method and the
qgraph package can also work with lavaan objects
library(psych)
lavaan.diagram(fit, title = "Psych package lavaan.diagram")
library(qgraph)
# Plot standardized model (numerical):
qgraph.lavaan(fit, layout = "tree", filetype = "", include = 2)
Have fun exploring new R packages!