---
title: "From ACE models to GCTA - continuous"
author: "[Baptiste Couvy-Duchesne](https://github.com/baptisteCD)"
date: "`r format(Sys.time(), '%d %B %Y')`"
mail: "baptiste.couvyduchesne@gmail.com"
github: "baptisteCD"
output:
  epuRate::ISGW:
    toc: TRUE
    number_sections: FALSE
    code_folding: "show"
---

```{r, echo=FALSE, message=FALSE, warning=FALSE}

# You need these libraries to run this template:
library(rmarkdown)    # install.packages("rmarkdown") 
# Get template
library(epuRate)      # devtools::install_github("baptisteCD/epuRate", force=TRUE)
   
# Soft-wrap code in knitting
knitr::opts_chunk$set(tidy.opts = list(width.cutoff = 60), tidy = TRUE)

```

<br><br><br><br>

# Environment

```{R, message=FALSE, warning=FALSE, eval=T}

library(OpenMx)

# Create table to store all ACE estimates
res=NULL

```

# ACE model - twins only


```{R, message=FALSE, warning=FALSE, eval=T}
# ----------------------------------------------------------------------------------------------------------------------
#  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)
colMeans(dzData[,c(selVars,covVars)],na.rm=T)
cov(mzData,use="pairwise.complete.obs")
cov(dzData,use="pairwise.complete.obs")

# Set Starting Values
sMu      <- 0                        # start value for means
sVa      <- .3                       # start value for A
sVc      <- .3                       # start value for C
sVe      <- .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

# Some different ways to extract output
fitACE$output$estimate
fitACE$VarC$result
fitACE$output$confidenceIntervals


```

## store results


```{R, message=FALSE, warning=FALSE, eval=T}

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


```

# ACE from twins and siblings

```{R, message=FALSE, warning=FALSE, eval=T}
# ----------------------------------------------------------------------------------------------------------------------
#  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")
cov(dzData[,selVars],use="pairwise.complete.obs")

# Set Starting Values
sMu      <- 0                        # start value for means
sVa      <- .3                       # start value for A
sVc      <- .3                       # start value for C
sVe      <- .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


```

## store results


```{R, message=FALSE, warning=FALSE, eval=T}

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


```


# Twin and siblings - alternate parametrisation

```{R, message=FALSE, warning=FALSE, eval=T}

# ----------------------------------------------------------------------------------------------------------------------
#  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")
cov(dzData[,selVars],use="pairwise.complete.obs")

# Set Starting Values
sMu      <- 0                        # start value for means
sVa      <- .3                       # start value for A
sVc      <- .3                       # start value for C
sVe      <- .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,.5,1,.5,1),  name="rAmz" ) 
relDZ      <- mxMatrix( type="Symm", nrow=nt, ncol=nt, free=FALSE, values=c(1,.5,.5,1,.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


```

## store results


```{R, message=FALSE, warning=FALSE, eval=T}

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


```

# Zygosity as definition variable

```{R, message=FALSE, warning=FALSE, eval=T}

# ----------------------------------------------------------------------------------------------------------------------
#  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")
cov(df[df$zygT==.5,selVars],use="pairwise.complete.obs")

# Set Starting Values
sMu      <- 0                        # start value for means
sVa      <- .3                       # start value for A
sVc      <- .3                       # start value for C
sVe      <- .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


```


```{R, message=FALSE, warning=FALSE, eval=T}

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


```

# ACE with relatedness 

```{R, message=FALSE, warning=FALSE, eval=T}

# ----------------------------------------------------------------------------------------------------------------------
#  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>.9,selVars],use="pairwise.complete.obs")
cov(df[df$s1<=.9,selVars],use="pairwise.complete.obs")

# Set Starting Values
sMu      <- 0                        # start value for means
sVa      <- .3                       # start value for A
sVc      <- .3                       # start value for C
sVe      <- .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

```

## store results

```{R, message=FALSE, warning=FALSE, eval=T}

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


```


# ACE - DZ only 

```{R, message=FALSE, warning=FALSE, eval=T}

# ----------------------------------------------------------------------------------------------------------------------
#  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")

# Set Starting Values
sMu      <- 0                        # start value for means
sVa      <- .3                       # start value for A
sVc      <- .3                       # start value for C
sVe      <- .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

```

## store results


```{R, message=FALSE, warning=FALSE, eval=T}

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


```

# AE on unrelated - GREML

```{R, message=FALSE, warning=FALSE, eval=T}
# 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)
#^^^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)])

# 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))
sum(is.na(phenfile$x))
#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.
#And check vice versa:
all(rownames(phenfile) %in% rownames(GRM)) #<--Should also be TRUE.
#In fact, the rows of phenfile are matched with those of the GRM:
all(rownames(phenfile)==rownames(GRM)) #<--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)
#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)
#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)
mxEval(h2,testrun,T)

#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)

# 95% confidence intervals
mxEval(h2,testrun,T)-1.96*h2se
mxEval(h2,testrun,T)+1.96*h2se

```

## store results


```{R, message=FALSE, warning=FALSE, eval=T}

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

```{R, message=FALSE, warning=FALSE, eval=T}

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


```{R, message=FALSE, warning=FALSE, eval=T}

colnames(res)=c("modelName","A_lb", "A", "A_ub", "C_lb", "C", "C_ub", "E_lb", "E", "E_ub")
print(res)

```
