# TABLE OF CONTENTS # Data simulation # Dual change score - Nathan Gillespie & Brad Verhulst 2018 Twin Workshop in Boulder # install.packages("devtools") # require("devtools") # devtools::install_github("tbates/umx") # install.packages("umx") #source("https://openmx.ssri.psu.edu/software/getOpenMx.R") # Load Libraries & Options rm(list=ls()) library(OpenMx); library(umx) library(psych); library(polycor); require(MASS) umx_set_optimizer(opt="NPSOL") #mxOption(NULL, 'Default optimizer', 'NPSOL') #mxOption(NULL, 'Default optimizer', 'SLSQP') umxVersion() # Data simulation set.seed(12345) ##$#$#$#$#$#$#$#$ # Latent Growth Simplex Model (i.e. Dual Change Score) on 5 Repeatedly Measured Latent Traits nVariables <- 5 nFactors <- 5 nSubjects <- 2000 nGrow <- 3 # The covariance of y is lambda.y %*% solve(I-B)%*% (gamma%*%phi%*%gamma + psi) %*% t(solve(I-B))%*%t(lambda.y) + theta.epsilon # Specify the lambda.y matrix - Factor loadings on the latent factors lambda.y <- diag(5) # Specify the I-B^-1 - The causal pathways between the latent factors I5 <- diag(5) B <- matrix(c(0,.9,0,0,0, # You can change these values 0,0,.9,0,0, 0,0,0,.9,0, 0,0,0,0,.9, 0,0,0,0,0),5) # Specify gamma - Factor loadings on the exogenous variables (in this case growth parameters) gamma <- matrix(c( 1, 1, 1, 1, 1, 0, 1, 2, 3, 4, 0, 1, 4, 9,16),5,3) # Specify phi - covariance of the exogenous variables phi <- matrix(c(.9,.4,.1, .4,.7,.2, .1,.2,.4),3 ) # You can change these values #eigen(phi)$values cov2cor(phi) #.9,.4,.1, #.4,.7,.2, #.1,.2,.4 #Specify psi - residuals of the latent factors psi <- .5 * diag(5) # You can change these values psi[1,1] <- 0 # Specify theta.epsilon - the residuals of the manifest variables theta.epsilon <- .1 * diag(nVariables) # You can change these values impCov <- lambda.y %*% solve(I5-B) %*% (gamma%*%phi%*%t(gamma) + psi) %*% t(solve(I5-B)) %*% t(lambda.y) + theta.epsilon h2 <- .5 c2 <- .25 MZ <- matrix(c(1, (h2 + c2), (h2+c2),1 ), 2) DZ <- matrix(c(1, (.5*h2 + c2), (.5*h2+c2),1 ), 2) impCovMZ <- MZ %x% impCov impCovDZ <- DZ %x% impCov LatGrowMean <- matrix(c(2, .8, -.1),3,1) # You can change these values LatFacMean <- solve(I5-B)%*%gamma %*% LatGrowMean OBSMeans <- lambda.y %*% LatFacMean set.seed(2016) MZdata <- mvrnorm(nSubjects, mu = c(OBSMeans, OBSMeans), Sigma = impCovMZ, empirical = T) colnames(MZdata) <- c(paste("t1", 1:5, sep = "_"), paste("t2", 1:5, sep = "_")) set.seed(2016) DZdata <- mvrnorm(nSubjects, mu = c(OBSMeans, OBSMeans), Sigma = impCovDZ, empirical = T) colnames(DZdata) <- c(paste("t1", 1:5, sep = "_"), paste("t2", 1:5, sep = "_")) # write.table(DZdata, "dz_longitudinal_dcs_data.txt") # write.table(MZdata, "mz_longitudinal_dcs_data.txt") # Dual change score - Genetically informative #Read data nvars <- 5 # Number of manifest/observed variables ntvars <- nvars*2 # nFactors <- 5 # Number of latent factor true scores nGrow <- 3 # Number of growth parameters: intercept; linear; & quadratic rates of change # Select Data # See DualChangeSim.R for simulated longitudinal twin pair data selVars <- c("t1_1", "t1_2", "t1_3", "t1_4", "t1_5", "t2_1", "t2_2", "t2_3", "t2_4", "t2_5") # mzData <- subset(data, zyg==1, selVars) # dzData <- subset(data, zyg==3, selVars) head(MZdata) # Based on simulation above head(DZdata) # Based on simulation above # Other Twin Data functions: # ?umx # umx_long2wide - Take a long twin-data file and make it wide (one family per row) # umx_wide2long - Change data 2-twin family data from wide to long format # umx_make_TwinData - Simulate twin data with control over A, C, and E parameters # umx_residualize - Easily residualize variables in long or wide dataframes # Print Descriptive Statistics round(colMeans(MZdata,na.rm=TRUE),2) round(colMeans(DZdata,na.rm=TRUE),2) round(cor(MZdata,use="complete"),2)[1:5,1:5] round(cor(DZdata,use="complete"),2)[1:5,1:5] round(cov(MZdata,use="complete"),2)[1:5,1:5] round(cov(DZdata,use="complete"),2)[1:5,1:5] # Convert a covariance matrix into a correlation matrix # umxCov2cor(cov(MZdata)) # umxCov2cor(cov(MZdata))[1:5,1:5] # round(umxCov2cor(cov(MZdata))[1:5,1:5],2) # Matrices in the full DCS model # ϕ phi = Cholesky Decomposition pathway coefficients to calculate correlations between Intercept, Slope and Quadratic # Γ gamma = Factor loadings from growth parameters to true scores / latent factors. Pre-defined/Set by you. # Ψ psi = Innovations, residuals for the latent factor true scores (variance not accounted by growth factors) # = Captures the age-specific / new variance 'emerging' at each time point. # Β beta = Auto-regression pathway. The contribution from one time point to the next. # I identity # Λ lamba = Factor loadings from true scores / latent factors to observed variable. Set to 1. # ϵ epsilon = Variance not captured by the growth or autoregressive processes. # ϕ phi # covariance between the exogenous variables i.e. between Intercept, Slope, & Quadractic latent factors phiA_labs <- c("VA_Int","Acov_IntS","Acov_IntQ","VA_S","Acov_SQ","VA_Q") phiA <- mxMatrix(type = "Symm", nrow = nGrow, ncol = nGrow, free = T, labels=phiA_labs, name = "phiA", values = c( 0.9, 0.5,0.3, 0.5,0.5,0.2)) # SymmMatrix 'phiA' # $labels # [,1] [,2] [,3] # [1,] "VA_I" "Acov_IS" "Acov_IQ" # [2,] "Acov_IS" "VA_S" Acov_SQ" # [3,] "Acov_IQ" "Acov_SQ" "VA_Q" # $values # [,1] [,2] [,3] # [1,] 0.9 0.5 0.3 # [2,] 0.5 0.5 0.5 # [3,] 0.3 0.5 0.2 phiC_labs <- c("VC_Int","Ccov_IntS","Ccov_IntQ","VC_S","Ccov_SQ","VC_Q") phiE_labs <- c("VE_Int","Ecov_IntS","Ecov_IntQ","VE_S","Ecov_SQ","VE_Q") phiC <- mxMatrix(type = "Symm", nrow = nGrow, ncol = nGrow, free = T, labels=phiC_labs, name = "phiC", values = c( 0.01, 0.01,0.01, 0.01,0.01,0.01)) phiE <- mxMatrix(type = "Symm", nrow = nGrow, ncol = nGrow, free = T, labels=phiE_labs, name = "phiE", values = c( 0.60, 0.21,0.60, 0.21,0.21,0.60)) # correlations between Intercept, Slope, & Quadractic latent factors corrA <- mxAlgebra( expression=cov2cor(phiA), name="corrA" ) # A latent factor correlations corrC <- mxAlgebra( expression=cov2cor(phiC), name="corrC" ) # C latent factor correlations corrE <- mxAlgebra( expression=cov2cor(phiE), name="corrE" ) # E latent factor correlations # Γ gamma # Factor loadings on the exogenous variables (intercept & growth parameters [I,S,Q]) onto the latent true scores (Eta 1-5) # Factor loadings from Intercept, Slope, & Quadractic to latent true scores # 5 (time points) x 3 (growth parameters) matrix Glabs <- c(paste("i",1:nvars,sep=""), paste("s",1:nvars,sep=""), paste("q",1:nvars,sep="")) gamma <- mxMatrix( type="Full", nrow=nFactors, ncol=nGrow, free=F, labels=Glabs, name="gamma", values=c(1,1,1,1,1, 0,1,2,3,4, 0,1,4,9,16)) # FullMatrix 'gamma' # [,1] [,2] [,3] # [1,] "i1" "s1" "q1" # [2,] "i2" "s2" "q2" # [3,] "i3" "s3" "q3" # [4,] "i4" "s4" "q4" # [5,] "i5" "s5" "q5" # $values # [,1] [,2] [,3] # [1,] 1 0 0 # [2,] 1 1 1 # [3,] 1 2 4 # [4,] 1 3 9 # [5,] 1 4 16 # Ψ psi # Innovations in the auto-regression model # Variance components to capture 'residuals' on the latent factor true scores (Eta 1-5) psiA_labs <- c("VA_i11","VA_i22","VA_i33","VA_i44","VA_i55") #psiA_labs <- c("VA_i11","VA_i22","VA_i22","VA_i22","VA_i22") psiA <- mxMatrix(type = "Diag", nrow = nFactors, ncol = nFactors, free = c(F,T,T,T,T), labels = psiA_labs, name = "psiA", values = c( 0,0.5,0.5,0.5,0.5)) # DiagMatrix 'psi_a' # $labels # [,1] [,2] [,3] [,4] [,5] # [1,] "VA_i11" NA NA NA NA # [2,] NA "VA_i22" NA NA NA # [3,] NA NA "VA_i33" NA NA # [4,] NA NA NA "VA_i44" NA # [5,] NA NA NA NA "VA_i55" # $values # [,1] [,2] [,3] [,4] [,5] # [1,] 0 0.0 0.0 0.0 0.0 # [2,] 0 0.9 0.0 0.0 0.0 # [3,] 0 0.0 0.1 0.0 0.0 # [4,] 0 0.0 0.0 0.1 0.0 # [5,] 0 0.0 0.0 0.0 0.1 psiC_labs <- c("VC_i11","VC_i22","VC_i33","VC_i44","VC_i55") psiE_labs <- c("VE_i11","VE_i22","VE_i33","VE_i44","VE_i55") #psiC_labs <- c("VC_i11","VC_i22","VC_i22","VC_i22","VC_i22") #psiE_labs <- c("VE_i11","VE_i22","VE_i22","VE_i22","VE_i22") psiC <- mxMatrix(type = "Diag", nrow = nFactors, ncol = nFactors, free = c(F,T,T,T,T), labels = psiC_labs, name = "psiC", values = c(0,0.5,0.5,0.5,0.5)) psiE <- mxMatrix(type = "Diag", nrow = nFactors, ncol = nFactors, free = c(F,T,T,T,T), labels = psiE_labs, name = "psiE", values = c(0,0.5,0.5,0.5,0.5)) # Β beta # Auto-regression pathways. The contribution from one time point to the next. # "The direction always goes down the column & out the row" betaF <- matrix(c(F,T,F,F,F, F,F,T,F,F, F,F,F,T,F, F,F,F,F,T, F,F,F,F,F),5) # Free elements betaS <- matrix(c(0,0.9,0,0,0, 0,0,0.9,0,0, 0,0,0,0.9,0, 0,0,0,0,0.9, 0,0,0,0,0),5) # Start values betalab <- matrix(c(NA,"b",NA,NA,NA, NA,NA,"b",NA,NA, NA,NA,NA,"b",NA, NA,NA,NA,NA,"b", NA,NA,NA,NA,NA),5) beta <- mxMatrix(type="Full", nrow = nFactors, ncol = nFactors, free = betaF, labels = betalab, values = betaS,name = "beta") # FullMatrix 'beta' # $labels # [,1] [,2] [,3] [,4] [,5] # [1,] NA NA NA NA NA # [2,] "b" NA NA NA NA # [3,] NA "b" NA NA NA # [4,] NA NA "b" NA NA # [5,] NA NA NA "b" NA # $values # [,1] [,2] [,3] [,4] [,5] # [1,] 0.0 0.0 0.0 0.0 0 # [2,] 0.9 0.0 0.0 0.0 0 # [3,] 0.0 0.9 0.0 0.0 0 # [4,] 0.0 0.0 0.9 0.0 0 # [5,] 0.0 0.0 0.0 0.9 0 # # $free # [,1] [,2] [,3] [,4] [,5] # [1,] FALSE FALSE FALSE FALSE FALSE # [2,] TRUE FALSE FALSE FALSE FALSE # [3,] FALSE TRUE FALSE FALSE FALSE # [4,] FALSE FALSE TRUE FALSE FALSE # [5,] FALSE FALSE FALSE TRUE FALSE ci <- mxCI(c("b")) I <- mxMatrix(type="Iden", nrow = nFactors, ncol = nFactors, name = "I") # Λ lamba # Factor loadings from true scores / latent factors to observed variable. Set to 1. # Matrix of factor loadings from the latent true scores (Etas 1-5) to the observed variables (Time 1-5) loadS <- diag(nFactors) loadF <- F lamba <- mxMatrix(type="Full", nrow = nvars, ncol = nFactors, free = loadF, values = loadS, name = "lamba") # FullMatrix 'lamba' # $labels: No labels assigned. # $values # [,1] [,2] [,3] [,4] [,5] # [1,] 1 0 0 0 0 # [2,] 0 1 0 0 0 # [3,] 0 0 1 0 0 # [4,] 0 0 0 1 0 # [5,] 0 0 0 0 1 # ϵ epsilon # Variance not captured by the growth or autoregressive processes # Residual A, C and E variance components on the observed phenotypic observations (Epsilon 1-5) epsilon_a <- mxMatrix(type = "Diag", nrow = nvars, ncol = nvars, free = T, labels = "res_a", values = 1, name = "epsilon_a") # DiagMatrix 'epsilon_a' # $labels # [,1] [,2] [,3] [,4] [,5] # [1,] "res_a" NA NA NA NA # [2,] NA "res_a" NA NA NA # [3,] NA NA "res_a" NA NA # [4,] NA NA NA "res_a" NA # [5,] NA NA NA NA "res_a" # $values # [,1] [,2] [,3] [,4] [,5] # [1,] 1 0 0 0 0 # [2,] 0 1 0 0 0 # [3,] 0 0 1 0 0 # [4,] 0 0 0 1 0 # [5,] 0 0 0 0 1 epsilon_c <- mxMatrix(type = "Diag", nrow = nvars, ncol = nvars, free = T, labels = "res_c" ,values = 1, name = "epsilon_c") epsilon_e <- mxMatrix(type = "Diag", nrow = nvars, ncol = nvars, free = T, labels = "res_e" ,values = 1, name = "epsilon_e") # Expected variance/covariance matrix # Λ * (I-Β)~ * (Γ * ϕ ' Γ' + Ψ) * (I-Β)~' * Λ' + ϵ # expCov <- mxAlgebra(Λ %*% solve(I-Β)%*% (Γ %*% ϕ %*% t(Γ) + Ψ) %*% t(solve(I-Β))%*%t(Λ) + ϵ, name = "expCov") A <- mxAlgebra(lamba %*% solve(I-beta)%*% (gamma %*% phiA %*% t(gamma) + psiA) %*% t(solve(I-beta))%*%t(lamba) + epsilon_a, name = "A") # C <- mxAlgebra(lamba %*% solve(I-beta)%*% (gamma %*% phiC %*% t(gamma) + psiC) %*% t(solve(I-beta))%*%t(lamba) + epsilon_c, name = "C") # E <- mxAlgebra(lamba %*% solve(I-beta)%*% (gamma %*% phiE %*% t(gamma) + psiE) %*% t(solve(I-beta))%*%t(lamba) + epsilon_e, name = "E") V <- mxAlgebra( expression= A+C+E, name="V" ) 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" ) # Standardize the A, C and E variance components on the latent growth processes totalVg <- mxAlgebra( expression= phiA+phiC+phiE, name="totalVg" ) # "VA_I" "Acov_IS" "Acov_IQ" # "Acov_IS" "VA_S" "Acov_SQ" # "Acov_IQ" "Acov_SQ" "VA_Q" # + # "VC_I" "Ccov_IS" "Ccov_IQ" # "Ccov_IS" "VC_S" "Ccov_SQ" # "Ccov_IQ" "Ccov_SQ" "VE_Q" # + # "VE_I" "Ecov_IS" "Ecov_IQ" # "Ecov_IS" "VE_S" "Ecov_SQ" # "Ecov_IQ" "Ecov_SQ" "VE_Q" stanVg <- mxAlgebra( expression=cbind(phiA/totalVg,phiC/totalVg, phiE/totalVg), name="stanVg") # Expected means # (Λ * ( (I-Β)~ * Γ * μ ) )' GroMean <- mxMatrix(type="Full", nrow = nGrow, ncol = 1, free = T,values = c(2, .8, -.1), labels=c("ui","us","uq"), name = "GroMean") FacMean <- mxAlgebra(solve(I-beta) %*% gamma %*% GroMean, name = "FacMean", dimnames = list(paste("F",1:5, sep=""),"FacMeans")) ManMean <- mxAlgebra(t(lamba %*% FacMean), "ManMean") expMean <- mxAlgebra( expression= cbind(ManMean, ManMean), name="expMean" ) # Expected means for twin pair # Data Objects for Multiple Groups # Read in the simulated data dataMZ <- mxData( observed=MZdata, type="raw" ) dataDZ <- mxData( observed=DZdata, type="raw" ) # Objective Objects for Multiple Groups # Define how the model expectations are calculated # The 'mxExpectationNormal' function uses the algebra defined by the 'covariance' and 'means' arguments to define # the expected covariance and means under the assumption of multivariate normality objMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selVars ) objDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars ) # Function to compute -2*(log likelihood) of the data given the current values of the free parameters and the expectation function fiML <- mxFitFunctionML() # Standardized variance Components rowVC <- rep('VC',nvars) # Create row labels colVC <- rep(c('A','C','E','SA','SC','SE'),each=nvars) # Create column labels VC <- mxAlgebra( expression=cbind(A,C,E,A/V,C/V,E/V), name="VC", dimnames=list(c("Time1","Time2","Time3","Time4","Time5"),colVC)) # Final steps pars <- list( phiA,phiC,phiE,corrA,corrC,corrE,gamma,psiA,psiC,psiE,beta,ci,I,lamba, epsilon_a,epsilon_c,epsilon_e,A,C,E,V,VC,totalVg,stanVg,GroMean,FacMean,ManMean,expMean) modelMZ <- mxModel( pars, dataMZ, covMZ, objMZ, fiML, name="MZ" ) # Function to create MxModel objects modelDZ <- mxModel( pars, dataDZ, covDZ, objDZ, fiML, name="DZ" ) # Function to create MxModel objects obj <- mxFitFunctionMultigroup(c("MZ","DZ")) # Function used to fit a multiple group model dcs_model <- mxModel( "dcs_model", pars, modelMZ, modelDZ, obj, ci ) # Check model identifitication (remove constraints) # mxCheckIdentification(dcs_model, details=TRUE) # Is the model locally identified? # Fetch model parameters omxGetParameters(dcs_model) # Number of p arameters # ϕ phi = 18 # Γ gamma = fixed # Ψ psi = 12 # Β beta = 1 # I identity = fixed # Λ lamba = fixed # ϵ epsilon = 3 # μ means = 3 # Total = ??? # Full (A+C+E) dual-change score model # NB: mxTryHard() uses (pseudo-)random-number generation ∴ set seed # Make multiple attempts to run a model set.seed(2016) summary(dcs_fit <- mxTryHard( dcs_model, extraTries=15, intervals = T, greenOK=T, checkHess=TRUE, fit2beat=Inf, exhaustive=T)) #summary(dcs_fit,verbose=T) # Model parameters contingent upon UPPER and LOWER 'b' parameters # Fetch ALL estimated model parameters round(dcs_fit$output$estimate,4) # Fetch unstandardized & standardized variance components dcs_fit$algebras$VC$result # Fetch just the standardized variance components # estVC <- mxAlgebra( expression=cbind(A,C,E,A/V,C/V,E/V), name="VC", dimnames=list(c("Time1","Time2","Time3","Time4","Time5"),colVC)) round(dcs_fit$algebras$VC$result[,16:30],2) # Examine the variance components on the growth processes # Fetch standardized A, C and E variance components on the intercept, slope (linear) and quadratric # stanVg <- mxAlgebra( expression=cbind(phiA/totalVg,phiC/totalVg, phiE/totalVg), name="stanVg") # SymmMatrix 'phiA' # [,1] [,2] [,3] # [1,] "VA_I" "Acov_IS" "Acov_IQ" # [2,] "Acov_IS" "VA_S" "Acov_SQ" / totalVg = total A, C and E variance/covariance in I,S & Q # [3,] "Acov_IQ" "Acov_SQ" "VC_Q" round(dcs_fit$algebras$stanVg$result,2) temp <- round(dcs_fit$algebras$stanVg$result,2) colnames(temp, do.NULL = FALSE) colnames(temp) <- c("A","A","A", "C","C","C", "E","E","E") rownames(temp) <- c("Intercetp","Linear slope","Quadratic") temp # A A A C C C E E E # Intercetp 0.52 0.55 0.88 0.22 0.19 -0.12 0.27 0.27 0.24 # Linear slope 0.55 0.54 0.66 0.19 0.20 0.08 0.27 0.27 0.26 # Quadratic 0.88 0.66 0.53 -0.12 0.08 0.21 0.24 0.26 0.26 # Fetch correlations for the A, C and E intercept, slope (linear) and quadratric latent factors dcs_fit$algebras$corrA dcs_fit$algebras$corrC dcs_fit$algebras$corrE corrA_isq<-dcs_fit$algebras$corrA$result colnames(corrA_isq, do.NULL = FALSE) colnames(corrA_isq) <- c("A_i","A_s","A_q") rownames(corrA_isq) <- c("A_i","A_s","A_q") round(corrA_isq,2) # AE Dual Change Score. Remove C phi (growth variance/covariance), remove C psi (auto-regression variance/covariance) AE_DCS <- mxModel( dcs_fit, name="AE_DCS" ) AE_DCS <- omxSetParameters( AE_DCS, label=phiC_labs, free=F, values=0) # Remove 'C' growth parameters (hint=phi) AE_DCS <- omxSetParameters( AE_DCS, label=psiC_labs, free=F, values=0) # Remove 'C' residuals (auto-regression) (hint=psi) AE_DCS <- omxSetParameters( AE_DCS, label=c("res_c"), free=F, values=0) # Remove 'C' residual errors on observed variables omxGetParameters(AE_DCS) set.seed(2016) summary(AE_DCS_fit <- mxTryHard(AE_DCS, extraTries=15, intervals = F, greenOK=T, checkHess=TRUE, fit2beat=Inf, exhaustive=T)) # Fetch ALL estimated model parameters round(AE_DCS_fit$output$estimate,4) # Fetch unstandardized & standardized variance components AE_DCS_fit$algebras$VC$result # Fetch just the standardized variance components round(AE_DCS_fit$algebras$VC$result[,16:30],2) # Fetch correlations for the A, C and E intercept, slope (linear) and quadratric latent factors ???$algebras$corrA ???$algebras$corrC ???$algebras$corrE # CE Dual Change Score. Remove A phi (growth variance/covariance), remove A psi (auto-regression variance/covariance) CE_DCS <- mxModel( dcs_fit, name="CE_DCS" ) CE_DCS <- omxSetParameters( CE_DCS, label=???, free=F, values=0) # Remove the 'A' latent growth process CE_DCS <- omxSetParameters( CE_DCS, label=???, free=F, values=0) # Remove the 'A' autoregression process CE_DCS <- omxSetParameters( CE_DCS, label=c("res_a"), free=F, values=0) # Remove the 'A' residual errors on observed variables omxGetParameters(CE_DCS) set.seed(2016) summary(CE_DCS_fit <- mxTryHard(CE_DCS, extraTries=15, intervals = F, greenOK=T, checkHess=TRUE, fit2beat=Inf, exhaustive=T)) # Fetch ALL estimated model parameters round(CE_DCS_fit$output$estimate,4) # E Dual Change Score. Remove A phi (growth variance/covariance), remove A psi (auto-regression variance/covariance) E_DCS <- mxModel( dcs_fit, name="E_DCS" ) E_DCS <- omxSetParameters( E_DCS, label=???, free=F, values=0) # Remove the 'A' latent growth process E_DCS <- omxSetParameters( E_DCS, label=???, free=F, values=0) # Remove the 'A' autoregression process E_DCS <- omxSetParameters( E_DCS, label=???, free=F, values=0) # Remove the 'C' latent growth process E_DCS <- omxSetParameters( E_DCS, label=???, free=F, values=0) # Remove the 'C' autoregression process E_DCS <- omxSetParameters( E_DCS, label=c("res_c"), free=F, values=0) E_DCS <- omxSetParameters( E_DCS, label=c("res_a"), free=F, values=0) omxGetParameters(E_DCS) set.seed(2016) summary(E_DCS_fit <- mxTryHard(E_DCS, extraTries=15, intervals = F, greenOK=T, checkHess=TRUE, fit2beat=Inf, exhaustive=T)) # Compare models subs <- c( AE_DCS_fit,CE_DCS_fit,E_DCS_fit ) comparisons<-mxCompare( dcs_fit, subs );comparisons # AE model with (i) an additive genetic latent growth process & (ii) a non-shared environmental autoregression 'simplex' process # A phi (growth variance/covariance) + E psi (auto-regression variance/covariance) # Remove all 'C' processes including residual errors Aphi_Epsi <- mxModel( dcs_fit, name="Aphi_Epsi" ) Aphi_Epsi <- omxSetParameters( Aphi_Epsi, label=???_labs, free=F, values=0) # Remove C latent growth process Aphi_Epsi <- omxSetParameters( Aphi_Epsi, label=???_labs, free=F, values=0) # Remove C autoregression process Aphi_Epsi <- omxSetParameters( Aphi_Epsi, label=c("res_c"), free=F, values=0) # Remove C residual errors Aphi_Epsi <- omxSetParameters( Aphi_Epsi, label=???_labs, free=F, values=0) # Remove A autoregression process Aphi_Epsi <- omxSetParameters( Aphi_Epsi, label=???_labs, free=F, values=0) # Remove E latent growth process Aphi_Epsi <- omxSetParameters( Aphi_Epsi, label=c("VE_i???"), free=T, values=.5) # Add first E innovation i.e. at Time 1 omxGetParameters(Aphi_Epsi) set.seed(2016) summary( Aphi_Epsi_fit <- mxTryHard( Aphi_Epsi, extraTries=15, intervals = F, greenOK=T, checkHess=TRUE, fit2beat=Inf, exhaustive=T)) # Compare models comparisons <- mxCompare( dcs_fit, Aphi_Epsi_fit );comparisons # umxCP(name="CP",selVars,MZdata,DZdata,nFac=1,freeLowerA=T,correlatedA=T,equateMeans=T,dzAr=0.5,dzCr=1,addStd=T)