From ACE models to GCTA - continuous
logo





Environment

library(OpenMx)

# Create table to store all ACE estimates
res = NULL

ACE model - twins only

# ----------------------------------------------------------------------------------------------------------------------
# Script: 00_ACEvc.R Author: Sarah Medland Date: 29 02 2020
# Twin Univariate ACE model to estimate causes of variation
# across multiple groups Matrix style model - Raw data -
# Continuous data
# -------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------|

# Set Working Directory


# Clear Workspace rm(list=ls())

# Load Libraries & Options
library(OpenMx)

# ----------------------------------------------------------------------------------------------------------------------
# PREPARE DATA

# Load Data
df <- read.table("twinsibData.txt", header = T)

# Select Variables for Analysis
nv <- 1  # number of traits
nt <- 2  # number of individuals
ntv <- nv * nt  # number of total variables
selVars <- c("Twin1", "Twin2")  # name of traits
covVars <- c("age1", "age2", "sex1", "sex2")  # list of covariates

# Select Data for Analysis
mzData <- subset(df, zygosity == 1, c(selVars, covVars))
dzData <- subset(df, zygosity == 2, c(selVars, covVars))

# Check Descriptives
colMeans(mzData[, c(selVars, covVars)], na.rm = T)
##      Twin1      Twin2       age1       age2       sex1       sex2 
##  0.2063185  0.1923213 19.9514978 19.9514978  0.4950000  0.4950000
colMeans(dzData[, c(selVars, covVars)], na.rm = T)
##      Twin1      Twin2       age1       age2       sex1       sex2 
##  0.2239943  0.3360167 20.0053464 20.0053464  0.5000000  0.5030000
cov(mzData, use = "pairwise.complete.obs")
##           Twin1     Twin2         age1         age2         sex1         sex2
## Twin1 8.8515304 6.5239037 0.1261025811 0.1261025811 0.1094463581 0.1094463581
## Twin2 6.5239037 8.6506656 0.1163547599 0.1163547599 0.1070644625 0.1070644625
## age1  0.1261026 0.1163548 0.4845264527 0.4845264527 0.0005952069 0.0005952069
## age2  0.1261026 0.1163548 0.4845264527 0.4845264527 0.0005952069 0.0005952069
## sex1  0.1094464 0.1070645 0.0005952069 0.0005952069 0.2502252252 0.2502252252
## sex2  0.1094464 0.1070645 0.0005952069 0.0005952069 0.2502252252 0.2502252252
cov(dzData, use = "pairwise.complete.obs")
##             Twin1      Twin2         age1         age2          sex1
## Twin1  8.07787519 4.15272695  0.330321792  0.330321792  0.0436080682
## Twin2  4.15272695 8.79148695  0.333837419  0.333837419  0.0157968951
## age1   0.33032179 0.33383742  0.513803535  0.513803535 -0.0148824163
## age2   0.33032179 0.33383742  0.513803535  0.513803535 -0.0148824163
## sex1   0.04360807 0.01579690 -0.014882416 -0.014882416  0.2502502503
## sex2  -0.04706772 0.04208049  0.005280718  0.005280718 -0.0005005005
##                sex2
## Twin1 -0.0470677249
## Twin2  0.0420804866
## age1   0.0052807180
## age2   0.0052807180
## sex1  -0.0005005005
## sex2   0.2502412412
# Set Starting Values
sMu <- 0  # start value for means
sVa <- 0.3  # start value for A
sVc <- 0.3  # start value for C
sVe <- 0.3  # start value for E

# ----------------------------------------------------------------------------------------------------------------------
# PREPARE MODEL

# Create Algebra for expected Mean Matrices
intercept <- mxMatrix(type = "Full", nrow = 1, ncol = ntv, free = TRUE,
    values = sMu, labels = "interC", name = "intercept")
betaS <- mxMatrix(type = "Full", nrow = 1, ncol = nv, free = TRUE,
    values = 0, labels = "betaS", name = "bS")
betaA <- mxMatrix(type = "Full", nrow = 1, ncol = nv, free = TRUE,
    values = 0, labels = "betaA", name = "bA")
defSex <- mxMatrix(type = "Full", nrow = 1, ncol = nt, free = FALSE,
    labels = c("data.sex1", "data.sex2"), name = "Sex")
defAge <- mxMatrix(type = "Full", nrow = 1, ncol = nt, free = FALSE,
    labels = c("data.age1", "data.age2"), name = "Age")
expMean <- mxAlgebra(expression = intercept + Sex %x% bS + Age %x%
    bA, name = "expMean")

# Create Matrices for Variance Components
covA <- mxMatrix(type = "Symm", nrow = nv, ncol = nv, free = TRUE,
    values = sVa, label = "VA11", name = "VA")
covC <- mxMatrix(type = "Symm", nrow = nv, ncol = nv, free = TRUE,
    values = sVc, label = "VC11", name = "VC")
covE <- mxMatrix(type = "Symm", nrow = nv, ncol = nv, free = TRUE,
    values = sVe, label = "VE11", name = "VE")

# Create Algebra for expected Variance/Covariance Matrices
# in MZ & DZ twins
covP <- mxAlgebra(expression = VA + VC + VE, name = "V")
covMZ <- mxAlgebra(expression = VA + VC, name = "cMZ")
covDZ <- mxAlgebra(expression = 0.5 %x% VA + VC, name = "cDZ")
expCovMZ <- mxAlgebra(expression = rbind(cbind(V, cMZ), cbind(t(cMZ),
    V)), name = "expCovMZ")
expCovDZ <- mxAlgebra(expression = rbind(cbind(V, cDZ), cbind(t(cDZ),
    V)), name = "expCovDZ")

# Create Data Objects for Multiple Groups
dataMZ <- mxData(observed = mzData, type = "raw")
dataDZ <- mxData(observed = dzData, type = "raw")

# Create Expectation Objects for Multiple Groups
expMZ <- mxExpectationNormal(covariance = "expCovMZ", means = "expMean",
    dimnames = selVars)
expDZ <- mxExpectationNormal(covariance = "expCovDZ", means = "expMean",
    dimnames = selVars)
funML <- mxFitFunctionML()

# Create Model Objects for Multiple Groups
defs <- list(defAge, defSex)
pars <- list(intercept, betaS, betaA, covA, covC, covE, covP)
modelMZ <- mxModel(defs, pars, expMean, covMZ, expCovMZ, dataMZ,
    expMZ, funML, name = "MZ")
modelDZ <- mxModel(defs, pars, expMean, covDZ, expCovDZ, dataDZ,
    expDZ, funML, name = "DZ")
multi <- mxFitFunctionMultigroup(c("MZ", "DZ"))

# Create Algebra for Variance Components
rowVC <- rep("VC", nv)
colVC <- rep(c("VA", "VC", "VE", "SA", "SC", "SE"), each = nv)
estVC <- mxAlgebra(expression = cbind(VA, VC, VE, VA/V, VC/V,
    VE/V), name = "VarC", dimnames = list(rowVC, colVC))

# Create Confidence Interval Objects
ciACE <- mxCI("VarC[1,4:6]")

# Build Model with Confidence Intervals
modelACE <- mxModel("ACEvc", modelMZ, modelDZ, multi, pars, estVC,
    ciACE)

# ----------------------------------------------------------------------------------------------------------------------
# RUN MODEL

