Twin analysis using TwinAnalysis

Ivan Voronin

January 12, 2021

Contents


The R code from this presentation is here: https://raw.githubusercontent.com/IvanVoronin/TwinAnalysis/master/slides/intro/intro.R

Twin method: Fundamentals

ACE

ADE

Univariate ACE

Univariate ACE

Univariate ACE

Univariate ACE

Univariate ACE

TwinAnalysis

Univariate toy data

# Making univariate toy data
library(dplyr,
        quietly = TRUE)

mz_data <- MASS::mvrnorm(
  n = 100,
  mu = c(0, 0),
  Sigma = matrix(c(1.0, 0.8,
                   0.8, 1.0),
                 nrow = 2)) %>%
  as.data.frame() %>%
  setNames(c('X1', 'X2'))

dz_data <- MASS::mvrnorm(
  n = 100,
  mu = c(0, 0),
  Sigma = matrix(c(1.0, 0.5,
                   0.5, 1.0),
                 nrow = 2)) %>%
  as.data.frame() %>%
  setNames(c('X1', 'X2'))

data <- rbind(data.frame(zyg = 'MZ', mz_data),
              data.frame(zyg = 'DZ', dz_data)) %>%
  mutate(
    zyg = factor(zyg, levels = c('MZ', 'DZ'))
  )

head(data)
##   zyg          X1          X2
## 1  MZ -0.72849176 -0.71698078
## 2  MZ  1.58692213 -0.09734755
## 3  MZ -0.46553953 -0.53837313
## 4  MZ  1.85518320  1.29294187
## 5  MZ  0.01350124 -0.60606016
## 6  MZ -1.26056323 -0.79239271

Descriptives

# Installing the packages
# install.packages('devtools')
# library(devtools)
# install_github('ivanvoronin/mlth.data.frame')
# install_github('ivanvoronin/TwinAnalysis')

# Loading the packages
library(OpenMx)
library(TwinAnalysis)

# Descriptive statistics
get_univ_descriptives(data, zyg = 'zyg', vars = 'X')
## $X
##    Twin1----------------- Twin2------------------- Cor--------------------------
##    n   M         SD       n   M          SD        r         lowerCI   upperCI  
## MZ 100 0.1148917 1.014334 100 0.07874427 0.9598355 0.7910906 0.7040748 0.8547013
## DZ 100 0.1994193 1.056857 100 0.07840418 1.1446936 0.6063534 0.4653590 0.7173375

Univariate ACE

# Univariate ACE model
model <- univariate_ace(data, zyg = 'zyg', ph = 'X')

model <- mxRun(model)

summary(model)
## Summary of ACE 
##  
## free parameters:
##            name matrix row col  Estimate  Std.Error A
## 1 ACE.mean[1,1]   mean   1   1 0.1193452 0.06804066  
## 2    ACE.a[1,1]      a   X   X 0.7235690 0.09621293  
## 3    ACE.c[1,1]      c   X   X 0.6076328 0.11776999  
## 4    ACE.e[1,1]      e   X   X 0.4542069 0.03247749  
## 
## Model Statistics: 
##                |  Parameters  |  Degrees of Freedom
##        Model:              4                    396
##    Saturated:             NA                     NA
## Independence:             NA                     NA
##                |  Fit (-2lnL units)
##        Model:              1025.121
##    Saturated:                    NA
## Independence:                    NA
## Number of observations/statistics: 200/400
## 
## Information Criteria: 
##       |  df Penalty  |  Parameters Penalty
## AIC:       233.1207               1033.121
## BIC:     -1073.0129               1046.314
##       |  Sample-Size Adjusted
## AIC:                 1033.326
## BIC:                 1033.642
## CFI: NA 
## TLI: 1   (also known as NNFI) 
## RMSEA:  0  [95% CI (NA, NA)]
## Prob(RMSEA <= 0.05): NA
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2021-01-12 04:22:33 
## Wall clock time: 0.07568693 secs 
## optimizer:  SLSQP 
## OpenMx version number: 2.18.1 
## Need help?  See help(mxSummary)

Output - parameters

