isolation_mhealth_overlap_ACE_sexdiffs.Rmd

---
title: "The overlap between social isolation and mental health problems: ACE herogeniety model (sex differences)" 
output:  
  html_document:
    df_print: paged
    toc: yes
    toc_depth: 2
    toc_float:
      collapsed: false
    number_sections: false
    highlight: monochrome
    theme: flatly
    code_folding: show
    includes:
      after_body: footer.html
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE,
                      comment = NA,
                      prompt = FALSE,
                      cache = FALSE,
                      message = FALSE,
                      warning = FALSE,
                      results = 'markup')

options(bitmapType = 'quartz') # to render fonts better
```

```{r Clear global environment, include=FALSE}
remove(list = ls())
```

```{r Load packages, include=FALSE}
library(knitr)
library(haven)
library(psych)   
library(bestNormalize)
library(OpenMx)
library(tidyr)
library(tidyverse)
library(dplyr)     # conflicts with tidyverse for e.g. rename and row_number
```

# Source functions

```{r source functions, include=FALSE}
source("isolation_mhealth_functions.R")
```

# Read in data

```{r source the data file path, include=FALSE}
# source raw data directory
source("../isolation_mentalhealth_data_path.R")
```

```{r read in dta data file, include=FALSE}
dat.raw <- read_dta(paste0(data_path_raw, "Katie_23Sep22.dta"))
colnames(dat.raw)
```

### Column names

```{r select variables needed}
dat <- dat.raw %>%
  dplyr::select(
         atwinid,
         btwinid,
         familyid,
         rorderp5,
         torder,
         zygosity,
         sampsex,
         sisoe12,
         sisoy12,
         masce12,     # anxiety
         mascy12,   
         cdie12,      # depression
         cdiy12,
         conec12,     # antisocial behaviour / conduct disorder
         conyc12,
         psysympe12,  # psychosis - why is this not the tot scale variable - to give more variation?
         psysympy12,
         socisoe18,
         socisoy18,
         gadsxe18,    # anxiety
         gadsxy18,
         mdesxe18,    # depression
         mdesxy18,
         cdsxe18,     # antisocial behaviour / conduct disorder
         cdsxy18,
         psyexpe18,   # psychotic experiences
         psyexpy18
  )

colnames(dat)
```

### Recode variables into factors {.tabset .tabset-fade}

#### Sex

```{r recode sex}
dat <- dat %>%
  mutate(
    sex = 
      recode_factor(as_factor(sampsex),
        "1" = "Male",
        "2" = "Female"))

table(dat$sex)
```

#### Zygosity

```{r recode zygosity}
dat <- dat %>%
  mutate(
    zygosity = 
      recode_factor(as_factor(zygosity),
        "1" = "MZ",
        "2" = "DZ"))
table(dat$zygosity)
```

### Convert variables to numeric

```{r create numeric isolation variables}
dat <- dat %>%
  mutate(
    sisoe12 = as.numeric(sisoe12),
    sisoy12 = as.numeric(sisoy12),
    masce12 = as.numeric(masce12),       # anxiety
    mascy12 = as.numeric(mascy12),   
    cdie12 = as.numeric(cdie12),         # depression
    cdiy12 = as.numeric(cdiy12),
    conec12 = as.numeric(conec12),       # antisocial behaviour / conduct disorder
    conyc12 = as.numeric(conyc12),
    psysympe12 = as.numeric(psysympe12), # psychosis - why is this not the tot scale variable - to give more variation?
    psysympy12 = as.numeric(psysympy12),
    socisoe18 = as.numeric(socisoe18),
    socisoy18 = as.numeric(socisoy18),
    gadsxe18 = as.numeric(gadsxe18),     # anxiety
    gadsxy18 = as.numeric(gadsxy18),
    mdesxe18 = as.numeric(mdesxe18),     # depression
    mdesxy18 = as.numeric(mdesxy18),
    cdsxe18 = as.numeric(cdsxe18),       # antisocial behaviour / conduct disorder
    cdsxy18 = as.numeric(cdsxy18),
    psyexpe18 = as.numeric(psyexpe18),   # psychotic experiences
    psyexpy18 = as.numeric(psyexpy18)
  ) %>%
  select(                                # remove variables not needed
    -c(sampsex)
  )
```

# Variable lists

```{r select variables - raw}
# social isolation
selvars_si12 <- c("sisoe12", "sisoy12")
selvars_si18 <- c("socisoe18", "socisoy18")

# anxiety
selvars_anx12 <- c("masce12", "mascy12")
selvars_anx18 <- c("gadsxe18", "gadsxy18")

# depression
selvars_dep12 <- c("cdie12", "cdiy12")
selvars_dep18 <- c("mdesxe18", "mdesxy18")

# conduct
selvars_con12 <- c("conec12", "conyc12")
selvars_con18 <- c("cdsxe18", "cdsxy18")

# psychosis
selvars_psy12 <- c("psysympe12", "psysympy12")
selvars_psy18 <- c("psyexpe18", "psyexpy18")

# all 
selvars <- c("sisoe12", "masce12", "cdie12", "conec12",  "psysympe12", "socisoe18", "gadsxe18", "mdesxe18", "cdsxe18", "psyexpe18",
             "sisoy12", "mascy12", "cdiy12", "conyc12", "psysympy12",  "socisoy18", "gadsxy18", "mdesxy18", "cdsxy18","psyexpy18")

selvars_e <- c("sisoe12", "masce12", "cdie12", "conec12",  "psysympe12", "socisoe18", "gadsxe18", "mdesxe18", "cdsxe18", "psyexpe18")
```

```{r select variables - normalised}
# social isolation
selvars_si12norm <- c("sisoe12norm", "sisoy12norm")
selvars_si18norm <- c("socisoe18norm", "socisoy18norm")

# anxiety
selvars_anx12norm <- c("masce12norm", "mascy12norm")
selvars_anx18norm <- c("gadsxe18norm", "gadsxy18norm")

# depression
selvars_dep12norm <- c("cdie12norm", "cdiy12norm")
selvars_dep18norm <- c("mdesxe18norm", "mdesxy18norm")

# conduct
selvars_con12norm <- c("conec12norm", "conyc12norm")
selvars_con18norm <- c("cdsxe18norm", "cdsxy18norm")

# psychosis
selvars_psy12norm <- c("psysympe12norm", "psysympy12norm")
selvars_psy18norm <- c("psyexpe18norm", "psyexpy18norm")

# all 
selvars_norm <- c("sisoe12norm", "masce12norm", "cdie12norm", "conec12norm", "psysympe12norm", "socisoe18norm", "gadsxe18norm", "mdesxe18norm", "cdsxe18norm", "psyexpe18norm",
                  "sisoy12norm", "mascy12norm", "cdiy12norm", "conyc12norm", "psysympy12norm", "socisoy18norm", "gadsxy18norm", "mdesxy18norm", "cdsxy18norm", "psyexpy18norm")

# all with non-normalised - not in the twin modelling specific order here only used to standardise all variables at once
selvars_norm_all <- c("sisoe12", "sisoy12", "masce12", "mascy12", "cdie12", "cdiy12", "conec12", "conyc12", "psysympe12", "psysympy12", 
                      "socisoe18", "socisoy18", "gadsxe18", "gadsxy18", "mdesxe18", "mdesxy18", "cdsxe18", "cdsxy18", "psyexpe18", "psyexpy18",
                      "sisoe12norm", "sisoy12norm", "masce12norm", "mascy12norm", "cdie12norm", "cdiy12norm", "conec12norm", "conyc12norm", "psysympe12norm", "psysympy12norm", 
                      "socisoe18norm", "socisoy18norm", "gadsxe18norm", "gadsxy18norm", "mdesxe18norm", "mdesxy18norm", "cdsxe18norm", "cdsxy18norm", "psyexpe18norm", "psyexpy18norm")
```

```{r select variables - normalised and regressed sex}
# social isolation
selvars_si12norm_reg <- c("sisoe12norm_reg", "sisoy12norm_reg")
selvars_si18norm_reg <- c("socisoe18norm_reg", "socisoy18norm_reg")

# anxiety
selvars_anx12norm_reg <- c("masce12norm_reg", "mascy12norm_reg")
selvars_anx18norm_reg <- c("gadsxe18norm_reg", "gadsxy18norm_reg")

# depression
selvars_dep12norm_reg <- c("cdie12norm_reg", "cdiy12norm_reg")
selvars_dep18norm_reg <- c("mdesxe18norm_reg", "mdesxy18norm_reg")

# conduct
selvars_con12norm_reg <- c("conec12norm_reg", "conyc12norm_reg")
selvars_con18norm_reg <- c("cdsxe18norm_reg", "cdsxy18norm_reg")

# psychosis
selvars_psy12norm_reg <- c("psysympe12norm_reg", "psysympy12norm_reg")
selvars_psy18norm_reg <- c("psyexpe18norm_reg", "psyexpy18norm_reg")

# all 
selvars_norm_reg <- c("sisoe12norm_reg", "masce12norm_reg", "cdie12norm_reg", "conec12norm_reg", "psysympe12norm_reg", "socisoe18norm_reg", "gadsxe18norm_reg", "mdesxe18norm_reg", "cdsxe18norm_reg", "psyexpe18norm_reg",
                  "sisoy12norm_reg", "mascy12norm_reg", "cdiy12norm_reg", "conyc12norm_reg", "psysympy12norm_reg", "socisoy18norm_reg", "gadsxy18norm_reg", "mdesxy18norm_reg", "cdsxy18norm_reg", "psyexpy18norm_reg")

# all with non-norm_regalised - not in the twin modelling specific order here only used to standardise all variables at once
selvars_norm_reg_all <- c("sisoe12", "sisoy12", "masce12", "mascy12", "cdie12", "cdiy12", "conec12", "conyc12", "psysympe12", "psysympy12", 
                      "socisoe18", "socisoy18", "gadsxe18", "gadsxy18", "mdesxe18", "mdesxy18", "cdsxe18", "cdsxy18", "psyexpe18", "psyexpy18",
                      "sisoe12norm_reg", "sisoy12norm_reg", "masce12norm_reg", "mascy12norm_reg", "cdie12norm_reg", "cdiy12norm_reg", "conec12norm_reg", "conyc12norm_reg", "psysympe12norm_reg", "psysympy12norm_reg", 
                      "socisoe18norm_reg", "socisoy18norm_reg", "gadsxe18norm_reg", "gadsxy18norm_reg", "mdesxe18norm_reg", "mdesxy18norm_reg", "cdsxe18norm_reg", "cdsxy18norm_reg", "psyexpe18norm_reg", "psyexpy18norm_reg")
```

# Data prep

## Skewness

```{r histograms}
# isolation
hist(dat$sisoe12)    # not normal
hist(dat$socisoe18)  # not normal
# anxiety
hist(dat$masce12)    # normal
hist(dat$gadsxe18)   # not normal
# depression
hist(dat$cdie12)     # not normal
hist(dat$mdesxe18)   # not normal
# conduct
hist(dat$conec12)    # not normal
hist(dat$cdsxe18)    # not normal
# psychosis
hist(dat$psysympe12) # not normal
hist(dat$psyexpe18)  # not normal
```

## Rank transformation

Almost all variables are non-normal. We will use the van der Waerden's rank-based transformation as used in [Rimfeld et al 2021](https://acamh.onlinelibrary.wiley.com/doi/full/10.1002/jcv2.12053). For analyses using transformed data, they conducted the van der Waerden transformation prior to residualizing for age and sex as recommended by [Pain et al. 2018](https://www.nature.com/articles/s41431-018-0159-6).  

I will transform all variables to get the normalised estimate.

```{r rank transform variables elder variables}
# isolation age 12
sisoe12_n <- bestNormalize(dat$sisoe12)    # select the type of transformation needed
dat$sisoe12norm <- predict(sisoe12_n)      # create normalised variable
hist(dat$sisoe12norm)                  
summary(dat$sisoe12norm)

# isolation age 18
socisoe18_n <- bestNormalize(dat$socisoe18)    
dat$socisoe18norm <- predict(socisoe18_n)  
hist(dat$socisoe18norm)                  
summary(dat$socisoe18norm)

# anxiety age 12 - wont actually use this as it's already normally distributed
masce12_n <- bestNormalize(dat$masce12)    
dat$masce12norm <- predict(masce12_n)  
hist(dat$masce12norm)                  
summary(dat$masce12norm)

# anxiety age 18
gadsxe18_n <- bestNormalize(dat$gadsxe18)    
dat$gadsxe18norm <- predict(gadsxe18_n)  
hist(dat$gadsxe18norm)                  
summary(dat$gadsxe18norm)

# depression age 12
cdie12_n <- bestNormalize(dat$cdie12)    
dat$cdie12norm <- predict(cdie12_n)  
hist(dat$cdie12norm)                  
summary(dat$cdie12norm)

# depression age 18
mdesxe18_n <- bestNormalize(dat$mdesxe18)    
dat$mdesxe18norm <- predict(mdesxe18_n)  
hist(dat$mdesxe18norm)                  
summary(dat$mdesxe18norm)

# conduct age 12
conec12_n <- bestNormalize(dat$conec12)    
dat$conec12norm <- predict(conec12_n)  
hist(dat$conec12norm)                  
summary(dat$conec12norm)

# conduct age 18
cdsxe18_n <- bestNormalize(dat$cdsxe18)    
dat$cdsxe18norm <- predict(cdsxe18_n)  
hist(dat$cdsxe18norm)                  
summary(dat$cdsxe18norm)

# psychosis age 12
psysympe12_n <- bestNormalize(dat$psysympe12)    
dat$psysympe12norm <- predict(psysympe12_n)  
hist(dat$psysympe12norm)                  
summary(dat$psysympe12norm)

# psychosis age 18
psyexpe18_n <- bestNormalize(dat$psyexpe18)    
dat$psyexpe18norm <- predict(psyexpe18_n)  
hist(dat$psyexpe18norm)                  
summary(dat$psyexpe18norm)
```

```{r rank transform variables younger variables}
# isolation age 12
sisoy12_n <- bestNormalize(dat$sisoy12)    # select the type of transformation needed
dat$sisoy12norm <- predict(sisoy12_n)  # create normalised variable
hist(dat$sisoy12norm)                  
summary(dat$sisoy12norm)

# isolation age 18
socisoy18_n <- bestNormalize(dat$socisoy18)    
dat$socisoy18norm <- predict(socisoy18_n)  
hist(dat$socisoy18norm)                  
summary(dat$socisoy18norm)

# anxiety age 12 - wont actually use this as it's already normally distributed
mascy12_n <- bestNormalize(dat$mascy12)    
dat$mascy12norm <- predict(mascy12_n)  
hist(dat$mascy12norm)                  
summary(dat$mascy12norm)

# anxiety age 18
gadsxy18_n <- bestNormalize(dat$gadsxy18)    
dat$gadsxy18norm <- predict(gadsxy18_n)  
hist(dat$gadsxy18norm)                  
summary(dat$gadsxy18norm)

# depression age 12
cdiy12_n <- bestNormalize(dat$cdiy12)    
dat$cdiy12norm <- predict(cdiy12_n)  
hist(dat$cdiy12norm)                  
summary(dat$cdiy12norm)

# depression age 18
mdesxy18_n <- bestNormalize(dat$mdesxy18)    
dat$mdesxy18norm <- predict(mdesxy18_n)  
hist(dat$mdesxy18norm)                  
summary(dat$mdesxy18norm)