# Run ACE Model
fitACE <- mxRun(modelACE, intervals = T)
sumACE <- summary(fitACE)
sumACE
## Summary of ACEvc 
##  
## free parameters:
##     name    matrix row col   Estimate  Std.Error A
## 1 interC intercept   1   1 -9.3256687 1.63631862  
## 2  betaS        bS   1   1  0.2753476 0.09439997  
## 3  betaA        bA   1   1  0.4720194 0.08180316  
## 4   VA11        VA   1   1  4.2031247 0.42144351  
## 5   VC11        VC   1   1  2.0183948 0.40546428  
## 6   VE11        VE   1   1  2.2134076 0.09780964  
## 
## confidence intervals:
##                    lbound  estimate    ubound note
## ACEvc.VarC[1,4] 0.4031709 0.4983001 0.5990074     
## ACEvc.VarC[1,5] 0.1443103 0.2392901 0.3269656     
## ACEvc.VarC[1,6] 0.2384267 0.2624098 0.2889128     
## 
## Model Statistics: 
##                |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
##        Model:              6                   3994              18823.12
##    Saturated:             NA                     NA                    NA
## Independence:             NA                     NA                    NA
## Number of observations/statistics: 2000/4000
## 
## Information Criteria: 
##       |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
## AIC:       10835.12               18835.12                 18835.16
## BIC:      -11534.88               18868.73                 18849.66
## 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: 2024-03-05 09:28:38 
## Wall clock time: 6.07047 secs 
## optimizer:  SLSQP 
## OpenMx version number: 2.21.11 
## Need help?  See help(mxSummary)
# Some different ways to extract output
fitACE$output$estimate
##     interC      betaS      betaA       VA11       VC11       VE11 
## -9.3256687  0.2753476  0.4720194  4.2031247  2.0183948  2.2134076
fitACE$VarC$result
##          VA       VC       VE        SA        SC        SE
## VC 4.203125 2.018395 2.213408 0.4983001 0.2392901 0.2624098
fitACE$output$confidenceIntervals
##                    lbound  estimate    ubound
## ACEvc.VarC[1,4] 0.4031709 0.4983001 0.5990074
## ACEvc.VarC[1,5] 0.1443103 0.2392901 0.3269656
## ACEvc.VarC[1,6] 0.2384267 0.2624098 0.2889128

store results

res = rbind(res, c("Twin", fitACE$output$confidenceIntervals[1,
    ], fitACE$output$confidenceIntervals[2, ], fitACE$output$confidenceIntervals[3,
    ]))

ACE from twins and siblings

# ----------------------------------------------------------------------------------------------------------------------
# Script: 01_extrasib.R Author: Sarah Medland Date: 29 02
# 2020 Twin Univariate ACE model to estimate causes of
# variation across multiple groups Matrix style model - Raw
# data - Continuous data
# -------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------|

# Set Working Directory

# Clear Workspace rm(list=ls())

# Load Libraries & Options
library(OpenMx)

# ----------------------------------------------------------------------------------------------------------------------
# PREPARE DATA

# Load Data
df <- read.table("twinsibData.txt", header = T)

# Select Variables for Analysis
nv <- 1  # number of traits
nt <- 3  # number of individuals
ntv <- nv * nt  # number of total variables
selVars <- c("Twin1", "Twin2", "Sib")  # name of traits
covVars <- c("age1", "age2", "age3", "sex1", "sex2", "sex3")  # list of covariates

# Select Data for Analysis
mzData <- subset(df, zygosity == 1, c(selVars, covVars))
dzData <- subset(df, zygosity == 2, c(selVars, covVars))
cov(mzData[, selVars], use = "pairwise.complete.obs")
##          Twin1    Twin2      Sib
## Twin1 8.851530 6.523904 5.745641
## Twin2 6.523904 8.650666 5.558020
## Sib   5.745641 5.558020 9.304954
cov(dzData[, selVars], use = "pairwise.complete.obs")
##          Twin1    Twin2      Sib
## Twin1 8.077875 4.152727 3.695605
## Twin2 4.152727 8.791487 4.301196
## Sib   3.695605 4.301196 8.325628
# Set Starting Values
sMu <- 0  # start value for means
sVa <- 0.3  # start value for A
sVc <- 0.3  # start value for C
sVe <- 0.3  # start value for E

# ----------------------------------------------------------------------------------------------------------------------
# PREPARE MODEL

# Create Algebra for expected Mean Matrices
intercept <- mxMatrix(type = "Full", nrow = 1, ncol = ntv, free = TRUE,
    values = sMu, labels = "interC", name = "intercept")
betaS <- mxMatrix(type = "Full", nrow = 1, ncol = nv, free = TRUE,
    values = 0, labels = "betaS", name = "bS")
betaA <- mxMatrix(type = "Full", nrow = 1, ncol = nv, free = TRUE,
    values = 0, labels = "betaA", name = "bA")
defSex <- mxMatrix(type = "Full", nrow = 1, ncol = nt, free = FALSE,
    labels = c("data.sex1", "data.sex2", "data.sex3"), name = "Sex")
defAge <- mxMatrix(type = "Full", nrow = 1, ncol = nt, free = FALSE,
    labels = c("data.age1", "data.age2", "data.age3"), name = "Age")
expMean <- mxAlgebra(expression = intercept + Sex %x% bS + Age %x%
    bA, name = "expMean")

# Create Matrices for Variance Components
covA <- mxMatrix(type = "Symm", nrow = nv, ncol = nv, free = TRUE,
    values = sVa, label = "VA11", name = "VA")
covC <- mxMatrix(type = "Symm", nrow = nv, ncol = nv, free = TRUE,
    values = sVc, label = "VC11", name = "VC")
covE <- mxMatrix(type = "Symm", nrow = nv, ncol = nv, free = TRUE,
    values = sVe, label = "VE11", name = "VE")

# Create Algebra for expected Variance/Covariance Matrices
# in MZ & DZ twins
covP <- mxAlgebra(expression = VA + VC + VE, name = "V")
covMZ <- mxAlgebra(expression = VA + VC, name = "cMZ")
covDZ <- mxAlgebra(expression = 0.5 %x% VA + VC, name = "cDZ")
expCovMZ <- mxAlgebra(expression = rbind(cbind(V, cMZ, cDZ),
    cbind(t(cMZ), V, cDZ), cbind(t(cDZ), t(cDZ), V)), name = "expCovMZ")

expCovDZ <- mxAlgebra(expression = rbind(cbind(V, cDZ, cDZ),
    cbind(t(cDZ), V, cDZ), cbind(t(cDZ), t(cDZ), V)), name = "expCovDZ")

# Create Data Objects for Multiple Groups
dataMZ <- mxData(observed = mzData, type = "raw")
dataDZ <- mxData(observed = dzData, type = "raw")

# Create Expectation Objects for Multiple Groups
expMZ <- mxExpectationNormal(covariance = "expCovMZ", means = "expMean",
    dimnames = selVars)
expDZ <- mxExpectationNormal(covariance = "expCovDZ", means = "expMean",
    dimnames = selVars)
funML <- mxFitFunctionML()

# Create Model Objects for Multiple Groups
defs <- list(defAge, defSex)
pars <- list(intercept, betaS, betaA, covA, covC, covE, covP)
modelMZ <- mxModel(defs, pars, expMean, covMZ, covDZ, expCovMZ,
    dataMZ, expMZ, funML, name = "MZ")
modelDZ <- mxModel(defs, pars, expMean, covDZ, expCovDZ, dataDZ,
    expDZ, funML, name = "DZ")
multi <- mxFitFunctionMultigroup(c("MZ", "DZ"))

# Create Algebra for Variance Components
rowVC <- rep("VC", nv)
colVC <- rep(c("VA", "VC", "VE", "SA", "SC", "SE"), each = nv)
estVC <- mxAlgebra(expression = cbind(VA, VC, VE, VA/V, VC/V,
    VE/V), name = "VarC", dimnames = list(rowVC, colVC))

# Create Confidence Interval Objects
ciACE <- mxCI("VarC[1,4:6]")

# Build Model with Confidence Intervals
modelACE <- mxModel("ACEsib", modelMZ, modelDZ, multi, pars,
    estVC, ciACE)


# ----------------------------------------------------------------------------------------------------------------------
# RUN MODEL