# Output tables from the univariate ACE
get_output_tables(model)
## $Variance_components
##        A(%)      C(%)      E(%)
## X 0.4763576 0.3359353 0.1877071
## 
## $Raw_variance
##      Total        A         C         E
## X 1.099073 0.523552 0.3692176 0.2063039

Output - model fit

# Model fit
model <- ref_models(model)

fit_stats(model)
##   model      base ep minus2LL  df      AIC       BIC
## 1   ACE Saturated  4 1025.121 396 233.1207 -1073.013
##        CFI      TLI RMSEA   diffLL diffdf         p
## 1 1.004281 1.001427     0 5.391483      6 0.4946683

Constrained models

# Constrained models
# This is a list of MxModels
nested_models <- def_nested_models(
  model,
  AE = mxConstraint(c[1, 1] == 0),
  CE = mxConstraint(a[1, 1] == 0),
  E = list(mxConstraint(c[1, 1] == 0),
           mxConstraint(a[1, 1] == 0)),
  run = TRUE
)

fit_stats(model, nested_models = nested_models)
## Warning in computeFitStatistics(likelihood, DoF, chi,
## chiDoF, retval[["numObs"]], : Your model may be mis-
## specified (and fit worse than an independence model),
## or you may be using the wrong independence model, see ?
## mxRefModels
##   model      base ep minus2LL  df      AIC        BIC
## 1   ACE Saturated  4 1025.121 396 233.1207 -1073.0129
## 2    AE       ACE  4 1030.892 397 236.8918 -1072.5402
## 3    CE       ACE  4 1041.887 397 247.8875 -1061.5445
## 4     E       ACE  4 1168.096 398 372.0962  -940.6342
##          CFI       TLI      RMSEA     diffLL diffdf
## 1 1.00428122 1.0014271 0.00000000   5.391483      6
## 2 0.97071468 0.9916328 0.05452727   5.771030      1
## 3 0.89335470 0.9695299 0.10405420  16.766705      1
## 4 0.01244993 0.7531125 0.29619134 142.975415      2
##              p
## 1 4.946683e-01
## 2 1.629248e-02
## 3 4.226849e-05
## 4 8.980141e-32

Bivariate toy data

# Making bivariate toy data
mz_data2 <- MASS::mvrnorm(
  n = 100,
  mu = c(0, 0, 0, 0),
  Sigma = matrix(
    c(1.0, 0.7, 0.6, 0.5,
      0.7, 1.0, 0.5, 0.8,
      0.6, 0.5, 1.0, 0.7,
      0.5, 0.8, 0.7, 1.0),
    nrow = 4)) %>%
  as.data.frame() %>%
  setNames(c('X1', 'Y1', 'X2', 'Y2'))

dz_data2 <- MASS::mvrnorm(
  n = 100,
  mu = c(0, 0, 0, 0),
  Sigma = matrix(
    c(1.0, 0.7, 0.4, 0.2,
      0.7, 1.0, 0.2, 0.4,
      0.4, 0.2, 1.0, 0.7,
      0.2, 0.4, 0.7, 1.0),
    nrow = 4)) %>%
  as.data.frame() %>%
  setNames(c('X1', 'Y1', 'X2', 'Y2'))

data2 <- rbind(data.frame(zyg = 'MZ', mz_data2),
               data.frame(zyg = 'DZ', dz_data2)) %>%
  mutate(
    zyg = factor(zyg, levels = c('MZ', 'DZ'))
  )

head(data2)
##   zyg         X1         Y1         X2          Y2
## 1  MZ -2.6916574 -2.8505475 -1.4463712 -1.76899631
## 2  MZ  0.7523445  0.8565155  1.0967468  0.62762898
## 3  MZ  0.7106782  0.6803035 -0.1330686  0.08203872
## 4  MZ -0.8245882 -0.7115066 -0.5609679 -0.01247062
## 5  MZ -1.3944491 -0.8541269 -1.4051872 -0.65260362
## 6  MZ  1.6806667  1.1647597 -0.2466635 -0.17531906

Descriptives

# Descriptive statistics
get_univ_descriptives(data2, 
                      zyg = 'zyg', 
                      vars = c('X', 'Y'))