# conduct age 12
conyc12_n <- bestNormalize(dat$conyc12)    
dat$conyc12norm <- predict(conyc12_n)  
hist(dat$conyc12norm)                  
summary(dat$conyc12norm)

# conduct age 18
cdsxy18_n <- bestNormalize(dat$cdsxy18)    
dat$cdsxy18norm <- predict(cdsxy18_n)  
hist(dat$cdsxy18norm)                  
summary(dat$cdsxy18norm)

# psychosis age 12
psysympy12_n <- bestNormalize(dat$psysympy12)    
dat$psysympy12norm <- predict(psysympy12_n)  
hist(dat$psysympy12norm)                  
summary(dat$psysympy12norm)

# psychosis age 18
psyexpy18_n <- bestNormalize(dat$psyexpy18)    
dat$psyexpy18norm <- predict(psyexpy18_n)  
hist(dat$psyexpy18norm)                  
summary(dat$psyexpy18norm)
```

## Regress out sex

We first normalised the twin variables, then regress out sex. We don't regress out age here as all twins were measured as close to thir birthday as possible. 

```{r Regress out age and sex}
# twin 1 - elder
## isolation
dat$sisoe12norm_reg <- (resid(lm(data = dat, sisoe12norm ~ sex, na.action = na.exclude))) + 1
dat$socisoe18norm_reg <- (resid(lm(data = dat, socisoe18norm ~ sex, na.action = na.exclude))) + 1
## anxiety
dat$masce12norm_reg <- (resid(lm(data = dat, masce12norm ~ sex, na.action = na.exclude))) + 1
dat$gadsxe18norm_reg <- (resid(lm(data = dat, gadsxe18norm ~ sex, na.action = na.exclude))) + 1
## depression
dat$cdie12norm_reg <- (resid(lm(data = dat, cdie12norm ~ sex, na.action = na.exclude))) + 1
dat$mdesxe18norm_reg <- (resid(lm(data = dat, mdesxe18norm ~ sex, na.action = na.exclude))) + 1
## conduct
dat$conec12norm_reg <- (resid(lm(data = dat, conec12norm ~ sex, na.action = na.exclude))) + 1
dat$cdsxe18norm_reg <- (resid(lm(data = dat, cdsxe18norm ~ sex, na.action = na.exclude))) + 1
## psychosis
dat$psysympe12norm_reg <- (resid(lm(data = dat, psysympe12norm ~ sex, na.action = na.exclude))) + 1
dat$psyexpe18norm_reg <- (resid(lm(data = dat, psyexpe18norm ~ sex, na.action = na.exclude))) + 1

# twin 2 - younger
## isolation
dat$sisoy12norm_reg <- (resid(lm(data = dat, sisoy12norm ~ sex, na.action = na.exclude)))  + 1
dat$socisoy18norm_reg <- (resid(lm(data = dat, socisoy18norm ~ sex, na.action = na.exclude))) + 1
## anxiety
dat$mascy12norm_reg <- (resid(lm(data = dat, mascy12norm ~ sex, na.action = na.exclude))) + 1
dat$gadsxy18norm_reg <- (resid(lm(data = dat, gadsxy18norm ~ sex, na.action = na.exclude))) + 1
## depression
dat$cdiy12norm_reg <- (resid(lm(data = dat, cdiy12norm ~ sex, na.action = na.exclude))) + 1
dat$mdesxy18norm_reg <- (resid(lm(data = dat, mdesxy18norm ~ sex, na.action = na.exclude))) + 1
## conduct
dat$conyc12norm_reg <- (resid(lm(data = dat, conyc12norm ~ sex, na.action = na.exclude))) + 1
dat$cdsxy18norm_reg <- (resid(lm(data = dat, cdsxy18norm ~ sex, na.action = na.exclude))) + 1
## psychosis
dat$psysympy12norm_reg <- (resid(lm(data = dat, psysympy12norm ~ sex, na.action = na.exclude))) + 1
dat$psyexpy18norm_reg <- (resid(lm(data = dat, psyexpy18norm ~ sex, na.action = na.exclude))) + 1
```

## Create twin dataset

To remove the double entry in the data, we will remove everyone who has a "random twin order" variable of 0. This will then remove any birth order effects. 

```{r remove one twin pair row}
dat.twin <- dat %>% filter(rorderp5 == "1")
```

```{r datasets for MZ and DZ}
dat.twin.MZ <- dat.twin %>% filter(zygosity == "MZ")
dat.twin.DZ <- dat.twin %>% filter(zygosity == "DZ")

# male only
dat.twin.MZm <- dat.twin.MZ %>% filter(sex == "Male")
dat.twin.DZm <- dat.twin.DZ %>% filter(sex == "Male")
# female only
dat.twin.MZf <- dat.twin.MZ %>% filter(sex == "Female")
dat.twin.DZf <- dat.twin.DZ %>% filter(sex == "Female")
```

## Summary of MZ and DZ data

### Overall

```{r describe MZ and DZ data}
MZ_summary <- describe(dat.twin.MZ, 
                       skew = FALSE, 
                       range = FALSE)

DZ_summary <- describe(dat.twin.DZ, 
                       skew = FALSE, 
                       range = FALSE)
```

### MZ correlations

```{r MZ matrices}
# covariance matrix
covar.mz <- cov(dat.twin.MZ[, selvars], use = "complete")
# correlation matrix (standardized covariance)
cor.mz <- cor(dat.twin.MZ[, selvars], use = "complete")
# round(cor.mz, 3)
```

### DZ correlations

```{r DZ matrices}
# covariance matrix
covar.dz <- cov(dat.twin.DZ[, selvars], use = "complete")
# correlation matrix (standardized covariance)
cor.dz <- cor(dat.twin.DZ[, selvars], use = "complete")
# round(cor.dz, 3)
```

# Heterogeneity ACE model

This model tests whether genes and environment influence a trait to a different degree in males and females - i.e. do they show quantitative sex differences. 

Sex-limited expression of genetic or environmental factors occurs in two basic forms. First, the effects of a factor may be larger on one sex than on another, which is known as *scalar sex limitation*. Second, some factors may have an effect on one sex but not on the other, which is called nonscalar sex limitation. In the classical twin study, scalar sex- limited effects cause same-sex male and same-sex female twin correlations to differ (Neale, Røysamb, and Jacobson, 2006). 

## Specify

This model will be applied to all variables at age 12 and 18 - social isolation, anxiety, depression, conduct disorder, psychosis. 

```{r number of variables (phenotypes)}
# number of variables
nv <- 1    					      # number of variables 
ntv <- nv*2 				      # number of twin variables

# default optimizer  
mxOption(NULL, "Default optimizer", "CSOLNP")
```

```{r start values}
# start values
svM	  <- c(rep(0.7, ntv))	  # means               
svSD  <- c(rep(1.1, ntv))   # standard deviations
svRmz <- c(0.3)             # correlations for MZ
svRdz <- c(0.15)            # correlations for DZ 

rowVars	<- rep('Vars', nv)
colVars	<- rep(c('h2', 'c2', 'e2'), each = nv)
```

```{r specify quantitative heterogeneity ACE model}
# Matrix & Algebra for expected means vector
MeanM		<- mxMatrix(type = "Full", nrow = 1, ncol = nv, free = TRUE, values = 0, label = "meanm", name = "Mm")
MeanMM	<- mxAlgebra(expression = cbind(Mm, Mm), name = "expMeanm")
MeanF		<- mxMatrix(type = "Full", nrow = 1, ncol = nv, free = TRUE, values = 0, label = "meanf", name = "Mf" )
MeanFF	<- mxAlgebra(expression = cbind(Mf, Mf), name = "expMeanf")

# Define Matrices to estimate a, c, and e path coefficients for Males And Females
pathAm <- mxMatrix(type = "Full", nrow = nv, ncol = nv, free = TRUE, values = 0.6, label = "am11", name = "am")
pathCm <- mxMatrix(type = "Full", nrow = nv, ncol = nv, free = TRUE, values = 0.6, label = "cm11", name = "cm")
pathEm <- mxMatrix(type = "Full", nrow = nv, ncol = nv, free = TRUE, values = 0.6, label = "em11", name = "em")
pathAf <- mxMatrix(type = "Full", nrow = nv, ncol = nv, free = TRUE, values = 0.6, label = "af11", name = "af")
pathCf <- mxMatrix(type = "Full", nrow = nv, ncol = nv, free = TRUE, values = 0.6, label = "cf11", name = "cf")
pathEf <- mxMatrix(type = "Full", nrow = nv, ncol = nv, free = TRUE, values = 0.6, label = "ef11", name = "ef")

# Matrices A, C, and E compute variance components
covAm	<- mxAlgebra(expression = am %*% t(am), name = "Am")
covCm	<- mxAlgebra(expression = cm %*% t(cm), name = "Cm")
covEm	<- mxAlgebra(expression = em %*% t(em), name = "Em")
covAf	<- mxAlgebra(expression = af %*% t(af), name = "Af")
covCf	<- mxAlgebra(expression = cf %*% t(cf), name = "Cf")
covEf	<- mxAlgebra(expression = ef %*% t(ef), name = "Ef")

# Algebra to compute total variances and standard deviations (diagonal only)
Varm <- mxAlgebra(expression = Am + Cm + Em, name = "Vm")
Varf <- mxAlgebra(expression = Af + Cf + Ef, name = "Vf")

# Algebra to compute standardized variance components
h2m	<- mxAlgebra(expression = Am/Vm, name = "hm2")
c2m	<- mxAlgebra(expression = Cm/Vm, name = "cm2")
e2m	<- mxAlgebra(expression = Em/Vm, name = "em2")
h2f	<- mxAlgebra(expression = Af/Vf, name = "hf2")
c2f	<- mxAlgebra(expression = Cf/Vf, name = "cf2")
e2f	<- mxAlgebra(expression = Ef/Vf, name = "ef2")
    
# Algebra for expected variance/covariance matrix in the 4 groups
covMZM <- mxAlgebra(expression = rbind(cbind(Am + Cm + Em, Am+Cm), cbind(Am + Cm, Am + Cm + Em)), name = "expCovMZM")
covMZF <- mxAlgebra(expression = rbind(cbind(Af + Cf + Ef, Af+Cf), cbind(Af + Cf, Af + Cf + Ef)), name = "expCovMZF")
covDZM <- mxAlgebra(expression = rbind(cbind(Am + Cm + Em, 0.5%x%Am + Cm), cbind(0.5%x%Am + Cm, Am + Cm + Em)), name = "expCovDZM")
covDZF <- mxAlgebra(expression = rbind(cbind(Af + Cf + Ef, 0.5%x%Af + Cf), cbind(0.5%x%Af + Cf, Af + Cf + Ef)), name = "expCovDZF")
    
# Data objects for Multiple Groups
dataMZM	<- mxData(observed = dat.twin.MZm, type = "raw")
dataDZM	<- mxData(observed = dat.twin.DZm, type = "raw")
dataMZF	<- mxData(observed = dat.twin.MZf, type = "raw")
dataDZF	<- mxData(observed = dat.twin.DZf, type = "raw")

# fit function
fitFunction <- mxFitFunctionML()
obj	<-mxFitFunctionAlgebra("m2LL")

# confidence intervals
ciM <- mxCI(c('hm2', 'cm2', 'em2'))
ciF <- mxCI(c('hf2', 'cf2', 'ef2'))
```

## Social isolation 

### Age 12

```{r select which variables to use}
# selvars_si12_chosen <- selvars_si12                 # raw
# selvars_si12_chosen <- selvars_si12_z_score         # standardised
# selvars_si12_chosen <- selvars_si12norm             # normalised
# selvars_si12_chosen <- selvars_si12_z_score_norm    # standardised and normalised
  selvars_si12_chosen <- selvars_si12norm_reg         # normalised and regressed sex
```

```{r fit het model for social isolation 12}
# Objective objects for Multiple Groups     
objMZM <- mxExpectationNormal(covariance = "expCovMZM", means = "expMeanm", dimnames = selvars_si12_chosen)
objDZM <- mxExpectationNormal(covariance = "expCovDZM", means = "expMeanm", dimnames = selvars_si12_chosen)
objMZF <- mxExpectationNormal(covariance = "expCovMZF", means = "expMeanf", dimnames = selvars_si12_chosen)
objDZF <- mxExpectationNormal(covariance = "expCovDZF", means = "expMeanf", dimnames = selvars_si12_chosen)

# Combine Groups
parsm <- list(pathAm, pathCm, pathEm, covAm, covCm, covEm, Varm, h2m, c2m, e2m)
parsf <- list(pathAf, pathCf, pathEf, covAf, covCf, covEf, Varf, h2f, c2f, e2f)

modelMZM <- mxModel(parsm, MeanM, MeanMM, covMZM, dataMZM, objMZM, fitFunction, name = "MZM")
modelDZM <- mxModel(parsm, MeanM, MeanMM, covDZM, dataDZM, objDZM, fitFunction, name = "DZM")
modelMZF <- mxModel(parsf, MeanF, MeanFF, covMZF, dataMZF, objMZF, fitFunction, name = "MZF")
modelDZF <- mxModel(parsf, MeanF, MeanFF, covDZF, dataDZF, objDZF, fitFunction, name = "DZF")

minus2ll <-mxAlgebra(expression = MZM.objective + DZM.objective + MZF.objective + DZF.objective, name = "m2LL")

hetace_isolation12	<-mxModel("univHetACE", parsm, parsf, modelMZM, modelDZM, modelMZF, modelDZF, minus2ll, obj, ciM, ciF) 
```

```{r run het model for social isolation 12}
# Run
hetace_fit_isolation12 <- mxTryHard(hetace_isolation12, intervals = TRUE)
hetace_isolation12_summ <- summary(hetace_fit_isolation12)

# Table
hetace_isolation12_table <- ACEhet_esitmates_table(data = hetace_isolation12_summ,
                                                   variable = "Social isolation",
                                                   age = "12")
hetace_isolation12_table
```

```{r specify and run the sub-model homogeneity ACE model isolation 12}
homace_isolation12 <- mxModel(hetace_fit_isolation12, name = "homace")

homace_isolation12	<- omxSetParameters(homace_isolation12, 
                                labels = c("am11", "cm11", "em11", "af11", "cf11", "ef11"), free = TRUE, values = 0.6,
                                newlabels = c("a11", "c11", "e11", "a11", "c11", "e11"))

homace_fit_isolation12	<- mxTryHard(homace_isolation12, intervals = TRUE)
homace_fit_isolation12_summ	<- summary(homace_fit_isolation12)

# parameter estimates (check that hf2,cf2 and ef2 are the same)
homace_isolation12_table <- ACEhet_esitmates_table(data = homace_fit_isolation12_summ,
                                                   variable = "Social isolation",
                                                   age = "12")
homace_isolation12_table

# compare het and hom
mxCompare(hetace_fit_isolation12, homace_fit_isolation12)
```

### Age 18

```{r select which variables to use}
# selvars_si18_chosen <- selvars_si18                 # raw
# selvars_si18_chosen <- selvars_si18_z_score         # standardised
# selvars_si18_chosen <- selvars_si18norm             # normalised
# selvars_si18_chosen <- selvars_si18_z_score_norm    # standardised and normalised
  selvars_si18_chosen <- selvars_si18norm_reg         # normalised and regressed sex
