Four-Factor Cholesky Decomposition

4 replies [Last post]
Pyang5's picture
Offline
Joined: 08/07/2012

Does the new OpenMX version allow me to model a Four-factor Cholesky Decomposition? I would like to use 2 continuous, 1 ordinal, and 1 dichotomous variable in the model? I was wondering if anyone might have code that is similar to post on the forum for 3 or 4 factor models? Below is my attempt but I get NaN for the standard error.

#vars v1 and v2 are continuous and v3 is ordinal (5 categories), but v4 is dichotomous
______________________________________________________________
nvar <- 4 # number of variables
tnvar <-8 # number of variables*max family size
nlower <-tnvar*(tnvar+1)/2 # number of free elements in a lower matrix tnvar*tnvar

Vars <- c("v1","v2","v3","v4")
selVars <- paste(Vars,c(rep(1,nvar),rep(2,nvar)),sep="")

# This yields a list of the 6 variables for analysis
#v11 v21 v31 v41 v12 v22 v32 v42

mzData <- as.matrix(subset(octotwin_data, zygosity==1, selVars))
dzData <- as.matrix(subset(octotwin_data, zygosity==2, selVars))

#Print Descriptive Statistics
#-------------------------------------------------------------------------
summary(mzData)
summary(dzData)

### Fit Multivariate ACE Model with RawData and Matrices Input ###

ACE_Cholesky_Model <- mxModel("ACE_Cholesky",
mxModel("ACE",
# Matrices a, c, and e to store a, c, and e path coefficients
mxMatrix( type="Lower", nrow=nvar, ncol=nvar, free=TRUE, values=.6, name="a" ),
mxMatrix( type="Lower", nrow=nvar, ncol=nvar, free=TRUE, values=.6, name="c" ),
mxMatrix( type="Lower", nrow=nvar, ncol=nvar, free=TRUE, values=.6, name="e" ),
# Matrices A, C, and E compute variance components
mxAlgebra( a %*% t(a), name="A" ),
mxAlgebra( c %*% t(c), name="C" ),
mxAlgebra( e %*% t(e), name="E" ),
# Algebra to compute total variances and standard deviations (diagonal only)
mxAlgebra( A+C+E, name="V" ),
mxMatrix( type="Iden", nrow=nvar, ncol=nvar, name="I"),
mxAlgebra( solve(sqrt(I*V)), name="iSD"),

# Matrix & Algebra for expected means vector
mxMatrix( type="Full", nrow=1, ncol=nvar, free=TRUE, values= 0, name="Mean" ),
mxAlgebra( cbind(Mean,Mean), name="expMean"),
# Algebra for expected variance/covariance matrix in MZ
mxAlgebra( rbind ( cbind(A+C+E , A+C),
cbind(A+C , A+C+E)), name="expCovMZ" ),
# Algebra for expected variance/covariance matrix in DZ
mxAlgebra( rbind( cbind(A+C+E , 0.5%x%A+C),
cbind(0.5%x%A+C , A+C+E)), name="expCovDZ" )
),
mxModel("MZ",
mxData( observed=mzData, type="raw" ),
mxFIMLObjective( covariance="ACE.expCovMZ", means="ACE.expMean", dimnames=selVars )
),
mxModel("DZ",
mxData( observed=dzData, type="raw" ),
mxFIMLObjective( covariance="ACE.expCovDZ", means="ACE.expMean", dimnames=selVars )
),
mxAlgebra( MZ.objective + DZ.objective, name="modelfit" ),
mxAlgebraObjective("modelfit")
)
ACE_Cholesky_Fit <- mxRun(ACE_Cholesky_Model)
ACE_Cholesky_Summ <- summary(ACE_Cholesky_Fit)
ACE_Cholesky_Summ

### Generate Multivariate Cholesky ACE Output ###

parameterSpecifications(ACE_Cholesky_Fit)
expectedMeansCovariances(ACE_Cholesky_Fit)
tableFitStatistics(ACE_Cholesky_Fit)