## $X
##    Twin1------------------ Twin2-------------------- Cor---------------------------
##    n   M          SD       n   M           SD        r         lowerCI    upperCI  
## MZ 100 0.07887346 1.081373 100 -0.04687941 0.9530093 0.6648379 0.53877652 0.7617783
## DZ 100 0.20394257 1.084314 100  0.04904954 0.9167004 0.2622338 0.06938975 0.4361851
## 
## $Y
##    Twin1----------------- Twin2------------------- Cor--------------------------
##    n   M         SD       n   M          SD        r         lowerCI   upperCI  
## MZ 100 0.1097118 1.073960 100 0.06303546 0.9887703 0.7824894 0.6925061 0.8484976
## DZ 100 0.1461730 1.026642 100 0.01767537 0.9693137 0.3916528 0.2115053 0.5460635

Twin correlations

# Cross-twin cross-trait correlations
get_cross_trait_cors(data2,
                     zyg = 'zyg',
                     vars = c('X', 'Y'))
## zyg: MZ
##           X1        Y1        X2        Y2
## X1 1.0000000 0.7758757 0.6648379 0.6373276
## Y1 0.7758757 1.0000000 0.5397699 0.7824894
## X2 0.6648379 0.5397699 1.0000000 0.7067741
## Y2 0.6373276 0.7824894 0.7067741 1.0000000
## -------------------------------------------- 
## zyg: DZ
##           X1         Y1         X2        Y2
## X1 1.0000000 0.69615564 0.26223377 0.1499928
## Y1 0.6961556 1.00000000 0.07471547 0.3916528
## X2 0.2622338 0.07471547 1.00000000 0.6412521
## Y2 0.1499928 0.39165278 0.64125211 1.0000000

Bivariate ACE

# Bivariate ACE model
model2 <- multivariate_ace(
  data = data2,
  zyg = 'zyg',
  ph = c('X', 'Y')
)

model2 <- mxRun(model2)

summary(model2)
## Summary of ACE 
##  
## free parameters:
##             name matrix row col      Estimate  Std.Error A
## 1  ACE.mean[1,1]   mean   1   1  7.463798e-02 0.06165076  
## 2  ACE.mean[1,2]   mean   1   2  8.757698e-02 0.06409393  
## 3     ACE.a[1,1]      a   X   X  7.808309e-07 0.23678477  
## 4     ACE.a[2,1]      a   Y   X  7.902705e-01 0.06220530  
## 5     ACE.a[2,2]      a   Y   Y  7.988041e-01 0.06628711  
## 6     ACE.c[1,1]      c   X   X -1.635437e-06 0.38064352  
## 7     ACE.c[2,1]      c   Y   X -2.007247e-01 0.10728556  
## 8     ACE.c[2,2]      c   Y   Y  3.921293e-01 0.10784317  
## 9     ACE.e[1,1]      e   X   X  4.801766e-01 0.02513380  
## 10    ACE.e[2,1]      e   Y   X  3.657950e-01 0.05190825  
## 11    ACE.e[2,2]      e   Y   Y  4.870613e-01 0.03402539  
## 
## Model Statistics: 
##                |  Parameters  |  Degrees of Freedom
##        Model:             11                    789
##    Saturated:             NA                     NA
## Independence:             NA                     NA
##                |  Fit (-2lnL units)
##        Model:              1826.955
##    Saturated:                    NA
## Independence:                    NA
## Number of observations/statistics: 200/800
## 
## Information Criteria: 
##       |  df Penalty  |  Parameters Penalty
## AIC:       248.9549               1848.955
## BIC:     -2353.4175               1885.236
##       |  Sample-Size Adjusted
## AIC:                 1850.359
## BIC:                 1850.387
## CFI: NA 
## TLI: 1   (also known as NNFI) 
## RMSEA:  0  [95% CI (NA, NA)]
## Prob(RMSEA <= 0.05): NA
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2021-01-12 04:22:34 
## Wall clock time: 0.08468008 secs 
## optimizer:  SLSQP 
## OpenMx version number: 2.18.1 
## Need help?  See help(mxSummary)

Output - parameters