# Run ACE Model
fitACE <- mxRun(modelACE, intervals = T)
sumACE <- summary(fitACE)
sumACE
## Summary of ACEsib 
##  
## free parameters:
##     name    matrix row col   Estimate  Std.Error A
## 1 interC intercept   1   1 -5.5920306 0.83974410  
## 2  betaS        bS   1   1  0.3613256 0.06485551  
## 3  betaA        bA   1   1  0.2814054 0.04186418  
## 4   VA11        VA   1   1  3.2101782 0.27289307  
## 5   VC11        VC   1   1  3.0337892 0.25834607  
## 6   VE11        VE   1   1  2.3026430 0.10473313  
## 
## confidence intervals:
##                     lbound  estimate    ubound note
## ACEsib.VarC[1,4] 0.3121488 0.3756084 0.4378840     
## ACEsib.VarC[1,5] 0.3028703 0.3549699 0.4054639     
## ACEsib.VarC[1,6] 0.2443557 0.2694218 0.2972491     
## 
## Model Statistics: 
##                |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
##        Model:              6                   5994               27891.1
##    Saturated:             NA                     NA                    NA
## Independence:             NA                     NA                    NA
## Number of observations/statistics: 2000/6000
## 
## Information Criteria: 
##       |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
## AIC:       15903.10               27903.10                 27903.14
## BIC:      -17668.71               27936.71                 27917.64
## 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: 2024-03-05 09:28:44 
## Wall clock time: 6.56052 secs 
## optimizer:  SLSQP 
## OpenMx version number: 2.21.11 
## Need help?  See help(mxSummary)

store results

res = rbind(res, c("Twin_sib", fitACE$output$confidenceIntervals[1,
    ], fitACE$output$confidenceIntervals[2, ], fitACE$output$confidenceIntervals[3,
    ]))

Twin and siblings - alternate parametrisation

# ----------------------------------------------------------------------------------------------------------------------
# Script: extrasib2.R Author: Sarah Medland Date: 29 02
# 2020 Twin Univariate ACE model to estimate causes of
# variation across multiple groups Matrix style model - Raw
# data - Continuous data
# -------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------|

# Set Working Directory

# Clear Workspace rm(list=ls())

# Load Libraries & Options
library(OpenMx)

# ----------------------------------------------------------------------------------------------------------------------
# PREPARE DATA

# Load Data
df <- read.table("twinsibData.txt", header = T)

# Select Variables for Analysis
nv <- 1  # number of traits
nt <- 3  # number of individuals
ntv <- nv * nt  # number of total variables
selVars <- c("Twin1", "Twin2", "Sib")  # name of traits
covVars <- c("age1", "age2", "age3", "sex1", "sex2", "sex3")  # list of covariates

# Select Data for Analysis
mzData <- subset(df, zygosity == 1, c(selVars, covVars))
dzData <- subset(df, zygosity == 2, c(selVars, covVars))
cov(mzData[, selVars], use = "pairwise.complete.obs")
##          Twin1    Twin2      Sib
## Twin1 8.851530 6.523904 5.745641
## Twin2 6.523904 8.650666 5.558020
## Sib   5.745641 5.558020 9.304954
cov(dzData[, selVars], use = "pairwise.complete.obs")
##          Twin1    Twin2      Sib
## Twin1 8.077875 4.152727 3.695605
## Twin2 4.152727 8.791487 4.301196
## Sib   3.695605 4.301196 8.325628
# Set Starting Values
sMu <- 0  # start value for means
sVa <- 0.3  # start value for A
sVc <- 0.3  # start value for C
sVe <- 0.3  # start value for E

# ----------------------------------------------------------------------------------------------------------------------
# PREPARE MODEL

# Create Algebra for expected Mean Matrices
intercept <- mxMatrix(type = "Full", nrow = 1, ncol = ntv, free = TRUE,
    values = sMu, labels = "interC", name = "intercept")
betaS <- mxMatrix(type = "Full", nrow = 1, ncol = nv, free = TRUE,
    values = 0, labels = "betaS", name = "bS")
betaA <- mxMatrix(type = "Full", nrow = 1, ncol = nv, free = TRUE,
    values = 0, labels = "betaA", name = "bA")
defSex <- mxMatrix(type = "Full", nrow = 1, ncol = nt, free = FALSE,
    labels = c("data.sex1", "data.sex2", "data.sex3"), name = "Sex")
defAge <- mxMatrix(type = "Full", nrow = 1, ncol = nt, free = FALSE,
    labels = c("data.age1", "data.age2", "data.age3"), name = "Age")
expMean <- mxAlgebra(expression = intercept + Sex %x% bS + Age %x%
    bA, name = "expMean")

# Create Matrices for Variance Components
covA <- mxMatrix(type = "Symm", nrow = nv, ncol = nv, free = TRUE,
    values = sVa, label = "VA11", name = "VA")
covC <- mxMatrix(type = "Symm", nrow = nv, ncol = nv, free = TRUE,
    values = sVc, label = "VC11", name = "VC")
covE <- mxMatrix(type = "Symm", nrow = nv, ncol = nv, free = TRUE,
    values = sVe, label = "VE11", name = "VE")

# Create Algebra for expected Variance/Covariance Matrices
# in MZ & DZ twins
covP <- mxAlgebra(expression = VA + VC + VE, name = "V")
relMZ <- mxMatrix(type = "Symm", nrow = nt, ncol = nt, free = FALSE,
    values = c(1, 1, 0.5, 1, 0.5, 1), name = "rAmz")
relDZ <- mxMatrix(type = "Symm", nrow = nt, ncol = nt, free = FALSE,
    values = c(1, 0.5, 0.5, 1, 0.5, 1), name = "rAdz")
relC <- mxMatrix(type = "Unit", nrow = nt, ncol = nt, free = FALSE,
    name = "rC")
relE <- mxMatrix(type = "Iden", nrow = nt, ncol = nt, free = FALSE,
    name = "rE")

expCovMZ <- mxAlgebra(expression = rAmz %x% VA + rC %x% VC +
    rE %x% VE, name = "expCovMZ")
expCovDZ <- mxAlgebra(expression = rAdz %x% VA + rC %x% VC +
    rE %x% VE, name = "expCovDZ")

# Create Data Objects for Multiple Groups
dataMZ <- mxData(observed = mzData, type = "raw")
dataDZ <- mxData(observed = dzData, type = "raw")

# Create Expectation Objects for Multiple Groups
expMZ <- mxExpectationNormal(covariance = "expCovMZ", means = "expMean",
    dimnames = selVars)
expDZ <- mxExpectationNormal(covariance = "expCovDZ", means = "expMean",
    dimnames = selVars)
funML <- mxFitFunctionML()

# Create Model Objects for Multiple Groups
defs <- list(defAge, defSex)
pars <- list(intercept, betaS, betaA, covA, covC, covE, covP,
    relMZ, relDZ, relC, relE)
modelMZ <- mxModel(defs, pars, expMean, expCovMZ, dataMZ, expMZ,
    funML, name = "MZ")
modelDZ <- mxModel(defs, pars, expMean, expCovDZ, dataDZ, expDZ,
    funML, name = "DZ")
multi <- mxFitFunctionMultigroup(c("MZ", "DZ"))

# Create Algebra for Variance Components
rowVC <- rep("VC", nv)
colVC <- rep(c("VA", "VC", "VE", "SA", "SC", "SE"), each = nv)
estVC <- mxAlgebra(expression = cbind(VA, VC, VE, VA/V, VC/V,
    VE/V), name = "VarC", dimnames = list(rowVC, colVC))

# Create Confidence Interval Objects
ciACE <- mxCI("VarC[1,4:6]")

# Build Model with Confidence Intervals
modelACE <- mxModel("ACEsib_alt", modelMZ, modelDZ, multi, pars,
    estVC, ciACE)

