# Latent Growth Model and distal outcome

Posted on

massimiliano.orri
Joined: 06/11/2018

Forums

Hi,

I would like to fit the twin model in the enclosed figure, but I am not sure if this is feasible.

Can someone help me understanding if the model is identified?

In case, are there any exemples of R code for this model?

Thanks a lot

M

I would like to fit the twin model in the enclosed figure, but I am not sure if this is feasible.

Can someone help me understanding if the model is identified?

In case, are there any exemples of R code for this model?

Thanks a lot

M

## seems OK, but no example script

Log in or register to post comments

In reply to seems OK, but no example script by AdminRobK

## Thank you Rob for your reply.

I tried to write the OpenMx syntax for this model. Can I ask you if this seems ok for you?

I also have 2 problems:

1) I found it hard to specify that my distal outcome is binary, so I wrote the script as it was continuous. How can I specify this?

2) My code doesn't work, and I obtain the following Error:

Error: The following error occurred while evaluating the subexpression 'MZ.gamma %&% MZ.phiA + MZ.psiA' during the evaluation of 'MZ.A' in model 'my_model' : Matrices must have same dimensions in as(e1, "dMatrix") + as(e2, "dMatrix")

Here my script:

#Read data

nvars <- 6 # Number of manifest/observed variables

ntvars <- nvars*2 #

nGrow <- 3 # Number of growth parameters: intercept, slope, and F for the distal outcome

selVars <- c("t1_1", "t1_2", "t1_3", "t1_4", "t1_5", "t1_6", # note: t1_6 and t2_6 are the distal oucome

"t2_1", "t2_2", "t2_3", "t2_4", "t2_5", "t2_6")

# ϕ phi

# covariance between the exogenous variables i.e. between Intercept, Slope, & F factors

phiA_labs <- c("VA_Int","Acov_IntS","Acov_IntF","VA_S","Acov_SF","VA_F")

phiC_labs <- c("VC_Int","Ccov_IntS","Ccov_IntF","VC_S","Ccov_SF","VC_F")

phiE_labs <- c("VE_Int","Ecov_IntS","Ecov_IntF","VE_S","Ecov_SF","VE_F")

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

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, & F 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 from Intercept, Slope, & F to observed scores

# 6 (time points) x 3 factors (I, S, F) matrix

Glabs <- c(paste("i",1:nvars,sep=""), paste("s",1:nvars,sep=""), paste("q",1:nvars,sep=""))

gamma <- mxMatrix( type="Full", nrow=nvars, ncol=nGrow,

free=F, labels=Glabs, name="gamma",

values=c(1,1,1,1,1,0, 0,1,2,3,4,0, 0,0,0,0,0,1) )

# Ψ psi

# Innovations

psiA_labs <- c("VA_i11","VA_i22","VA_i33","VA_i44","VA_i55","VA_i66")

psiC_labs <- c("VC_i11","VC_i22","VC_i33","VC_i44","VC_i55","VC_i66")

psiE_labs <- c("VE_i11","VE_i22","VE_i33","VE_i44","VE_i55","VE_i66")

psiA <- mxMatrix(type = "Diag", nrow = nvars, ncol = nvars, free = c(F,T,T,T,T,F), labels = psiA_labs, name = "psiA",

values = c( 0,0.5,0.5,0.5,0.5,0))

psiC <- mxMatrix(type = "Diag", nrow = nvars,ncol = nvars, free = c(F,T,T,T,T,F), labels = psiC_labs, name = "psiC",

values = c(0,0.5,0.5,0.5,0.5,0))

psiE <- mxMatrix(type = "Diag", nrow = nvars, ncol = nvars, free = c(F,T,T,T,T,F), labels = psiE_labs, name = "psiE",

values = c(0,0.5,0.5,0.5,0.5,0))

# Expected variance/covariance matrix

A <- mxAlgebra(gamma %&% phiA + psiA, name = "A")