# Output tables from the bivariate ACE
get_output_tables(model2)
## $Variance_components
##        A(%)       C(%)      E(%)
## X 0.6068125 0.03914755 0.3540399
## Y 0.6200555 0.14941993 0.2305245
## 
## $Covariation_components
##             Total         A           C         E
## X <-> X 1.0291933 0.6245274  0.04029040 0.3643755
## X <-> Y 0.7307259 0.6312713 -0.07871003 0.1781646
## Y <-> Y 1.0290821 0.6380881  0.15376538 0.2372287
##              A(%)       C(%)      E(%)
## X <-> X 0.6068125 0.03914755 0.3540399
## X <-> Y 0.7107743 0.08862286 0.2006028
## Y <-> Y 0.6200555 0.14941993 0.2305245
## 
## $All_correlations
##         r_A r_C       r_E    Total
## X <-> X   1   1 1.0000000 1.000000
## X <-> Y   1  -1 0.6059868 0.710037
## Y <-> Y   1   1 1.0000000 1.000000

Output - model fit

# Model fit
model2 <- ref_models(model2)

fit_stats(model2)
##   model      base ep minus2LL  df      AIC       BIC
## 1   ACE Saturated 11 1826.955 789 248.9549 -2353.418
##         CFI       TLI      RMSEA   diffLL diffdf         p
## 1 0.9933821 0.9953285 0.02989659 20.03894     17 0.2722424

Vignette

OpenMx

A model (mxModel):

  • matrices (mxMatrix)
  • matrix algebra (mxAlgebra)
  • data (mxData)
  • fit function (mxFitFunction…)
  • expectation (mxExpectation…)
  • sub-models (mxModel)
  • CIs (mxCI), constraints (mxConstraint)

Univatiate ACE

# Twin models in OpenMx
library(OpenMx)

# Univariate ACE model
model <- mxModel(
  name = 'ACE',
  
  # Means
  mxMatrix(type = 'Full',
           nrow = 1, ncol = 1,
           free = TRUE,
           name = 'mean'),
  mxAlgebra(cbind(mean, mean),
            name = 'expMeans'),
  
  # Variance components
  mxMatrix(type = 'Lower',
           nrow = 1, ncol = 1,
           free = TRUE,
           name = 'a'),
  mxMatrix(type = 'Lower',
           nrow = 1, ncol = 1,
           free = TRUE,
           name = 'c'),
  mxMatrix(type = 'Lower',
           nrow = 1, ncol = 1,
           free = TRUE,
           name = 'e'),
  mxAlgebra(t(a) %*% a, name = 'A'),
  mxAlgebra(t(c) %*% c, name = 'C'),
  mxAlgebra(t(e) %*% e, name = 'E'),
  
  # Expected covariance
  mxAlgebra(rbind(cbind(A + C + E, A + C    ),
                  cbind(A + C,     A + C + E)),
            name = 'expCovMZ'),
  mxAlgebra(rbind(cbind(A + C + E, .5%x%A + C),
                  cbind(.5%x%A + C, A + C + E)),
            name = 'expCovDZ'),
  
  # MZ submodel
  mxModel(name = 'MZ',
          mxData(observed = mz_data,
                 type = 'raw'),
          mxExpectationNormal(covariance = 'ACE.expCovMZ',
                              means = 'ACE.expMeans',
                              dimnames = c('X1', 'X2')),
          mxFitFunctionML()),
  
  # DZ submodel
  mxModel(name = 'DZ',
          mxData(observed = dz_data,
                 type = 'raw'),
          mxExpectationNormal(covariance = 'ACE.expCovDZ',
                              means = 'ACE.expMeans',
                              dimnames = c('X1', 'X2')),
          mxFitFunctionML()),
  
  mxFitFunctionMultigroup(c('MZ','DZ')),
  
  # Output tables
  mxAlgebra(A + C + E, name = 'V'),
  mxAlgebra(cbind(V, A, C, E),
            dimnames = list('X', c('Total', 'A', 'C', 'E')),
            name = 'Raw_variance'),
  mxAlgebra(cbind(A / V, C / V, E / V),
            dimnames = list('X', c('A(%)', 'C(%)', 'E(%)')),
            name = 'Variance_components')
)