ACEpathMatrices <- c("ACE.a","ACE.c","ACE.e","ACE.iSD","ACE.iSD %*% ACE.a","ACE.iSD %*% ACE.c","ACE.iSD %*% ACE.e")
ACEpathLabels <- c("pathEst_a","pathEst_c","pathEst_e","sd","stPathEst_a","stPathEst_c","stPathEst_e")
formatOutputMatrices(ACE_Cholesky_Fit,ACEpathMatrices,ACEpathLabels,Vars,3)

ACEcovMatrices <- c("ACE.A","ACE.C","ACE.E","ACE.V","ACE.A/ACE.V","ACE.C/ACE.V","ACE.E/ACE.V")
ACEcovLabels <- c("covComp_A","covComp_C","covComp_E","Var","stCovComp_A","stCovComp_C","stCovComp_E")
formatOutputMatrices(ACE_Cholesky_Fit,ACEcovMatrices,ACEcovLabels,Vars,3)

# Genetic and Environmental Correlations
ACEcorMatrices <- c("solve(sqrt(ACE.I*ACE.A)) %*% ACE.A %*% solve(sqrt(ACE.I*ACE.A))",
"solve(sqrt(ACE.I*ACE.C)) %*% ACE.C %*% solve(sqrt(ACE.I*ACE.C))",
"solve(sqrt(ACE.I*ACE.E)) %*% ACE.E %*% solve(sqrt(ACE.I*ACE.E))",
"solve(sqrt(ACE.I*ACE.V)) %*% ACE.V %*% solve(sqrt(ACE.I*ACE.V))")
ACEcorLabels <- c("corA","corC","corE","cor")
formatOutputMatrices(ACE_Cholesky_Fit,ACEcorMatrices,ACEcorLabels, Vars, 4)

Ryne's picture
Offline
Joined: 07/31/2009
This looks pretty good. NaN

This looks pretty good. NaN standard error usually means that the Hessian isn't invertible. Did you get any warnings when you ran the model?

Pyang5's picture
Offline
Joined: 08/07/2012
Warning Message

Warning message:
In model 'ACE_Cholesky' NPSOL returned a non-zero status code 1. The final iterate satisfies the optimality conditions to the accuracy requested, but the sequence of iterates has not yet converged. NPSOL was terminated because no further improvement could be made in the merit function (Mx status GREEN).

I'm posting the output shortly after this one.

Pyang5's picture
Offline
Joined: 08/07/2012
Output