```

```{r fit het model for social isolation 18}
# Objective objects for Multiple Groups     
objMZM <- mxExpectationNormal(covariance = "expCovMZM", means = "expMeanm", dimnames = selvars_si18_chosen)
objDZM <- mxExpectationNormal(covariance = "expCovDZM", means = "expMeanm", dimnames = selvars_si18_chosen)
objMZF <- mxExpectationNormal(covariance = "expCovMZF", means = "expMeanf", dimnames = selvars_si18_chosen)
objDZF <- mxExpectationNormal(covariance = "expCovDZF", means = "expMeanf", dimnames = selvars_si18_chosen)

# Combine Groups
parsm <- list(pathAm, pathCm, pathEm, covAm, covCm, covEm, Varm, h2m, c2m, e2m)
parsf <- list(pathAf, pathCf, pathEf, covAf, covCf, covEf, Varf, h2f, c2f, e2f)

modelMZM <- mxModel(parsm, MeanM, MeanMM, covMZM, dataMZM, objMZM, fitFunction, name = "MZM")
modelDZM <- mxModel(parsm, MeanM, MeanMM, covDZM, dataDZM, objDZM, fitFunction, name = "DZM")
modelMZF <- mxModel(parsf, MeanF, MeanFF, covMZF, dataMZF, objMZF, fitFunction, name = "MZF")
modelDZF <- mxModel(parsf, MeanF, MeanFF, covDZF, dataDZF, objDZF, fitFunction, name = "DZF")

minus2ll <-mxAlgebra(expression = MZM.objective + DZM.objective + MZF.objective + DZF.objective, name = "m2LL")

hetace_isolation18	<-mxModel("univHetACE", parsm, parsf, modelMZM, modelDZM, modelMZF, modelDZF, minus2ll, obj, ciM, ciF) 
```

```{r run het model for social isolation 18}
# Run
hetace_fit_isolation18 <- mxTryHard(hetace_isolation18, intervals = TRUE)
hetace_isolation18_summ <- summary(hetace_fit_isolation18)

# Table
hetace_isolation18_table <- ACEhet_esitmates_table(data = hetace_isolation18_summ,
                                                   variable = "Social isolation",
                                                   age = "18")
hetace_isolation18_table
```

```{r specify and run the sub-model homogeneity ACE model isolation 18}
homace_isolation18 <- mxModel(hetace_fit_isolation18, name = "homace")

homace_isolation18	<- omxSetParameters(homace_isolation18, 
                                labels = c("am11", "cm11", "em11", "af11", "cf11", "ef11"), free = TRUE, values = 0.6,
                                newlabels = c("a11", "c11", "e11", "a11", "c11", "e11"))

homace_fit_isolation18	<- mxTryHard(homace_isolation18, intervals = TRUE)
homace_fit_isolation18_summ	<- summary(homace_fit_isolation18)

# parameter estimates (check that hf2,cf2 and ef2 are the same)
homace_isolation18_table <- ACEhet_esitmates_table(data = homace_fit_isolation18_summ,
                                                   variable = "Social isolation",
                                                   age = "18")
homace_isolation18_table

# compare het and hom
mxCompare(hetace_fit_isolation18, homace_fit_isolation18)
```

## Anxiety

### Age 12

```{r select which variables to use}
# selvars_anx12_chosen <- selvars_anx12                 # raw
# selvars_anx12_chosen <- selvars_anx12_z_score         # standardised
# selvars_anx12_chosen <- selvars_anx12norm             # normalised
# selvars_anx12_chosen <- selvars_anx12_z_score_norm    # standardised and normalised
  selvars_anx12_chosen <- selvars_anx12norm_reg         # normalised and regressed sex
```

```{r fit het model for anxiety 12}
# Objective objects for Multiple Groups     
objMZM <- mxExpectationNormal(covariance = "expCovMZM", means = "expMeanm", dimnames = selvars_anx12_chosen)
objDZM <- mxExpectationNormal(covariance = "expCovDZM", means = "expMeanm", dimnames = selvars_anx12_chosen)
objMZF <- mxExpectationNormal(covariance = "expCovMZF", means = "expMeanf", dimnames = selvars_anx12_chosen)
objDZF <- mxExpectationNormal(covariance = "expCovDZF", means = "expMeanf", dimnames = selvars_anx12_chosen)

# Combine Groups
parsm <- list(pathAm, pathCm, pathEm, covAm, covCm, covEm, Varm, h2m, c2m, e2m)
parsf <- list(pathAf, pathCf, pathEf, covAf, covCf, covEf, Varf, h2f, c2f, e2f)

modelMZM <- mxModel(parsm, MeanM, MeanMM, covMZM, dataMZM, objMZM, fitFunction, name = "MZM")
modelDZM <- mxModel(parsm, MeanM, MeanMM, covDZM, dataDZM, objDZM, fitFunction, name = "DZM")
modelMZF <- mxModel(parsf, MeanF, MeanFF, covMZF, dataMZF, objMZF, fitFunction, name = "MZF")
modelDZF <- mxModel(parsf, MeanF, MeanFF, covDZF, dataDZF, objDZF, fitFunction, name = "DZF")

minus2ll <-mxAlgebra(expression = MZM.objective + DZM.objective + MZF.objective + DZF.objective, name = "m2LL")

hetace_anxiety12	<-mxModel("univHetACE", parsm, parsf, modelMZM, modelDZM, modelMZF, modelDZF, minus2ll, obj, ciM, ciF) 
```

```{r run het model for anxiety 12}
# Run
hetace_fit_anxiety12 <- mxTryHard(hetace_anxiety12, intervals = TRUE)
hetace_anxiety12_summ <- summary(hetace_fit_anxiety12)

# Table
hetace_anxiety12_table <- ACEhet_esitmates_table(data = hetace_anxiety12_summ,
                                                   variable = "Anxiety",
                                                   age = "12")
hetace_anxiety12_table
```

```{r specify and run the sub-model homogeneity ACE model anxiety 12}
homace_anxiety12 <- mxModel(hetace_fit_anxiety12, name = "homace")

homace_anxiety12	<- omxSetParameters(homace_anxiety12, 
                                labels = c("am11", "cm11", "em11", "af11", "cf11", "ef11"), free = TRUE, values = 0.6,
                                newlabels = c("a11", "c11", "e11", "a11", "c11", "e11"))

homace_fit_anxiety12	<- mxTryHard(homace_anxiety12, intervals = TRUE)
homace_fit_anxiety12_summ	<- summary(homace_fit_anxiety12)

# parameter estimates (check that hf2,cf2 and ef2 are the same)
homace_anxiety12_table <- ACEhet_esitmates_table(data = homace_fit_anxiety12_summ,
                                                   variable = "Anxiety",
                                                   age = "12")
homace_anxiety12_table

# compare het and hom
mxCompare(hetace_fit_anxiety12, homace_fit_anxiety12)
```

### Age 18

```{r select which variables to use}
# selvars_anx18_chosen <- selvars_anx18                 # raw
# selvars_anx18_chosen <- selvars_anx18_z_score         # standardised
# selvars_anx18_chosen <- selvars_anx18norm             # normalised
# selvars_anx18_chosen <- selvars_anx18_z_score_norm    # standardised and normalised
  selvars_anx18_chosen <- selvars_anx18norm_reg         # normalised and regressed sex
```

```{r fit het model for anxiety 18}
# Objective objects for Multiple Groups     
objMZM <- mxExpectationNormal(covariance = "expCovMZM", means = "expMeanm", dimnames = selvars_anx18_chosen)
objDZM <- mxExpectationNormal(covariance = "expCovDZM", means = "expMeanm", dimnames = selvars_anx18_chosen)
objMZF <- mxExpectationNormal(covariance = "expCovMZF", means = "expMeanf", dimnames = selvars_anx18_chosen)
objDZF <- mxExpectationNormal(covariance = "expCovDZF", means = "expMeanf", dimnames = selvars_anx18_chosen)

# Combine Groups
parsm <- list(pathAm, pathCm, pathEm, covAm, covCm, covEm, Varm, h2m, c2m, e2m)
parsf <- list(pathAf, pathCf, pathEf, covAf, covCf, covEf, Varf, h2f, c2f, e2f)

modelMZM <- mxModel(parsm, MeanM, MeanMM, covMZM, dataMZM, objMZM, fitFunction, name = "MZM")
modelDZM <- mxModel(parsm, MeanM, MeanMM, covDZM, dataDZM, objDZM, fitFunction, name = "DZM")
modelMZF <- mxModel(parsf, MeanF, MeanFF, covMZF, dataMZF, objMZF, fitFunction, name = "MZF")
modelDZF <- mxModel(parsf, MeanF, MeanFF, covDZF, dataDZF, objDZF, fitFunction, name = "DZF")

minus2ll <-mxAlgebra(expression = MZM.objective + DZM.objective + MZF.objective + DZF.objective, name = "m2LL")

hetace_anxiety18	<-mxModel("univHetACE", parsm, parsf, modelMZM, modelDZM, modelMZF, modelDZF, minus2ll, obj, ciM, ciF) 
```

```{r run het model for anxiety 18}
# Run
hetace_fit_anxiety18 <- mxTryHard(hetace_anxiety18, intervals = TRUE)
hetace_anxiety18_summ <- summary(hetace_fit_anxiety18)

# Table
hetace_anxiety18_table <- ACEhet_esitmates_table(data = hetace_anxiety18_summ,
                                                   variable = "Anxiety",
                                                   age = "18")
hetace_anxiety18_table
```

```{r specify and run the sub-model homogeneity ACE model anxiety 18}
homace_anxiety18 <- mxModel(hetace_fit_anxiety18, name = "homace")

homace_anxiety18	<- omxSetParameters(homace_anxiety18, 
                                labels = c("am11", "cm11", "em11", "af11", "cf11", "ef11"), free = TRUE, values = 0.6,
                                newlabels = c("a11", "c11", "e11", "a11", "c11", "e11"))

homace_fit_anxiety18	<- mxTryHard(homace_anxiety18, intervals = TRUE)
homace_fit_anxiety18_summ	<- summary(homace_fit_anxiety18)

# parameter estimates (check that hf2,cf2 and ef2 are the same)
homace_anxiety18_table <- ACEhet_esitmates_table(data = homace_fit_anxiety18_summ,
                                                   variable = "Anxiety",
                                                   age = "18")
homace_anxiety18_table

# compare het and hom
mxCompare(hetace_fit_anxiety18, homace_fit_anxiety18)
```

## Depression

### Age 12

```{r select which variables to use}
# selvars_dep12_chosen <- selvars_dep12                 # raw
# selvars_dep12_chosen <- selvars_dep12_z_score         # standardised
# selvars_dep12_chosen <- selvars_dep12norm             # normalised
# selvars_dep12_chosen <- selvars_dep12_z_score_norm    # standardised and normalised
  selvars_dep12_chosen <- selvars_dep12norm_reg         # normalised and regressed sex
```

```{r fit het model for depression 12}
# Objective objects for Multiple Groups     
objMZM <- mxExpectationNormal(covariance = "expCovMZM", means = "expMeanm", dimnames = selvars_dep12_chosen)
objDZM <- mxExpectationNormal(covariance = "expCovDZM", means = "expMeanm", dimnames = selvars_dep12_chosen)
objMZF <- mxExpectationNormal(covariance = "expCovMZF", means = "expMeanf", dimnames = selvars_dep12_chosen)
objDZF <- mxExpectationNormal(covariance = "expCovDZF", means = "expMeanf", dimnames = selvars_dep12_chosen)

# Combine Groups
parsm <- list(pathAm, pathCm, pathEm, covAm, covCm, covEm, Varm, h2m, c2m, e2m)
parsf <- list(pathAf, pathCf, pathEf, covAf, covCf, covEf, Varf, h2f, c2f, e2f)

modelMZM <- mxModel(parsm, MeanM, MeanMM, covMZM, dataMZM, objMZM, fitFunction, name = "MZM")
modelDZM <- mxModel(parsm, MeanM, MeanMM, covDZM, dataDZM, objDZM, fitFunction, name = "DZM")
modelMZF <- mxModel(parsf, MeanF, MeanFF, covMZF, dataMZF, objMZF, fitFunction, name = "MZF")
modelDZF <- mxModel(parsf, MeanF, MeanFF, covDZF, dataDZF, objDZF, fitFunction, name = "DZF")

minus2ll <-mxAlgebra(expression = MZM.objective + DZM.objective + MZF.objective + DZF.objective, name = "m2LL")

hetace_depression12	<-mxModel("univHetACE", parsm, parsf, modelMZM, modelDZM, modelMZF, modelDZF, minus2ll, obj, ciM, ciF) 
```

```{r run het model for depression 12}
# Run
hetace_fit_depression12 <- mxTryHard(hetace_depression12, intervals = TRUE)
hetace_depression12_summ <- summary(hetace_fit_depression12)

# Table
hetace_depression12_table <- ACEhet_esitmates_table(data = hetace_depression12_summ,
                                                   variable = "Depression",
                                                   age = "12")
hetace_depression12_table
```

```{r specify and run the sub-model homogeneity ACE model depression 12}
homace_depression12 <- mxModel(hetace_fit_depression12, name = "homace")

homace_depression12	<- omxSetParameters(homace_depression12, 
                                labels = c("am11", "cm11", "em11", "af11", "cf11", "ef11"), free = TRUE, values = 0.6,
                                newlabels = c("a11", "c11", "e11", "a11", "c11", "e11"))

homace_fit_depression12	<- mxTryHard(homace_depression12, intervals = TRUE)
homace_fit_depression12_summ	<- summary(homace_fit_depression12)

# parameter estimates (check that hf2,cf2 and ef2 are the same)
homace_depression12_table <- ACEhet_esitmates_table(data = homace_fit_depression12_summ,
                                                   variable = "Depression",
                                                   age = "12")
homace_depression12_table

# compare het and hom
mxCompare(hetace_fit_depression12, homace_fit_depression12)
```

### Age 18

```{r select which variables to use}
# selvars_dep18_chosen <- selvars_dep18                 # raw
# selvars_dep18_chosen <- selvars_dep18_z_score         # standardised
# selvars_dep18_chosen <- selvars_dep18norm             # normalised
# selvars_dep18_chosen <- selvars_dep18_z_score_norm    # standardised and normalised
  selvars_dep18_chosen <- selvars_dep18norm_reg         # normalised and regressed sex
```

```{r fit het model for depression 18}
# Objective objects for Multiple Groups     
objMZM <- mxExpectationNormal(covariance = "expCovMZM", means = "expMeanm", dimnames = selvars_dep18_chosen)
objDZM <- mxExpectationNormal(covariance = "expCovDZM", means = "expMeanm", dimnames = selvars_dep18_chosen)
objMZF <- mxExpectationNormal(covariance = "expCovMZF", means = "expMeanf", dimnames = selvars_dep18_chosen)
objDZF <- mxExpectationNormal(covariance = "expCovDZF", means = "expMeanf", dimnames = selvars_dep18_chosen)

# Combine Groups
parsm <- list(pathAm, pathCm, pathEm, covAm, covCm, covEm, Varm, h2m, c2m, e2m)
parsf <- list(pathAf, pathCf, pathEf, covAf, covCf, covEf, Varf, h2f, c2f, e2f)