model <- mxRun(model)

summary(model)
## Summary of ACE 
##  
## free parameters:
##            name matrix row col  Estimate  Std.Error A
## 1 ACE.mean[1,1]   mean   1   1 0.1193452 0.06804098  
## 2    ACE.a[1,1]      a   1   1 0.7235689 0.09621386  
## 3    ACE.c[1,1]      c   1   1 0.6076328 0.11777050  
## 4    ACE.e[1,1]      e   1   1 0.4542069 0.03247763  
## 
## Model Statistics: 
##                |  Parameters  |  Degrees of Freedom
##        Model:              4                    396
##    Saturated:             NA                     NA
## Independence:             NA                     NA
##                |  Fit (-2lnL units)
##        Model:              1025.121
##    Saturated:                    NA
## Independence:                    NA
## Number of observations/statistics: 200/400
## 
## Information Criteria: 
##       |  df Penalty  |  Parameters Penalty
## AIC:       233.1207               1033.121
## BIC:     -1073.0129               1046.314
##       |  Sample-Size Adjusted
## AIC:                 1033.326
## BIC:                 1033.642
## CFI: NA 
## TLI: 1   (also known as NNFI) 
## RMSEA:  0  [95% CI (NA, NA)]
## Prob(RMSEA <= 0.05): NA
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2021-01-12 04:22:34 
## Wall clock time: 0.04203105 secs 
## optimizer:  SLSQP 
## OpenMx version number: 2.18.1 
## Need help?  See help(mxSummary)

Univariate output

# Univariate output
mxEval(Raw_variance, model)
##      Total        A         C         E
## X 1.099073 0.523552 0.3692176 0.2063039
mxEval(Variance_components, model)
##        A(%)      C(%)      E(%)
## X 0.4763576 0.3359353 0.1877071

Bivariate ACE

# Bivariate ACE model
model2 <- mxModel(
  name = 'ACE',
  
  # Means
  mxMatrix(type = 'Full',
           nrow = 1, ncol = 2,
           free = TRUE,
           name = 'mean'),
  mxAlgebra(cbind(mean, mean),
            name = 'expMeans'),
  
  # Variance components
  mxMatrix(type = 'Lower',
           nrow = 2, ncol = 2,
           free = TRUE,
           name = 'a'),
  mxMatrix(type = 'Lower',
           nrow = 2, ncol = 2,
           free = TRUE,
           name = 'c'),
  mxMatrix(type = 'Lower',
           nrow = 2, ncol = 2,
           free = TRUE,
           name = 'e'),
  mxAlgebra(t(a) %*% a, name = 'A'),
  mxAlgebra(t(c) %*% c, name = 'C'),
  mxAlgebra(t(e) %*% e, name = 'E'),
  
  # Covariance
  mxAlgebra(rbind(cbind(A + C + E, A + C    ),
                  cbind(A + C,     A + C + E)),
            name = 'expCovMZ'),
  mxAlgebra(rbind(cbind(A + C + E, .5%x%A + C),
                  cbind(.5%x%A + C, A + C + E)),
            name = 'expCovDZ'),
  
  # MZ submodel
  mxModel(name = 'MZ',
          mxData(observed = mz_data2,
                 type = 'raw'),
          mxExpectationNormal(covariance = 'ACE.expCovMZ',
                              means = 'ACE.expMeans',
                              dimnames = c('X1', 'Y1', 'X2', 'Y2')),
          mxFitFunctionML()),
  
  # DZ submodel
  mxModel(name = 'DZ',
          mxData(observed = dz_data2,
                 type = 'raw'),
          mxExpectationNormal(covariance = 'ACE.expCovDZ',
                              means = 'ACE.expMeans',
                              dimnames = c('X1', 'Y1', 'X2', 'Y2')),
          mxFitFunctionML()),
  
  mxFitFunctionMultigroup(c('MZ','DZ')),
  
  # Output tables
  mxAlgebra(A + C + E, name = 'V'),
  mxAlgebra(abs(A) + abs(C) + abs(E), name = 'Vabs'),
  mxAlgebra(abs(A) / Vabs, name = 'Aperc'),
  mxAlgebra(abs(C) / Vabs, name = 'Cperc'),
  mxAlgebra(abs(E) / Vabs, name = 'Eperc'),
  mxAlgebra(cbind(diag2vec(Aperc),
                  diag2vec(Cperc),
                  diag2vec(Eperc)),
            dimnames = list(c('X', 'Y'), c('A(%)', 'C(%)', 'E(%)')),
            name = 'Variance_components'),
  
  mxAlgebra(cbind(vech(V),
                  vech(A), vech(C), vech(E),
                  vech(Aperc), vech(Cperc), vech(Eperc)),
            dimnames = list(c('X <-> X', 'X <-> Y', 'Y <-> Y'),
                            c('Total', 'A', 'C', 'E', 'A(%)', 'C(%)', 'E(%)')),
            name = 'Covariation_components'),
  
  mxMatrix(type = "Iden",
           nrow = 2, ncol = 2,
           name = "I"),
  mxAlgebra(sqrt(I * V), name = "SD_total"),
  mxAlgebra(sqrt(I * A), name = 'SD_A'),
  mxAlgebra(sqrt(I * C), name = 'SD_C'),
  mxAlgebra(sqrt(I * E), name = 'SD_E'),
  
  mxAlgebra(solve(SD_total) %&% V, name='r_total'),
  mxAlgebra(solve(SD_A) %&% A, name = 'r_A'),
  mxAlgebra(solve(SD_C) %&% C, name = 'r_C'),
  mxAlgebra(solve(SD_E) %&% E, name = 'r_E'),
  
  mxAlgebra(cbind(vech(r_A), vech(r_C), vech(r_E),
                  vech(r_total)),
            dimnames = list(c('X <-> X', 'X <-> Y', 'Y <-> Y'),
                            c('r_A', 'r_C', 'r_E', 'Total')),
            name='All_correlations')
)

