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)

dfBin <- df
dfBin$Twin1 <- ifelse(dfBin$Twin1 > 0, 1, 0)
dfBin$Twin2 <- ifelse(dfBin$Twin2 > 0, 1, 0)
dfBin$Sib <- ifelse(dfBin$Sib > 0, 1, 0)

# Check the data
table(dfBin$Twin1)
## 
##    0    1 
##  921 1079
table(dfBin$Twin2)
## 
##    0    1 
##  940 1060
table(dfBin$Sib)
## 
##    0    1 
##  949 1051
# Ensure that the data are encoded as an ordered factor
dfBin$Twin1 <- mxFactor(dfBin$Twin1, levels = 0:1)
dfBin$Twin2 <- mxFactor(dfBin$Twin2, levels = 0:1)
dfBin$Sib <- mxFactor(dfBin$Sib, levels = 0:1)

mzData <- subset(dfBin, zygosity == 1)
dzData <- subset(dfBin, zygosity == 2)

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

# Set Starting Values
sVa <- 0.3  # start value for A
sVc <- 0.3  # start value for C
sVe <- 0.3  # start value for E
ThrVals <- 0  # start for the threshold

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

# Create Algebra for expected Mean Matrices
intercept <- mxMatrix(type = "Full", nrow = 1, ncol = ntv, free = FALSE,
    values = 0, 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")
expThr <- mxMatrix(type = "Full", nrow = nTH, ncol = ntv, free = TRUE,
    values = ThrVals, labels = paste("th", 1:nTH, sep = ""),
    name = "expThr")

# 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")
cons <- mxConstraint(VA + VC + VE == 1, name = "cons")  # Constraint so that the variance equals 1 (necessary for identification)

# 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",
    thresholds = "expThr", dimnames = selVars)
expDZ <- mxExpectationNormal(covariance = "expCovDZ", means = "expMean",
    thresholds = "expThr", dimnames = selVars)
funML <- mxFitFunctionML()

# Create Model Objects for Multiple Groups
defs <- list(defAge, defSex)
pars <- list(intercept, betaS, betaA, expThr, 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, cons)

# ----------------------------------------------------------------------------------------------------------------------
# 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 betaS     bS   1   1 0.08950213 0.04254163  
## 2 betaA     bA   1   1 0.19623496 0.03348751  
## 3   th1 expThr   1   1 3.87646753 0.66993395  
## 4  VA11     VA   1   1 0.53811726 0.10498124  
## 5  VC11     VC   1   1 0.17754657 0.09121863  
## 6  VE11     VE   1   1 0.28433617 0.02999880  
## 
## confidence intervals:
##                       lbound  estimate    ubound note
## ACEvc.VarC[1,4]  0.333397704 0.5381173 0.7456443     
## ACEvc.VarC[1,5] -0.006892831 0.1775466 0.3512009     
## ACEvc.VarC[1,6]  0.230272307 0.2843362 0.3471491     
## 
## Model Statistics: 
##                |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
##        Model:              6                   3995              5119.589
##    Saturated:             NA                     NA                    NA
## Independence:             NA                     NA                    NA
## Number of observations/statistics: 2000/4001
## 
## Constraint 'cons' contributes 1 observed statistic. 
## 
## Information Criteria: 
##       |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
## AIC:      -2870.411               5131.589                 5131.631
## BIC:     -25246.016               5165.195                 5146.132
## 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:32:22 
## Wall clock time: 13.05658 secs 
## optimizer:  SLSQP 
## OpenMx version number: 2.21.11 
## Need help?  See help(mxSummary)
# Some different ways to extract output
fitACE$output$estimate
##      betaS      betaA        th1       VA11       VC11       VE11 
## 0.08950213 0.19623496 3.87646753 0.53811726 0.17754657 0.28433617
fitACE$VarC$result
##           VA        VC        VE        SA        SC        SE
## VC 0.5381173 0.1775466 0.2843362 0.5381173 0.1775466 0.2843362
fitACE$output$confidenceIntervals
##                       lbound  estimate    ubound
## ACEvc.VarC[1,4]  0.333397704 0.5381173 0.7456443
## ACEvc.VarC[1,5] -0.006892831 0.1775466 0.3512009
## ACEvc.VarC[1,6]  0.230272307 0.2843362 0.3471491

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)

dfBin <- df
dfBin$Twin1 <- ifelse(dfBin$Twin1 > 0, 1, 0)
dfBin$Twin2 <- ifelse(dfBin$Twin2 > 0, 1, 0)
dfBin$Sib <- ifelse(dfBin$Sib > 0, 1, 0)

