You are here

Latent Growth Model and distal outcome

15 posts / 0 new
Last post
massimiliano.orri's picture
Joined: 06/11/2018 - 16:13
Latent Growth Model and distal outcome
AttachmentSize
Image icon model.png330.63 KB

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

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
seems OK, but no example script

I don't see any reason why it wouldn't work or wouldn't be identified. I don't know of any example script for such a model, though.

massimiliano.orri's picture
Joined: 06/11/2018 - 16:13
Thank you Rob for your reply.

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

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
some remarks
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.

massimiliano.orri's picture
Joined: 06/11/2018 - 16:13
Thank you again for your

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

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
phenotypic means
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.

massimiliano.orri's picture
Joined: 06/11/2018 - 16:13
I really appreciate your help

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

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
put objects into models; restore `selVars`
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.

massimiliano.orri's picture
Joined: 06/11/2018 - 16:13
Optimizer returned a non-zero status code 5

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?

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
status RED (code 5 or 6)
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().

massimiliano.orri's picture
Joined: 06/11/2018 - 16:13
Error in mxCheckIdentification()

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: ....

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
probably OK
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.

massimiliano.orri's picture
Joined: 06/11/2018 - 16:13
Strange parameters

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?

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
you are mistaken

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?

massimiliano.orri's picture
Joined: 06/11/2018 - 16:13
Thanks!

ok, I see. Thank you for your explanations, I really appreciated