Login

News

Latent Growth Model and distal outcome

Posted on
No user picture. massimiliano.orri Joined: 06/11/2018
Forums
Behavioral Genetics Models
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
Replied on Thu, 09/13/2018 - 14:28
No user picture. massimiliano.orri Joined: 06/11/2018

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
Replied on Fri, 09/14/2018 - 12:35
Picture of user. AdminRobK Joined: 01/24/2014

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

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?

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?

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
Replied on Fri, 09/14/2018 - 18:31
No user picture. massimiliano.orri Joined: 06/11/2018

Thank you again for your reply.
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
Replied on Sat, 09/15/2018 - 16:22
Picture of user. AdminRobK Joined: 01/24/2014

Thank you again for your reply.
For the Error, I solved the problem updating openMx.

OK, good!

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

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

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 don't know. Perhaps.

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

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.

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 ?

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
Replied on Mon, 09/17/2018 - 17:20
No user picture. massimiliano.orri Joined: 06/11/2018

I really appreciate your help, thank you so much!
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
Replied on Tue, 09/18/2018 - 10:56
Picture of user. AdminRobK Joined: 01/24/2014

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?

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
Replied on Tue, 09/18/2018 - 13:52
No user picture. massimiliano.orri Joined: 06/11/2018

Thank you Rob - and sorry for my very basic mistakes!
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
Replied on Tue, 09/18/2018 - 15:16
Picture of user. AdminRobK Joined: 01/24/2014

Now the code works.

Good to hear.

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?

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
Replied on Tue, 09/18/2018 - 17:00
No user picture. massimiliano.orri Joined: 06/11/2018

I tried with *mxCheckIdentification()* but it gave me this error:
*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
Replied on Wed, 09/19/2018 - 10:38
Picture of user. AdminRobK Joined: 01/24/2014

I tried with mxCheckIdentification() but it gave me this error:
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.

That might be a bug.

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

I think it's pretty likely to be due to 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: ....

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
Replied on Thu, 09/20/2018 - 16:14
No user picture. massimiliano.orri Joined: 06/11/2018

Thank you very much for your help.
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
Replied on Fri, 09/21/2018 - 15:01
Picture of user. AdminRobK Joined: 01/24/2014

You are mistaken. Those three "standardized" matrices are respectively the common-factor _A_, _C_, and _E_ covariance matrices, elementwise-divided by the total common-factor covariance matrix. Their off-diagonals can certainly be greater than one, or less than zero. More or less this same topic has been discussed before on these forums here, here, here, here, and here.

Perhaps you want to use, e.g., cov2cor(phiA) instead?

Log in or register to post comments
  Login

News