modelMZM <- mxModel(parsm, MeanM, MeanMM, covMZM, dataMZM, objMZM, fitFunction, name = "MZM")
modelDZM <- mxModel(parsm, MeanM, MeanMM, covDZM, dataDZM, objDZM, fitFunction, name = "DZM")
modelMZF <- mxModel(parsf, MeanF, MeanFF, covMZF, dataMZF, objMZF, fitFunction, name = "MZF")
modelDZF <- mxModel(parsf, MeanF, MeanFF, covDZF, dataDZF, objDZF, fitFunction, name = "DZF")

minus2ll <-mxAlgebra(expression = MZM.objective + DZM.objective + MZF.objective + DZF.objective, name = "m2LL")

hetace_depression18	<-mxModel("univHetACE", parsm, parsf, modelMZM, modelDZM, modelMZF, modelDZF, minus2ll, obj, ciM, ciF) 
```

```{r run het model for depression 18}
# Run
hetace_fit_depression18 <- mxTryHard(hetace_depression18, intervals = TRUE)
hetace_depression18_summ <- summary(hetace_fit_depression18)

# Table
hetace_depression18_table <- ACEhet_esitmates_table(data = hetace_depression18_summ,
                                                   variable = "Depression",
                                                   age = "18")
hetace_depression18_table
```

```{r specify and run the sub-model homogeneity ACE model depression 18}
homace_depression18 <- mxModel(hetace_fit_depression18, name = "homace")

homace_depression18	<- omxSetParameters(homace_depression18, 
                                labels = c("am11", "cm11", "em11", "af11", "cf11", "ef11"), free = TRUE, values = 0.6,
                                newlabels = c("a11", "c11", "e11", "a11", "c11", "e11"))

homace_fit_depression18	<- mxTryHard(homace_depression18, intervals = TRUE)
homace_fit_depression18_summ	<- summary(homace_fit_depression18)

# parameter estimates (check that hf2,cf2 and ef2 are the same)
homace_depression18_table <- ACEhet_esitmates_table(data = homace_fit_depression18_summ,
                                                   variable = "Depression",
                                                   age = "18")
homace_depression18_table

# compare het and hom
mxCompare(hetace_fit_depression18, homace_fit_depression18)
```

## Conduct disorder

### Age 12

```{r select which variables to use}
# selvars_con12_chosen <- selvars_con12                 # raw
# selvars_con12_chosen <- selvars_con12_z_score         # standardised
# selvars_con12_chosen <- selvars_con12norm             # normalised
# selvars_con12_chosen <- selvars_con12_z_score_norm    # standardised and normalised
  selvars_con12_chosen <- selvars_con12norm_reg         # normalised and regressed sex
```

```{r fit het model for conduct 12}
# Objective objects for Multiple Groups     
objMZM <- mxExpectationNormal(covariance = "expCovMZM", means = "expMeanm", dimnames = selvars_con12_chosen)
objDZM <- mxExpectationNormal(covariance = "expCovDZM", means = "expMeanm", dimnames = selvars_con12_chosen)
objMZF <- mxExpectationNormal(covariance = "expCovMZF", means = "expMeanf", dimnames = selvars_con12_chosen)
objDZF <- mxExpectationNormal(covariance = "expCovDZF", means = "expMeanf", dimnames = selvars_con12_chosen)

# Combine Groups
parsm <- list(pathAm, pathCm, pathEm, covAm, covCm, covEm, Varm, h2m, c2m, e2m)
parsf <- list(pathAf, pathCf, pathEf, covAf, covCf, covEf, Varf, h2f, c2f, e2f)

modelMZM <- mxModel(parsm, MeanM, MeanMM, covMZM, dataMZM, objMZM, fitFunction, name = "MZM")
modelDZM <- mxModel(parsm, MeanM, MeanMM, covDZM, dataDZM, objDZM, fitFunction, name = "DZM")
modelMZF <- mxModel(parsf, MeanF, MeanFF, covMZF, dataMZF, objMZF, fitFunction, name = "MZF")
modelDZF <- mxModel(parsf, MeanF, MeanFF, covDZF, dataDZF, objDZF, fitFunction, name = "DZF")

minus2ll <-mxAlgebra(expression = MZM.objective + DZM.objective + MZF.objective + DZF.objective, name = "m2LL")

hetace_conduct12	<-mxModel("univHetACE", parsm, parsf, modelMZM, modelDZM, modelMZF, modelDZF, minus2ll, obj, ciM, ciF) 
```

```{r run het model for conduct 12}
# Run
hetace_fit_conduct12 <- mxTryHard(hetace_conduct12, intervals = TRUE)
hetace_conduct12_summ <- summary(hetace_fit_conduct12)

# Table
hetace_conduct12_table <- ACEhet_esitmates_table(data = hetace_conduct12_summ,
                                                   variable = "Conduct disorder",
                                                   age = "12")
hetace_conduct12_table
```

```{r specify and run the sub-model homogeneity ACE model conduct 12}
homace_conduct12 <- mxModel(hetace_fit_conduct12, name = "homace")

homace_conduct12	<- omxSetParameters(homace_conduct12, 
                                labels = c("am11", "cm11", "em11", "af11", "cf11", "ef11"), free = TRUE, values = 0.6,
                                newlabels = c("a11", "c11", "e11", "a11", "c11", "e11"))

homace_fit_conduct12	<- mxTryHard(homace_conduct12, intervals = TRUE)
homace_fit_conduct12_summ	<- summary(homace_fit_conduct12)

# parameter estimates (check that hf2,cf2 and ef2 are the same)
homace_conduct12_table <- ACEhet_esitmates_table(data = homace_fit_conduct12_summ,
                                                   variable = "Conduct disorder",
                                                   age = "12")
homace_conduct12_table

# compare het and hom
mxCompare(hetace_fit_conduct12, homace_fit_conduct12)
```

### Age 18

```{r select which variables to use}
# selvars_con18_chosen <- selvars_con18                 # raw
# selvars_con18_chosen <- selvars_con18_z_score         # standardised
# selvars_con18_chosen <- selvars_con18norm             # normalised
# selvars_con18_chosen <- selvars_con18_z_score_norm    # standardised and normalised
  selvars_con18_chosen <- selvars_con18norm_reg         # normalised and regressed sex
```

```{r fit het model for conduct 18}
# Objective objects for Multiple Groups     
objMZM <- mxExpectationNormal(covariance = "expCovMZM", means = "expMeanm", dimnames = selvars_con18_chosen)
objDZM <- mxExpectationNormal(covariance = "expCovDZM", means = "expMeanm", dimnames = selvars_con18_chosen)
objMZF <- mxExpectationNormal(covariance = "expCovMZF", means = "expMeanf", dimnames = selvars_con18_chosen)
objDZF <- mxExpectationNormal(covariance = "expCovDZF", means = "expMeanf", dimnames = selvars_con18_chosen)

# Combine Groups
parsm <- list(pathAm, pathCm, pathEm, covAm, covCm, covEm, Varm, h2m, c2m, e2m)
parsf <- list(pathAf, pathCf, pathEf, covAf, covCf, covEf, Varf, h2f, c2f, e2f)

modelMZM <- mxModel(parsm, MeanM, MeanMM, covMZM, dataMZM, objMZM, fitFunction, name = "MZM")
modelDZM <- mxModel(parsm, MeanM, MeanMM, covDZM, dataDZM, objDZM, fitFunction, name = "DZM")
modelMZF <- mxModel(parsf, MeanF, MeanFF, covMZF, dataMZF, objMZF, fitFunction, name = "MZF")
modelDZF <- mxModel(parsf, MeanF, MeanFF, covDZF, dataDZF, objDZF, fitFunction, name = "DZF")

minus2ll <-mxAlgebra(expression = MZM.objective + DZM.objective + MZF.objective + DZF.objective, name = "m2LL")

hetace_conduct18	<-mxModel("univHetACE", parsm, parsf, modelMZM, modelDZM, modelMZF, modelDZF, minus2ll, obj, ciM, ciF) 
```

```{r run het model for conduct 18}
# Run
hetace_fit_conduct18 <- mxTryHard(hetace_conduct18, intervals = TRUE)
hetace_conduct18_summ <- summary(hetace_fit_conduct18)

# Table
hetace_conduct18_table <- ACEhet_esitmates_table(data = hetace_conduct18_summ,
                                                   variable = "Conduct disorder",
                                                   age = "18")
hetace_conduct18_table
```

```{r specify and run the sub-model homogeneity ACE model conduct 18}
homace_conduct18 <- mxModel(hetace_fit_conduct18, name = "homace")

homace_conduct18	<- omxSetParameters(homace_conduct18, 
                                labels = c("am11", "cm11", "em11", "af11", "cf11", "ef11"), free = TRUE, values = 0.6,
                                newlabels = c("a11", "c11", "e11", "a11", "c11", "e11"))

homace_fit_conduct18	<- mxTryHard(homace_conduct18, intervals = TRUE)
homace_fit_conduct18_summ	<- summary(homace_fit_conduct18)

# parameter estimates (check that hf2,cf2 and ef2 are the same)
homace_conduct18_table <- ACEhet_esitmates_table(data = homace_fit_conduct18_summ,
                                                   variable = "Conduct disorder",
                                                   age = "18")
homace_conduct18_table

# compare het and hom
mxCompare(hetace_fit_conduct18, homace_fit_conduct18)
```

## Psychotic experiences

### Age 12

```{r select which variables to use}
# selvars_psy12_chosen <- selvars_psy12                 # raw
# selvars_psy12_chosen <- selvars_psy12_z_score         # standardised
# selvars_psy12_chosen <- selvars_psy12norm             # normalised
# selvars_psy12_chosen <- selvars_psy12_z_score_norm    # standardised and normalised
  selvars_psy12_chosen <- selvars_psy12norm_reg         # normalised and regressed sex
```

```{r fit het model for psychosis 12}
# Objective objects for Multiple Groups     
objMZM <- mxExpectationNormal(covariance = "expCovMZM", means = "expMeanm", dimnames = selvars_psy12_chosen)
objDZM <- mxExpectationNormal(covariance = "expCovDZM", means = "expMeanm", dimnames = selvars_psy12_chosen)
objMZF <- mxExpectationNormal(covariance = "expCovMZF", means = "expMeanf", dimnames = selvars_psy12_chosen)
objDZF <- mxExpectationNormal(covariance = "expCovDZF", means = "expMeanf", dimnames = selvars_psy12_chosen)

# Combine Groups
parsm <- list(pathAm, pathCm, pathEm, covAm, covCm, covEm, Varm, h2m, c2m, e2m)
parsf <- list(pathAf, pathCf, pathEf, covAf, covCf, covEf, Varf, h2f, c2f, e2f)

modelMZM <- mxModel(parsm, MeanM, MeanMM, covMZM, dataMZM, objMZM, fitFunction, name = "MZM")
modelDZM <- mxModel(parsm, MeanM, MeanMM, covDZM, dataDZM, objDZM, fitFunction, name = "DZM")
modelMZF <- mxModel(parsf, MeanF, MeanFF, covMZF, dataMZF, objMZF, fitFunction, name = "MZF")
modelDZF <- mxModel(parsf, MeanF, MeanFF, covDZF, dataDZF, objDZF, fitFunction, name = "DZF")

minus2ll <-mxAlgebra(expression = MZM.objective + DZM.objective + MZF.objective + DZF.objective, name = "m2LL")

hetace_psychosis12	<-mxModel("univHetACE", parsm, parsf, modelMZM, modelDZM, modelMZF, modelDZF, minus2ll, obj, ciM, ciF) 
```

```{r run het model for psychosis 12}
# Run
hetace_fit_psychosis12 <- mxTryHard(hetace_psychosis12, intervals = TRUE)
hetace_psychosis12_summ <- summary(hetace_fit_psychosis12)

# Table
hetace_psychosis12_table <- ACEhet_esitmates_table(data = hetace_psychosis12_summ,
                                                   variable = "Psychosis",
                                                   age = "12")
hetace_psychosis12_table
```

```{r specify and run the sub-model homogeneity ACE model psychosis 12}
homace_psychosis12 <- mxModel(hetace_fit_psychosis12, name = "homace")

homace_psychosis12	<- omxSetParameters(homace_psychosis12, 
                                labels = c("am11", "cm11", "em11", "af11", "cf11", "ef11"), free = TRUE, values = 0.6,
                                newlabels = c("a11", "c11", "e11", "a11", "c11", "e11"))

homace_fit_psychosis12	<- mxTryHard(homace_psychosis12, intervals = TRUE)
homace_fit_psychosis12_summ	<- summary(homace_fit_psychosis12)

# parameter estimates (check that hf2,cf2 and ef2 are the same)
homace_psychosis12_table <- ACEhet_esitmates_table(data = homace_fit_psychosis12_summ,
                                                   variable = "Psychosis",
                                                   age = "12")
homace_psychosis12_table

# compare het and hom
mxCompare(hetace_fit_psychosis12, homace_fit_psychosis12)
```

### Age 18

```{r select which variables to use}
# selvars_psy18_chosen <- selvars_psy18                 # raw
# selvars_psy18_chosen <- selvars_psy18_z_score         # standardised
# selvars_psy18_chosen <- selvars_psy18norm             # normalised
# selvars_psy18_chosen <- selvars_psy18_z_score_norm    # standardised and normalised
  selvars_psy18_chosen <- selvars_psy18norm_reg         # normalised and regressed sex
```

```{r fit het model for psychosis 18}
# Objective objects for Multiple Groups     
objMZM <- mxExpectationNormal(covariance = "expCovMZM", means = "expMeanm", dimnames = selvars_psy18_chosen)
objDZM <- mxExpectationNormal(covariance = "expCovDZM", means = "expMeanm", dimnames = selvars_psy18_chosen)
objMZF <- mxExpectationNormal(covariance = "expCovMZF", means = "expMeanf", dimnames = selvars_psy18_chosen)
objDZF <- mxExpectationNormal(covariance = "expCovDZF", means = "expMeanf", dimnames = selvars_psy18_chosen)

# Combine Groups
parsm <- list(pathAm, pathCm, pathEm, covAm, covCm, covEm, Varm, h2m, c2m, e2m)
parsf <- list(pathAf, pathCf, pathEf, covAf, covCf, covEf, Varf, h2f, c2f, e2f)

modelMZM <- mxModel(parsm, MeanM, MeanMM, covMZM, dataMZM, objMZM, fitFunction, name = "MZM")
modelDZM <- mxModel(parsm, MeanM, MeanMM, covDZM, dataDZM, objDZM, fitFunction, name = "DZM")
modelMZF <- mxModel(parsf, MeanF, MeanFF, covMZF, dataMZF, objMZF, fitFunction, name = "MZF")
modelDZF <- mxModel(parsf, MeanF, MeanFF, covDZF, dataDZF, objDZF, fitFunction, name = "DZF")

minus2ll <-mxAlgebra(expression = MZM.objective + DZM.objective + MZF.objective + DZF.objective, name = "m2LL")

hetace_psychosis18	<-mxModel("univHetACE", parsm, parsf, modelMZM, modelDZM, modelMZF, modelDZF, minus2ll, obj, ciM, ciF) 
```

```{r run het model for psychosis 18}
# Run
hetace_fit_psychosis18 <- mxTryHard(hetace_psychosis18, intervals = TRUE)
hetace_psychosis18_summ <- summary(hetace_fit_psychosis18)