>
> nvar <- 4 # number of variables
> tnvar <-8 # number of variables*max family size
> nlower <-tnvar*(tnvar+1)/2 # number of free elements in a lower matrix tnvar*tnvar
>
> Vars <- c("v1","v2","v3","v4")
> selVars <- paste(Vars,c(rep(1,nvar),rep(2,nvar)),sep="")
>
> # This yields a list of the 6 variables for analysis
> #v11v21v31 v41 v12v22v32 v42
>
>
> mzData <- as.matrix(subset(octotwin_data, zygosity==1, selVars))
> dzData <- as.matrix(subset(octotwin_data, zygosity==2, selVars))
>
> #Print Descriptive Statistics
> #-------------------------------------------------------------------------
> summary(mzData)
v11 v21 v31 v41
Min. :-2.2443 Min. :-2.302 Min. :-2.2443 Min. :0.0000
1st Qu.:-0.4045 1st Qu.:-2.266 1st Qu.:-0.9063 1st Qu.:0.0000
Median : 0.4318 Median :-2.233 Median :-0.2372 Median :0.0000
Mean : 0.1889 Mean :-2.221 Mean :-0.2562 Mean :0.1812
3rd Qu.: 0.9336 3rd Qu.:-2.167 3rd Qu.: 0.5154 3rd Qu.:0.0000
Max. : 1.4353 Max. :-2.076 Max. : 1.7698 Max. :1.0000
NA's :24.0000 NA's :21.000 NA's :52.0000
v12 v22 v32 v42
Min. :-2.2443 Min. :-2.302 Min. :-2.24427 Min. :0.0000
1st Qu.:-0.5718 1st Qu.:-2.275 1st Qu.:-0.82263 1st Qu.:0.0000
Median : 0.4318 Median :-2.236 Median :-0.07000 Median :0.0000
Mean : 0.1134 Mean :-2.223 Mean :-0.08286 Mean :0.1477
3rd Qu.: 0.9336 3rd Qu.:-2.173 3rd Qu.: 0.68264 3rd Qu.:0.0000
Max. : 1.4353 Max. :-2.076 Max. : 2.68966 Max. :1.0000
NA's :19.0000 NA's :21.000 NA's :45.00000
> summary(dzData)
v11 v21 v31 v41
Min. :-2.24433 Min. :-2.302 Min. :-2.2443 Min. :0.0000
1st Qu.:-0.82264 1st Qu.:-2.302 1st Qu.:-1.1571 1st Qu.:0.0000
Median : 0.01365 Median :-2.236 Median :-0.3209 Median :0.0000
Mean :-0.10698 Mean :-2.231 Mean :-0.4165 Mean :0.2129
3rd Qu.: 0.68268 3rd Qu.:-2.192 3rd Qu.: 0.2645 3rd Qu.:0.0000
Max. : 1.43534 Max. :-2.052 Max. : 2.1043 Max. :1.0000
NA's :37.00000 NA's :32.000 NA's :70.0000
v12 v22 v32 v42
Min. :-2.2443 Min. :-2.302 Min. :-2.2443 Min. :0.0000
1st Qu.:-0.9899 1st Qu.:-2.266 1st Qu.:-1.0317 1st Qu.:0.0000
Median : 0.1391 Median :-2.236 Median :-0.4881 Median :0.0000
Mean :-0.1278 Mean :-2.230 Mean :-0.4008 Mean :0.1733
3rd Qu.: 0.6827 3rd Qu.:-2.192 3rd Qu.: 0.2645 3rd Qu.:0.0000
Max. : 1.4353 Max. :-2.099 Max. : 1.6862 Max. :1.0000
NA's :40.0000 NA's :33.000 NA's :67.0000
>
> ### Fit Multivariate ACE Model with RawData and Matrices Input ###
>
> ACE_Cholesky_Model <- mxModel("ACE_Cholesky",
+ mxModel("ACE",
+ # Matrices a, c, and e to store a, c, and e path coefficients
+ mxMatrix( type="Lower", nrow=nvar, ncol=nvar, free=TRUE, values=.6, name="a" ),
+ mxMatrix( type="Lower", nrow=nvar, ncol=nvar, free=TRUE, values=.6, name="c" ),
+ mxMatrix( type="Lower", nrow=nvar, ncol=nvar, free=TRUE, values=.6, name="e" ),
+ # Matrices A, C, and E compute variance components
+ mxAlgebra( a %*% t(a), name="A" ),
+ mxAlgebra( c %*% t(c), name="C" ),
+ mxAlgebra( e %*% t(e), name="E" ),
+ # Algebra to compute total variances and standard deviations (diagonal only)
+ mxAlgebra( A+C+E, name="V" ),
+ mxMatrix( type="Iden", nrow=nvar, ncol=nvar, name="I"),
+ mxAlgebra( solve(sqrt(I*V)), name="iSD"),
+
+ # Matrix & Algebra for expected means vector
+ mxMatrix( type="Full", nrow=1, ncol=nvar, free=TRUE, values= 0, name="Mean" ),
+ mxAlgebra( cbind(Mean,Mean), name="expMean"),
+ # Algebra for expected variance/covariance matrix in MZ
+ mxAlgebra( rbind ( cbind(A+C+E , A+C),
+ cbind(A+C , A+C+E)), name="expCovMZ" ),
+ # Algebra for expected variance/covariance matrix in DZ
+ mxAlgebra( rbind( cbind(A+C+E , 0.5%x%A+C),
+ cbind(0.5%x%A+C , A+C+E)), name="expCovDZ" )
+ ),
+ mxModel("MZ",
+ mxData( observed=mzData, type="raw" ),
+ mxFIMLObjective( covariance="ACE.expCovMZ", means="ACE.expMean", dimnames=selVars )
+ ),
+ mxModel("DZ",
+ mxData( observed=dzData, type="raw" ),
+ mxFIMLObjective( covariance="ACE.expCovDZ", means="ACE.expMean", dimnames=selVars )
+ ),
+ mxAlgebra( MZ.objective + DZ.objective, name="modelfit" ),
+ mxAlgebraObjective("modelfit")
+ )
> ACE_Cholesky_Fit <- mxRun(ACE_Cholesky_Model)
Running ACE_Cholesky
Warning message:
In model 'ACE_Cholesky' NPSOL returned a non-zero status code 1. The final iterate satisfies the optimality conditions to the accuracy requested, but the sequence of iterates has not yet converged. NPSOL was terminated because no further improvement could be made in the merit function (Mx status GREEN).
> ACE_Cholesky_Summ <- summary(ACE_Cholesky_Fit)
> ACE_Cholesky_Summ
data:
$MZ.data
v11 v21 v31 v41 v12 v22 v32 v42
Min. :-2.2443 Min. :-2.302 Min. :-2.2443 Min. :0.0000 Min. :-2.2443 Min. :-2.302 Min. :-2.24427 Min. :0.0000
1st Qu.:-0.4045 1st Qu.:-2.266 1st Qu.:-0.9063 1st Qu.:0.0000 1st Qu.:-0.5718 1st Qu.:-2.275 1st Qu.:-0.82263 1st Qu.:0.0000
Median : 0.4318 Median :-2.233 Median :-0.2372 Median :0.0000 Median : 0.4318 Median :-2.236 Median :-0.07000 Median :0.0000
Mean : 0.1889 Mean :-2.221 Mean :-0.2562 Mean :0.1812 Mean : 0.1134 Mean :-2.223 Mean :-0.08286 Mean :0.1477
3rd Qu.: 0.9336 3rd Qu.:-2.167 3rd Qu.: 0.5154 3rd Qu.:0.0000 3rd Qu.: 0.9336 3rd Qu.:-2.173 3rd Qu.: 0.68264 3rd Qu.:0.0000
Max. : 1.4353 Max. :-2.076 Max. : 1.7698 Max. :1.0000 Max. : 1.4353 Max. :-2.076 Max. : 2.68966 Max. :1.0000
NA's :24.0000 NA's :21.000 NA's :52.0000 NA's :19.0000 NA's :21.000 NA's :45.00000