# ----------------------------------------------------------------------------------------------------------------------
# RUN MODEL

# Run ACE Model
fitACE <- mxRun(modelACE, intervals = T)
sumACE <- summary(fitACE)
sumACE
## Summary of ACEsib_alt 
##  
## free parameters:
##     name    matrix row col   Estimate  Std.Error A
## 1 interC intercept   1   1 -5.5920306 0.83974410  
## 2  betaS        bS   1   1  0.3613256 0.06485551  
## 3  betaA        bA   1   1  0.2814054 0.04186418  
## 4   VA11        VA   1   1  3.2101782 0.27289307  
## 5   VC11        VC   1   1  3.0337892 0.25834607  
## 6   VE11        VE   1   1  2.3026430 0.10473313  
## 
## confidence intervals:
##                         lbound  estimate    ubound note
## ACEsib_alt.VarC[1,4] 0.3121488 0.3756084 0.4378840     
## ACEsib_alt.VarC[1,5] 0.3028703 0.3549699 0.4054639     
## ACEsib_alt.VarC[1,6] 0.2443557 0.2694218 0.2972491     
## 
## Model Statistics: 
##                |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
##        Model:              6                   5994               27891.1
##    Saturated:             NA                     NA                    NA
## Independence:             NA                     NA                    NA
## Number of observations/statistics: 2000/6000
## 
## Information Criteria: 
##       |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
## AIC:       15903.10               27903.10                 27903.14
## BIC:      -17668.71               27936.71                 27917.64
## 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: 2024-03-05 09:28:51 
## Wall clock time: 6.55051 secs 
## optimizer:  SLSQP 
## OpenMx version number: 2.21.11 
## Need help?  See help(mxSummary)

store results

res = rbind(res, c("Twin_sib2", fitACE$output$confidenceIntervals[1,
    ], fitACE$output$confidenceIntervals[2, ], fitACE$output$confidenceIntervals[3,
    ]))

Zygosity as definition variable

# ----------------------------------------------------------------------------------------------------------------------
# Script: zygdef.R Author: Sarah Medland Date: 29 02 2020
# Twin Univariate ACE model to estimate causes of variation
# across multiple groups Matrix style model - Raw data -
# Continuous data
# -------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------|

# Set Working Directory

# Clear Workspace rm(list=ls())

# Load Libraries & Options
library(OpenMx)

# ----------------------------------------------------------------------------------------------------------------------
# PREPARE DATA

# Load Data
df <- read.table("twinsibData.txt", header = T)

# Select Variables for Analysis
nv <- 1  # number of traits
nt <- 3  # number of individuals
ntv <- nv * nt  # number of total variables
selVars <- c("Twin1", "Twin2", "Sib")  # name of traits
covVars <- c("age1", "age2", "age3", "sex1", "sex2", "sex3")  # list of covariates

# Check Correlations of Data for Analysis
cov(df[df$zygT == 1, selVars], use = "pairwise.complete.obs")
##          Twin1    Twin2      Sib
## Twin1 8.851530 6.523904 5.745641
## Twin2 6.523904 8.650666 5.558020
## Sib   5.745641 5.558020 9.304954
cov(df[df$zygT == 0.5, selVars], use = "pairwise.complete.obs")
##          Twin1    Twin2      Sib
## Twin1 8.077875 4.152727 3.695605
## Twin2 4.152727 8.791487 4.301196
## Sib   3.695605 4.301196 8.325628
# Set Starting Values
sMu <- 0  # start value for means
sVa <- 0.3  # start value for A
sVc <- 0.3  # start value for C
sVe <- 0.3  # start value for E

# ----------------------------------------------------------------------------------------------------------------------
# PREPARE MODEL

# Create Algebra for expected Mean Matrices
intercept <- mxMatrix(type = "Full", nrow = 1, ncol = ntv, free = TRUE,
    values = sMu, labels = "interC", name = "intercept")
betaS <- mxMatrix(type = "Full", nrow = 1, ncol = nv, free = TRUE,
    values = 0, labels = "betaS", name = "bS")
betaA <- mxMatrix(type = "Full", nrow = 1, ncol = nv, free = TRUE,
    values = 0, labels = "betaA", name = "bA")
defSex <- mxMatrix(type = "Full", nrow = 1, ncol = nt, free = FALSE,
    labels = c("data.sex1", "data.sex2", "data.sex3"), name = "Sex")
defAge <- mxMatrix(type = "Full", nrow = 1, ncol = nt, free = FALSE,
    labels = c("data.age1", "data.age2", "data.age3"), name = "Age")
expMean <- mxAlgebra(expression = intercept + Sex %x% bS + Age %x%
    bA, name = "expMean")

# Create Matrices for Variance Components
covA <- mxMatrix(type = "Symm", nrow = nv, ncol = nv, free = TRUE,
    values = sVa, label = "VA11", name = "VA")
covC <- mxMatrix(type = "Symm", nrow = nv, ncol = nv, free = TRUE,
    values = sVc, label = "VC11", name = "VC")
covE <- mxMatrix(type = "Symm", nrow = nv, ncol = nv, free = TRUE,
    values = sVe, label = "VE11", name = "VE")

# Create Algebra for expected Variance/Covariance Matrices
# in MZ & DZ twins
covP <- mxAlgebra(expression = VA + VC + VE, name = "V")
relA <- mxMatrix(type = "Stand", nrow = ntv, ncol = ntv, free = FALSE,
    labels = c("data.zygT", "data.zygS", "data.zygS"), name = "rA")
relC <- mxMatrix(type = "Unit", nrow = ntv, ncol = ntv, free = FALSE,
    name = "rC")
relE <- mxMatrix(type = "Iden", nrow = ntv, ncol = ntv, free = FALSE,
    name = "rE")

expCov <- mxAlgebra(expression = rA %x% VA + rC %x% VC + rE %x%
    VE, name = "expCov")

# Create Data Object
dataTW <- mxData(observed = df, type = "raw")

# Create Expectation Objects
expTW <- mxExpectationNormal(covariance = "expCov", means = "expMean",
    dimnames = selVars)
funML <- mxFitFunctionML()

# Create Model Objects
defs <- list(defAge, defSex, relA)
pars <- list(intercept, betaS, betaA, covA, covC, covE, covP,
    relC, relE)

# Create Algebra for Variance Components
rowVC <- rep("VC", nv)
colVC <- rep(c("VA", "VC", "VE", "SA", "SC", "SE"), each = nv)
estVC <- mxAlgebra(expression = cbind(VA, VC, VE, VA/V, VC/V,
    VE/V), name = "VarC", dimnames = list(rowVC, colVC))

# Create Confidence Interval Objects
ciACE <- mxCI("VarC[1,4:6]")

# Build Model with Confidence Intervals
modelACE <- mxModel("zygdef", defs, pars, expMean, expCov, dataTW,
    expTW, funML, estVC, ciACE)

# ----------------------------------------------------------------------------------------------------------------------
# RUN MODEL

# Run ACE Model
fitACE <- mxRun(modelACE, intervals = T)
sumACE <- summary(fitACE)
sumACE
## Summary of zygdef 
##  
## free parameters:
##     name    matrix row col   Estimate  Std.Error A
## 1 interC intercept   1   1 -5.5920305 0.83979380  
## 2  betaS        bS   1   1  0.3613256 0.06493829  
## 3  betaA        bA   1   1  0.2814054 0.04185934  
## 4   VA11        VA   1   1  3.2101782 0.27126856  
## 5   VC11        VC   1   1  3.0337892 0.25695191  
## 6   VE11        VE   1   1  2.3026429 0.10436147  
## 
## confidence intervals:
##                     lbound  estimate    ubound note
## zygdef.VarC[1,4] 0.3121488 0.3756084 0.4378840     
## zygdef.VarC[1,5] 0.3028703 0.3549699 0.4054639     
## zygdef.VarC[1,6] 0.2443557 0.2694218 0.2972491     
## 
## Model Statistics: 
##                |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
##        Model:              6                   5994               27891.1
##    Saturated:              9                   5991                    NA
## Independence:              6                   5994                    NA
## Number of observations/statistics: 2000/6000
## 
## Information Criteria: 
##       |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
## AIC:       15903.10               27903.10                 27903.14
## BIC:      -17668.71               27936.71                 27917.64
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2024-03-05 09:29:00 
## Wall clock time: 9.101765 secs 
## optimizer:  SLSQP 
## OpenMx version number: 2.21.11 
## Need help?  See help(mxSummary)
res = rbind(res, c("ZygVar", fitACE$output$confidenceIntervals[1,
    ], fitACE$output$confidenceIntervals[2, ], fitACE$output$confidenceIntervals[3,
    ]))

