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 )
| dsClean | The base::data.frame containing complete cases for the |
|---|---|
| estimateA | Should the A variance component be estimated? A^2 represents the proportion of variability due to a shared genetic influence. |
| estimateC | Should the C variance component be estimated? C^2 represents the proportion of variability due to a shared environmental influence. |
| printOutput | 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.
Will Beasley
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.6219253 0.2097338 0.1683408 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.#> #>#> lavaan 0.6-9 ended normally after 54 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 df #> 4.000 0.027 447.241 16.000 #> pvalue baseline.chisq baseline.df baseline.pvalue #> 0.000 2106.324 4.000 0.000 #> cfi tli nnfi rfi #> 0.795 0.949 0.949 NA #> nfi pnfi ifi rni #> NA 3.151 0.794 0.795 #> logl unrestricted.logl aic bic #> -65103.779 -64880.158 130215.557 130243.671 #> ntotal bic2 rmsea rmsea.ci.lower #> 8338.000 130230.960 0.114 0.105 #> rmsea.ci.upper rmsea.pvalue rmr rmr_nomean #> 0.123 0.000 9.992 12.765 #> srmr srmr_bentler srmr_bentler_nomean crmr #> 0.130 0.130 0.089 0.178 #> crmr_nomean srmr_mplus srmr_mplus_nomean cn_05 #> 0.026 0.153 0.083 491.246 #> cn_01 gfi agfi pgfi #> 597.581 0.999 0.999 0.799 #> mfi #> 0.974#Examine low-level details like each group's individual parameter estimates and standard errors. summary(GetDetails(ace))#> lavaan 0.6-9 ended normally after 54 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 #> E1 0.000 #> E2 0.000 #> A1 0.000 #> A2 0.000 #> C1 0.000 #> C2 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 #> E1 0.000 #> E2 0.000 #> A1 0.000 #> A2 0.000 #> C1 0.000 #> C2 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 #> E1 0.000 #> E2 0.000 #> A1 0.000 #> A2 0.000 #> C1 0.000 #> C2 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 #> E1 0.000 #> E2 0.000 #> A1 0.000 #> A2 0.000 #> C1 0.000 #> C2 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#> $`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 #> #>