$DZ.data
v11 v21 v31 v41 v12 v22 v32 v42
Min. :-2.24433 Min. :-2.302 Min. :-2.2443 Min. :0.0000 Min. :-2.2443 Min. :-2.302 Min. :-2.2443 Min. :0.0000
1st Qu.:-0.82264 1st Qu.:-2.302 1st Qu.:-1.1571 1st Qu.:0.0000 1st Qu.:-0.9899 1st Qu.:-2.266 1st Qu.:-1.0317 1st Qu.:0.0000
Median : 0.01365 Median :-2.236 Median :-0.3209 Median :0.0000 Median : 0.1391 Median :-2.236 Median :-0.4881 Median :0.0000
Mean :-0.10698 Mean :-2.231 Mean :-0.4165 Mean :0.2129 Mean :-0.1278 Mean :-2.230 Mean :-0.4008 Mean :0.1733
3rd Qu.: 0.68268 3rd Qu.:-2.192 3rd Qu.: 0.2645 3rd Qu.:0.0000 3rd Qu.: 0.6827 3rd Qu.:-2.192 3rd Qu.: 0.2645 3rd Qu.:0.0000
Max. : 1.43534 Max. :-2.052 Max. : 2.1043 Max. :1.0000 Max. : 1.4353 Max. :-2.099 Max. : 1.6862 Max. :1.0000
NA's :37.00000 NA's :32.000 NA's :70.0000 NA's :40.0000 NA's :33.000 NA's :67.0000