C <- mxAlgebra(gamma %&% phiC + psiC, name = "C")

E <- mxAlgebra(gamma %&% phiE + psiE, 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" )

stanVg <- mxAlgebra( expression=cbind(phiA/totalVg,phiC/totalVg, phiE/totalVg), name="stanVg")

# Expected means

ManMean <- mxMatrix(type="Full", nrow = 1, ncol = ntvars, free = T,values = .1, name = "ManMean")

# Data Objects for Multiple Groups

dataMZ <- mxData( observed=MZdata, type="raw" )

dataDZ <- mxData( observed=DZdata, type="raw" )

# Objects for Multiple Groups

objMZ <- mxExpectationNormal( covariance="expCovMZ", means="ManMean", dimnames=selVars )

objDZ <- mxExpectationNormal( covariance="expCovDZ", means="ManMean", 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","Time6"),colVC))

# Final steps

pars <- list( phiA,phiC,phiE,corrA,corrC,corrE,gamma,psiA,psiC,psiE,

A,C,E,V,VC,totalVg,stanVg,ManMean)

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

my_model <- mxModel( "my_model", pars, modelMZ, modelDZ, obj)

# Fit model

fit <- mxRun( my_model, intervals=F )

summary(my_fit)

Thank you very much!

M

Log in or register to post comments

In reply to Thank you Rob for your reply. by massimiliano.orri

## some remarks

That's very strange, because I don't see from your syntax why that error would be occurring. What is your output from `mxVersion()`? Maybe you need to encase

`gamma %&% phiA`

in parentheses?I'll assume you want to treat it as a threshold variable. There are several changes you'll need to make to your script. First, run the distal outcome through

`mxFactor()`

. Then, you'll need to make an MxMatrix containing the thresholds for the binary variable. That matrix of thresholds needs to go into the MZ and DZ models, and you'll need to provide its name as the value of argument `thresholds` to `mxExpectationNormal()`. Either the intercept or the threshold for the binary variable must be fixed, to identify the location of the latent trait presumed to underlie the observed dichotomy. The scale of the latent trait also needs to be identified. Numerically speaking, the easiest way to do that is to fix its nonshared environmental variance to 1.0. You could instead use an MxConstraint to fix its total variance to 1.0, which provides results that are easier to interpret, but makes optimization more difficult.One comment I have about your script is that you should model the phenotypic mean for the five timepoints in terms of the intercept and slope. That is, there should be a mean intercept and a mean slope, and the model-expected phenotypic mean at a given timepoint would be the mean intercept plus the product of the mean slope and the timepoint's loading on the latent slope.

Log in or register to post comments

In reply to some remarks by AdminRobK

## Thank you again for your

For the Error, I solved the problem updating openMx.

That's exact, I want to treat the binary distal outcome as a threshold variable. I'll try to follow your indications.

For your comment (modeling the phenotypic mean for the five timepoints) I see what you are suggesting, but unfortunately I could not implement that successfully...can you please help me understand how can I modify my code to model these phenotypic means ?

Relatedly, I used *mxCheckIdentification()* and I obtained that the model is not locally identified, and that the non-identified parameters are all the means of the manifest variables (ManMean in my script). Is it because I didn't model the phenotypic mean as you suggested?

I also would like to show you what I obtained when I ran the model (using *mxTryHard*, as not possible otherwise):

*Not all eigenvalues of Hessian are greater than 0*

What does this message mean? Is there something I can do about it ?

When I tried to estimate the following AE model instead of ACE to simplify (AE is what I expect from prior analysis), I obtained the same warning and bizarre results

Log in or register to post comments

In reply to Thank you again for your by massimiliano.orri

## phenotypic means

OK, good!

OK. You'll probably find it helpful to look up an example script that involves a threshold variable in a twin analysis.

I don't know. Perhaps.

That all sounds like it could be symptomatic of unidentification. "Not all eigenvalues of Hessian are greater than 0" suggests that the optimizer didn't find a local minimum of the objective function, and is fairly common when the model is unidentified.

Sure...and I'll also try to sort out some other potential issues for you as well. First, redefine `gamma` as:

gamma <- mxMatrix( type="Full", nrow=nvars, ncol=nGrow, free=F, labels=Glabs, name="gamma",

values=c(1,1,1,1,1,0, -2,-1,0,1,2,0, 0,0,0,0,0,1) )

This way, the loadings on the latent intercept will be orthogonal to the loadings onto the latent slope. The orthogonality will improve the numerical stability of the model and make optimization easier. Naturally, it will also change the interpretation of the mean intercept, which will now be the model-expected mean phenotype score at the third, not first, timepoint. With that in mind, define these new objects:

iLoadings <- mxMatrix(type="Full",nrow=1,ncol=6,free=F,values=c(1,1,1,1,1,0),name="iLoads")

sLoadings <- mxMatrix(type="Full",nrow=1,ncol=6,free=F,values=c(-2,-1,0,1,2,0),name="sLoads")

#You could get good empirical start values for the next two free parameters

#from lm() in base R, or from lmer() in package 'lme4':

meanIntrcpt <- mxMatrix(type="Full",nrow=1,ncol=1,free=T,values=0.1,labels="muIntrcpt",name="MuIntrcpt")

meanSlope <- mxMatrix(type="Full",nrow=1,ncol=1,free=T,values=0.1,labels="muSlope",name="MuSlope")

manMean <- mxAlgebra(

cbind(MuIntrcpt*iLoads + MuSlope*sLoads, MuIntrcpt*iLoads + MuSlope*sLoads),

name="ManMean"

)

Note that this syntax fixes the mean of the distal outcome to zero, meaning that its threshold (when you get around to that) will need to be a free parameter. You will quite likely also want to adjust the phenotypic mean for age and sex. Now suppose the age of twin #1, sex of twin #1, age of twin #2, and sex of twin #2 have column names in your dataset "age1", "sex1", "age2", "sex2" (and that the sex variables are dummy codes). Then, you would define these additional objects:

defAge <- mxMatrix(type="Full",nrow=1,ncol=2,free=F,labels=c("data.age1","data.age2"),name="xAge")

defSex <- mxMatrix(type="Full",nrow=1,ncol=2,free=F,labels=c("data.sex1","data.sex2"),name="xSex")

#In the below, the "t" suffix refers to the time-varying trait, and the "d" suffix, to the distal outcome:

betaAge <- mxMatrix(

type="Full",nrow=1,ncol=6,free=T,values=0.1,

labels=c(rep("bAget"),"bAged"),name="BAge")

betaSex <- mxMatrix(

type="Full",nrow=1,ncol=6,free=T,values=0.1,

labels=c(rep("bSext"),"bSexd"),name="BSex")

This syntax models the effects of age and sex on the time-varying trait as constant across time, which I think is reasonable. You would then redefine the phenotypic mean vector as:

manMean <- mxAlgebra(

cbind(MuIntrcpt*iLoads + MuSlope*sLoads, MuIntrcpt*iLoads + MuSlope*sLoads) +

(xAge%x%BAge) + (xSex%x%BSex)

name="ManMean"

)

Note that this syntax adjust the distal outcome for the effects of age and sex, but fixes its regression intercept at zero.

Log in or register to post comments

In reply to phenotypic means by AdminRobK

## I really appreciate your help

Following your indication I changed my model, and now it has a binary distal outcome specified as threshold (hope I did that in the right way...)

It also has an adjustment for sex (not for age since the twins have the same age).

This is the model:

#Read data

nvars <- 6 # Number of manifest/observed variables

ntvars <- nvars*2 #

nGrow <- 3 # Number of growth parameters: intercept, slope, and F for the distal outcome

selVars <- c("t1_1", "t1_2", "t1_3", "t1_4", "t1_5",

"t2_1", "t2_2", "t2_3", "t2_4", "t2_5")

# note: t1_6 and t2_6 are the distal outcome

nvarsBin <- 1

ntvarsBin <- nvarsBin*2

selVarsBin <- c("t1_6", "t2_6")

# Definition variables

defSex <- mxMatrix(type="Full",nrow=1,ncol=2,free=F,labels=c("data.sex1","data.sex2"),name="xSex")

betaSex <- mxMatrix(type="Full",nrow=1,ncol=6,free=T,values=0.1,labels=c(rep("bSext"),"bSexd"),name="BSex")

# Distal outcome as categorical

MZdata$t1_6 <- mxFactor(MZdata$t1_6, levels=c(0:1) )

DZdata$t1_6 <- mxFactor(DZdata$t1_6, levels=c(0:1) )

MZdata$t2_6 <- mxFactor(MZdata$t2_6, levels=c(0:1) )

DZdata$t2_6 <- mxFactor(DZdata$t2_6, levels=c(0:1) )

threshold <- mxMatrix( type="Full", nrow=1, ncol=ntvarsBin, free=TRUE, values=.5, label="thre", name="threshold" )

expThresh <- mxAlgebra( expression= threshold, name="expThresh" )

# ϕ phi

# covariance between the exogenous variables i.e. between Intercept, Slope, & F factors

phiA_labs <- c("VA_Int","Acov_IntS","Acov_IntF","VA_S","Acov_SF","VA_F")

phiC_labs <- c("VC_Int","Ccov_IntS","Ccov_IntF","VC_S","Ccov_SF","VC_F")

phiE_labs <- c("VE_Int","Ecov_IntS","Ecov_IntF","VE_S","Ecov_SF","VE_F")

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

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 = c(T,T,T,T,T,F), labels=phiE_labs, name = "phiE",

values = c( 0.60,

0.21,0.60,

0.21,0.21,1))

# correlations between Intercept, Slope, & F 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 from Intercept, Slope, & F to observed scores

# 6 (time points) x 3 factors (I, S, F) matrix

Glabs <- c(paste("i",1:nvars,sep=""), paste("s",1:nvars,sep=""), paste("q",1:nvars,sep=""))

gamma <- mxMatrix( type="Full", nrow=nvars, ncol=nGrow, free=F, labels=Glabs, name="gamma",

values=c(1,1,1,1,1,0, -2,-1,0,1,2,0, 0,0,0,0,0,1) )

# Ψ psi

# Innovations

psiA_labs <- c("VA_i11","VA_i22","VA_i33","VA_i44","VA_i55","VA_i66")

psiC_labs <- c("VC_i11","VC_i22","VC_i33","VC_i44","VC_i55","VC_i66")

psiE_labs <- c("VE_i11","VE_i22","VE_i33","VE_i44","VE_i55","VE_i66")

psiA <- mxMatrix(type = "Diag", nrow = nvars, ncol = nvars, free = c(F,T,T,T,T,F), labels = psiA_labs, name = "psiA",

values = c( 0,0.5,0.5,0.5,0.5,0))

psiC <- mxMatrix(type = "Diag", nrow = nvars,ncol = nvars, free = c(F,T,T,T,T,F), labels = psiC_labs, name = "psiC",

values = c(0,0.5,0.5,0.5,0.5,0))

psiE <- mxMatrix(type = "Diag", nrow = nvars, ncol = nvars, free = c(F,T,T,T,T,F), labels = psiE_labs, name = "psiE",

values = c(0,0.5,0.5,0.5,0.5,0))

# Expected variance/covariance matrix

A <- mxAlgebra(gamma %&% phiA + psiA, name = "A")

C <- mxAlgebra(gamma %&% phiC + psiC, name = "C")

E <- mxAlgebra(gamma %&% phiE + psiE, 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" )

stanVg <- mxAlgebra( expression=cbind(phiA/totalVg,phiC/totalVg, phiE/totalVg), name="stanVg")

# Expected means

iLoadings <- mxMatrix(type="Full",nrow=1,ncol=6,free=F,values=c( 1, 1,1,1,1,0),name="iLoads")

sLoadings <- mxMatrix(type="Full",nrow=1,ncol=6,free=F,values=c(-2,-1,0,1,2,0),name="sLoads")

meanIntrcpt <- mxMatrix(type="Full",nrow=1,ncol=1,free=T,values=0.1,labels="muIntrcpt",name="MuIntrcpt")

meanSlope <- mxMatrix(type="Full",nrow=1,ncol=1,free=T,values=0.1,labels="muSlope",name="MuSlope")

ManMean <- mxAlgebra(cbind(MuIntrcpt*iLoads + MuSlope*sLoads, MuIntrcpt*iLoads + MuSlope*sLoads) + (xSex%x%BSex), name="ManMean" )

# Data Objects for Multiple Groups

dataMZ <- mxData( observed=MZdata, type="raw" )

dataDZ <- mxData( observed=DZdata, type="raw" )

# Objects for Multiple Groups

objMZ <- mxExpectationNormal( covariance="expCovMZ", means="ManMean", dimnames=selVars,

thresholds="expThresh", threshnames=c("t1_6","t2_6"))

objDZ <- mxExpectationNormal( covariance="expCovDZ", means="ManMean", dimnames=selVars, thresholds="expThresh",

threshnames=c("t1_6","t2_6"))

# 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","Time6"),colVC))

# Final steps

pars <- list( phiA,phiC,phiE,corrA,corrC,corrE,gamma,psiA,psiC,psiE,

A,C,E,V,VC,totalVg,stanVg,ManMean,threshold)

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

my_model <- mxModel( "my_model", pars, modelMZ, modelDZ, obj)

`# Fit model`

summary(my_fit <- mxRun( my_model, intervals=F ) )

However, it doesn't work because of the following error:

*Error: Unknown reference 'muIntrcpt' detected in the entity 'ManMean' in model 'my_model'*

But it seems to me the the muIntrcpt objet is correctly specified and referenced in the model, isn't it?

An additional question: I think that in this model I am not adjusting my distal outcome (t1_6 and t2_6) for the effect of sex...if I wanted to adjust my outcome for sex, I should replace this:

threshold <- mxMatrix( type="Full", nrow=1, ncol=ntvarsBin, free=TRUE, values=.5, label="thre", name="threshold" )

expThresh <- mxAlgebra( expression= threshold, name="expThresh" )

by this:

threshold <- mxMatrix( type="Full", nrow=1, ncol=ntvarsBin, free=TRUE, values=.5, label="thre", name="threshold" )

expThresh <- mxAlgebra( expression= threshold + (xSex%x%BSex), name="expThresh" )

right?

Thanks again

M

Log in or register to post comments

In reply to I really appreciate your help by massimiliano.orri

## put objects into models; restore `selVars`

Objects `iLoadings`, `sLoadings`, `meanIntercept`, and `meanSlope` have to go into the `mxModel()` statement somehow as well (perhaps as part of list `pars`). The objects have to be part of the MxModel's namespace for all this to work properly.

Once you fix that, though, your model still isn't going to run. The syntax I gave you for the phenotypic mean assumed that `selVars` would be defined as it was in your first attempt at a script. But in your newest script, you dropped the two distal outcomes from `selVars`. Since you provided `selVars` as the `dimnames` argument to `mxExpectationNormal()`, OpenMx is only going to model the ten data columns containing the time-varying trait. But the model-expected, twin-pair phenotypic mean vector and covariance matrix are too big, i.e. 12x12 and 1x12, respectively. That's going to cause an error.

Your specification of the threshold looks OK to me (although you could just directly use the 'threshold' MxMatrix, that is, the 'expThresh' MxAlgebra isn't necessary).

The syntax I gave you adjusts the "liability-scale" mean of the distal outcome for sex. You could adjust its threshold for sex instead, similarly to what you propose, except you'd need to define a new 1x1 MxMatrix to hold the free parameter for the effect of sex on the threshold. If you do that, then fix parameter 'bSexd' to zero in `betaSex`.

Log in or register to post comments

In reply to put objects into models; restore `selVars` by AdminRobK

## Optimizer returned a non-zero status code 5

Now the code works.

I obtained the following Warning:

*Warning message:

In model 'my model' Optimizer returned a non-zero status code 5. The Hessian at the solution does not appear to be convex. See ?mxCheckIdentification for possible diagnosis (Mx status RED).*

You said that it may me symptomatic of non-identification. Any suggestion on how to solve the problem?

Log in or register to post comments

In reply to Optimizer returned a non-zero status code 5 by massimiliano.orri

## status RED (code 5 or 6)

Good to hear.

Status RED can occur even when the model's identified. You could again try `mxCheckIdentification()` and see what it says. If it say the model isn't locally identified, then lingering unidentification is likely the cause.

Otherwise, see here for a thread about what to do about status RED in general. Status RED is also quite common in analyses of threshold traits, sometimes unavoidably so, even when the optimizer has indeed found a local minimum. In your particular case, once you're pretty sure the model is identified, then be sure you're using the CSOLNP optimizer, and replace `mxRun()` in your script with `mxTryHardOrdinal()`.

Log in or register to post comments

In reply to status RED (code 5 or 6) by AdminRobK

## Error in mxCheckIdentification()

*Error in dimnames(x) <- dn :

length of 'dimnames' [2] not equal to array extent*

Of course, I don't understand what it is about, a part that 2 objects are not of the same size while they are expected to be of the same size.

Is there a way to fix it or to check identification in another way? To understand if my model is really not identified of it is just an issue with the threshold variable

The model with CSOLNP optimizer (or SLSQP, I tried this too) and *mxTryHardOrdinal()* gives pretty much the same results as the previous one, ie always the same problem of Hessian. For instance, at each iteration of *mxTryHard()* the result is:

*Not all eigenvalues of Hessian are greater than 0: ....*

Log in or register to post comments

In reply to Error in mxCheckIdentification() by massimiliano.orri

## probably OK

That might be a bug.

I think it's pretty likely to be due to the threshold variable.

If you're getting pretty much the same results--as in, similar parameter estimates and similar -2logL values--then that's a pretty good sign that the model is indeed identified, and that the optimizer is finding a local minimum. I think you are probably OK.

Log in or register to post comments

In reply to probably OK by AdminRobK

## Strange parameters

Indeed, I obtained the same -2LogL and parameter estimates with all the optimizers.

However, I am not sure that everything is right when I look at the standardized variance components...

Here the results (AE model):

mxAlgebra 'stanVg'

$formula: cbind(phiA/totalVg, phiC/totalVg, phiE/totalVg)

$result:

[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]

[1,] 0.8159866 0.6890866 0.3237748 0 0 0 0.1840134 0.3109134 0.6762252

[2,] 0.6890866 0.4352550 -28.6371619 0 0 0 0.3109134 0.5647450 29.6371619

[3,] 0.3237748 -28.6371619 0.5483035 0 0 0 0.6762252 29.6371619 0.4516965

The are 2 values >1, which should not be possible....

Any idea of why it is the case? Or potential solutions?

Log in or register to post comments

In reply to Strange parameters by massimiliano.orri

## you are mistaken

Perhaps you want to use, e.g.,

`cov2cor(phiA)`

instead?Log in or register to post comments

In reply to you are mistaken by AdminRobK

## Thanks!

Log in or register to post comments