R/ace-lavaan-group.R
AceLavaanGroup.RdThis function uses the lavaan package to estimate a univariate ACE model, using multiple groups.
Each group has a unique value of R (i.e., the Relatedness coefficient).
AceLavaanGroup(
dsClean,
estimateA = TRUE,
estimateC = TRUE,
printOutput = FALSE
)The base::data.frame containing complete cases for the R groups to be included in the estimation.
Should the A variance component be estimated? A^2 represents the proportion of variability due to a shared genetic influence.
Should the C variance component be estimated? C^2 represents the proportion of variability due to a shared environmental influence.
Indicates if the estimated parameters and fit statistics are printed to the console.
An AceEstimate object.
The variance component for E is always estimated, while the A and C estimates can be fixed to zero (when estimateA and/or estimateC are set to FALSE).
Currently, the variables in dsClean must be named O1, O2 and R; the letter 'O' stands for Outcome. This may not be as restrictive as it initially seems, because dsClean is intended to be produced by CleanSemAceDataset(). If this is too restrictive for your uses, we'd like to here about it (please email wibeasley at hotmail period com).
The lavaan package is developed by Yves Rosseel at Ghent University. Three good starting points are the package website (http://lavaan.ugent.be/), the package documentation (https://cran.r-project.org/package=lavaan) and the JSS paper.
Rosseel, Yves (2012), lavaan: An R Package for Structural Equation Modeling. Journal of Statistical Software, 48, (2), 1-36.
CleanSemAceDataset(). Further ACE model details are discussed in our package's vignettes.
library(NlsyLinks) # Load the package into the current R session.
dsLinks <- Links79PairExpanded # Start with the built-in data.frame in NlsyLinks
dsLinks <- dsLinks[dsLinks$RelationshipPath == "Gen2Siblings", ] # Use only Gen2 Siblings (NLSY79-C)
oName_S1 <- "MathStandardized_S1" # Stands for Outcome1
oName_S2 <- "MathStandardized_S2" # Stands for Outcome2
dsGroupSummary <- RGroupSummary(dsLinks, oName_S1, oName_S2)
dsClean <- CleanSemAceDataset(dsDirty = dsLinks, dsGroupSummary, oName_S1, oName_S2)
ace <- AceLavaanGroup(dsClean)
ace
#> [1] "Results of ACE estimation: [show]"
#> ASquared CSquared ESquared CaseCount
#> 0.6219254 0.2097338 0.1683407 8338.0000000
# Should produce:
# [1] "Results of ACE estimation: [show]"
# ASquared CSquared ESquared CaseCount
# 0.6681874 0.1181227 0.2136900 8390.0000000
library(lavaan) # Load the package to access methods of the lavaan class.
#> This is lavaan 0.6-19
#> lavaan is FREE software! Please report any bugs.
GetDetails(ace)
#> lavaan 0.6-19 ended normally after 39 iterations
#>
#> Estimator ML
#> Optimization method NLMINB
#> Number of model parameters 32
#> Number of equality constraints 28
#>
#> Number of observations per group:
#> 1 2689
#> 2 137
#> 3 5491
#> 4 21
#>
#> Model Test User Model:
#>
#> Test statistic 447.241
#> Degrees of freedom 16
#> P-value (Chi-square) 0.000
#> Test statistic for each group:
#> 1 281.866
#> 2 30.277
#> 3 127.671
#> 4 7.428
# Exmaine fit stats like Chi-Squared, RMSEA, CFI, etc.
fitMeasures(GetDetails(ace)) # The function 'fitMeasures' is defined in the lavaan package.
#> npar fmin chisq
#> 4.000 0.027 447.241
#> df pvalue baseline.chisq
#> 16.000 0.000 2106.324
#> baseline.df baseline.pvalue cfi
#> 4.000 0.000 0.795
#> tli nnfi rfi
#> 0.949 0.949 NA
#> nfi pnfi ifi
#> NA 3.151 0.794
#> rni logl unrestricted.logl
#> 0.795 -65103.779 -64880.158
#> aic bic ntotal
#> 130215.557 130243.671 8338.000
#> bic2 rmsea rmsea.ci.lower
#> 130230.960 0.114 0.105
#> rmsea.ci.upper rmsea.ci.level rmsea.pvalue
#> 0.123 0.900 0.000
#> rmsea.close.h0 rmsea.notclose.pvalue rmsea.notclose.h0
#> 0.050 1.000 0.080
#> rmr rmr_nomean srmr
#> 9.992 12.765 0.130
#> srmr_bentler srmr_bentler_nomean crmr
#> 0.130 0.089 0.138
#> crmr_nomean srmr_mplus srmr_mplus_nomean
#> 0.026 0.153 0.083
#> cn_05 cn_01 gfi
#> 491.246 597.581 0.999
#> agfi pgfi mfi
#> 0.999 0.799 0.974
#> ecvi
#> 0.055
# Examine low-level details like each group's individual parameter estimates and standard errors.
summary(GetDetails(ace))
#> lavaan 0.6-19 ended normally after 39 iterations
#>
#> Estimator ML
#> Optimization method NLMINB
#> Number of model parameters 32
#> Number of equality constraints 28
#>
#> Number of observations per group:
#> 1 2689
#> 2 137
#> 3 5491
#> 4 21
#>
#> Model Test User Model:
#>
#> Test statistic 447.241
#> Degrees of freedom 16
#> P-value (Chi-square) 0.000
#> Test statistic for each group:
#> 1 281.866
#> 2 30.277
#> 3 127.671
#> 4 7.428
#>
#> Parameter Estimates:
#>
#> Standard errors Standard
#> Information Observed
#> Observed information based on Hessian
#>
#>
#> Group 1 [1]:
#>
#> Latent Variables:
#> Estimate Std.Err z-value P(>|z|)
#> E1 =~
#> O1 (e) 5.257 0.539 9.754 0.000
#> E2 =~
#> O2 (e) 5.257 0.539 9.754 0.000
#> A1 =~
#> O1 (a) 10.103 0.507 19.915 0.000
#> A2 =~
#> O2 (a) 10.103 0.507 19.915 0.000
#> C1 =~
#> O1 (c) 5.867 0.425 13.820 0.000
#> C2 =~
#> O2 (c) 5.867 0.425 13.820 0.000
#>
#> Covariances:
#> Estimate Std.Err z-value P(>|z|)
#> .O1 ~~
#> .O2 0.000
#> E1 ~~
#> E2 0.000
#> A1 ~~
#> A2 0.250
#> C1 ~~
#> C2 1.000
#>
#> Intercepts:
#> Estimate Std.Err z-value P(>|z|)
#> .O1 (int) 98.317 0.121 814.898 0.000
#> .O2 (int) 98.317 0.121 814.898 0.000
#>
#> Variances:
#> Estimate Std.Err z-value P(>|z|)
#> E1 1.000
#> E2 1.000
#> A1 1.000
#> A2 1.000
#> C1 1.000
#> C2 1.000
#> .O1 0.000
#> .O2 0.000
#>
#>
#> Group 2 [2]:
#>
#> Latent Variables:
#> Estimate Std.Err z-value P(>|z|)
#> E1 =~
#> O1 (e) 5.257 0.539 9.754 0.000
#> E2 =~
#> O2 (e) 5.257 0.539 9.754 0.000
#> A1 =~
#> O1 (a) 10.103 0.507 19.915 0.000
#> A2 =~
#> O2 (a) 10.103 0.507 19.915 0.000
#> C1 =~
#> O1 (c) 5.867 0.425 13.820 0.000
#> C2 =~
#> O2 (c) 5.867 0.425 13.820 0.000
#>
#> Covariances:
#> Estimate Std.Err z-value P(>|z|)
#> .O1 ~~
#> .O2 0.000
#> E1 ~~
#> E2 0.000
#> A1 ~~
#> A2 0.375
#> C1 ~~
#> C2 1.000
#>
#> Intercepts:
#> Estimate Std.Err z-value P(>|z|)
#> .O1 (int) 98.317 0.121 814.898 0.000
#> .O2 (int) 98.317 0.121 814.898 0.000
#>
#> Variances:
#> Estimate Std.Err z-value P(>|z|)
#> E1 1.000
#> E2 1.000
#> A1 1.000
#> A2 1.000
#> C1 1.000
#> C2 1.000
#> .O1 0.000
#> .O2 0.000
#>
#>
#> Group 3 [3]:
#>
#> Latent Variables:
#> Estimate Std.Err z-value P(>|z|)
#> E1 =~
#> O1 (e) 5.257 0.539 9.754 0.000
#> E2 =~
#> O2 (e) 5.257 0.539 9.754 0.000
#> A1 =~
#> O1 (a) 10.103 0.507 19.915 0.000
#> A2 =~
#> O2 (a) 10.103 0.507 19.915 0.000
#> C1 =~
#> O1 (c) 5.867 0.425 13.820 0.000
#> C2 =~
#> O2 (c) 5.867 0.425 13.820 0.000
#>
#> Covariances:
#> Estimate Std.Err z-value P(>|z|)
#> .O1 ~~
#> .O2 0.000
#> E1 ~~
#> E2 0.000
#> A1 ~~
#> A2 0.500
#> C1 ~~
#> C2 1.000
#>
#> Intercepts:
#> Estimate Std.Err z-value P(>|z|)
#> .O1 (int) 98.317 0.121 814.898 0.000
#> .O2 (int) 98.317 0.121 814.898 0.000
#>
#> Variances:
#> Estimate Std.Err z-value P(>|z|)
#> E1 1.000
#> E2 1.000
#> A1 1.000
#> A2 1.000
#> C1 1.000
#> C2 1.000
#> .O1 0.000
#> .O2 0.000
#>
#>
#> Group 4 [4]:
#>
#> Latent Variables:
#> Estimate Std.Err z-value P(>|z|)
#> E1 =~
#> O1 (e) 5.257 0.539 9.754 0.000
#> E2 =~
#> O2 (e) 5.257 0.539 9.754 0.000
#> A1 =~
#> O1 (a) 10.103 0.507 19.915 0.000
#> A2 =~
#> O2 (a) 10.103 0.507 19.915 0.000
#> C1 =~
#> O1 (c) 5.867 0.425 13.820 0.000
#> C2 =~
#> O2 (c) 5.867 0.425 13.820 0.000
#>
#> Covariances:
#> Estimate Std.Err z-value P(>|z|)
#> .O1 ~~
#> .O2 0.000
#> E1 ~~
#> E2 0.000
#> A1 ~~
#> A2 1.000
#> C1 ~~
#> C2 1.000
#>
#> Intercepts:
#> Estimate Std.Err z-value P(>|z|)
#> .O1 (int) 98.317 0.121 814.898 0.000
#> .O2 (int) 98.317 0.121 814.898 0.000
#>
#> Variances:
#> Estimate Std.Err z-value P(>|z|)
#> E1 1.000
#> E2 1.000
#> A1 1.000
#> A2 1.000
#> C1 1.000
#> C2 1.000
#> .O1 0.000
#> .O2 0.000
#>
#> Defined Parameters:
#> Estimate Std.Err z-value P(>|z|)
#> e2 27.631 5.665 4.877 0.000
#> a2 102.081 10.252 9.958 0.000
#> c2 34.425 4.982 6.910 0.000
#>
# Extract low-level details. This may be useful when programming simulations.
inspect(GetDetails(ace), what = "converged") # The lavaan package defines 'inspect'.
#> [1] TRUE
inspect(GetDetails(ace), what = "coef")
#> $`1`
#> $`1`$lambda
#> E1 E2 A1 A2 C1 C2
#> O1 5.257 0.000 10.103 0.000 5.867 0.000
#> O2 0.000 5.257 0.000 10.103 0.000 5.867
#>
#> $`1`$theta
#> O1 O2
#> O1 0
#> O2 0 0
#>
#> $`1`$psi
#> E1 E2 A1 A2 C1 C2
#> E1 1.00
#> E2 0.00 1.00
#> A1 0.00 0.00 1.00
#> A2 0.00 0.00 0.25 1.00
#> C1 0.00 0.00 0.00 0.00 1.00
#> C2 0.00 0.00 0.00 0.00 1.00 1.00
#>
#> $`1`$nu
#> intrcp
#> O1 98.317
#> O2 98.317
#>
#> $`1`$alpha
#> intrcp
#> E1 0
#> E2 0
#> A1 0
#> A2 0
#> C1 0
#> C2 0
#>
#>
#> $`2`
#> $`2`$lambda
#> E1 E2 A1 A2 C1 C2
#> O1 5.257 0.000 10.103 0.000 5.867 0.000
#> O2 0.000 5.257 0.000 10.103 0.000 5.867
#>
#> $`2`$theta
#> O1 O2
#> O1 0
#> O2 0 0
#>
#> $`2`$psi
#> E1 E2 A1 A2 C1 C2
#> E1 1.000
#> E2 0.000 1.000
#> A1 0.000 0.000 1.000
#> A2 0.000 0.000 0.375 1.000
#> C1 0.000 0.000 0.000 0.000 1.000
#> C2 0.000 0.000 0.000 0.000 1.000 1.000
#>
#> $`2`$nu
#> intrcp
#> O1 98.317
#> O2 98.317
#>
#> $`2`$alpha
#> intrcp
#> E1 0
#> E2 0
#> A1 0
#> A2 0
#> C1 0
#> C2 0
#>
#>
#> $`3`
#> $`3`$lambda
#> E1 E2 A1 A2 C1 C2
#> O1 5.257 0.000 10.103 0.000 5.867 0.000
#> O2 0.000 5.257 0.000 10.103 0.000 5.867
#>
#> $`3`$theta
#> O1 O2
#> O1 0
#> O2 0 0
#>
#> $`3`$psi
#> E1 E2 A1 A2 C1 C2
#> E1 1.0
#> E2 0.0 1.0
#> A1 0.0 0.0 1.0
#> A2 0.0 0.0 0.5 1.0
#> C1 0.0 0.0 0.0 0.0 1.0
#> C2 0.0 0.0 0.0 0.0 1.0 1.0
#>
#> $`3`$nu
#> intrcp
#> O1 98.317
#> O2 98.317
#>
#> $`3`$alpha
#> intrcp
#> E1 0
#> E2 0
#> A1 0
#> A2 0
#> C1 0
#> C2 0
#>
#>
#> $`4`
#> $`4`$lambda
#> E1 E2 A1 A2 C1 C2
#> O1 5.257 0.000 10.103 0.000 5.867 0.000
#> O2 0.000 5.257 0.000 10.103 0.000 5.867
#>
#> $`4`$theta
#> O1 O2
#> O1 0
#> O2 0 0
#>
#> $`4`$psi
#> E1 E2 A1 A2 C1 C2
#> E1 1
#> E2 0 1
#> A1 0 0 1
#> A2 0 0 1 1
#> C1 0 0 0 0 1
#> C2 0 0 0 0 1 1
#>
#> $`4`$nu
#> intrcp
#> O1 98.317
#> O2 98.317
#>
#> $`4`$alpha
#> intrcp
#> E1 0
#> E2 0
#> A1 0
#> A2 0
#> C1 0
#> C2 0
#>
#>