The final iterate satisfies the optimality conditions to the accuracy requested, but the sequence of iterates has not yet converged. NPSOL was terminated because no further improvement could be made in the merit function (Mx status GREEN).

free parameters:
name matrix row col Estimate Std.Error
1 ACE.a 1 1 7.915107e-01 0.085341221
2 ACE.a 2 1 2.890662e-02 0.006778621
3 ACE.a 3 1 6.396641e-01 0.125868806
4 ACE.a 4 1 -7.514610e-02 0.046775965
5 ACE.a 2 2 -2.810453e-02 0.008087888
6 ACE.a 3 2 -1.649436e-01 0.203130382
7 ACE.a 4 2 3.368846e-02 0.069319677
8 ACE.a 3 3 3.455544e-01 0.190306485
9 ACE.a 4 3 -1.214520e-01 0.142814567
10 ACE.a 4 4 -1.427293e-01 0.151420792
11 ACE.c 1 1 -4.091653e-01 0.150965339
12 ACE.c 2 1 2.152481e-03 0.010640423
13 ACE.c 3 1 -1.621195e-02 0.260958391
14 ACE.c 4 1 8.829372e-02 0.062912267
15 ACE.c 2 2 -2.304735e-04 0.016423243
16 ACE.c 3 2 2.115071e-01 0.317464791
17 ACE.c 4 2 4.381787e-02 0.088886860
18 ACE.c 3 3 2.290132e-05 NaN
19 ACE.c 4 3 4.619492e-06 NaN
20 ACE.c 4 4 -8.934258e-09 0.091540369
21 ACE.e 1 1 4.953322e-01 0.033523570
22 ACE.e 2 1 1.000396e-02 0.004415674
23 ACE.e 3 1 1.882684e-01 0.061370898
24 ACE.e 4 1 -5.383125e-02 0.028801489
25 ACE.e 2 2 4.323006e-02 0.002409342
26 ACE.e 3 2 3.853718e-02 0.057773223
27 ACE.e 4 2 2.086212e-02 0.026412341
28 ACE.e 3 3 5.826048e-01 0.043130385
29 ACE.e 4 3 -2.634869e-03 0.030165499
30 ACE.e 4 4 -3.053772e-01 0.017078148
31 ACE.Mean 1 1 -1.017620e-01 0.050736267
32 ACE.Mean 1 2 -2.231421e+00 0.002790341
33 ACE.Mean 1 3 -4.447926e-01 0.050626622
34 ACE.Mean 1 4 1.823697e-01 0.016334605

observed statistics: 2347
estimated parameters: 34
degrees of freedom: 2313
-2 log likelihood: 1355.343
saturated -2 log likelihood: NA
number of observations: 351
chi-square: NA
p: NA
AIC (Mx): -3270.657
BIC (Mx): -6100.328
adjusted BIC:
RMSEA: NA
timestamp: 2012-07-20 15:20:16
frontend time: 3.152 secs
backend time: 12.452 secs
independent submodels time: 0 secs
wall clock time: 15.604 secs
cpu time: 15.604 secs
openmx version number: 1.0.6-1581

>
>
>
> ### Generate Multivariate Cholesky ACE Output ###
>
> parameterSpecifications(ACE_Cholesky_Fit)
model:ACE, matrix:a
[,1] [,2] [,3] [,4]
[1,] [NA] 0 0 0
[2,] [NA] [NA] 0 0
[3,] [NA] [NA] [NA] 0
[4,] [NA] [NA] [NA] [NA]

model:ACE, matrix:c
[,1] [,2] [,3] [,4]
[1,] [NA] 0 0 0
[2,] [NA] [NA] 0 0
[3,] [NA] [NA] [NA] 0
[4,] [NA] [NA] [NA] [NA]

model:ACE, matrix:e
[,1] [,2] [,3] [,4]
[1,] [NA] 0 0 0
[2,] [NA] [NA] 0 0
[3,] [NA] [NA] [NA] 0
[4,] [NA] [NA] [NA] [NA]