ACE with relatedness

# ----------------------------------------------------------------------------------------------------------------------
# Script: 04_relatedness.R Author: Sarah Medland Date: 29
# 02 2020 Twin Univariate ACE model to estimate causes of
# variation across multiple groups Matrix style model - Raw
# data - Continuous data
# -------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------|

# Set Working Directory

# Clear Workspace rm(list=ls())

# Load Libraries & Options
library(OpenMx)

# ----------------------------------------------------------------------------------------------------------------------
# PREPARE DATA Load Data
df <- read.table("twinsibData.txt", header = T)

# Select Variables for Analysis
nv <- 1  # number of traits
nt <- 3  # number of individuals
ntv <- nv * nt  # number of total variables
selVars <- c("Twin1", "Twin2", "Sib")  # name of traits
covVars <- c("age1", "age2", "age3", "sex1", "sex2", "sex3")  # list of covariates

# Check distribution of genetic relatedness data
# s1=Twin1-Twin2; s2=Twin1-Sib; s3=Twin2-Sib
hist(df$s1)

hist(df$s2)

hist(df$s3)

# Check Correlations of Data for Analysis
cov(df[df$s1 > 0.9, selVars], use = "pairwise.complete.obs")
##          Twin1    Twin2      Sib
## Twin1 8.851530 6.523904 5.745641
## Twin2 6.523904 8.650666 5.558020
## Sib   5.745641 5.558020 9.304954
cov(df[df$s1 <= 0.9, selVars], use = "pairwise.complete.obs")
##          Twin1    Twin2      Sib
## Twin1 8.077875 4.152727 3.695605
## Twin2 4.152727 8.791487 4.301196
## Sib   3.695605 4.301196 8.325628
# Set Starting Values
sMu <- 0  # start value for means
sVa <- 0.3  # start value for A
sVc <- 0.3  # start value for C
sVe <- 0.3  # start value for E

# ----------------------------------------------------------------------------------------------------------------------
# PREPARE MODEL

# Create Algebra for expected Mean Matrices
intercept <- mxMatrix(type = "Full", nrow = 1, ncol = ntv, free = TRUE,
    values = sMu, labels = "interC", name = "intercept")
betaS <- mxMatrix(type = "Full", nrow = 1, ncol = nv, free = TRUE,
    values = 0, labels = "betaS", name = "bS")
betaA <- mxMatrix(type = "Full", nrow = 1, ncol = nv, free = TRUE,
    values = 0, labels = "betaA", name = "bA")
defSex <- mxMatrix(type = "Full", nrow = 1, ncol = nt, free = FALSE,
    labels = c("data.sex1", "data.sex2", "data.sex3"), name = "Sex")
defAge <- mxMatrix(type = "Full", nrow = 1, ncol = nt, free = FALSE,
    labels = c("data.age1", "data.age2", "data.age3"), name = "Age")
expMean <- mxAlgebra(expression = intercept + Sex %x% bS + Age %x%
    bA, name = "expMean")

# Create Matrices for Variance Components
covA <- mxMatrix(type = "Symm", nrow = nv, ncol = nv, free = TRUE,
    values = sVa, label = "VA11", name = "VA")
covC <- mxMatrix(type = "Symm", nrow = nv, ncol = nv, free = TRUE,
    values = sVc, label = "VC11", name = "VC")
covE <- mxMatrix(type = "Symm", nrow = nv, ncol = nv, free = TRUE,
    values = sVe, label = "VE11", name = "VE")

# Create Algebra for expected Variance/Covariance Matrices
# in MZ & DZ twins
covP <- mxAlgebra(expression = VA + VC + VE, name = "V")
relA <- mxMatrix(type = "Stand", nrow = nt, ncol = nt, free = FALSE,
    labels = c("data.s1", "data.s2", "data.s3"), name = "rA")
relC <- mxMatrix(type = "Unit", nrow = nt, ncol = nt, free = FALSE,
    name = "rC")
relE <- mxMatrix(type = "Iden", nrow = nt, ncol = nt, free = FALSE,
    name = "rE")

expCov <- mxAlgebra(expression = rA %x% VA + rC %x% VC + rE %x%
    VE, name = "expCov")

# Create Data Objects for Multiple Groups
dataTW <- mxData(observed = df, type = "raw")

# Create Expectation Objects for Multiple Groups
expTW <- mxExpectationNormal(covariance = "expCov", means = "expMean",
    dimnames = selVars)
funML <- mxFitFunctionML()

# Create Model Objects for Multiple Groups
defs <- list(defAge, defSex, relA)
pars <- list(intercept, betaS, betaA, covA, covC, covE, covP,
    relC, relE)

# Create Algebra for Variance Components
rowVC <- rep("VC", nv)
colVC <- rep(c("VA", "VC", "VE", "SA", "SC", "SE"), each = nv)
estVC <- mxAlgebra(expression = cbind(VA, VC, VE, VA/V, VC/V,
    VE/V), name = "VarC", dimnames = list(rowVC, colVC))

# Create Confidence Interval Objects
ciACE <- mxCI("VarC[1,4:6]")

# Build Model with Confidence Intervals
modelACE <- mxModel("grm", defs, pars, expMean, expCov, dataTW,
    expTW, funML, estVC, ciACE)

# ----------------------------------------------------------------------------------------------------------------------
# RUN MODEL

# Run ACE Model
fitACE <- mxRun(modelACE, intervals = T)
sumACE <- summary(fitACE)
sumACE
## Summary of grm 
##  
## free parameters:
##     name    matrix row col   Estimate  Std.Error A
## 1 interC intercept   1   1 -5.6308479 0.84586675  
## 2  betaS        bS   1   1  0.3613433 0.06481819  
## 3  betaA        bA   1   1  0.2833848 0.04216737  
## 4   VA11        VA   1   1  3.2116426 0.27318054  
## 5   VC11        VC   1   1  3.0372916 0.25838644  
## 6   VE11        VE   1   1  2.3022458 0.10462711  
## 
## confidence intervals:
##                  lbound  estimate    ubound note
## grm.VarC[1,4] 0.3123296 0.3755789 0.4376419     
## grm.VarC[1,5] 0.3032235 0.3551898 0.4055261     
## grm.VarC[1,6] 0.2443023 0.2692314 0.2969891     
## 
## Model Statistics: 
##                |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
##        Model:              6                   5994               27889.4
##    Saturated:              9                   5991                    NA
## Independence:              6                   5994                    NA
## Number of observations/statistics: 2000/6000
## 
## Information Criteria: 
##       |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
## AIC:       15901.40               27901.40                 27901.44
## BIC:      -17670.41               27935.01                 27915.94
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2024-03-05 09:29:13 
## Wall clock time: 12.10256 secs 
## optimizer:  SLSQP 
## OpenMx version number: 2.21.11 
## Need help?  See help(mxSummary)

store results

res = rbind(res, c("Twin_GRM", fitACE$output$confidenceIntervals[1,
    ], fitACE$output$confidenceIntervals[2, ], fitACE$output$confidenceIntervals[3,
    ]))