# Table
hetace_psychosis18_table <- ACEhet_esitmates_table(data = hetace_psychosis18_summ,
                                                   variable = "Psychosis",
                                                   age = "18")
hetace_psychosis18_table
```

```{r specify and run the sub-model homogeneity ACE model psychosis 18}
homace_psychosis18 <- mxModel(hetace_fit_psychosis18, name = "homace")

homace_psychosis18	<- omxSetParameters(homace_psychosis18, 
                                labels = c("am11", "cm11", "em11", "af11", "cf11", "ef11"), free = TRUE, values = 0.6,
                                newlabels = c("a11", "c11", "e11", "a11", "c11", "e11"))

homace_fit_psychosis18	<- mxTryHard(homace_psychosis18, intervals = TRUE)
homace_fit_psychosis18_summ	<- summary(homace_fit_psychosis18)

# parameter estimates (check that hf2,cf2 and ef2 are the same)
homace_psychosis18_table <- ACEhet_esitmates_table(data = homace_fit_psychosis18_summ,
                                                   variable = "Psychosis",
                                                   age = "18")
homace_psychosis18_table

# compare het and hom
mxCompare(hetace_fit_psychosis18, homace_fit_psychosis18)
```


***

# Scalar ACE model 

The 'non-scalar' sex limitation refers to qualitative sex differences - the actual components having an effect differ between sexes - e.g., different A and C influences. Whereas the first one - the scalar sex limitation can be thought about in 2 ways: A 'true ' quantitative sex difference - where it is the same components in males and females but having different amounts (quantities/magnitudes of influences). In the 'false' quantitative sex difference - the total variance is larger in females compared to males or vice versa, thus, the 'raw' amount of A, C and E components (i.e., unstandardised) will be naturally larger in females - even if the standardised variances (i.e., relative proportions) are comparable across sex. This 'false' quantitative sex difference is what we test in the scalar model.

This will model one set of ACE parameters for males and females, but specifying a scalar to mutiply e.g. the female variances. 

Note: This only makes sense if you know from the univariate results that both variables show a variance difference between males and females. If only one shows variance differences, then more complicated hybrid models are needed. 

## Specify

```{r check out the male and female datasets}
# overall
psych::describe(dat.twin[selvars])
# males
psych::describe(dat.twin.MZm[selvars])
psych::describe(dat.twin.DZm[selvars])
# females
psych::describe(dat.twin.MZf[selvars])
psych::describe(dat.twin.DZf[selvars])
```

```{r number of variables}
nv			<- 1			        	# number of variables for a twin = 1 in Univariate
ntv			<- 2*nv		        	# number of variables for a pair = 2* 1 for Univariate
nlower	<- ntv*(ntv+1)/2 		# number of free elements in a lower matrix ntv*ntv
ncor		<- (nv*(nv+1)/2)-nv	# number of free elements in a correlation matrix nv*nv
```

```{r labels for means, sds, and scalar}
# labels for the means
LabMM	  <- paste("mm", 1:nv, sep = "")		# Male singleton
LabMF	  <- paste("mf", 1:nv, sep = "")		# Female singleton
LabMMM	<- c(LabMM, LabMM)			          # male-male pairs
LabMFF	<- c(LabMF, LabMF)			          # female-female pairs

# labels for the A, C and E standard deviation matrices (paths from latent factors to variables) no sex dif
LabAm	  <- paste("apm", 1:nv, sep = "")		
LabCm	  <- paste("cpm", 1:nv, sep = "")		
LabEm	  <- paste("epm", 1:nv, sep = "")		

# labels for the scalar matrix
LabS		<- paste("s", 1:nv, sep = "")		
LabSS  	<- c(LabS, LabS)
LabNA	  <- c("NA")
# Labdos	<- c(LabNA, LabS) # dont think we need this one
PatF		<- c(F)
PatT		<- c(T)
Pat	  	<- c(PatF, PatT)
```

```{r start values}
Stpath	  <- c(.6) 			# change here according to the variances of the variables
Stcor		  <- c(.6)
Stmean	  <- 4.6
StScalar	<-c(.8)				# change here to >1 if female variance is larger than males or if it is mixed
DumOnes	  <-c(1)
StSdos	  <- c(DumOnes,StScalar)		# Scalar for Males is just 1 - do we need this??
```

```{r specify}
# Matrices to estimate Means for Male/Female pairs
MeanMM	<- mxMatrix("Full", 1, ntv, free = TRUE, values = c(Stmean, Stmean), labels = LabMMM, name = "expMeanMM")
MeanM		<- mxMatrix("Full", 1, nv,  free = TRUE, values = c(Stmean), 		     labels = LabMM, 	name = "expMeanM")
MeanFF	<- mxMatrix("Full", 1, ntv, free = TRUE, values = c(Stmean, Stmean), labels = LabMFF,	name = "expMeanFF")
MeanF		<- mxMatrix("Full", 1, nv,  free = TRUE, values = c(Stmean), 		     labels = LabMF, 	name = "expMeanF")
# MeanMF	<- mxMatrix("Full", 1, ntv, free = TRUE, values = c(Stmean,Stmean), 	labels=LabMMF, 	name="expMeanMF" )

# ACE paths
pathAm	<- mxMatrix(type = "Diag", nrow = nv, ncol = nv, free = TRUE, values = Stpath, label = LabAm, name = "am")
pathCm	<- mxMatrix(type = "Diag", nrow = nv, ncol = nv, free = TRUE, values = Stpath, label = LabCm, name = "cm")
pathEm	<- mxMatrix(type = "Diag", nrow = nv, ncol = nv, free = TRUE, values = Stpath, label = LabEm, name = "em")

# Algebra to compute the variance components
covAm		<- mxAlgebra(expression = am %*% t(am), name = "Am")
covCm		<- mxAlgebra(expression = cm %*% t(cm), name = "Cm")
covEm		<- mxAlgebra(expression = em %*% t(em), name = "Em")

# Algebra to compute standardized variance components
covm		<- mxAlgebra(expression = Am + Cm + Em, name = "Vm")
StAm		<- mxAlgebra(expression = Am/Vm, name = "h2m")
StCm		<- mxAlgebra(expression = Cm/Vm, name = "c2m")
StEm		<- mxAlgebra(expression = Em/Vm, name = "e2m")

# scalar
ScalarF	<- mxMatrix(type = "Diag", nrow = 2, ncol = 2, free = TRUE, values = c(StScalar, StScalar), label = LabSS, name = "ScF")
# ScalarOS	<- mxMatrix( type="Diag", nrow=2, ncol=2, free=Pat, values=StSdos, label=Labdos, name="ScOS")

# Algebra for expected variance/covariance matrix in the 4 groups
covMZM	<- mxAlgebra(expression = rbind(cbind(Am + Cm + Em, Am+Cm),	cbind(Am + Cm, Am + Cm + Em)), name = "ExpCovMZM")
covDZM	<- mxAlgebra(expression = rbind(cbind(Am + Cm + Em, 0.5%x%Am + Cm),	cbind(0.5%x%Am + Cm, Am + Cm + Em)), name = "ExpCovDZM")

covMZF	<- mxAlgebra(expression = ScF%&%(rbind(cbind(Am + Cm + Em, Am + Cm), cbind(Am + Cm, Am + Cm + Em))), name = "ExpCovMZF")
covDZF	<- mxAlgebra(expression = ScF%&%(rbind(cbind(Am + Cm + Em, 0.5%x%Am + Cm), cbind(0.5%x%Am + Cm, Am + Cm + Em))), name = "ExpCovDZF")

# covDOS	<-mxAlgebra( expression= ScOS %&% (rbind  ( cbind(Am+Cm+Em, 0.5%x%Am+Cm),	cbind(0.5%x%Am+Cm, Am+Cm+Em)	)),	name="ExpCovDOS")

# Data objects for Multiple Groups
dataMZM	<- mxData(dat.twin.MZm, type = "raw")
dataMZF	<- mxData(dat.twin.MZf, type = "raw")
dataDZM	<- mxData(dat.twin.DZm, type = "raw")
dataDZF	<- mxData(dat.twin.DZf, type = "raw")

fitFunction <- mxFitFunctionML()
```

## Univariate ACE model 

```{r specity ACE}
# Matrix & Algebra for expected means vector (free)
MeanMZ <- mxMatrix("Full", 1, ntv, free = T, values = c(90,90), labels = c("MZm1", "MZm2"), name = "expMeanMZ")
MeanDZ <- mxMatrix("Full", 1, ntv, free = T, values = c(90,90), labels = c("DZm1", "DZm2"), name = "expMeanDZ")

# Matrices to estimate free a, c, and e path coefficients
pathA	<- mxMatrix(type = "Full", nrow = nv, ncol = nv, free = T, values = 10, label = "a11", name = "a")
pathC	<- mxMatrix(type = "Full", nrow = nv, ncol = nv, free = T, values = 10, label = "c11", name = "c")
pathE	<- mxMatrix(type = "Full", nrow = nv, ncol = nv, free = T, values = 10, label = "e11", name = "e")

# Algebra to generate Matrices to hold A, C, and E computed Variance Components
covA <- mxAlgebra(expression = a %*% t(a), name = "A" )
covC <- mxAlgebra(expression = c %*% t(c), name = "C" )
covE <- mxAlgebra(expression = e %*% t(e), name = "E" )

# Algebra to compute Total Variance
covP <- mxAlgebra(expression = A+C+E, name = "V" )
sd   <- mxAlgebra(expression = sqrt(V), name = "SD")

# Algebra to compute standardized path estimates 
stpathA	<- mxAlgebra(expression = a/SD, name = "sta")
stpathC <- mxAlgebra(expression = c/SD, name = "stc")
stpathE <- mxAlgebra(expression = e/SD, name = "ste")

# Algebra to compute standardized variance components 
stcovA <- mxAlgebra(expression = A/V, name = "h2")
stcovC <- mxAlgebra(expression = C/V, name = "c2")
stcovE <- mxAlgebra(expression = E/V, name = "e2")

# Algebra for Expected Variance/Covariance Matrices in MZ & DZ twins
covMZ <- mxAlgebra(expression = rbind(cbind(A+C+E,A+C), cbind(A+C,A+C+E)), name = "expCovMZ")
covDZ <- mxAlgebra(expression = rbind(cbind(A+C+E,0.5%x%A+C),cbind(0.5%x%A+C,A+C+E)),	name = "expCovDZ")

# Data objects for Multiple Groups
dataMZ <-mxData(observed = dat.twin.MZ, type = "raw")
dataDZ <-mxData(observed = dat.twin.DZ, type = "raw")

fitFunction <- mxFitFunctionML()
```

## Social isolation

### Age 12

```{r select which variables to use}
# selvars_si12_chosen <- selvars_si12                 # raw
# selvars_si12_chosen <- selvars_si12_z_score         # standardised
# selvars_si12_chosen <- selvars_si12norm             # normalised
# selvars_si12_chosen <- selvars_si12_z_score_norm    # standardised and normalised
  selvars_si12_chosen <- selvars_si12norm_reg         # normalised and regressed sex
```

```{r social isolation scalar}
# objectives   
objMZM	  <- mxExpectationNormal(covariance = "ExpCovMZM", means = "expMeanMM", dimnames = selvars_si12_chosen)
objDZM	  <- mxExpectationNormal(covariance = "ExpCovDZM", means = "expMeanMM", dimnames = selvars_si12_chosen)
objMZF	  <- mxExpectationNormal(covariance = "ExpCovMZF", means = "expMeanFF", dimnames = selvars_si12_chosen)
objDZF	  <- mxExpectationNormal(covariance = "ExpCovDZF", means = "expMeanFF", dimnames = selvars_si12_chosen)

# Combine Groups
parsm	  	<- list(pathAm, pathCm, pathEm, covAm, covCm, covEm, covm, StAm, StCm, StEm, fitFunction)
modelMZM 	<- mxModel(parsm, MeanM, MeanMM, covMZM, dataMZM, objMZM,  name = "MZM")
modelDZM 	<- mxModel(parsm, MeanM, MeanMM, covDZM, dataDZM, objDZM,  name = "DZM")
modelMZF 	<- mxModel(parsm, MeanF, MeanFF, covMZF, dataMZF, objMZF, ScalarF, name = "MZF")
modelDZF 	<- mxModel(parsm, MeanF, MeanFF, covDZF, dataDZF, objDZF, ScalarF, name = "DZF")

minus2ll  <- mxAlgebra(expression =	MZM.objective + DZM.objective + MZF.objective + DZF.objective, name = "m2LL")
obj	    	<- mxFitFunctionAlgebra("m2LL")
ci	    	<-mxCI(c('MZM.h2m[1,1]','MZM.c2m[1,1]','MZM.e2m[1,1]') )

ScACE2Model	<-mxModel("ScACE2", parsm, modelMZM, modelDZM, modelMZF, modelDZF, minus2ll, obj, ci) 

# Run
ScACE_isolation12_fit	<- mxTryHard(ScACE2Model, intervals = TRUE)
ScACE_isolation12_sum <- summary(ScACE_isolation12_fit)
```

```{r output social isolation 12 scalar}
# ACE estimates
scalar_isolation12_table <- ACE_esitmates_table(
  data = ScACE_isolation12_sum, 
  variable = "Social isolation",
  age = "12")

scalar_isolation12_table

# scalar estimate
mxEval(MZF.ScF, ScACE_isolation12_fit)
```

```{r run age 12 isolation ACE}
# Objective objects for Multiple Groups
objMZ <- mxExpectationNormal(covariance = "expCovMZ", means = "expMeanMZ", dimnames = selvars_si12_chosen)
objDZ <- mxExpectationNormal(covariance = "expCovDZ", means = "expMeanDZ", dimnames = selvars_si12_chosen)

# Combine Groups
pars <- list(pathA, pathC, pathE, covA, covC, covE, covP, sd, stpathA, stpathC, stpathE, stcovA, stcovC, stcovE, fitFunction)

# MZ and DZ models
modelMZ	<- mxModel(pars, MeanMZ, covMZ, dataMZ, objMZ, name = "MZ")
modelDZ	<- mxModel(pars, MeanDZ, covDZ, dataDZ, objDZ, name = "DZ")

# Fitting
minus2ll	<- mxAlgebra(expression = MZ.objective + DZ.objective, name = "m2LL")
obj	    	<- mxFitFunctionAlgebra("m2LL")
Conf	  	<- mxCI(c('ACE.h2', 'ACE.c2', 'ACE.e2'))

# Model
ACEModel <- mxModel("ACE", pars, modelMZ, modelDZ, minus2ll, obj, Conf) 

# Run
ACEFit_isolation12 <- mxTryHard(ACEModel, intervals = TRUE)
isolation12_estimates <- summary(ACEFit_isolation12)

# Table
isolation12_table <- ACE_esitmates_table(
  data = isolation12_estimates, 
  variable = "Social isolation",
  age = "12")