model:ACE, matrix:Mean
[,1] [,2] [,3] [,4]
[1,] [NA] [NA] [NA] [NA]

> expectedMeansCovariances(ACE_Cholesky_Fit)
model:MZ, covariance:ACE.expCovMZ
v11 v21 v31 v41 v12 v22 v32 v42
v11 1.03925954 0.026954461 0.60618982 -0.122270026 0.79390555 0.021999179 0.51293440 -0.095605674
v21 0.02695446 0.003599061 0.02659194 -0.002575721 0.02199918 0.001630144 0.02304255 -0.002939066
v31 0.60618982 0.026591943 0.97714101 -0.098622716 0.51293440 0.023042549 0.60078254 -0.087756870
v41 -0.12227003 -0.002575721 -0.09862272 0.148215079 -0.09560567 -0.002939066 -0.08775687 0.051619898
v12 0.79390555 0.021999179 0.51293440 -0.095605674 1.03925954 0.026954461 0.60618982 -0.122270026
v22 0.02199918 0.001630144 0.02304255 -0.002939066 0.02695446 0.003599061 0.02659194 -0.002575721
v32 0.51293440 0.023042549 0.60078254 -0.087756870 0.60618982 0.026591943 0.97714101 -0.098622716
v42 -0.09560567 -0.002939066 -0.08775687 0.051619898 -0.12227003 -0.002575721 -0.09862272 0.148215079

model:MZ, means:ACE.expMean
v11 v21 v31 v41 v12 v22 v32 v42
[1,] -0.1017620 -2.231421 -0.4447926 0.1823697 -0.1017620 -2.231421 -0.4447926 0.1823697

model:DZ, covariance:ACE.expCovDZ
v11 v21 v31 v41 v12 v22 v32 v42
v11 1.03925954 0.026954461 0.60618982 -0.122270026 0.48066091 0.010559229 0.25978388 -0.065866202
v21 0.02695446 0.003599061 0.02659194 -0.002575721 0.01055923 0.000817415 0.01147945 -0.001379557
v31 0.60618982 0.026591943 0.97714101 -0.098622716 0.25978388 0.011479453 0.32289032 -0.039960245
v41 -0.12227003 -0.002575721 -0.09862272 0.148215079 -0.06586620 -0.001379557 -0.03996024 0.030667842
v12 0.48066091 0.010559229 0.25978388 -0.065866202 1.03925954 0.026954461 0.60618982 -0.122270026
v22 0.01055923 0.000817415 0.01147945 -0.001379557 0.02695446 0.003599061 0.02659194 -0.002575721
v32 0.25978388 0.011479453 0.32289032 -0.039960245 0.60618982 0.026591943 0.97714101 -0.098622716
v42 -0.06586620 -0.001379557 -0.03996024 0.030667842 -0.12227003 -0.002575721 -0.09862272 0.148215079

model:DZ, means:ACE.expMean
v11 v21 v31 v41 v12 v22 v32 v42
[1,] -0.1017620 -2.231421 -0.4447926 0.1823697 -0.1017620 -2.231421 -0.4447926 0.1823697

> tableFitStatistics(ACE_Cholesky_Fit)
Name ep -2LL df AIC
Model 1 : ACE_Cholesky 34 1355.34 2313 -3270.66

>
> ACEpathMatrices <- c("ACE.a","ACE.c","ACE.e","ACE.iSD","ACE.iSD %*% ACE.a","ACE.iSD %*% ACE.c","ACE.iSD %*% ACE.e")
> ACEpathLabels <- c("pathEst_a","pathEst_c","pathEst_e","sd","stPathEst_a","stPathEst_c","stPathEst_e")
> formatOutputMatrices(ACE_Cholesky_Fit,ACEpathMatrices,ACEpathLabels,Vars,3)
[1] "Matrix ACE.a"
pathEst_a1 pathEst_a2 pathEst_a3 pathEst_a4
v1 0.792 0.000 0.000 0.000
v2 0.029 -0.028 0.000 0.000
v3 0.640 -0.165 0.346 0.000
v4 -0.075 0.034 -0.121 -0.143