ACE - DZ only

# ----------------------------------------------------------------------------------------------------------------------
# Script: 05_relatednessDZonlu.R Author: Sarah Medland
# Date: 29 02 2020 Twin Univariate ACE model to estimate
# causes of variation across multiple groups Matrix style
# model - Raw data - Continuous data
# -------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------|

# Set Working Directory

# Clear Workspace rm(list=ls())

# Load Libraries & Options
library(OpenMx)

# ----------------------------------------------------------------------------------------------------------------------
# PREPARE DATA Load Data
df <- read.table("DZonly10k.txt", header = T)

# Select Variables for Analysis
nv <- 1  # number of traits
nt <- 3  # number of individuals
ntv <- nv * nt  # number of total variables
selVars <- c("Twin1", "Twin2", "Sib")  # name of traits
covVars <- c("age1", "age2", "age3", "sex1", "sex2", "sex3")  # list of covariates

# Check Correlations of Data for Analysis
cov(df[, c(selVars)], use = "pairwise.complete.obs")
##          Twin1    Twin2      Sib
## Twin1 8.534289 4.319766 4.292122
## Twin2 4.319766 8.595571 4.365393
## Sib   4.292122 4.365393 8.704640
# Set Starting Values
sMu <- 0  # start value for means
sVa <- 0.3  # start value for A
sVc <- 0.3  # start value for C
sVe <- 0.3  # start value for E

# ----------------------------------------------------------------------------------------------------------------------
# PREPARE MODEL

# Create Algebra for expected Mean Matrices
intercept <- mxMatrix(type = "Full", nrow = 1, ncol = ntv, free = TRUE,
    values = sMu, labels = "interC", name = "intercept")
betaS <- mxMatrix(type = "Full", nrow = 1, ncol = nv, free = TRUE,
    values = 0, labels = "betaS", name = "bS")
betaA <- mxMatrix(type = "Full", nrow = 1, ncol = nv, free = TRUE,
    values = 0, labels = "betaA", name = "bA")
defSex <- mxMatrix(type = "Full", nrow = 1, ncol = nt, free = FALSE,
    labels = c("data.sex1", "data.sex2", "data.sex3"), name = "Sex")
defAge <- mxMatrix(type = "Full", nrow = 1, ncol = nt, free = FALSE,
    labels = c("data.age1", "data.age2", "data.age3"), name = "Age")
expMean <- mxAlgebra(expression = intercept + Sex %x% bS + Age %x%
    bA, name = "expMean")

# Create Matrices for Variance Components
covA <- mxMatrix(type = "Symm", nrow = nv, ncol = nv, free = TRUE,
    values = sVa, label = "VA11", name = "VA")
covC <- mxMatrix(type = "Symm", nrow = nv, ncol = nv, free = TRUE,
    values = sVc, label = "VC11", name = "VC")
covE <- mxMatrix(type = "Symm", nrow = nv, ncol = nv, free = TRUE,
    values = sVe, label = "VE11", name = "VE")

# Create Algebra for expected Variance/Covariance Matrices
# in MZ & DZ twins
covP <- mxAlgebra(expression = VA + VC + VE, name = "V")
relA <- mxMatrix(type = "Stand", nrow = nt, ncol = nt, free = FALSE,
    labels = c("data.s1", "data.s2", "data.s3"), name = "rA")
relC <- mxMatrix(type = "Unit", nrow = nt, ncol = nt, free = FALSE,
    name = "rC")
relE <- mxMatrix(type = "Iden", nrow = nt, ncol = nt, free = FALSE,
    name = "rE")

expCov <- mxAlgebra(expression = rA %x% VA + rC %x% VC + rE %x%
    VE, name = "expCov")

# Create Data Objects for Multiple Groups
dataTW <- mxData(observed = df, type = "raw")

# Create Expectation Objects for Multiple Groups
expTW <- mxExpectationNormal(covariance = "expCov", means = "expMean",
    dimnames = selVars)
funML <- mxFitFunctionML()

# Create Model Objects for Multiple Groups
defs <- list(defAge, defSex, relA)
pars <- list(intercept, betaS, betaA, covA, covC, covE, covP,
    relC, relE)

# Create Algebra for Variance Components
rowVC <- rep("VC", nv)
colVC <- rep(c("VA", "VC", "VE", "SA", "SC", "SE"), each = nv)
estVC <- mxAlgebra(expression = cbind(VA, VC, VE, VA/V, VC/V,
    VE/V), name = "VarC", dimnames = list(rowVC, colVC))

# Create Confidence Interval Objects
ciACE <- mxCI("VarC[1,4:6]")

# Build Model with Confidence Intervals
modelACE <- mxModel("DZonly", pars, defs, expMean, expCov, expTW,
    dataTW, funML, estVC, ciACE)

# ----------------------------------------------------------------------------------------------------------------------
# RUN MODEL

# Run ACE Model
fitACE <- mxRun(modelACE, intervals = T)
sumACE <- summary(fitACE)
sumACE
## Summary of DZonly 
##  
## The Hessian at the solution does not appear to be convex. See ?mxCheckIdentification for possible diagnosis (Mx status RED). 
##  
## free parameters:
##     name    matrix row col   Estimate  Std.Error A
## 1 interC intercept   1   1 -9.2148352 0.36918169  
## 2  betaS        bS   1   1  0.1247789 0.02677999  
## 3  betaA        bA   1   1  0.4588013 0.01837662  
## 4   VA11        VA   1   1  0.3358636         NA  
## 5   VC11        VC   1   1  3.9003608         NA  
## 6   VE11        VE   1   1  4.0855831         NA  
## 
## confidence intervals:
##                      lbound   estimate    ubound note
## DZonly.VarC[1,4] -0.2001905 0.04035946 0.2814802     
## DZonly.VarC[1,5]  0.3472009 0.46869154 0.5891498     
## DZonly.VarC[1,6]  0.3701761 0.49094900 0.6127423     
## 
## Model Statistics: 
##                |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
##        Model:              6                  29994              142099.8
##    Saturated:              9                  29991                    NA
## Independence:              6                  29994                    NA
## Number of observations/statistics: 10000/30000
## 
## 
## ** Information matrix is not positive definite (not at a candidate optimum).
##   Be suspicious of these results. At minimum, do not trust the standard errors.
## 
## Information Criteria: 
##       |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
## AIC:       82111.82               142111.8                 142111.8
## BIC:     -134155.13               142155.1                 142136.0
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2024-03-05 09:30:37 
## Wall clock time: 83.06414 secs 
## optimizer:  SLSQP 
## OpenMx version number: 2.21.11 
## Need help?  See help(mxSummary)

store results

res = rbind(res, c("DZonly_GRM", fitACE$output$confidenceIntervals[1,
    ], fitACE$output$confidenceIntervals[2, ], fitACE$output$confidenceIntervals[3,
    ]))

AE on unrelated - GREML

# Copyright 2019-2022 by Robert M. Kirkpatrick
# Licensed under CC BY 4.0 <http://creativecommons.org/licenses/by/4.0/>

# This is an adaptation of a script previously used at the 2017 AGES Workshop at Virginia Commonwealth University.

#This script demonstrates the use of the GREML feature in a simple but realistic example.
#It first loads a genomic-relatedness matrix (GRM) and a file containing the phenotype and a covariate.  Then, it
#fits a simple GREML model to estimate additive-genetic variance, residual variance, and heritability.
#The GRM is 2000x2000 and was simulated from the phenotype
#Note that this analysis is simple enough that it could be done in GCTA (and GCTA would do it more quickly).

require(OpenMx)

options(mxCondenseMatrixSlots=TRUE)  #<--Saves memory
#The MxModel in this script has only 2 explicit free parameters, so there's almost nothing to be gained by
#setting the number of threads above 2:
mxOption(NULL,"Number of threads",2)