dfBin$Twin1 <- mxFactor(dfBin$Twin1, levels = 0:1)
dfBin$Twin2 <- mxFactor(dfBin$Twin2, levels = 0:1)
dfBin$Sib <- mxFactor(dfBin$Sib, levels = 0:1)

mzData <- subset(dfBin, zygosity == 1)
dzData <- subset(dfBin, zygosity == 2)

# Select Variables for Analysis
nv <- 1  # number of traits
nt <- 3  # number of individuals
nTH <- 1  # number of thresholds
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


# Set Starting Values
sVa <- 0.3  # start value for A
sVc <- 0.3  # start value for C
sVe <- 0.3  # start value for E
ThrVals <- 0  # start for the threshold

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

# Create Algebra for expected Mean Matrices
intercept <- mxMatrix(type = "Full", nrow = 1, ncol = ntv, free = FALSE,
    values = 0, 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")
expThr <- mxMatrix("Full", nTH, ntv, free = TRUE, values = ThrVals,
    labels = paste("th", 1:nTH, sep = ""), name = "expThr")

# 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")
cons <- mxConstraint(VA + VC + VE == 1, name = "cons")  # Constraint so that the variance equals 1 (necessary for identification)

# 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",
    thresholds = "expThr", dimnames = selVars)
expDZ <- mxExpectationNormal(covariance = "expCovDZ", means = "expMean",
    thresholds = "expThr", dimnames = selVars)
funML <- mxFitFunctionML()

# Create Model Objects for Multiple Groups
defs <- list(defAge, defSex)
pars <- list(intercept, betaS, betaA, expThr, 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, cons)

# ----------------------------------------------------------------------------------------------------------------------
# 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 betaS     bS   1   1 0.1179450 0.030981258  
## 2 betaA     bA   1   1 0.1049937 0.007750931  
## 3   th1 expThr   1   1 2.0778973 0.156340424  
## 4  VA11     VA   1   1 0.2543580 0.071152596 !
## 5  VC11     VC   1   1 0.4331444 0.049126864  
## 6  VE11     VE   1   1 0.3124976 0.031818950  
## 
## confidence intervals:
##                     lbound  estimate    ubound note
## ACEsib.VarC[1,4] 0.1130752 0.2543580 0.3913480     
## ACEsib.VarC[1,5] 0.3349206 0.4331444 0.5283482     
## ACEsib.VarC[1,6] 0.2537166 0.3124976 0.3784697     
## 
## Model Statistics: 
##                |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
##        Model:              6                   5995              7428.427
##    Saturated:             NA                     NA                    NA
## Independence:             NA                     NA                    NA
## Number of observations/statistics: 2000/6001
## 
## Constraint 'cons' contributes 1 observed statistic. 
## 
## Information Criteria: 
##       |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
## AIC:      -4561.573               7440.427                 7440.469
## BIC:     -38138.983               7474.032                 7454.970
## 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:32:57 
## Wall clock time: 35.15035 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)

dfBin <- df
dfBin$Twin1 <- ifelse(dfBin$Twin1 > 0, 1, 0)
dfBin$Twin2 <- ifelse(dfBin$Twin2 > 0, 1, 0)
dfBin$Sib <- ifelse(dfBin$Sib > 0, 1, 0)

dfBin$Twin1 <- mxFactor(dfBin$Twin1, levels = 0:1)
dfBin$Twin2 <- mxFactor(dfBin$Twin2, levels = 0:1)
dfBin$Sib <- mxFactor(dfBin$Sib, levels = 0:1)

mzData <- subset(dfBin, zygosity == 1)
dzData <- subset(dfBin, zygosity == 2)

# Select Variables for Analysis
nv <- 1  # number of traits
nt <- 3  # number of individuals
nTH <- 1  # number of thresholds
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

# Set Starting Values
sVa <- 0.3  # start value for A
sVc <- 0.3  # start value for C
sVe <- 0.3  # start value for E
ThrVals <- 0  # start for the threshold

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

# Create Algebra for expected Mean Matrices
intercept <- mxMatrix(type = "Full", nrow = 1, ncol = ntv, free = FALSE,
    values = 0, 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")
expThr <- mxMatrix("Full", nrow = nTH, ncol = ntv, free = TRUE,
    values = ThrVals, labels = paste("th", 1:nTH, sep = ""),
    name = "expThr")

# 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")
cons <- mxConstraint(VA + VC + VE == 1, name = "cons")  # Constraint so that the variance equals 1 (necessary for identification)

# 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",
    thresholds = "expThr", dimnames = selVars)