[1] "Matrix ACE.c"
pathEst_c1 pathEst_c2 pathEst_c3 pathEst_c4
v1 -0.409 0.000 0.000 0.000
v2 0.002 0.000 0.000 0.000
v3 -0.016 0.212 0.000 0.000
v4 0.088 0.044 0.000 0.000

[1] "Matrix ACE.e"
pathEst_e1 pathEst_e2 pathEst_e3 pathEst_e4
v1 0.495 0.000 0.000 0.000
v2 0.010 0.043 0.000 0.000
v3 0.188 0.039 0.583 0.000
v4 -0.054 0.021 -0.003 -0.305

[1] "Matrix ACE.iSD"
sd1 sd2 sd3 sd4
v1 0.981 0.000 0.000 0.000
v2 0.000 16.669 0.000 0.000
v3 0.000 0.000 1.012 0.000
v4 0.000 0.000 0.000 2.597

[1] "Matrix ACE.iSD %*% ACE.a"
stPathEst_a1 stPathEst_a2 stPathEst_a3 stPathEst_a4
v1 0.776 0.000 0.000 0.000
v2 0.482 -0.468 0.000 0.000
v3 0.647 -0.167 0.350 0.000
v4 -0.195 0.088 -0.315 -0.371

[1] "Matrix ACE.iSD %*% ACE.c"
stPathEst_c1 stPathEst_c2 stPathEst_c3 stPathEst_c4
v1 -0.401 0.000 0.000 0.000
v2 0.036 -0.004 0.000 0.000
v3 -0.016 0.214 0.000 0.000
v4 0.229 0.114 0.000 0.000

[1] "Matrix ACE.iSD %*% ACE.e"
stPathEst_e1 stPathEst_e2 stPathEst_e3 stPathEst_e4
v1 0.486 0.000 0.000 0.000
v2 0.167 0.721 0.000 0.000
v3 0.190 0.039 0.589 0.000
v4 -0.140 0.054 -0.007 -0.793

>
> ACEcovMatrices <- c("ACE.A","ACE.C","ACE.E","ACE.V","ACE.A/ACE.V","ACE.C/ACE.V","ACE.E/ACE.V")
> ACEcovLabels <- c("covComp_A","covComp_C","covComp_E","Var","stCovComp_A","stCovComp_C","stCovComp_E")
> formatOutputMatrices(ACE_Cholesky_Fit,ACEcovMatrices,ACEcovLabels,Vars,3)
[1] "Matrix ACE.A"
covComp_A1 covComp_A2 covComp_A3 covComp_A4
v1 0.626 0.023 0.506 -0.059
v2 0.023 0.002 0.023 -0.003
v3 0.506 0.023 0.556 -0.096
v4 -0.059 -0.003 -0.096 0.042

[1] "Matrix ACE.C"
covComp_C1 covComp_C2 covComp_C3 covComp_C4
v1 0.167 -0.001 0.007 -0.036
v2 -0.001 0.000 0.000 0.000
v3 0.007 0.000 0.045 0.008
v4 -0.036 0.000 0.008 0.010

[1] "Matrix ACE.E"
covComp_E1 covComp_E2 covComp_E3 covComp_E4
v1 0.245 0.005 0.093 -0.027
v2 0.005 0.002 0.004 0.000
v3 0.093 0.004 0.376 -0.011
v4 -0.027 0.000 -0.011 0.097

[1] "Matrix ACE.V"
Var1 Var2 Var3 Var4
v1 1.039 0.027 0.606 -0.122
v2 0.027 0.004 0.027 -0.003
v3 0.606 0.027 0.977 -0.099
v4 -0.122 -0.003 -0.099 0.148

[1] "Matrix ACE.A/ACE.V"
stCovComp_A1 stCovComp_A2 stCovComp_A3 stCovComp_A4
v1 0.603 0.849 0.835 0.486
v2 0.849 0.452 0.870 1.211
v3 0.835 0.870 0.569 0.969
v4 0.486 1.211 0.969 0.283