model2 <- mxRun(model2)

summary(model2)
## Summary of ACE 
##  
## free parameters:
##             name matrix row col      Estimate  Std.Error A
## 1  ACE.mean[1,1]   mean   1   1  7.463802e-02 0.06165076  
## 2  ACE.mean[1,2]   mean   1   2  8.757693e-02 0.06409414  
## 3     ACE.a[1,1]      a   1   1  5.971474e-07 0.23679816  
## 4     ACE.a[2,1]      a   2   1  7.902704e-01 0.06220489  
## 5     ACE.a[2,2]      a   2   2  7.988042e-01 0.06628864  
## 6     ACE.c[1,1]      c   1   1  7.628998e-07 0.38072582  
## 7     ACE.c[2,1]      c   2   1 -2.007245e-01 0.10728655  
## 8     ACE.c[2,2]      c   2   2  3.921294e-01 0.10784464  
## 9     ACE.e[1,1]      e   1   1  4.801766e-01 0.02513386  
## 10    ACE.e[2,1]      e   2   1  3.657949e-01 0.05190817  
## 11    ACE.e[2,2]      e   2   2  4.870612e-01 0.03402557  
## 
## Model Statistics: 
##                |  Parameters  |  Degrees of Freedom
##        Model:             11                    789
##    Saturated:             NA                     NA
## Independence:             NA                     NA
##                |  Fit (-2lnL units)
##        Model:              1826.955
##    Saturated:                    NA
## Independence:                    NA
## Number of observations/statistics: 200/800
## 
## Information Criteria: 
##       |  df Penalty  |  Parameters Penalty
## AIC:       248.9549               1848.955
## BIC:     -2353.4175               1885.236
##       |  Sample-Size Adjusted
## AIC:                 1850.359
## BIC:                 1850.387
## CFI: NA 
## TLI: 1   (also known as NNFI) 
## RMSEA:  0  [95% CI (NA, NA)]
## Prob(RMSEA <= 0.05): NA
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2021-01-12 04:22:34 
## Wall clock time: 0.09264636 secs 
## optimizer:  SLSQP 
## OpenMx version number: 2.18.1 
## Need help?  See help(mxSummary)

Bivariate output