isolation12_table
```

Now to compare the univariate ACE model and the scalar model:

```{r compare scalar and uni isolation12}
# compare scalar and univariate model using the AIC
mxCompare(ACEFit_isolation12, ScACE_isolation12_fit)
```

Here, the AIC is smaller for the scalar model, therefore the Scalar model has a better fit to the data. 

Compare to heterogeneity model:

```{r compare scalar and uni isolation12}
# compare scalar and univariate model using the AIC
mxCompare(hetace_fit_isolation12, ScACE_isolation12_fit)
```

**The scalar model is still the best fit for isolation age 12** - compared to the heterogeneity model. 

### Age 18

```{r select which variables to use}
# selvars_si18_chosen <- selvars_si18                 # raw
# selvars_si18_chosen <- selvars_si18_z_score         # standardised
# selvars_si18_chosen <- selvars_si18norm             # normalised
# selvars_si18_chosen <- selvars_si18_z_score_norm    # standardised and normalised
  selvars_si18_chosen <- selvars_si18norm_reg         # normalised and regressed sex
```

```{r social isolation scalar}
# objectives   
objMZM	  <- mxExpectationNormal(covariance = "ExpCovMZM", means = "expMeanMM", dimnames = selvars_si18_chosen)
objDZM	  <- mxExpectationNormal(covariance = "ExpCovDZM", means = "expMeanMM", dimnames = selvars_si18_chosen)
objMZF	  <- mxExpectationNormal(covariance = "ExpCovMZF", means = "expMeanFF", dimnames = selvars_si18_chosen)
objDZF	  <- mxExpectationNormal(covariance = "ExpCovDZF", means = "expMeanFF", dimnames = selvars_si18_chosen)

# Combine Groups
parsm	  	<- list(pathAm, pathCm, pathEm, covAm, covCm, covEm, covm, StAm, StCm, StEm, fitFunction)
modelMZM 	<- mxModel(parsm, MeanM, MeanMM, covMZM, dataMZM, objMZM,  name = "MZM")
modelDZM 	<- mxModel(parsm, MeanM, MeanMM, covDZM, dataDZM, objDZM,  name = "DZM")
modelMZF 	<- mxModel(parsm, MeanF, MeanFF, covMZF, dataMZF, objMZF, ScalarF, name = "MZF")
modelDZF 	<- mxModel(parsm, MeanF, MeanFF, covDZF, dataDZF, objDZF, ScalarF, name = "DZF")

minus2ll  <- mxAlgebra(expression =	MZM.objective + DZM.objective + MZF.objective + DZF.objective, name = "m2LL")
obj	    	<- mxFitFunctionAlgebra("m2LL")
ci	    	<-mxCI(c('MZM.h2m[1,1]','MZM.c2m[1,1]','MZM.e2m[1,1]') )

ScACE2Model	<-mxModel("ScACE2", parsm, modelMZM, modelDZM, modelMZF, modelDZF, minus2ll, obj, ci) 

# Run
ScACE_isolation18_fit	<- mxTryHard(ScACE2Model, intervals = TRUE)
ScACE_isolation18_sum <- summary(ScACE_isolation18_fit)
```

```{r output social isolation 18 scalar}
# ACE estimates
scalar_isolation18_table <- ACE_esitmates_table(
  data = ScACE_isolation18_sum, 
  variable = "Social isolation",
  age = "18")

scalar_isolation18_table

# scalar estimate? 
mxEval(MZF.ScF, ScACE_isolation18_fit)
```

```{r run age 18 isolation ACE}
# Objective objects for Multiple Groups
objMZ <- mxExpectationNormal(covariance = "expCovMZ", means = "expMeanMZ", dimnames = selvars_si18_chosen)
objDZ <- mxExpectationNormal(covariance = "expCovDZ", means = "expMeanDZ", dimnames = selvars_si18_chosen)

# Combine Groups
pars <- list(pathA, pathC, pathE, covA, covC, covE, covP, sd, stpathA, stpathC, stpathE, stcovA, stcovC, stcovE, fitFunction)

# MZ and DZ models
modelMZ	<- mxModel(pars, MeanMZ, covMZ, dataMZ, objMZ, name = "MZ")
modelDZ	<- mxModel(pars, MeanDZ, covDZ, dataDZ, objDZ, name = "DZ")

# Fitting
minus2ll	<- mxAlgebra(expression = MZ.objective + DZ.objective, name = "m2LL")
obj	    	<- mxFitFunctionAlgebra("m2LL")
Conf	  	<- mxCI(c('ACE.h2', 'ACE.c2', 'ACE.e2'))

# Model
ACEModel <- mxModel("ACE", pars, modelMZ, modelDZ, minus2ll, obj, Conf) 

# Run
ACEFit_isolation18 <- mxTryHard(ACEModel, intervals = TRUE)
isolation18_estimates <- summary(ACEFit_isolation18)

# Table
isolation18_table <- ACE_esitmates_table(
  data = isolation18_estimates, 
  variable = "Social isolation",
  age = "18")

isolation18_table
```

Now to compare the univariate ACE model and the scalar model:

```{r compare scalar and uni isolation18}
# compare scalar and univariate model using the AIC
mxCompare(ACEFit_isolation18, ScACE_isolation18_fit)
```

Here, the AIC is smaller for the scalar model, therefore the Scalar model has the best fit for the data. 

Compare to heterogeneity model:

```{r compare scalar and uni isolation18}
# compare scalar and univariate model using the AIC
mxCompare(hetace_fit_isolation18, ScACE_isolation18_fit)
```

**The scalar model is still the best fit for isolation age 18** - compared to the heterogeneity model. 

## Anxiety

### Age 12

```{r select which variables to use}
# selvars_anx12_chosen <- selvars_anx12                 # raw
# selvars_anx12_chosen <- selvars_anx12_z_score         # standardised
# selvars_anx12_chosen <- selvars_anx12norm             # normalised
# selvars_anx12_chosen <- selvars_anx12_z_score_norm    # standardised and normalised
  selvars_anx12_chosen <- selvars_anx12norm_reg         # normalised and regressed sex
```

```{r social isolation scalar}
# objectives   
objMZM	  <- mxExpectationNormal(covariance = "ExpCovMZM", means = "expMeanMM", dimnames = selvars_anx12_chosen)
objDZM	  <- mxExpectationNormal(covariance = "ExpCovDZM", means = "expMeanMM", dimnames = selvars_anx12_chosen)
objMZF	  <- mxExpectationNormal(covariance = "ExpCovMZF", means = "expMeanFF", dimnames = selvars_anx12_chosen)
objDZF	  <- mxExpectationNormal(covariance = "ExpCovDZF", means = "expMeanFF", dimnames = selvars_anx12_chosen)

# Combine Groups
parsm	  	<- list(pathAm, pathCm, pathEm, covAm, covCm, covEm, covm, StAm, StCm, StEm, fitFunction)
modelMZM 	<- mxModel(parsm, MeanM, MeanMM, covMZM, dataMZM, objMZM,  name = "MZM")
modelDZM 	<- mxModel(parsm, MeanM, MeanMM, covDZM, dataDZM, objDZM,  name = "DZM")
modelMZF 	<- mxModel(parsm, MeanF, MeanFF, covMZF, dataMZF, objMZF, ScalarF, name = "MZF")
modelDZF 	<- mxModel(parsm, MeanF, MeanFF, covDZF, dataDZF, objDZF, ScalarF, name = "DZF")

minus2ll  <- mxAlgebra(expression =	MZM.objective + DZM.objective + MZF.objective + DZF.objective, name = "m2LL")
obj	    	<- mxFitFunctionAlgebra("m2LL")
ci	    	<-mxCI(c('MZM.h2m[1,1]','MZM.c2m[1,1]','MZM.e2m[1,1]') )

ScACE2Model	<-mxModel("ScACE2", parsm, modelMZM, modelDZM, modelMZF, modelDZF, minus2ll, obj, ci) 

# Run
ScACE_anxiety12_fit	<- mxTryHard(ScACE2Model, intervals = TRUE)
ScACE_anxiety12_sum <- summary(ScACE_anxiety12_fit)
```

```{r output social isolation 18 scalar}
# ACE estimates
ACE_esitmates_table(
  data = ScACE_anxiety12_sum, 
  variable = "Anxiety",
  age = "12")

# scalar estimate? 
mxEval(MZF.ScF, ScACE_anxiety12_fit)
```

```{r run age 12 anxiety ACE}
# Objective objects for Multiple Groups
objMZ <- mxExpectationNormal(covariance = "expCovMZ", means = "expMeanMZ", dimnames = selvars_anx12_chosen)
objDZ <- mxExpectationNormal(covariance = "expCovDZ", means = "expMeanDZ", dimnames = selvars_anx12_chosen)

# Combine Groups
pars <- list(pathA, pathC, pathE, covA, covC, covE, covP, sd, stpathA, stpathC, stpathE, stcovA, stcovC, stcovE, fitFunction)

# MZ and DZ models
modelMZ	<- mxModel(pars, MeanMZ, covMZ, dataMZ, objMZ, name = "MZ")
modelDZ	<- mxModel(pars, MeanDZ, covDZ, dataDZ, objDZ, name = "DZ")

# Fitting
minus2ll	<- mxAlgebra(expression = MZ.objective + DZ.objective, name = "m2LL")
obj	    	<- mxFitFunctionAlgebra("m2LL")
Conf	  	<- mxCI(c('ACE.h2', 'ACE.c2', 'ACE.e2'))

# Model
ACEModel <- mxModel("ACE", pars, modelMZ, modelDZ, minus2ll, obj, Conf) 

# Run
ACEFit_anxiety12 <- mxTryHard(ACEModel, intervals = TRUE)
anxiety12_estimates <- summary(ACEFit_anxiety12)

# Table
anxiety12_table <- ACE_esitmates_table(
  data = anxiety12_estimates, 
  variable = "Anxiety",
  age = "12")

anxiety12_table
```

Now to compare the univariate ACE model and the scalar model:

```{r compare scalar and uni anxiety12}
# compare scalar and univariate model using the AIC
mxCompare(ACEFit_anxiety12, ScACE_anxiety12_fit)
```

Here, the scalar model has the same AIC as the ACE model. Therefore, we would not take this as a better fit. 

Compare to heterogeneity model:

```{r compare scalar and uni anxiety12}
# compare scalar and univariate model using the AIC
mxCompare(hetace_fit_anxiety12, ScACE_anxiety12_fit)
```

**The homogeneity ACE model is still the best fit for anxiety age 12** - compared to the heterogeneity model. 

### Age 18

```{r select which variables to use}
# selvars_anx18_chosen <- selvars_anx18                 # raw
# selvars_anx18_chosen <- selvars_anx18_z_score         # standardised
# selvars_anx18_chosen <- selvars_anx18norm             # normalised
# selvars_anx18_chosen <- selvars_anx18_z_score_norm    # standardised and normalised
  selvars_anx18_chosen <- selvars_anx18norm_reg         # normalised and regressed sex
```

```{r social isolation scalar}
# objectives   
objMZM	  <- mxExpectationNormal(covariance = "ExpCovMZM", means = "expMeanMM", dimnames = selvars_anx18_chosen)
objDZM	  <- mxExpectationNormal(covariance = "ExpCovDZM", means = "expMeanMM", dimnames = selvars_anx18_chosen)
objMZF	  <- mxExpectationNormal(covariance = "ExpCovMZF", means = "expMeanFF", dimnames = selvars_anx18_chosen)
objDZF	  <- mxExpectationNormal(covariance = "ExpCovDZF", means = "expMeanFF", dimnames = selvars_anx18_chosen)

# Combine Groups
parsm	  	<- list(pathAm, pathCm, pathEm, covAm, covCm, covEm, covm, StAm, StCm, StEm, fitFunction)
modelMZM 	<- mxModel(parsm, MeanM, MeanMM, covMZM, dataMZM, objMZM,  name = "MZM")
modelDZM 	<- mxModel(parsm, MeanM, MeanMM, covDZM, dataDZM, objDZM,  name = "DZM")
modelMZF 	<- mxModel(parsm, MeanF, MeanFF, covMZF, dataMZF, objMZF, ScalarF, name = "MZF")
modelDZF 	<- mxModel(parsm, MeanF, MeanFF, covDZF, dataDZF, objDZF, ScalarF, name = "DZF")

minus2ll  <- mxAlgebra(expression =	MZM.objective + DZM.objective + MZF.objective + DZF.objective, name = "m2LL")
obj	    	<- mxFitFunctionAlgebra("m2LL")
ci	    	<- mxCI(c('MZM.h2m[1,1]','MZM.c2m[1,1]','MZM.e2m[1,1]') )

ScACE2Model	<-mxModel("ScACE2", parsm, modelMZM, modelDZM, modelMZF, modelDZF, minus2ll, obj, ci) 

# Run
ScACE_anxiety18_fit	<- mxTryHard(ScACE2Model, intervals = TRUE)
ScACE_anxiety18_sum <- summary(ScACE_anxiety18_fit)
```

```{r output social anxiety 18 scalar}
# ACE estimates
ACE_esitmates_table(
  data = ScACE_anxiety18_sum, 
  variable = "Anxiety",
  age = "18")

# scalar estimate? 
mxEval(MZF.ScF, ScACE_anxiety18_fit)
```

```{r run age 18 anxiety ACE}
# Objective objects for Multiple Groups
objMZ <- mxExpectationNormal(covariance = "expCovMZ", means = "expMeanMZ", dimnames = selvars_anx18_chosen)
objDZ <- mxExpectationNormal(covariance = "expCovDZ", means = "expMeanDZ", dimnames = selvars_anx18_chosen)

# Combine Groups
pars <- list(pathA, pathC, pathE, covA, covC, covE, covP, sd, stpathA, stpathC, stpathE, stcovA, stcovC, stcovE, fitFunction)

# MZ and DZ models
modelMZ	<- mxModel(pars, MeanMZ, covMZ, dataMZ, objMZ, name = "MZ")
modelDZ	<- mxModel(pars, MeanDZ, covDZ, dataDZ, objDZ, name = "DZ")

# Fitting
minus2ll	<- mxAlgebra(expression = MZ.objective + DZ.objective, name = "m2LL")
obj	    	<- mxFitFunctionAlgebra("m2LL")
Conf	  	<- mxCI(c('ACE.h2', 'ACE.c2', 'ACE.e2'))

# Model
ACEModel <- mxModel("ACE", pars, modelMZ, modelDZ, minus2ll, obj, Conf) 

# Run
ACEFit_anxiety18 <- mxTryHard(ACEModel, intervals = TRUE)
anxiety18_estimates <- summary(ACEFit_anxiety18)

# Table
anxiety18_table <- ACE_esitmates_table(
  data = anxiety18_estimates, 
  variable = "Anxiety",
  age = "18")

anxiety18_table
```

Now to compare the scalar model and the univariate mdoel

```{r compare scalar and uni anxiety18}
# compare scalar and univariate model using the AIC
mxCompare(ACEFit_anxiety18, ScACE_anxiety18_fit)
```

Here, the AIC is lower for the scalar model, therefore the Scalar model is the best fit to the data. 

Compare to heterogeneity model:

```{r compare scalar and uni anxiety18}
# compare scalar and univariate model using the AIC
mxCompare(hetace_fit_anxiety18, ScACE_anxiety18_fit)
```

**The scalar model is still the best fit for anxiety age 18** - compared to the heterogeneity model. 

## Depression

### Age 12

```{r select which variables to use}
# selvars_dep12_chosen <- selvars_dep12                 # raw
# selvars_dep12_chosen <- selvars_dep12_z_score         # standardised
# selvars_dep12_chosen <- selvars_dep12norm             # normalised
# selvars_dep12_chosen <- selvars_dep12_z_score_norm    # standardised and normalised
  selvars_dep12_chosen <- selvars_dep12norm_reg         # normalised and regressed sex