[1] "Matrix ACE.C/ACE.V"
stCovComp_C1 stCovComp_C2 stCovComp_C3 stCovComp_C4
v1 0.161 -0.033 0.011 0.295
v2 -0.033 0.001 -0.003 -0.070
v3 0.011 -0.003 0.046 -0.079
v4 0.295 -0.070 -0.079 0.066

[1] "Matrix ACE.E/ACE.V"
stCovComp_E1 stCovComp_E2 stCovComp_E3 stCovComp_E4
v1 0.236 0.184 0.154 0.218
v2 0.184 0.547 0.133 -0.141
v3 0.154 0.133 0.385 0.110
v4 0.218 -0.141 0.110 0.652

>
> # Genetic and Environmental Correlations
> ACEcorMatrices <- c("solve(sqrt(ACE.I*ACE.A)) %*% ACE.A %*% solve(sqrt(ACE.I*ACE.A))",
+ "solve(sqrt(ACE.I*ACE.C)) %*% ACE.C %*% solve(sqrt(ACE.I*ACE.C))",
+ "solve(sqrt(ACE.I*ACE.E)) %*% ACE.E %*% solve(sqrt(ACE.I*ACE.E))",
+ "solve(sqrt(ACE.I*ACE.V)) %*% ACE.V %*% solve(sqrt(ACE.I*ACE.V))")
> ACEcorLabels <- c("corA","corC","corE","cor")
> formatOutputMatrices(ACE_Cholesky_Fit,ACEcorMatrices,ACEcorLabels, Vars, 4)
[1] "Matrix solve(sqrt(ACE.I*ACE.A)) %*% ACE.A %*% solve(sqrt(ACE.I*ACE.A))"
corA1 corA2 corA3 corA4
v1 1.0000 0.7170 0.8580 -0.3671
v2 0.7170 1.0000 0.7694 -0.3779
v3 0.8580 0.7694 1.0000 -0.6264
v4 -0.3671 -0.3779 -0.6264 1.0000

[1] "Matrix solve(sqrt(ACE.I*ACE.C)) %*% ACE.C %*% solve(sqrt(ACE.I*ACE.C))"
corC1 corC2 corC3 corC4
v1 1.0000 -0.9943 0.0764 -0.8958
v2 -0.9943 1.0000 -0.1821 0.8433
v3 0.0764 -0.1821 1.0000 0.3748
v4 -0.8958 0.8433 0.3748 1.0000

[1] "Matrix solve(sqrt(ACE.I*ACE.E)) %*% ACE.E %*% solve(sqrt(ACE.I*ACE.E))"
corE1 corE2 corE3 corE4
v1 1.0000 0.2255 0.3069 -0.1732
v2 0.2255 1.0000 0.1304 0.0263
v3 0.3069 0.1304 1.0000 -0.0570
v4 -0.1732 0.0263 -0.0570 1.0000

[1] "Matrix solve(sqrt(ACE.I*ACE.V)) %*% ACE.V %*% solve(sqrt(ACE.I*ACE.V))"
cor1 cor2 cor3 cor4
v1 1.0000 0.4407 0.6015 -0.3115
v2 0.4407 1.0000 0.4484 -0.1115
v3 0.6015 0.4484 1.0000 -0.2592
v4 -0.3115 -0.1115 -0.2592 1.0000

>

Ryne's picture
Offline
Joined: 07/31/2009
Huh. I don't see any problems

Huh. I don't see any problems beyond the NaN standard error. If it is caused by a negative diagonal element in the inverse Hessian, then that might mean that your model hasn't properly converged. As this model runs very quickly, try a few different sets of starting values, especially for the c matrix. Hopefully you'll get a different solution with a better -2LL, no status 1 code and no NaN standard error.

If you always get the same answer regardless of starting values (-2LL=1355.343), then we'll have to look a little closer as to why those last C matrix parameters are so problematic.