#Neues R Skript - ACE Univariate

#Install packages
devtools::install_github('mirkoruks/twinflex')

install.packages("tidyverse")
install.packages("OpenMx")
install.packages("haven")

library(OpenMx)
library(twinflex)
library(tidyverse)
library(haven) #to read data

#install.packages("umx")
#install.packages("car")
#install.packages("sjlabelled")

#You can set you directory path here
#setwd("C:/Users/R") 

#Load data-set 
data_all <- read_dta("ZA6701_person_wid1_v5-0-0.dta")

data_all$zyg <- as.numeric(as.character(data_all$zyg0102))

#Select a subset
twinlife_data_reduced <- data_all[c("fid","pid","cgr","ptyp", "per0100", "per0200", "per0300"
                                    , "par0100t", "par0100u", "cer2200", "cer2201", "cer1400"
                                    , "cer1500", "zygot")]
#twinlife_data_reduced <- subset(twinlife_data_reduce, select=c(VARIABLEN))

#Select only certain cases, e.g. twins in cohort 2
twinlife_data_c2_twin <- subset(data_all, cgr == 2 & ptyp < 3, select=c("fid","cgr"
                              ,"ptyp", "per0100", "per0200", "per0300", "cer2200"
                              , "cer2201", "cer1400", "cer1500", "zygot"))

#Convert variable to numeric
twinlife_data_c2_twin$cer1500 <- as.numeric(as.character(twinlife_data_c2_twin$cer1500))
twinlife_data_c2_twin$cer2201 <- as.numeric(as.character(twinlife_data_c2_twin$cer2201))

#declare missing values
twinlife_data_c2_twin$cer1500[twinlife_data_c2_twin$cer1500<0] <- NA
twinlife_data_c2_twin$cer2201[twinlife_data_c2_twin$cer2201<0] <- NA
twinlife_data_c2_twin$zygot[twinlife_data_c2_twin$zygot<0] <- NA

twinlife_data_c2_twin$per0100[twinlife_data_c2_twin$per0100<0] <- NA


#Create variable containing the information
twinlife_data_c2_twin$german_grade <- twinlife_data_c2_twin$cer2201
twinlife_data_c2_twin$german_grade[is.na(twinlife_data_c2_twin$german_grade)
                                       ] <- twinlife_data_c2_twin$cer1500[is.na(twinlife_data_c2_twin$german_grade)]

#OR: within(twinlife_data_c2_twin, german_grade[is.na(german_grade)] <- cer1400[is.na(german_grade)])

#Create wide-format data-set
twinlife_data_wide <- stats::reshape(data=as.data.frame(twinlife_data_c2_twin), 
                     v.names=c("per0100", "per0200", "per0300", "cer2200", "cer2201", "cer1400", "cer1500", "german_grade")
                     , timevar="ptyp", idvar="fid", direction="wide", sep="_")


#*******************************************************************************
#Analysis unsing twinflex
help(twinflex)

dataACE <- twinlife_data_wide %>% filter(!is.na(zygot) & !is.na(german_grade_1) & !is.na(german_grade_2))

univ_cont_tf <- twinflex(acevars = "german_grade", data = dataACE, sep = "_")
mxTryHard(univ_cont_tf)
summary(univ_cont_tf)

modelAE <- omxSetParameters(univ_cont_tf, labels= c('c11'), free=FALSE, values=0)
modelfitAE <- mxRun(modelAE)
modelfitAE <- mxTryHard(modelAE)

summary(modelfitAE)

#-------------------------------------------------------------------------------
#*******************************************************************************
#Analysis unsing UMX
library(umx)
help(umxACE)

#Creates two seperate data-sets: one for MZ, one for DZ
dz = subset(twinlife_data_wide, zygot == 2)
mz = subset(twinlife_data_wide, zygot == 1)

#ACE-Model
ACE1 = umxACE(selDVs = "german_grade", sep = "_", dzData = dz, mzData = mz,  addStd = TRUE, addCI = TRUE,) 
plot(ACE1)

ACE1v = umxACEv(selDVs = "math_grade", sep = "_", dzData = dz, mzData = mz,  addStd = TRUE, addCI = TRUE,) 

umxCI(ACE1, run = "yes")
summary(ACE1)

#Compares possible models: ACE, ADE, AE, CE, (DE)
umxReduceACE(ACE1)

# Using TryHard
set.seed(102)
ACE1 = umxACE(selDVs = "math_grade", sep = "", dzData = dz, mzData = mz,  addStd = TRUE, addCI = TRUE, tryHard = "yes")