# Bivariate output
mxEval(Variance_components, model2)
##        A(%)      C(%)      E(%)
## X 0.6068126 0.0391475 0.3540399
## Y 0.6200556 0.1494200 0.2305244
mxEval(All_correlations, model2)
##         r_A r_C       r_E    Total
## X <-> X   1   1 1.0000000 1.000000
## X <-> Y   1  -1 0.6059866 0.710037
## Y <-> Y   1   1 1.0000000 1.000000

Bivariate output

mxEval(Covariation_components, model2)
##             Total         A           C         E
## X <-> X 1.0291932 0.6245273  0.04029034 0.3643755
## X <-> Y 0.7307258 0.6312714 -0.07871000 0.1781645
## Y <-> Y 1.0290823 0.6380882  0.15376550 0.2372286
##              A(%)       C(%)      E(%)
## X <-> X 0.6068126 0.03914750 0.3540399
## X <-> Y 0.7107744 0.08862283 0.2006027
## Y <-> Y 0.6200556 0.14942002 0.2305244

mlth.data.frame

library(mlth.data.frame)

# mlth.data.frame - multiheaded data.frame
model2 <- multivariate_ace(
  data = data2,
  zyg = 'zyg',
  ph = c('X', 'Y')
)
model2 <- def_ci(model2, ci = 'Variance_components')
model2 <- mxRun(model2, intervals = TRUE)

get_output_tables(model2)
## $Variance_components
##   A(%)------------------------- C(%)----------------------------- E(%)-------------------------
##   est       lCI       uCI       est        lCI          uCI       est       lCI       uCI      
## X 0.6068125 0.4472799 0.7191698 0.03914755 6.708449e-22 0.1590670 0.3540399 0.2648382 0.4683790
## Y 0.6200555 0.4191779 0.7638497 0.14941993 2.689729e-02 0.3204491 0.2305245 0.1692345 0.3152961
## 
## $Covariation_components
##             Total         A           C         E
## X <-> X 1.0291933 0.6245274  0.04029040 0.3643755
## X <-> Y 0.7307259 0.6312713 -0.07871003 0.1781646
## Y <-> Y 1.0290821 0.6380881  0.15376538 0.2372287
##              A(%)       C(%)      E(%)
## X <-> X 0.6068125 0.03914755 0.3540399
## X <-> Y 0.7107743 0.08862286 0.2006028
## Y <-> Y 0.6200555 0.14941993 0.2305245
## 
## $All_correlations
##         r_A r_C       r_E    Total
## X <-> X   1   1 1.0000000 1.000000
## X <-> Y   1  -1 0.6059868 0.710037
## Y <-> Y   1   1 1.0000000 1.000000

M <- get_output_tables(model2)$Variance_components

M[['A(%)']]
##   est       lCI       uCI      
## X 0.6068125 0.4472799 0.7191698
## Y 0.6200555 0.4191779 0.7638497

# Formatting the table for html or latex report
# Powered by knitr and kableExtra
library(kableExtra)

M %>% 
  behead() %>%
  kable(digits = 3,
        align = 'c') %>%
  add_complex_header_above(M)
A(%)
C(%)
E(%)
est lCI uCI est lCI uCI est lCI uCI
X 0.607 0.447 0.719 0.039 0.000 0.159 0.354 0.265 0.468
Y 0.620 0.419 0.764 0.149 0.027 0.320 0.231 0.169 0.315 

# Saving the table for the output
M %>%
  register_output(
    name = 'Table 1',
    caption = 'Table 1: Variance components')
##   A(%)------------------------- C(%)----------------------------- E(%)-------------------------
##   est       lCI       uCI       est        lCI          uCI       est       lCI       uCI      
## X 0.6068125 0.4472799 0.7191698 0.03914755 6.708449e-22 0.1590670 0.3540399 0.2648382 0.4683790
## Y 0.6200555 0.4191779 0.7638497 0.14941993 2.689729e-02 0.3204491 0.2305245 0.1692345 0.3152961

# Run at the end of the script
# to write all output tables into a single file
# openxlsx package required
write.xlsx.output(file = 'output.xlsx')

Vignettes