```

```{r social isolation scalar}
# objectives   
objMZM	  <- mxExpectationNormal(covariance = "ExpCovMZM", means = "expMeanMM", dimnames = selvars_dep12_chosen)
objDZM	  <- mxExpectationNormal(covariance = "ExpCovDZM", means = "expMeanMM", dimnames = selvars_dep12_chosen)
objMZF	  <- mxExpectationNormal(covariance = "ExpCovMZF", means = "expMeanFF", dimnames = selvars_dep12_chosen)
objDZF	  <- mxExpectationNormal(covariance = "ExpCovDZF", means = "expMeanFF", dimnames = selvars_dep12_chosen)

# Combine Groups
parsm	  	<- list(pathAm, pathCm, pathEm, covAm, covCm, covEm, covm, StAm, StCm, StEm, fitFunction)
modelMZM 	<- mxModel(parsm, MeanM, MeanMM, covMZM, dataMZM, objMZM,  name = "MZM")
modelDZM 	<- mxModel(parsm, MeanM, MeanMM, covDZM, dataDZM, objDZM,  name = "DZM")
modelMZF 	<- mxModel(parsm, MeanF, MeanFF, covMZF, dataMZF, objMZF, ScalarF, name = "MZF")
modelDZF 	<- mxModel(parsm, MeanF, MeanFF, covDZF, dataDZF, objDZF, ScalarF, name = "DZF")

minus2ll  <- mxAlgebra(expression =	MZM.objective + DZM.objective + MZF.objective + DZF.objective, name = "m2LL")
obj	    	<- mxFitFunctionAlgebra("m2LL")
ci	    	<- mxCI(c('MZM.h2m[1,1]','MZM.c2m[1,1]','MZM.e2m[1,1]') )

ScACE2Model	<-mxModel("ScACE2", parsm, modelMZM, modelDZM, modelMZF, modelDZF, minus2ll, obj, ci) 

# Run
ScACE_depression12_fit	<- mxTryHard(ScACE2Model, intervals = TRUE)
ScACE_depression12_sum <- summary(ScACE_depression12_fit)
```

```{r output social depression 12 scalar}
# ACE estimates
scalar_depression12_table <- ACE_esitmates_table(
  data = ScACE_depression12_sum, 
  variable = "Depression",
  age = "12")

scalar_depression12_table

# scalar estimate 
mxEval(MZF.ScF, ScACE_depression12_fit)
```

```{r run age 12 depression ACE}
# Objective objects for Multiple Groups
objMZ <- mxExpectationNormal(covariance = "expCovMZ", means = "expMeanMZ", dimnames = selvars_dep12_chosen)
objDZ <- mxExpectationNormal(covariance = "expCovDZ", means = "expMeanDZ", dimnames = selvars_dep12_chosen)

# Combine Groups
pars <- list(pathA, pathC, pathE, covA, covC, covE, covP, sd, stpathA, stpathC, stpathE, stcovA, stcovC, stcovE, fitFunction)

# MZ and DZ models
modelMZ	<- mxModel(pars, MeanMZ, covMZ, dataMZ, objMZ, name = "MZ")
modelDZ	<- mxModel(pars, MeanDZ, covDZ, dataDZ, objDZ, name = "DZ")

# Fitting
minus2ll	<- mxAlgebra(expression = MZ.objective + DZ.objective, name = "m2LL")
obj	    	<- mxFitFunctionAlgebra("m2LL")
Conf	  	<- mxCI(c('ACE.h2', 'ACE.c2', 'ACE.e2'))

# Model
ACEModel <- mxModel("ACE", pars, modelMZ, modelDZ, minus2ll, obj, Conf) 

# Run
ACEFit_depression12 <- mxTryHard(ACEModel, intervals = TRUE)
depression12_estimates <- summary(ACEFit_depression12)

# Table
depression12_table <- ACE_esitmates_table(
  data = depression12_estimates, 
  variable = "Depression",
  age = "12")
depression12_table
```

Compare the scalar and univariate models:

```{r compare scalar and uni anxiety18}
# compare scalar and uni
mxCompare(ACEFit_depression12, ScACE_depression12_fit)
```
Here, the AIC is not lower for the scalar model, therefore we would say the univariate ACE has the better fit. 

Compare to heterogeneity model:

```{r compare scalar and uni depression12}
# compare scalar and univariate model using the AIC
mxCompare(hetace_fit_depression12, ScACE_depression12_fit)
```

**The univariate ACE model is still the best fit for depression age 12** - compared to the heterogeneity model. 

### Age 18

```{r select which variables to use}
# selvars_dep18_chosen <- selvars_dep18                 # raw
# selvars_dep18_chosen <- selvars_dep18_z_score         # standardised
# selvars_dep18_chosen <- selvars_dep18norm             # normalised
# selvars_dep18_chosen <- selvars_dep18_z_score_norm    # standardised and normalised
  selvars_dep18_chosen <- selvars_dep18norm_reg         # normalised and regressed sex
```

```{r social isolation scalar}
# objectives   
objMZM	  <- mxExpectationNormal(covariance = "ExpCovMZM", means = "expMeanMM", dimnames = selvars_dep18_chosen)
objDZM	  <- mxExpectationNormal(covariance = "ExpCovDZM", means = "expMeanMM", dimnames = selvars_dep18_chosen)
objMZF	  <- mxExpectationNormal(covariance = "ExpCovMZF", means = "expMeanFF", dimnames = selvars_dep18_chosen)
objDZF	  <- mxExpectationNormal(covariance = "ExpCovDZF", means = "expMeanFF", dimnames = selvars_dep18_chosen)

# Combine Groups
parsm	  	<- list(pathAm, pathCm, pathEm, covAm, covCm, covEm, covm, StAm, StCm, StEm, fitFunction)
modelMZM 	<- mxModel(parsm, MeanM, MeanMM, covMZM, dataMZM, objMZM,  name = "MZM")
modelDZM 	<- mxModel(parsm, MeanM, MeanMM, covDZM, dataDZM, objDZM,  name = "DZM")
modelMZF 	<- mxModel(parsm, MeanF, MeanFF, covMZF, dataMZF, objMZF, ScalarF, name = "MZF")
modelDZF 	<- mxModel(parsm, MeanF, MeanFF, covDZF, dataDZF, objDZF, ScalarF, name = "DZF")

minus2ll  <- mxAlgebra(expression =	MZM.objective + DZM.objective + MZF.objective + DZF.objective, name = "m2LL")
obj	    	<- mxFitFunctionAlgebra("m2LL")
ci	    	<- mxCI(c('MZM.h2m[1,1]','MZM.c2m[1,1]','MZM.e2m[1,1]') )

ScACE2Model	<-mxModel("ScACE2", parsm, modelMZM, modelDZM, modelMZF, modelDZF, minus2ll, obj, ci) 

# Run
ScACE_depression18_fit	<- mxTryHard(ScACE2Model, intervals = TRUE)
ScACE_depression18_sum <- summary(ScACE_depression18_fit)
```

```{r output social depression 18 scalar}
# ACE estimates
scalar_depression18_table <- ACE_esitmates_table(
  data = ScACE_depression18_sum, 
  variable = "Depression",
  age = "18")

scalar_depression18_table

# scalar estimate 
mxEval(MZF.ScF, ScACE_depression18_fit)
```

```{r run age 18 depression ACE}
# Objective objects for Multiple Groups
objMZ <- mxExpectationNormal(covariance = "expCovMZ", means = "expMeanMZ", dimnames = selvars_dep18_chosen)
objDZ <- mxExpectationNormal(covariance = "expCovDZ", means = "expMeanDZ", dimnames = selvars_dep18_chosen)

# Combine Groups
pars <- list(pathA, pathC, pathE, covA, covC, covE, covP, sd, stpathA, stpathC, stpathE, stcovA, stcovC, stcovE, fitFunction)

# MZ and DZ models
modelMZ	<- mxModel(pars, MeanMZ, covMZ, dataMZ, objMZ, name = "MZ")
modelDZ	<- mxModel(pars, MeanDZ, covDZ, dataDZ, objDZ, name = "DZ")

# Fitting
minus2ll	<- mxAlgebra(expression = MZ.objective + DZ.objective, name = "m2LL")
obj	    	<- mxFitFunctionAlgebra("m2LL")
Conf	  	<- mxCI(c('ACE.h2', 'ACE.c2', 'ACE.e2'))

# Model
ACEModel <- mxModel("ACE", pars, modelMZ, modelDZ, minus2ll, obj, Conf) 

# Run
ACEFit_depression18 <- mxTryHard(ACEModel, intervals = TRUE)
depression18_estimates <- summary(ACEFit_depression18)

# Table
depression18_table <- ACE_esitmates_table(
  data = depression18_estimates, 
  variable = "Depression",
  age = "18")
depression18_table
```

Compare the scalar and univariate models:

```{r compare scalar and uni depression18}
# compare scalar and uni
mxCompare(ACEFit_depression18, ScACE_depression18_fit)
```

Here, the AIC is lower for the scalar model, therefore the scalar model is the best fit for the data.

Compare to heterogeneity model:

```{r compare scalar and uni depression18}
# compare scalar and univariate model using the AIC
mxCompare(hetace_fit_depression18, ScACE_depression18_fit)
```

**The scalar model is still the best fit for depression age 18** - compared to the heterogeneity model. 

## Conduct disorder

### Age 12

```{r select which variables to use}
# selvars_con12_chosen <- selvars_con12                 # raw
# selvars_con12_chosen <- selvars_con12_z_score         # standardised
# selvars_con12_chosen <- selvars_con12norm             # normalised
# selvars_con12_chosen <- selvars_con12_z_score_norm    # standardised and normalised
  selvars_con12_chosen <- selvars_con12norm_reg         # normalised and regressed sex
```

```{r social isolation scalar}
# objectives   
objMZM	  <- mxExpectationNormal(covariance = "ExpCovMZM", means = "expMeanMM", dimnames = selvars_con12_chosen)
objDZM	  <- mxExpectationNormal(covariance = "ExpCovDZM", means = "expMeanMM", dimnames = selvars_con12_chosen)
objMZF	  <- mxExpectationNormal(covariance = "ExpCovMZF", means = "expMeanFF", dimnames = selvars_con12_chosen)
objDZF	  <- mxExpectationNormal(covariance = "ExpCovDZF", means = "expMeanFF", dimnames = selvars_con12_chosen)

# Combine Groups
parsm	  	<- list(pathAm, pathCm, pathEm, covAm, covCm, covEm, covm, StAm, StCm, StEm, fitFunction)
modelMZM 	<- mxModel(parsm, MeanM, MeanMM, covMZM, dataMZM, objMZM,  name = "MZM")
modelDZM 	<- mxModel(parsm, MeanM, MeanMM, covDZM, dataDZM, objDZM,  name = "DZM")
modelMZF 	<- mxModel(parsm, MeanF, MeanFF, covMZF, dataMZF, objMZF, ScalarF, name = "MZF")
modelDZF 	<- mxModel(parsm, MeanF, MeanFF, covDZF, dataDZF, objDZF, ScalarF, name = "DZF")

minus2ll  <- mxAlgebra(expression =	MZM.objective + DZM.objective + MZF.objective + DZF.objective, name = "m2LL")
obj	    	<- mxFitFunctionAlgebra("m2LL")
ci	    	<- mxCI(c('MZM.h2m[1,1]','MZM.c2m[1,1]','MZM.e2m[1,1]') )

ScACE2Model	<-mxModel("ScACE2", parsm, modelMZM, modelDZM, modelMZF, modelDZF, minus2ll, obj, ci) 

# Run
ScACE_conduct12_fit	<- mxTryHard(ScACE2Model, intervals = TRUE)
ScACE_conduct12_sum <- summary(ScACE_conduct12_fit)
```

```{r output social conduct 12 scalar}
# ACE estimates
scalar_conduct12_table <- ACE_esitmates_table(
  data = ScACE_conduct12_sum, 
  variable = "Conduct disorder",
  age = "12")

scalar_conduct12_table

# scalar estimate 
mxEval(MZF.ScF, ScACE_conduct12_fit)
```

```{r run age 12 conduct ACE}
# Objective objects for Multiple Groups
objMZ <- mxExpectationNormal(covariance = "expCovMZ", means = "expMeanMZ", dimnames = selvars_con12_chosen)
objDZ <- mxExpectationNormal(covariance = "expCovDZ", means = "expMeanDZ", dimnames = selvars_con12_chosen)

# Combine Groups
pars <- list(pathA, pathC, pathE, covA, covC, covE, covP, sd, stpathA, stpathC, stpathE, stcovA, stcovC, stcovE, fitFunction)

# MZ and DZ models
modelMZ	<- mxModel(pars, MeanMZ, covMZ, dataMZ, objMZ, name = "MZ")
modelDZ	<- mxModel(pars, MeanDZ, covDZ, dataDZ, objDZ, name = "DZ")

# Fitting
minus2ll	<- mxAlgebra(expression = MZ.objective + DZ.objective, name = "m2LL")
obj	    	<- mxFitFunctionAlgebra("m2LL")
Conf	  	<- mxCI(c('ACE.h2', 'ACE.c2', 'ACE.e2'))

# Model
ACEModel <- mxModel("ACE", pars, modelMZ, modelDZ, minus2ll, obj, Conf) 

# Run
ACEFit_conduct12 <- mxTryHard(ACEModel, intervals = TRUE)
conduct12_estimates <- summary(ACEFit_conduct12)

# Table
conduct12_table <- ACE_esitmates_table(
  data = conduct12_estimates, 
  variable = "Conduct disorder",
  age = "12")
conduct12_table
```

Compare the scalar and univariate models:

```{r compare scalar and uni depression18}
# compare scalar and uni
mxCompare(ACEFit_conduct12, ScACE_conduct12_fit)
```

Here, the AIC is lower for the scalar model, therefore the scalar model is the best fit for the data. 

Compare to heterogeneity model:

```{r compare scalar and uni conduct12}
# compare scalar and univariate model using the AIC
mxCompare(hetace_fit_conduct12, ScACE_conduct12_fit)
```

**The heterogeneity model is still the best fit for conduct age 12** - compared to the scalar model. 

### Age 18

```{r select which variables to use}
# selvars_con18_chosen <- selvars_con18                 # raw
# selvars_con18_chosen <- selvars_con18_z_score         # standardised
# selvars_con18_chosen <- selvars_con18norm             # normalised
# selvars_con18_chosen <- selvars_con18_z_score_norm    # standardised and normalised
  selvars_con18_chosen <- selvars_con18norm_reg         # normalised and regressed sex
```

```{r social isolation scalar}
# objectives   
objMZM	  <- mxExpectationNormal(covariance = "ExpCovMZM", means = "expMeanMM", dimnames = selvars_con18_chosen)
objDZM	  <- mxExpectationNormal(covariance = "ExpCovDZM", means = "expMeanMM", dimnames = selvars_con18_chosen)
objMZF	  <- mxExpectationNormal(covariance = "ExpCovMZF", means = "expMeanFF", dimnames = selvars_con18_chosen)
objDZF	  <- mxExpectationNormal(covariance = "ExpCovDZF", means = "expMeanFF", dimnames = selvars_con18_chosen)