expDZ <- mxExpectationNormal(covariance = "expCovDZ", means = "expMean",
    thresholds = "expThr", dimnames = selVars)
funML <- mxFitFunctionML()

# Create Model Objects for Multiple Groups
defs <- list(defAge, defSex)
pars <- list(intercept, betaS, betaA, expThr, 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, cons)

# ----------------------------------------------------------------------------------------------------------------------
# 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 betaS     bS   1   1 0.1179450 0.030981258  
## 2 betaA     bA   1   1 0.1049937 0.007750931  
## 3   th1 expThr   1   1 2.0778973 0.156340424  
## 4  VA11     VA   1   1 0.2543580 0.071152596 !
## 5  VC11     VC   1   1 0.4331444 0.049126864  
## 6  VE11     VE   1   1 0.3124976 0.031818950  
## 
## confidence intervals:
##                         lbound  estimate    ubound note
## ACEsib_alt.VarC[1,4] 0.1130752 0.2543580 0.3913480     
## ACEsib_alt.VarC[1,5] 0.3349206 0.4331444 0.5283482     
## ACEsib_alt.VarC[1,6] 0.2537166 0.3124976 0.3784697     
## 
## Model Statistics: 
##                |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
##        Model:              6                   5995              7428.427
##    Saturated:             NA                     NA                    NA
## Independence:             NA                     NA                    NA
## Number of observations/statistics: 2000/6001
## 
## Constraint 'cons' contributes 1 observed statistic. 
## 
## Information Criteria: 
##       |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
## AIC:      -4561.573               7440.427                 7440.469
## BIC:     -38138.983               7474.032                 7454.970
## 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:33:33 
## Wall clock time: 35.12386 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)

dfBin <- df
dfBin$Twin1 <- ifelse(dfBin$Twin1 > 0, 1, 0)
dfBin$Twin2 <- ifelse(dfBin$Twin2 > 0, 1, 0)
dfBin$Sib <- ifelse(dfBin$Sib > 0, 1, 0)

dfBin$Twin1 <- mxFactor(dfBin$Twin1, levels = 0:1)
dfBin$Twin2 <- mxFactor(dfBin$Twin2, levels = 0:1)
dfBin$Sib <- mxFactor(dfBin$Sib, levels = 0:1)

# Select Variables for Analysis
nv <- 1  # number of traits
nt <- 3  # number of individuals
nTH <- 1  # number of thresholds
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

# Set Starting Values
sVa <- 0.3  # start value for A
sVc <- 0.3  # start value for C
sVe <- 0.3  # start value for E
ThrVals <- 0  # start for the threshold

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

# Create Algebra for expected Mean Matrices
intercept <- mxMatrix(type = "Full", nrow = 1, ncol = ntv, free = FALSE,
    values = 0, 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")
expThr <- mxMatrix("Full", nrow = nTH, ncol = ntv, free = TRUE,
    values = ThrVals, labels = paste("th", 1:nTH, sep = ""),
    name = "expThr")

# 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")
cons <- mxConstraint(VA + VC + VE == 1, name = "cons")  # Constraint so that the variance equals 1 (necessary for identification)

# 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 = nt, free = FALSE,
    labels = c("data.zygT", "data.zygS", "data.zygS"), name = "rA")
relC <- mxMatrix(type = "Unit", nrow = ntv, ncol = nt, free = FALSE,
    name = "rC")
relE <- mxMatrix(type = "Iden", nrow = ntv, ncol = nt, free = FALSE,
    name = "rE")

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

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

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

# Create Model Objects
defs <- list(defAge, defSex, relA)
pars <- list(intercept, betaS, betaA, expThr, 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, cons)

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