#Total number of participants:
N <- 2000
# Load and check data: ##################################################################

#Load the GRM's data.  Argument 'prefix' to omxReadGRMBin() should be everything in the filename and path of one of the GRM's
#that precedes the ".grm" :
grmstuff <- omxReadGRMBin(prefix="GRM",returnList=T)
#closeAllConnections()
str(grmstuff)
## List of 4
##  $ diag: num [1:2000] 1.001 1.01 0.99 1.005 0.997 ...
##  $ off : num [1:1999000] 0.00259 0.00302 0.00387 -0.00322 -0.00408 ...
##  $ id  :'data.frame':    2000 obs. of  2 variables:
##   ..$ V1: chr [1:2000] "ID1" "ID2" "ID3" "ID4" ...
##   ..$ V2: chr [1:2000] "ID1" "ID2" "ID3" "ID4" ...
##  $ N   : num 5e+05
#^^^grmstuff will be a list of 4 elements:
#     $diag is a numeric vector containing the GRM's diagonal elements.
#     $off is a numeric vector containing the GRM's off-diagonal elements.
#     $id is a dataframe that contains the family and individual IDs corresponding to the rows and columns of the GRM.
#     $N is the number of markers used to compute the GRM.

#Use the elements of grmstuff to make a proper GRM:
#Initialize:
GRM <- matrix(
    0.0,nrow=N,ncol=N,
    dimnames=list(grmstuff$id$V1,grmstuff$id$V1))
#^^^We're giving the GRM dimnames based on IDs.
#Populate the upper triangle of the GRM, excluding diagonal, with grmstuff$off:
GRM[!lower.tri(GRM,diag=T)] <- grmstuff$off
#Populate the lower triangle, excluding diagonal, of the GRM:
GRM <- GRM + t(GRM)
#Populate the diagonal of the GRM:
diag(GRM) <- grmstuff$diag

#Look at upper left corner of GRM:
hist(GRM[lower.tri(GRM)])

range(GRM[lower.tri(GRM)])
## [1] -0.02166385  0.02240862
# Look at the diagonal of GRM
hist(diag(GRM))

#Load phenotype file:
df <- read.table("twinsibData.txt", header=T)
# We add IDs for Twin1- so they match those from the GRM
df$fid<-paste0("ID", 1:2000)
df$id<-paste0("ID", 1:2000)

# We only keep the variables of interest in the analyses
phenfile <- df[,c("fid", "id", "Twin1",  "age1", "sex1")]

N <- 2000
#The columns are family ID, individual ID, phenotype, and covariate:
head(phenfile)
####################
# Some checks
#There is one missing value for the phenotype and for the covariate:
sum(is.na(phenfile$y))
## [1] 0
sum(is.na(phenfile$x))
## [1] 0
#Give phenfile rownames based on IDs:
#rownames(phenfile) <- phenfile$famid + phenfile$id
rownames(phenfile) <- phenfile$id

#Check to make sure that there are no individuals represented in the GRM but not in the phenotype file (there won't
#be, since the data are simulated, but this is a good check to make when working with real data):
all(rownames(GRM) %in% rownames(phenfile)) #<--Should return TRUE.
## [1] TRUE
#And check vice versa:
all(rownames(phenfile) %in% rownames(GRM)) #<--Should also be TRUE.
## [1] TRUE
#In fact, the rows of phenfile are matched with those of the GRM:
all(rownames(phenfile)==rownames(GRM)) #<--TRUE.
## [1] TRUE
# Use mxGREMLDataHandler(): ################################################

#Give the data-handler the dataset, and tell it that the one phenotype has column label 'y',
#that the one covariate has column label 'x', and that a lead column of ones needs to be appended to the 'X'
#matrix (for the intercept):
gremldat <- mxGREMLDataHandler(data=phenfile, yvars="Twin1", Xvars="age1", addOnes=T)
#We can see that participants #7 and #77 in phenfile were dropped due to missing data:
str(gremldat)
## List of 2
##  $ yX         : num [1:2000, 1:3] -3.666 3.297 -1.219 0.385 1.458 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:3] "Twin1" "Intrcpt" "age1"
##  $ casesToDrop: int(0)
#We can see that the first column of gremldata$yX is phenotype vector y, and that the remaining columns constitute
#the X matrix--a regression constant and the covariate:
head(gremldat$yX)
##           Twin1 Intrcpt     age1
## [1,] -3.6657747       1 19.96535
## [2,]  3.2971352       1 19.89334
## [3,] -1.2191284       1 19.87797
## [4,]  0.3845710       1 19.75907
## [5,]  1.4577869       1 20.32925
## [6,]  0.8721012       1 20.23061
#Number of datapoints, post-data-handling:
Np <- nrow(gremldat$yX)
#Drop rows and columns #7 and #77 from the GRM:

# Create MxData object, expectation, fitfunction, and compute plan: ##########################################

#The MxData object.  We will use gremldat$yX as the raw dataset to be used in our MxModel:
mxdat <- mxData(observed=gremldat$yX, type="raw", sort=FALSE)
#^^^N.B. the 'sort=FALSE' argument!  I CANNOT EMPHASIZE ENOUGH HOW IMPORTANT THAT IS!!!  By default, OpenMx
#automatically sorts the rows of a raw dataset at runtime, to try to achieve a perfomance boost.  However, we
#currently have the rows of y and X sorted exactly the same as in the GRM, so we do not want OpenMx to do any
#sorting!

#We  tell the GREML espectation that the name of the model-expected covariance matrix is "V", and that the dataset
#being input is y horizontally adhered to X (and therefore, it should not call the data-handler at runtime):
ge <- mxExpectationGREML(V="V",dataset.is.yX=TRUE)
#^^^Note: mxExpectationGREML() has argument 'casesToDropFromV', to which we could have passed
#gremldat$casesToDrop.  In that case, we could define our covariance matrix and its derivatives as
#full-size NxN matrices, and allow the backend to use the value of 'casesToDropFromV' to automatically
#trim rows and columns in those matrices that correspond to missing observations (i.e., participants #7 and #77).
#However, that backend  matrix-resizing process carries a performance cost.

#The GREML fitfunction tells OpenMx that the derivative of 'V' with respect to free parameter
# 'va'(the additive-genetic variance) is a matrix named 'A', and that the derivative of 'V' w/r/t free parameter
# 've' is a matrix named 'I'.  At runtime, the GREML fitfunction will use these derivatives to help with
#optimization and to compute standard errors:
gff <- mxFitFunctionGREML(dV=c(va="A",ve="I"))

#This is a custom compute plan.  It is necessary because we want to use the Newton-Raphson optimizer, which
#can use analytic first and second derivatives of the GREML fitfunction to speed up convergence:
plan <- mxComputeSequence(
    #List of steps.  The first is the most notable.  It tells OpenMx to use the Newton-Raphson optimizer, and to
    #display what it's doing in real time (via the verbose argument):
    steps=list(
        mxComputeNewtonRaphson(verbose=5L),
        mxComputeOnce("fitfunction", c("gradient","hessian")),
        mxComputeStandardError(),
        mxComputeHessianQuality(),
        mxComputeReportDeriv(),
        mxComputeReportExpectation()
))
#^^^Note:  If you are running the R GUI under Windows, delete the 'verbose=5L' argument in the above.


# Create the MxModel, and run it: ###############################################################################
#(We will create a lot of objects inside the mxModel() statement.  This helps to save memory.)