# Combine Groups
parsm	  	<- list(pathAm, pathCm, pathEm, covAm, covCm, covEm, covm, StAm, StCm, StEm, fitFunction)
modelMZM 	<- mxModel(parsm, MeanM, MeanMM, covMZM, dataMZM, objMZM,  name = "MZM")
modelDZM 	<- mxModel(parsm, MeanM, MeanMM, covDZM, dataDZM, objDZM,  name = "DZM")
modelMZF 	<- mxModel(parsm, MeanF, MeanFF, covMZF, dataMZF, objMZF, ScalarF, name = "MZF")
modelDZF 	<- mxModel(parsm, MeanF, MeanFF, covDZF, dataDZF, objDZF, ScalarF, name = "DZF")

minus2ll  <- mxAlgebra(expression =	MZM.objective + DZM.objective + MZF.objective + DZF.objective, name = "m2LL")
obj	    	<- mxFitFunctionAlgebra("m2LL")
ci	    	<- mxCI(c('MZM.h2m[1,1]','MZM.c2m[1,1]','MZM.e2m[1,1]') )

ScACE2Model	<-mxModel("ScACE2", parsm, modelMZM, modelDZM, modelMZF, modelDZF, minus2ll, obj, ci) 

# Run
ScACE_conduct18_fit	<- mxTryHard(ScACE2Model, intervals = TRUE)
ScACE_conduct18_sum <- summary(ScACE_conduct18_fit)
```

```{r output conduct 18 scalar}
# ACE estimates
scalar_conduct18_table <- ACE_esitmates_table(
  data = ScACE_conduct18_sum, 
  variable = "Conduct disorder",
  age = "18")

scalar_conduct18_table

# scalar estimate 
mxEval(MZF.ScF, ScACE_conduct18_fit)
```

```{r run age 18 conduct ACE}
# Objective objects for Multiple Groups
objMZ <- mxExpectationNormal(covariance = "expCovMZ", means = "expMeanMZ", dimnames = selvars_con18_chosen)
objDZ <- mxExpectationNormal(covariance = "expCovDZ", means = "expMeanDZ", dimnames = selvars_con18_chosen)

# Combine Groups
pars <- list(pathA, pathC, pathE, covA, covC, covE, covP, sd, stpathA, stpathC, stpathE, stcovA, stcovC, stcovE, fitFunction)

# MZ and DZ models
modelMZ	<- mxModel(pars, MeanMZ, covMZ, dataMZ, objMZ, name = "MZ")
modelDZ	<- mxModel(pars, MeanDZ, covDZ, dataDZ, objDZ, name = "DZ")

# Fitting
minus2ll	<- mxAlgebra(expression = MZ.objective + DZ.objective, name = "m2LL")
obj	    	<- mxFitFunctionAlgebra("m2LL")
Conf	  	<- mxCI(c('ACE.h2', 'ACE.c2', 'ACE.e2'))

# Model
ACEModel <- mxModel("ACE", pars, modelMZ, modelDZ, minus2ll, obj, Conf) 

# Run
ACEFit_conduct18 <- mxTryHard(ACEModel, intervals = TRUE)
conduct18_estimates <- summary(ACEFit_conduct18)

# Table
conduct18_table <- ACE_esitmates_table(
  data = conduct18_estimates, 
  variable = "Conduct disorder",
  age = "18")
conduct18_table
```

Compare the scalar and univariate models:

```{r compare scalar and uni conduct18}
# compare scalar and uni
mxCompare(ACEFit_conduct18, ScACE_conduct18_fit)
```

Here, the AIC is (only just) lower for the scalar model, therefore the scalar model is the best fit for the data. 

Compare to heterogeneity model:

```{r compare scalar and uni conduct18}
# compare scalar and univariate model using the AIC
mxCompare(hetace_fit_conduct18, ScACE_conduct18_fit)
```

**The scalar model is still the best fit for conduct age 18** - compared to the heterogeneity model. 


## Psychotic experiences

### Age 12

```{r select which variables to use}
# selvars_psy12_chosen <- selvars_psy12                 # raw
# selvars_psy12_chosen <- selvars_psy12_z_score         # standardised
# selvars_psy12_chosen <- selvars_psy12norm             # normalised
# selvars_psy12_chosen <- selvars_psy12_z_score_norm    # standardised and normalised
  selvars_psy12_chosen <- selvars_psy12norm_reg         # normalised and regressed sex
```

```{r social isolation scalar}
# objectives   
objMZM	  <- mxExpectationNormal(covariance = "ExpCovMZM", means = "expMeanMM", dimnames = selvars_psy12_chosen)
objDZM	  <- mxExpectationNormal(covariance = "ExpCovDZM", means = "expMeanMM", dimnames = selvars_psy12_chosen)
objMZF	  <- mxExpectationNormal(covariance = "ExpCovMZF", means = "expMeanFF", dimnames = selvars_psy12_chosen)
objDZF	  <- mxExpectationNormal(covariance = "ExpCovDZF", means = "expMeanFF", dimnames = selvars_psy12_chosen)

# Combine Groups
parsm	  	<- list(pathAm, pathCm, pathEm, covAm, covCm, covEm, covm, StAm, StCm, StEm, fitFunction)
modelMZM 	<- mxModel(parsm, MeanM, MeanMM, covMZM, dataMZM, objMZM,  name = "MZM")
modelDZM 	<- mxModel(parsm, MeanM, MeanMM, covDZM, dataDZM, objDZM,  name = "DZM")
modelMZF 	<- mxModel(parsm, MeanF, MeanFF, covMZF, dataMZF, objMZF, ScalarF, name = "MZF")
modelDZF 	<- mxModel(parsm, MeanF, MeanFF, covDZF, dataDZF, objDZF, ScalarF, name = "DZF")

minus2ll  <- mxAlgebra(expression =	MZM.objective + DZM.objective + MZF.objective + DZF.objective, name = "m2LL")
obj	    	<- mxFitFunctionAlgebra("m2LL")
ci	    	<- mxCI(c('MZM.h2m[1,1]','MZM.c2m[1,1]','MZM.e2m[1,1]') )

ScACE2Model	<-mxModel("ScACE2", parsm, modelMZM, modelDZM, modelMZF, modelDZF, minus2ll, obj, ci) 

# Run
ScACE_psychosis12_fit	<- mxTryHard(ScACE2Model, intervals = TRUE)
ScACE_psychosis12_sum <- summary(ScACE_psychosis12_fit)
```

```{r output psychosis 12 scalar}
# ACE estimates
scalar_psychosis12_table <- ACE_esitmates_table(
  data = ScACE_psychosis12_sum, 
  variable = "Psychosis",
  age = "12")

scalar_psychosis12_table

# scalar estimate 
mxEval(MZF.ScF, ScACE_psychosis12_fit)
```

```{r run age 12 psychosis ACE}
# Objective objects for Multiple Groups
objMZ <- mxExpectationNormal(covariance = "expCovMZ", means = "expMeanMZ", dimnames = selvars_psy12_chosen)
objDZ <- mxExpectationNormal(covariance = "expCovDZ", means = "expMeanDZ", dimnames = selvars_psy12_chosen)

# Combine Groups
pars <- list(pathA, pathC, pathE, covA, covC, covE, covP, sd, stpathA, stpathC, stpathE, stcovA, stcovC, stcovE, fitFunction)

# MZ and DZ models
modelMZ	<- mxModel(pars, MeanMZ, covMZ, dataMZ, objMZ, name = "MZ")
modelDZ	<- mxModel(pars, MeanDZ, covDZ, dataDZ, objDZ, name = "DZ")

# Fitting
minus2ll	<- mxAlgebra(expression = MZ.objective + DZ.objective, name = "m2LL")
obj	    	<- mxFitFunctionAlgebra("m2LL")
Conf	  	<- mxCI(c('ACE.h2', 'ACE.c2', 'ACE.e2'))

# Model
ACEModel <- mxModel("ACE", pars, modelMZ, modelDZ, minus2ll, obj, Conf) 

# Run
ACEFit_psychosis12 <- mxTryHard(ACEModel, intervals = TRUE)
psychosis12_estimates <- summary(ACEFit_psychosis12)

# Table
psychosis12_table <- ACE_esitmates_table(
  data = psychosis12_estimates, 
  variable = "Psychosis",
  age = "12")
psychosis12_table
```

Compare the scalar and univariate models:

```{r compare scalar and uni conduct18}
# compare scalar and uni
mxCompare(ACEFit_psychosis12, ScACE_psychosis12_fit)
```

The scalar has a lower AIC and therefore a better fit to the data. 

Compare to heterogeneity model:

```{r compare scalar and uni psychosis12}
# compare scalar and univariate model using the AIC
mxCompare(hetace_fit_psychosis12, ScACE_psychosis12_fit)
```

**The scalar model is still the best fit for psychosis age 12** - compared to the heterogeneity model. 

### Age 18

```{r select which variables to use}
# selvars_psy18_chosen <- selvars_psy18                 # raw
# selvars_psy18_chosen <- selvars_psy18_z_score         # standardised
# selvars_psy18_chosen <- selvars_psy18norm             # normalised
# selvars_psy18_chosen <- selvars_psy18_z_score_norm    # standardised and normalised
  selvars_psy18_chosen <- selvars_psy18norm_reg         # normalised and regressed sex
```

```{r social isolation scalar}
# objectives   
objMZM	  <- mxExpectationNormal(covariance = "ExpCovMZM", means = "expMeanMM", dimnames = selvars_psy18_chosen)
objDZM	  <- mxExpectationNormal(covariance = "ExpCovDZM", means = "expMeanMM", dimnames = selvars_psy18_chosen)
objMZF	  <- mxExpectationNormal(covariance = "ExpCovMZF", means = "expMeanFF", dimnames = selvars_psy18_chosen)
objDZF	  <- mxExpectationNormal(covariance = "ExpCovDZF", means = "expMeanFF", dimnames = selvars_psy18_chosen)

# Combine Groups
parsm	  	<- list(pathAm, pathCm, pathEm, covAm, covCm, covEm, covm, StAm, StCm, StEm, fitFunction)
modelMZM 	<- mxModel(parsm, MeanM, MeanMM, covMZM, dataMZM, objMZM,  name = "MZM")
modelDZM 	<- mxModel(parsm, MeanM, MeanMM, covDZM, dataDZM, objDZM,  name = "DZM")
modelMZF 	<- mxModel(parsm, MeanF, MeanFF, covMZF, dataMZF, objMZF, ScalarF, name = "MZF")
modelDZF 	<- mxModel(parsm, MeanF, MeanFF, covDZF, dataDZF, objDZF, ScalarF, name = "DZF")

minus2ll  <- mxAlgebra(expression =	MZM.objective + DZM.objective + MZF.objective + DZF.objective, name = "m2LL")
obj	    	<- mxFitFunctionAlgebra("m2LL")
ci	    	<- mxCI(c('MZM.h2m[1,1]','MZM.c2m[1,1]','MZM.e2m[1,1]') )

ScACE2Model	<-mxModel("ScACE2", parsm, modelMZM, modelDZM, modelMZF, modelDZF, minus2ll, obj, ci) 

# Run
ScACE_psychosis18_fit	<- mxTryHard(ScACE2Model, intervals = TRUE)
ScACE_psychosis18_sum <- summary(ScACE_psychosis18_fit)
```

```{r output psychosis 18 scalar}
# ACE estimates
scalar_psychosis18_table <- ACE_esitmates_table(
  data = ScACE_psychosis18_sum, 
  variable = "Psychosis",
  age = "18")

scalar_psychosis18_table

# scalar estimate 
mxEval(MZF.ScF, ScACE_psychosis18_fit)
```

```{r run age 18 psychosis ACE}
# Objective objects for Multiple Groups
objMZ <- mxExpectationNormal(covariance = "expCovMZ", means = "expMeanMZ", dimnames = selvars_psy18_chosen)
objDZ <- mxExpectationNormal(covariance = "expCovDZ", means = "expMeanDZ", dimnames = selvars_psy18_chosen)

# Combine Groups
pars <- list(pathA, pathC, pathE, covA, covC, covE, covP, sd, stpathA, stpathC, stpathE, stcovA, stcovC, stcovE, fitFunction)

# MZ and DZ models
modelMZ	<- mxModel(pars, MeanMZ, covMZ, dataMZ, objMZ, name = "MZ")
modelDZ	<- mxModel(pars, MeanDZ, covDZ, dataDZ, objDZ, name = "DZ")

# Fitting
minus2ll	<- mxAlgebra(expression = MZ.objective + DZ.objective, name = "m2LL")
obj	    	<- mxFitFunctionAlgebra("m2LL")
Conf	  	<- mxCI(c('ACE.h2', 'ACE.c2', 'ACE.e2'))

# Model
ACEModel <- mxModel("ACE", pars, modelMZ, modelDZ, minus2ll, obj, Conf) 

# Run
ACEFit_psychosis18 <- mxTryHard(ACEModel, intervals = TRUE)
psychosis18_estimates <- summary(ACEFit_psychosis18)

# Table 
psychosis18_table <- ACE_esitmates_table(
  data = psychosis18_estimates, 
  variable = "Psychosis",
  age = "18")
psychosis18_table
```

Compare the scalar and univariate models:

```{r compare scalar and uni conduct18}
# compare scalar and uni
mxCompare(ACEFit_psychosis18, ScACE_psychosis18_fit)
```

The scalar has a lower AIC and therefore a better fit to the data. 

Compare to heterogeneity model:

```{r compare scalar and uni psychosis18}
# compare scalar and univariate model using the AIC
mxCompare(hetace_fit_psychosis18, ScACE_psychosis18_fit)
```

**The scalar model is still the best fit for psychosis age 18** - compared to the heterogeneity model.

# Output tables

```{r output tables}
ACE_estimates_hetace <- rbind(hetace_isolation12_table, 
                          #   hetace_anxiety12_table,
                              hetace_depression12_table,
                              hetace_conduct12_table,
                              hetace_psychosis12_table, 
                              hetace_isolation18_table,
                          #   hetace_anxiety18_table,
                              hetace_depression18_table,  
                              hetace_conduct18_table,
                              hetace_psychosis18_table) 

ACE_estimates_hetace

ACE_estimates_homace <- rbind(homace_isolation12_table, 
                          #   homace_anxiety12_table,
                              homace_depression12_table,
                              homace_conduct12_table,
                              homace_psychosis12_table, 
                              homace_isolation18_table,
                          #   homace_anxiety18_table,
                              homace_depression18_table,  
                              homace_conduct18_table,
                              homace_psychosis18_table)  

ACE_estimates_homace

ACE_estimates <- rbind(isolation12_table, 
                          #   anxiety12_table,
                              depression12_table,
                              conduct12_table,
                              psychosis12_table, 
                              isolation18_table,
                          #   anxiety18_table,
                              depression18_table,  
                              conduct18_table,
                              psychosis18_table) 

ACE_estimates

ACE_estimates_scalar <- rbind(scalar_isolation12_table, 
                          #   scalar_anxiety12_table,
                              scalar_depression12_table,
                              scalar_conduct12_table,
                              scalar_psychosis12_table, 
                              scalar_isolation18_table,
                          #   scalar_anxiety18_table,
                              scalar_depression18_table,  
                              scalar_conduct18_table,
                              scalar_psychosis18_table) 

ACE_estimates_scalar
```