# Run ACE Model
fitACE <- mxRun(modelACE, intervals = T)
sumACE <- summary(fitACE)
sumACE
## Summary of zygdef 
##  
## The model does not satisfy the first-order optimality conditions to the required accuracy, and no improved point for the merit function could be found during the final linesearch (Mx status RED) 
##  
## Your ordinal model may converge if you reduce mvnRelEps
##    try: model <- mxOption(model, 'mvnRelEps', mxOption(model, 'mvnRelEps')/5)
## 
## free parameters:
##    name matrix row   col  Estimate   Std.Error A
## 1 betaS     bS   1     1 0.1175745 0.032170065  
## 2 betaA     bA   1     1 0.1038757 0.008782796  
## 3   th1 expThr   1 Twin1 2.0547080 0.175246885  
## 4  VA11     VA   1     1 0.2460526 0.084154822 !
## 5  VC11     VC   1     1 0.4379598 0.055726869  
## 6  VE11     VE   1     1 0.3159876 0.036761189 !
## 
## confidence intervals:
##                     lbound  estimate    ubound note
## zygdef.VarC[1,4] 0.1123027 0.2460526 0.3914871     
## zygdef.VarC[1,5] 0.3356056 0.4379598 0.5286203     
## zygdef.VarC[1,6] 0.2536727 0.3159876 0.3789370     
## 
## Model Statistics: 
##                |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
##        Model:              6                   5995              7428.444
##    Saturated:              9                   5992                    NA
## Independence:              6                   5995                    NA
## Number of observations/statistics: 2000/6001
## 
## Constraint 'cons' contributes 1 observed statistic. 
## 
## Information Criteria: 
##       |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
## AIC:      -4561.556               7440.444                 7440.486
## BIC:     -38138.966               7474.049                 7454.987
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2024-03-05 09:34:12 
## Wall clock time: 38.71569 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)

dfBin <- df
dfBin$T1 <- ifelse(dfBin$Twin1 > 0, 1, 0)
dfBin$T2 <- ifelse(dfBin$Twin2 > 0, 1, 0)
dfBin$T3 <- ifelse(dfBin$Sib > 0, 1, 0)

dfBin$T1 <- mxFactor(dfBin$T1, levels = 0:1)
dfBin$T2 <- mxFactor(dfBin$T2, levels = 0:1)
dfBin$T3 <- mxFactor(dfBin$T3, levels = 0:1)

# Select Variables for Analysis
nv <- 1  # number of traits
nt <- 3  # number of individuals
nTH <- 1  # number of thresholds
ntv <- nv * nt  # number of total variables
selVars <- c("T1", "T2", "T3")  # 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)

# Set Starting Values
sVa <- 0.3  # start value for A
sVc <- 0.3  # start value for C
sVe <- 0.3  # start value for E
ThrVals <- 0  # start for the threshold

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