testmod <- mxModel(
    "GREML_1GRM_1trait", #<--Model name
    #1x1 matrix containing nonshared-environmental (residual) variance component:
    mxMatrix(type = "Full", nrow = 1, ncol=1, free=T, values = var(gremldat$yX[,"Twin1"])/2, labels = "ve",
                     lbound = 0.0001, name = "Ve"),
    #1x1 matrix containing additive-genetic variance component:
    mxMatrix(type = "Full", nrow = 1, ncol=1, free=T, values = var(gremldat$yX[,"Twin1"])/2, labels = "va",
                     lbound=0, name = "Va"), #<--Note the lower bound on zero.
    #998x998 identity matrix--the "relatedness matrix" for the residuals:
    mxMatrix("Iden",nrow=Np,name="I"),
    #The GRM:
    mxMatrix("Symm",nrow=Np,free=F,values=GRM,name="A"),
    #The model-expected covariance matrix:
    mxAlgebra((A%x%Va) + (I%x%Ve), name="V"),
    #An MxAlgebra for the heritability:
    mxAlgebra(Va/(Va+Ve), name="h2"),
    mxdat, #<--MxData object
    ge, #<--GREML expectation
    gff, #<--GREML fitfunction
    plan #<--Custom compute plan
)


#Run and view summary:
testrun <- mxRun(testmod)
#testrun<-mxTryHard(testrun, extraTries = 10)
summary(testrun)
## Summary of GREML_1GRM_1trait 
##  
## free parameters:
##   name matrix row col Estimate Std.Error lbound ubound
## 1   ve     Ve   1   1 7.096399  1.294151  1e-04       
## 2   va     Va   1   1 1.263205  1.282592      0       
## 
## regression coefficients:
##      name     coeff         se
## 1 Intrcpt -8.825507 1.82815749
## 2    age1  0.452379 0.09144978
## 
## Model Statistics: 
##                |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
##        Model:              4                   1996              9929.998
##    Saturated:             NA                     NA                    NA
## Independence:             NA                     NA                    NA
## Number of observations/statistics: 1/2000
## 
## Information Criteria: 
##       |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
## AIC:       5937.998               9933.998                 9927.998
## BIC:       9929.998               9929.998                 9925.839
## 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: 2024-03-05 09:31:02 
## Wall clock time: 21.45609 secs 
## OpenMx version number: 2.21.11 
## Need help?  See help(mxSummary)
mxEval(h2,testrun,T)
##           [,1]
## [1,] 0.1511082
#If we had added the appropriate additional step to the compute plan, we could have requested
#confidence intervals for va, ve, and/or h2 (but not for the regression coefficients).
#We have asymptotic-ML-theory standard errors for va and ve.  It is possible to obtain an SE for h2
#using a delta-method approximation, though its accuracy is quesionable given the fairly small N:
scm <- testrun$output$vcov #<--Sampling covariance matrix for ve and va
pointest <- testrun$output$estimate #<--Point estimates of ve and va
h2se <- sqrt(
    (pointest[2]/(pointest[1]+pointest[2]))^2 * (
        (scm[2,2]/pointest[2]^2) - (2*scm[1,2]/pointest[1]/(pointest[1]+pointest[2])) +
            (sum(scm)*(pointest[1]+pointest[2])^-2)
    ))
print(h2se)
##        va 
## 0.1575237
# 95% confidence intervals
mxEval(h2,testrun,T)-1.96*h2se
##            [,1]
## [1,] -0.1576382
mxEval(h2,testrun,T)+1.96*h2se
##           [,1]
## [1,] 0.4598546

store results

res = rbind(res, c("GREML", mxEval(h2, testrun, T) - 1.96 * h2se,
    h2se, mxEval(h2, testrun, T) + 1.96 * h2se, NA, NA, NA, NA,
    1 - h2se, NA))

Plot results

error.bar <- function(x, y, upper, lower=upper, length=0.1,...){
  arrows(x,upper, x, lower, angle=90, code=3, length=length, ...)
}
 

#png("A_C_estimates_allMethods.png", width = 20, height = 15, units = "cm", res = 400 )
par(mfrow=c(1,2), mar=c(10,4,1,1))
#layout(matrix(c(1,2,2), 1, 2, byrow = TRUE))
plot( res[,3], pch=20, cex=2, ylim=c(-0.1,1), ylab="A", xaxt='n', xlab="", col="darkblue")
error.bar(1:7, y = as.numeric(res[,3]), upper = as.numeric(res[,4]), lower = as.numeric(res[,2]), col="darkblue" )
grid()
abline(h=0, lwd=2, col="black")
text(x = 1:7,
     y = par("usr")[3] - 0.1,
     labels = c("ACE twins", "Twins and Sibs", "Alternate twins and sibs", "Zygosity as definition variable", "Measured genetic relationship", "DZ only", "AE unrelated - GREML"),
     xpd = NA,
     ## Rotate the labels by 35 degrees.
     srt = 45,
     cex = 0.9, adj=1)

#axis(1, at = c(1:7), las=2, labels = c("ACE twins", "Twins and Sibs", "Alternate twins and sibs", "Zygosity as definition variable", "Measured genetic relationship", "DZ only", "AE unrelated - GREML"), )

plot( res[,6], pch=20, cex=2, ylim=c(-0.1,1), ylab="C", xaxt='n', xlab="", col="goldenrod4")
error.bar(1:7, y = as.numeric(res[,6]), upper = as.numeric(res[,7]), lower = as.numeric(res[,5]), col="goldenrod4" )
grid()
abline(h=0, lwd=2, col= "black")
text(x = 1:7,
     y = par("usr")[3] - 0.1,
     labels = c("ACE twins", "Twins and Sibs", "Alternate twins and sibs", "Zygosity as definition variable", "Measured genetic relationship", "DZ only", "AE unrelated - GREML"),
     xpd = NA,
     ## Rotate the labels by 35 degrees.
     srt = 45,
     cex = 0.9, adj=1)

#dev.off()

Show values

colnames(res) = c("modelName", "A_lb", "A", "A_ub", "C_lb", "C",
    "C_ub", "E_lb", "E", "E_ub")
print(res)
##      modelName    A_lb                 A                    A_ub               
## [1,] "Twin"       "0.403170946807557"  "0.498300061082194"  "0.599007353663766"
## [2,] "Twin_sib"   "0.312148768490362"  "0.375608350283173"  "0.437883993042671"
## [3,] "Twin_sib2"  "0.312148768490362"  "0.375608350283173"  "0.437883993042671"
## [4,] "ZygVar"     "0.312148762145735"  "0.375608351833094"  "0.437883993032954"
## [5,] "Twin_GRM"   "0.31232964872899"   "0.375578874300066"  "0.437641948870637"
## [6,] "DZonly_GRM" "-0.200190470850726" "0.0403594584376225" "0.281480163423856"
## [7,] "GREML"      "-0.157638196233984" "0.157523670059722"  "0.459854590400124"
##      C_lb                C                   C_ub               
## [1,] "0.144310315683289" "0.239290130920058" "0.326965639095612"
## [2,] "0.302870332983627" "0.354969872444485" "0.40546393501121" 
## [3,] "0.302870332983627" "0.354969872444485" "0.40546393501121" 
## [4,] "0.302870332546734" "0.354969871241579" "0.405463933962984"
## [5,] "0.303223504106537" "0.355189767536258" "0.405526105479709"
## [6,] "0.347200873643996" "0.468691540576436" "0.589149814869066"
## [7,] NA                  NA                  NA                 
##      E_lb                E                   E_ub               
## [1,] "0.238426708822721" "0.262409807997747" "0.288912818053576"
## [2,] "0.244355651981887" "0.269421777272342" "0.29724908252969" 
## [3,] "0.244355651981887" "0.269421777272342" "0.29724908252969" 
## [4,] "0.244355651400181" "0.269421776925327" "0.297249082347298"
## [5,] "0.244302306626096" "0.269231358163676" "0.29698910698171" 
## [6,] "0.370176113774172" "0.490949000985941" "0.612742289647207"
## [7,] NA                  "0.842476329940278" NA
 



A work by Baptiste Couvy-Duchesne