# Create Algebra for expected Mean Matrices
intercept <- mxMatrix(type = "Full", nrow = 1, ncol = ntv, free = FALSE,
    values = 0, 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")
expThr <- mxMatrix("Full", nrow = nTH, ncol = ntv, free = TRUE,
    values = ThrVals, labels = paste("th", 1:nTH, sep = ""),
    name = "expThr")

# 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")
cons <- mxConstraint(VA + VC + VE == 1, name = "cons")  # Constraint so that the variance equals 1 (necessary for identification)

# 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 = dfBin, type = "raw")

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

# Create Model Objects for Multiple Groups
defs <- list(defAge, defSex, relA)
pars <- list(intercept, betaS, betaA, expThr, 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, cons)

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

# Run ACE Model
fitACE <- mxRun(modelACE, intervals = T)
sumACE <- summary(fitACE)
sumACE
## Summary of grm 
##  
## The model does not satisfy the first-order optimality conditions to the required accuracy, and no improved point for the merit function could be found during the final linesearch (Mx status RED) 
##  
## Your ordinal model may converge if you reduce mvnRelEps
##    try: model <- mxOption(model, 'mvnRelEps', mxOption(model, 'mvnRelEps')/5)
## 
## free parameters:
##    name matrix row col  Estimate   Std.Error A
## 1 betaS     bS   1   1 0.1177752 0.032796217  
## 2 betaA     bA   1   1 0.1054690 0.008532588  
## 3   th1 expThr   1  T1 2.0872152 0.170968052  
## 4  VA11     VA   1   1 0.2660027 0.078241442 !
## 5  VC11     VC   1   1 0.4256882 0.052907834 !
## 6  VE11     VE   1   1 0.3083091 0.034290837  
## 
## confidence intervals:
##                  lbound  estimate    ubound note
## grm.VarC[1,4] 0.1260427 0.2660027 0.4015862     
## grm.VarC[1,5] 0.3281679 0.4256882 0.5200095     
## grm.VarC[1,6] 0.2501881 0.3083091 0.3729450     
## 
## Model Statistics: 
##                |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
##        Model:              6                   5995              7426.989
##    Saturated:              9                   5992                    NA
## Independence:              6                   5995                    NA
## Number of observations/statistics: 2000/6001
## 
## Constraint 'cons' contributes 1 observed statistic. 
## 
## Information Criteria: 
##       |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
## AIC:      -4563.011               7438.989                 7439.031
## BIC:     -38140.421               7472.594                 7453.532
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2024-03-05 09:34:55 
## Wall clock time: 42.71655 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)

dfBin <- df
dfBin$Twin1 <- ifelse(dfBin$Twin1 > 0, 1, 0)
dfBin$Twin2 <- ifelse(dfBin$Twin2 > 0, 1, 0)
dfBin$Sib <- ifelse(dfBin$Sib > 0, 1, 0)

dfBin$Twin1 <- mxFactor(dfBin$Twin1, levels = 0:1)
dfBin$Twin2 <- mxFactor(dfBin$Twin2, levels = 0:1)
dfBin$Sib <- mxFactor(dfBin$Sib, levels = 0:1)

table(dfBin$Twin1)
## 
##    0    1 
## 4993 5007
table(dfBin$Twin2)
## 
##    0    1 
## 4997 5003
table(dfBin$Sib)
## 
##    0    1 
## 4982 5018
# Select Variables for Analysis
nv <- 1  # number of traits
nt <- 3  # number of individuals
nTH <- 1  # number of thresholds
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

# Set Starting Values
sVa <- 0.4  # start value for A
sVc <- 0.2  # start value for C
sVe <- 0.3  # start value for E
ThrVals <- 4  # start for the threshold

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

# Create Algebra for expected Mean Matrices
intercept <- mxMatrix(type = "Full", nrow = 1, ncol = ntv, free = FALSE,
    values = 0, labels = "interC", name = "intercept")
betaS <- mxMatrix(type = "Full", nrow = 1, ncol = nv, free = TRUE,
    values = 0.01, labels = "betaS", name = "bS")
betaA <- mxMatrix(type = "Full", nrow = 1, ncol = nv, free = TRUE,
    values = 0.1, 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")
expThr <- mxMatrix("Full", nrow = nTH, ncol = ntv, free = TRUE,
    values = ThrVals, labels = paste("th", 1:nTH, sep = ""),
    name = "expThr")

# 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")
cons <- mxConstraint(VA + VC + VE == 1, name = "cons")  # Constraint so that the variance equals 1 (necessary for identification)

# 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 = dfBin, type = "raw")

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

# Create Model Objects for Multiple Groups
defs <- list(defAge, defSex, relA)
pars <- list(intercept, betaS, betaA, expThr, 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", defs, pars, expMean, expCov, dataTW,
    expTW, funML, estVC, ciACE, cons)

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

# Run ACE Model
fitACE <- mxRun(modelACE, intervals = T)  # This can take several minutes to run with CIs
sumACE <- summary(fitACE)
sumACE
## Summary of DZonly 
##  
## The model does not satisfy the first-order optimality conditions to the required accuracy, and no improved point for the merit function could be found during the final linesearch (Mx status RED) 
##  
## Your ordinal model may converge if you reduce mvnRelEps
##    try: model <- mxOption(model, 'mvnRelEps', mxOption(model, 'mvnRelEps')/5)
## 
## free parameters:
##    name matrix row   col   Estimate   Std.Error A
## 1 betaS     bS   1     1 0.03054091 0.013243369  
## 2 betaA     bA   1     1 0.18054482 0.005609056  
## 3   th1 expThr   1 Twin1 3.62451105 0.112415848  
## 4  VA11     VA   1     1 0.01074084 0.236192647  
## 5  VC11     VC   1     1 0.48503394 0.118483642  
## 6  VE11     VE   1     1 0.50422521 0.118418721  
## 
## confidence intervals:
##                      lbound   estimate    ubound note
## DZonly.VarC[1,4] -0.4515270 0.01074084 0.4775872     
## DZonly.VarC[1,5]  0.2502554 0.48503394 0.7158699     
## DZonly.VarC[1,6]  0.2699095 0.50422521 0.7376346     
## 
## Model Statistics: 
##                |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
##        Model:              6                  29995              38192.36
##    Saturated:              9                  29992                    NA
## Independence:              6                  29995                    NA
## Number of observations/statistics: 10000/30001
## 
## Constraint 'cons' contributes 1 observed statistic. 
## 
## Information Criteria: 
##       |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
## AIC:      -21797.64               38204.36                 38204.36
## BIC:     -238071.80               38247.62                 38228.55
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2024-03-05 09:39:30 
## Wall clock time: 275.1607 secs 
## optimizer:  SLSQP 
## OpenMx version number: 2.21.11 
## Need help?  See help(mxSummary)
# mxOption(NULL,'Default optimizer', 'CSOLNP') fitACE <-
# mxTryHardOrdinal( modelACE, intervals=T )

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 1000x1000, and was computed from 50,000 simulated SNPs in linkage equilibrium.
#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)
#You need to set R's working directory to the directory containing the data files for this demo.
#(i.e., YOU MUST CHANGE THE NEXT LINE TO REFLECT WHERE, ON YOUR COMPUTER, YOU'VE PLACED THE DATA FILES):
#setwd("./AGES2017/data")

#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)
df$Twin1 <- ifelse(df$Twin1 > 0, 1, 0)

# 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] 0 1 0 1 1 1 0 0 1 1 ...
##   ..- 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,]     0       1 19.96535
## [2,]     1       1 19.89334
## [3,]     0       1 19.87797
## [4,]     1       1 19.75907
## [5,]     1       1 20.32925
## [6,]     1       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 0.18610296 0.03796054  1e-04       
## 2   va     Va   1   1 0.06031352 0.03791159      0       
## 
## regression coefficients:
##      name       coeff         se
## 1 Intrcpt -0.76259592 0.31364843
## 2    age1  0.06518182 0.01568965
## 
## Model Statistics: 
##                |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
##        Model:              4                   1996              2887.319
##    Saturated:             NA                     NA                    NA
## Independence:             NA                     NA                    NA
## Number of observations/statistics: 1/2000
## 
## Information Criteria: 
##       |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
## AIC:      -1104.681               2891.319                 2885.319
## BIC:       2887.319               2887.319                 2883.160
## 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:39:56 
## Wall clock time: 21.54765 secs 
## OpenMx version number: 2.21.11 
## Need help?  See help(mxSummary)
mxEval(h2,testrun,T)
##           [,1]
## [1,] 0.2447625
#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.1655626
# 95% confidence intervals
mxEval(h2,testrun,T)-1.96*h2se
##             [,1]
## [1,] -0.07974026
mxEval(h2,testrun,T)+1.96*h2se
##           [,1]
## [1,] 0.5692653

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_binary.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                   
## [1,] "Twin"       "0.333397703530672"   "0.538117255044524" 
## [2,] "Twin_sib"   "0.113075218713683"   "0.254357985817382" 
## [3,] "Twin_sib2"  "0.113075218713683"   "0.254357985817382" 
## [4,] "ZygVar"     "0.112302738969462"   "0.24605261509722"  
## [5,] "Twin_GRM"   "0.126042709921443"   "0.266002675893346" 
## [6,] "DZonly_GRM" "-0.451527015889042"  "0.0107408443839886"
## [7,] "GREML"      "-0.0797402587557465" "0.165562638441246" 
##      A_ub                C_lb                   C                  
## [1,] "0.745644310199822" "-0.00689283099769655" "0.177546570650915"
## [2,] "0.391347994284373" "0.334920580378078"    "0.433144405860925"
## [3,] "0.391347994284373" "0.334920580378078"    "0.433144405860925"
## [4,] "0.391487111614681" "0.335605611403928"    "0.437959812834118"
## [5,] "0.401586171160483" "0.328167879544956"    "0.425688194317904"
## [6,] "0.477587200459368" "0.250255415046221"    "0.485033941163498"
## [7,] "0.569265283933939" NA                     NA                 
##      C_ub                E_lb                E                  
## [1,] "0.35120094557072"  "0.230272306513428" "0.284336174304561"
## [2,] "0.528348151293424" "0.253716590615379" "0.312497608321694"
## [3,] "0.528348151293424" "0.253716590615379" "0.312497608321694"
## [4,] "0.528620281482147" "0.253672740785762" "0.315987572068662"
## [5,] "0.520009506808218" "0.250188076195774" "0.30830912978875" 
## [6,] "0.715869909033963" "0.269909534739407" "0.504225214452513"
## [7,] NA                  NA                  "0.834437361558754"
##      E_ub               
## [1,] "0.347149090163806"
## [2,] "0.378469731714189"
## [3,] "0.378469731714189"
## [4,] "0.378936977026472"
## [5,] "0.372945020062481"
## [6,] "0.737634648712551"
## [7,] NA
 



A work by Baptiste Couvy-Duchesne