# the thresholds in SEM of ordinal variables

1) I found tow kinds of scripts of saturate SEM of one ordinal variable, one script constrains the free status of thresholds, and another one does not. I can't figure out which one to follow.

[the first one:](https://static1.squarespace.com/static/58b2481a9f7456906a3b9600/t/5bcf80379140b7614a07d8d1/1540325431670/oneSATma.pdf)

`nth <- 3`

frTh <- matrix(rep(c(F,F,(rep(T,nth-2)))),nrow=nth,ncol=nv) # free status for thresholds

thinMZ <- mxMatrix( type="Full", nrow=nth, ncol=ntv, free=frTh, values=svTh, lbound=lbTh, labels=labThMZ, name="thinMZ" )

[the second one:](http://ibg.colorado.edu/cdrom2016/maes/UnivariateAnalysis/onea/oneSAToa.R)

`nth <- 3`

thinMZ <- mxMatrix( type="Full", nrow=nth, ncol=ntv, free=TRUE, values=svTh, lbound=lbTh, labels=labThMZ, name="thinMZ" )

2)

If I do need to constrain the free status of thresholds, then how many thresholds should I fix?

In the first script, the number of thresholds is 3 and two shresholds are fixed. If my ordinal variable is categorized into 5 groups, and then the number of thresholds should be 4, then how many thresholds should I fix? three or two? In fact, I don't konw the reason for fixing thresholds. So, please don't mind my kind of stupid question.

I would be really appreciated for your reply!

Thanks!

## model identification

The choices are arbitrary, but sometimes one choice is preferable to others because it makes results easier to interpret. For instance, option (3) is often used with longitudinal data, because it's easier to interpret the variance of the latent continuum changing with time, as opposed to interpreting the spread of the thresholds changing with time.

Note that in biometrical variance-components analysis (as with twin data, etc.), fixing the variance to 1 will often require an MxConstraint. An alternative to using an MxConstraint is to fix the nonshared-environmental variance component to 1, which does change the interpretation of the genetic and shared-environmental variance components..

Log in or register to post comments

In reply to model identification by AdminRobK

## Model is not locally identified

vars <- 'tea_5g' # list of variables names

covars <- c("age","region_type")

nv <- 1 # number of variables

ncv <- 2 # number of covariates

ntv <- nv*2 # number of total variables

nth <- 4

selVars <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="")

covVars <- paste(covars,c(rep(1,ncv),rep(2,ncv)),sep="")

# Select Data for Analysis

mzData <- subset(tea_5g, zygosity1==0, c(selVars,covVars))

dzData <- subset(tea_5g, zygosity1==1, c(selVars,covVars))

# tell mx the data is ordinal

mzData[,1:2] <- mxFactor(mzData[,1:2], levels= c(0:nth))

dzData[,1:2] <- mxFactor(dzData[,1:2], levels= c(0:nth))

# Set Starting Values

svBe <- 0.01 # start value for beta regressions

svLTh <- -1.5 # start value for first threshold

svITh <- 1 # start value for increments

svTh <- matrix(rep(c(svLTh,(rep(svITh,nth-1)))),nrow=nth,ncol=nv) # start value for thresholds

lbTh <- matrix(rep(c(-3,(rep(0.001,nth-1))),nv),nrow=nth,ncol=nv) # lower bounds for thresholds

svCor <- .5 # start value for correlations

lbCor <- -0.99 # lower bounds for correlations

ubCor <- 0.99 # upper bounds for correlations

labThMZ <- c(paste("t",1:nth,"MZ1",sep=""),paste("t",1:nth,"MZ2",sep=""))

labThDZ <- c(paste("t",1:nth,"DZ1",sep=""),paste("t",1:nth,"DZ2",sep=""))

`modelMZ<- mxModel("MZ",`

# Define definition variables

mxMatrix( type="Full", nrow=ncv, ncol=2, free=F, label=c("data.age1","data.region_type1","data.age2","data.region_type2"), name="DefVars"),

mxMatrix( type="Full", nrow=1, ncol=ncv, free=TRUE, values=svBe,labels=c("Bage","Bregion_type"), name="BageTH"),

# expected mean Matrices

mxMatrix( type="Zero", nrow=1, ncol=ntv, name="meanG" ), #intercepts#

mxAlgebra( expression= meanG + BageTH%*%DefVars, name="expMeanMZ" ),

# Define threshold matrices : nth--number of thresholds

mxMatrix( type="Full", nrow=nth, ncol=ntv, free=T, values=svTh, lbound=lbTh, labels=labThMZ, name="thinMZ" ),

mxMatrix( type="Lower", nrow=nth, ncol=nth, free=FALSE, values=1, name="inc" ),

mxAlgebra( expression= inc %*% thinMZ, name="threMZ" ),

# expected variance–covariance Matrices

mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=TRUE, values=svCor, lbound=lbCor, ubound=ubCor, labels="rMZ",name="corMZ" ),

# Read in the data

mxData( observed=mzData, type="raw" ),

mxExpectationNormal( covariance="corMZ", means="expMeanMZ", dimnames=selVars, thresholds="threMZ"),

mxFitFunctionML()

)

modelDZ<- mxModel("DZ",

mxMatrix( type="Full", nrow=ncv, ncol=2, free=F, label=c("data.age1","data.region_type1","data.age2","data.region_type2"), name="DefVars"),

mxMatrix( type="Full", nrow=1, ncol=ncv, free=TRUE, values=svBe,labels=c("Bage","Bregion_type"), name="BageTH"),

mxMatrix( type="Zero", nrow=1, ncol=ntv, name="meanG" ), #two labels#

mxAlgebra( expression= meanG + BageTH%*%DefVars, name="expMeanDZ" ),

mxMatrix( type="Full", nrow=nth, ncol=ntv, free=frTh, values=svTh, lbound=lbTh, labels=labThDZ, name="thinDZ" ),

mxMatrix( type="Lower", nrow=nth, ncol=nth, free=FALSE, values=1, name="inc" ),

mxAlgebra( expression= inc %*% thinDZ, name="threDZ" ),

mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=TRUE, values=svCor, lbound=lbCor, ubound=ubCor, labels="rDZ",name="corDZ" ),

mxData( observed=dzData, type="raw"),

mxExpectationNormal( covariance="corDZ", means="expMeanDZ", dimnames=selVars, thresholds="threDZ" ),

mxFitFunctionML()

)

Conf <- mxCI (c ('MZ.corMZ','DZ.corDZ'))

SatModel <- mxModel( "Sat", modelMZ, modelDZ,mxFitFunctionMultigroup(c('MZ.fitfunction','DZ.fitfunction')), Conf )

# -----------------------------------------------------------------------------------------------

# RUN Saturated Model

SatFit <- mxTryHard(SatModel, intervals=F,extraTries = 41) #,bestInitsOutput=T,silent=F#

(SatSumm <- summary(SatFit))

mxCheckIdentification(SatFit, details=F)

# the result showed 'Model is not locally identified'

Thank you so much!

Log in or register to post comments

In reply to Model is not locally identified by diana

## probably false alarm

That's probably a false alarm. Try `mxCheckIdentification()` again, with `details=TRUE`, to see the function's guess as to which parameters are unidentified. Also, try running your model (perhaps with `mxTryHardOrdinal()`), and check to see if `mxCheckIdentification()` says it's locally unidentified at the solution, too.

Log in or register to post comments

In reply to probably false alarm by AdminRobK

## I have one more question, how

#select ordinal variable 1-------------------

table(tea5g_smk4g$tea_5g1,tea5g_smk4g$tea_5g2)

nth1 <- 4 # number of thresholds

vars1 <- c('tea_5g') # list of ordinal variables names

nvo1 <- 1 # number of ordinal variables

ntvo1 <- nvo1*2 # number of total ordinal variables

ordVars1 <- paste(vars1,c(rep(1,nvo1),rep(2,nvo1)),sep="")

#select ordinal variable 2-------------------

table(tea5g_smk4g$smk_num_4g1,tea5g_smk4g$smk_num_4g2)

nth2 <- 3 # number of thresholds

vars2 <- c('smk_num_4g') # list of ordinal variables names

nvo2 <- 1 # number of ordinal variables

ntvo2 <- nvo2*2 # number of total ordinal variables

ordVars2 <- paste(vars2,c(rep(1,nvo2),rep(2,nvo2)),sep="")

# Select Variables for Analysis---------------------

vars <- c('tea_5g','smk_num_4g') # list of variables names

covars <- c("age","region_type")

nv <- nvo1+nvo2 # number of variables

ncv <- 2 # number of covariates

ntv <- nv*2 # number of total variables

selVars <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="") #"tea_5g1" "smk_num_4g1" "tea_5g2" "smk_num_4g2"

covVars <- paste(covars,c(rep(1,ncv),rep(2,ncv)),sep="")

# Select Data for Analysis---------------------------

mzData <- subset(tea5g_smk4g, zygosity2==0, c(selVars,covVars))

dzData <- subset(tea5g_smk4g, zygosity2==1, c(selVars,covVars))

`mzDataF <- cbind(mxFactor( x=mzData[,ordVars1], levels=c(0:nth1)),mxFactor( x=mzData[,ordVars2], levels=c(0:nth2)),mzData[,covVars] )`

dzDataF <- cbind(mxFactor( x=dzData[,ordVars1], levels=c(0:nth1)),mxFactor( x=dzData[,ordVars2], levels=c(0:nth2)),dzData[,covVars] )

I don't know how to define threshold matrix.

Thanks so much!

Log in or register to post comments

In reply to I have one more question, how by diana

## see User Guide

`NA`

."Log in or register to post comments

In reply to see User Guide by AdminRobK

## some results and new problems

1) After I ran

`mxCheckIdentification(), details=TRUE`

, the result showed the coefficients of covatiates were not identified. Does that matters?2) I followed your advice by fixing corresponding thresholds with starting values of NA. However, I got the waring as follows:

**Computing Hessian and/or standard errors and/or confidence intervals from imperfect solutionError : In model 'Sat' I was expecting 3 thresholds in column 'smk_num_4g1' of matrix/algebra 'MZ.threMZ' but I hit NA values after only 0 thresholds. You need to increase the number of thresholds for 'smk_num_4g1' and give them values other than NA

Retry limit reached; solution not found. Best fit=25943569000 (started at 47142898000) (16 attempt(s): 16 valid, 0 errors)**

Here is my script and I'm sorry to bother you to check what's wrong with it. Thank you so much!!!!

#select ordinal variable 1-------------------

table(tea5g_smk4g$tea_5g1,tea5g_smk4g$tea_5g2)

nth1 <- 4 # number of thresholds

vars1 <- c('tea_5g') # list of ordinal variables names

nvo1 <- 1 # number of ordinal variables

ntvo1 <- nvo1*2 # number of total ordinal variables

ordVars1 <- paste(vars1,c(rep(1,nvo1),rep(2,nvo1)),sep="")

#select ordinal variable 2------------------

table(tea5g_smk4g$smk_num_4g1,tea5g_smk4g$smk_num_4g2)

nth2 <- 3 # number of thresholds

vars2 <- c('smk_num_4g') # list of ordinal variables names

nvo2 <- 1 # number of ordinal variables

ntvo2 <- nvo2*2 # number of total ordinal variables

ordVars2 <- paste(vars2,c(rep(1,nvo2),rep(2,nvo2)),sep="")

# Select Variables for Analysis---------------------

vars <- c('tea_5g','smk_num_4g') # list of variables names

covars <- c("age","region_type")

nv <- nvo1+nvo2 # number of variables

ncv <- 2 # number of covariates

ntv <- nv*2 # number of total variables

selVars <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="") #"tea_5g1" "smk_num_4g1" "tea_5g2" "smk_num_4g2"

covVars <- paste(covars,c(rep(1,ncv),rep(2,ncv)),sep="")

# Select Data for Analysis---------------------------

mzData <- subset(tea5g_smk4g, zygosity2==0, c(selVars,covVars))

dzData <- subset(tea5g_smk4g, zygosity2==1, c(selVars,covVars))

mzDataF <- cbind(mxFactor( x=mzData[,ordVars1], levels=c(0:nth1)),mxFactor( x=mzData[,ordVars2], levels=c(0:nth2)),mzData[,covVars] )

dzDataF <- cbind(mxFactor( x=dzData[,ordVars1], levels=c(0:nth1)),mxFactor( x=dzData[,ordVars2], levels=c(0:nth2)),dzData[,covVars] )

# Set Starting Values

frMV <- c(FALSE,FALSE) # free status for meanG of variables

frTH <- matrix(c(rep(T,nth1),c(rep(T,nth2),F)),4,4) # free status for thresholds

svMe <- c(0,0) # start value for means

svBe <- 0.01 # start value for beta regressions

svLTh <- 0 # start value for first threshold

svITh <- 1 # start value for increments

svTh <- matrix(c(svLTh,rep(svITh,nth1-1),svLTh,rep(svITh,nth2-1),NA),nrow=nth1,ncol=ntv) # start value for thresholds

lbTh <- matrix(c(-3,rep(0.001,nth1-1),-3,rep(0.001,nth2-1),NA),nrow=nth1,ncol=ntv) # lower bound for thresholds

svCor <- .5 # start value for correlations

# Create Labels

labdata <- cbind(paste("data.",covars[1],c(1,2),sep = ""),paste("data.",covars[2],c(1,2),sep = ""))

# [,1] [,2]

# [1,] "data.age1" "data.region_type1"

# [2,] "data.age2" "data.region_type2"

labeta <-cbind(paste("B",covars,"_",vars[1],sep = ""),paste("B",covars,"_",vars[2],sep = ""))

# [,1] [,2]

# [1,] "Bage_tea_5g" "Bage_smk_num_4g"

# [2,] "Bregion_type_tea_5g" "Bregion_type_smk_num_4g"

labTh <- function(lab,vars,nth) { paste(paste("t",1:nth,lab,sep=""),rep(vars,each=nth),sep="") }

labThMZ <- matrix(c(labTh("MZ","tea_5g1",nth1),labTh("MZ","smk_num_4g1",nth2),NA,labTh("MZ","tea_5g2",nth1),labTh("MZ","smk_num_4g2",nth2),NA),nth1,ntv)

labThDZ <- matrix(c(labTh("DZ","tea_5g1",nth1),labTh("DZ","smk_num_4g1",nth2),NA,labTh("DZ","tea_5g2",nth1),labTh("DZ","smk_num_4g2",nth2),NA),nth1,ntv)

labThZ <- matrix(c(labTh("Z","tea_5g1",nth1),labTh("Z","smk_num_4g1",nth2),NA,labTh("Z","tea_5g2",nth1),labTh("Z","smk_num_4g2",nth2),NA),nth1,ntv)

`labCvMZ <- c("MZCVT1", "MZWP1", "MZCTT1T2","MZCTT2T1", "MZWP2", "MZCVT2") # labels for (co)variances for MZ twins`

labCvDZ <- c("DZCVT1", "DZWP1", "DZCTT1T2","DZCTT2T1", "DZWP2", "DZCVT2") # labels for (co)variances for DZ twins

labCvZ <- c("ZCVT1", "ZWP1", "ZCTT1T2","ZCTT2T1", "ZWP2", "ZCVT2")

##########################################

# wc1 smk1 | wc2 smk2

# wc1 VP1T1 |

# smk1 CVT1* VP2T1 |

#-----------------------------------------

# wc2 WP1* CTT2T1 | VP1T2

# smk2 CTT1T2* WP2* | CVT2 VP2T2

##########################################

modelMZ <- mxModel( "MZ",

mxMatrix( type="Full", nrow=2, ncol=ncv, free=F, label=labdata, name="MZDefVars"),

mxMatrix( type="Full", nrow=ncv, ncol=nv, free=T, values=svBe,labels=labeta, name="Beta"),

mxMatrix( type="Full", nrow=1, ncol=nv, free=frMV, values=svMe, name="meanG" ),

mxAlgebra( expression= cbind(meanG + MZDefVars[1,] %*% Beta,meanG + MZDefVars[2,] %*% Beta), name="expMeanMZ" ),

mxMatrix( type="Full", nrow=nth1, ncol=ntv, free=frTH, values=svTh, lbound=lbTh, labels=labThMZ, name="thinMZ" ),

mxMatrix( type="Lower", nrow=nth1, ncol=nth1, free=FALSE, values=1, name="inc" ),

mxAlgebra( expression= inc %*% thinMZ, name="threMZ" ),

mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=svCor, labels=labCvMZ, name="expCovMZ"), #lbound=-.99, ubound=.99

mxData( observed=mzDataF, type="raw" ),

mxExpectationNormal( covariance="expCovMZ", means="expMeanMZ", dimnames=selVars, thresholds="threMZ" ),

mxFitFunctionML()

)

modelDZ <- mxModel( "DZ",

mxMatrix( type="Full", nrow=2, ncol=ncv, free=F, label=labdata, name="DZDefVars"),

mxMatrix( type="Full", nrow=ncv, ncol=nv, free=T, values=svBe,labels=labeta, name="Beta"),

mxMatrix( type="Full", nrow=1, ncol=nv, free=frMV, values=svMe, name="meanG" ),

mxAlgebra( expression= cbind(meanG + DZDefVars[1,] %*% Beta,meanG + DZDefVars[2,] %*% Beta), name="expMeanDZ" ),

mxMatrix( type="Full", nrow=nth1, ncol=ntv, free=frTH, values=svTh, lbound=lbTh, labels=labThDZ, name="thinDZ" ),

mxMatrix( type="Lower", nrow=nth1, ncol=nth1, free=FALSE, values=1, name="inc" ),

mxAlgebra( expression= inc %*% thinDZ, name="threDZ" ),

mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=svCor, labels=labCvDZ, name="expCovDZ"), # lbound=-.99, ubound=.99

mxData( observed=dzDataF, type="raw" ),

mxExpectationNormal( covariance="expCovDZ", means="expMeanDZ", dimnames=selVars, thresholds="threDZ" ),

mxFitFunctionML()

)

# Combine Groups

Conf <- mxCI (c ('MZ.expCovMZ[3,1]','DZ.expCovDZ[3,1]','MZ.expCovMZ[4,2]','DZ.expCovDZ[4,2]','MZ.expCovMZ[2,1]','DZ.expCovDZ[2,1]','MZ.expCovMZ[4,1]','DZ.expCovDZ[4,1]') )#'MZ.expCovMZ[3,1]','MZ.expCovMZ[4,2]'

SatModel <- mxModel( "Sat", modelMZ, modelDZ,mxFitFunctionMultigroup(c('MZ','DZ')), Conf )

SatFit <- mxTryHardOrdinal(SatModel, intervals=F,extraTries = 15)

**Wish you have a good day!**

Thanks♪(･ω･)ﾉ

Log in or register to post comments

In reply to some results and new problems by diana

## Data objects

Log in or register to post comments

In reply to Data objects by AdminRobK

## im sorry i dont have my data

please leave more information since when i get up you would possibly get off work.

thank you very much！！！

Log in or register to post comments

In reply to some results and new problems by diana

## unsure

Did `mxTryHardOrdinal()` return a fitted MxModel object? If so, does `mxCheckIdentification()` still say the model is unidentified?

Log in or register to post comments

In reply to unsure by AdminRobK

## Hi, Rob!

1）After I ran

`str(mzDataF)`

, it turned out that it was a dataframe. According to the warning, it seems that there is something wrong in the matrix of **MZ.threMZ**?2)

`mxTryHardOrdinal()`

can return a fitted solution, but`mxCheckIdentification()`

still say the model is not identified and the non-identified parameters are two coefficients of age and region_type.Thanks again!

Log in or register to post comments

In reply to Hi, Rob! by diana

## What's the actual output of

Log in or register to post comments

In reply to What's the actual output of by AdminRobK

## Hi, the result was in the

Thanks!

Log in or register to post comments

In reply to Hi, Rob! by diana

## definition variables

`mxCheckIdentification()`

doesn't necessarily give accurate results when the MxModel uses definition variables.Log in or register to post comments

## did you solve the issue?

I am wondering if you solved the issue with the error OpenMx gives when using this script? I have one ordinal variable with three thresholds and one dichotomous, so the matrix of thresholds have two NA values as it is written in the documentation. This gives the following error: "I was expecting 1 thresholds in column ... of matrix/algebra 'MZ.expThreshold' but I hit NA values after only 0 thresholds. You need to increase the number of thresholds for ... and give them values other than NA". I suspect this comes from multiplication of a lower matrix with the matrix that contains NA.

When I replace NA's with zero's for example, I don't get this error and the model runs. However, a new error message occurs: "...fit is not finite (Found 1 thresholds too close together in column 3.)"

So how can we specify a threshold matrix for several variables when variables have different number of thresholds, so that multiplication works?

Log in or register to post comments

In reply to did you solve the issue? by Julia

## Would you mind posting the

Log in or register to post comments

## Thank you for your reply!

nv = 4

nvo = 2

nth = 3

ninc <- nth-1

frMV = c(F,T,T,F)

frTh = c(T,T,T,T,F,F)

meanLabs <- paste(vars, "mean",sep="_")

threshLabs <- paste(rep(varso,each=nth),c("thresh",rep("inc",ninc)),c(1,1:ninc),sep='_')

svMe <- c(0,0,0,0)

svTh <- 0

defAge1 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels="data.Zage_s_1", name="Age1")

defAge2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels="data.Zage_s_2", name="Age2")

defSex1 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels="data.sex_s_1", name="Sex1")

defSex2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels="data.sex_s_2", name="Sex2")

# Regression effects

B_Age <- mxMatrix( type="Full", nrow=1, ncol=nv, free=frMV, values=c(0,.1,.1,0), labels=betaLabs_age, name="bAge" )

B_Age_Thre <-mxMatrix( type="Full", nrow=nth, ncol=nvo, free=frTh, values=c(rep(.2,nth),0.2,NA,NA), labels=betaLabsThre_age, name="bAgeThre" )

B_Sex <- mxMatrix( type="Full", nrow=1, ncol=nv, free=frMV, values=c(0,.1,.1,0), labels=betaLabs_sex, name="bSex" )

B_Sex_Thre <-mxMatrix( type="Full", nrow=nth, ncol=nvo, free=frTh, values=c(rep(.2,nth),0.2,NA,NA), labels=betaLabsThre_sex, name="bSexThre" )

#setting up the regression

intercept <- mxMatrix( type="Full", nrow=1, ncol=nv, free=frMV, values=svMe, labels=meanLabs, name="intercept" )

expMean <- mxAlgebra( expression = cbind(intercept + (bAge%x%Age1) + (bSex%x%Sex1), intercept + (bAge%x%Age2) + (bSex%x%Sex2)) , name="expMean")

`inc <-mxMatrix( type="Lower",nrow=nth, ncol=nth, free=F, values=1, name="Low")`

thre <-mxMatrix( type="Full", nrow=nth, ncol=nvo, free=frTh, values=svTh, lbound=c(-3,rep(.001,ninc)), ubound=3, labels=threshLabs, name="Threshold")

expThre <-mxAlgebra( expression= cbind(Low%*%Threshold + (bAgeThre%x%Age1) + (bSexThre%x%Sex1), Low%*%Threshold + (bAgeThre%x%Age2) + (bSexThre%x%Sex2) ), name="expThreshold")

I believe the product of Low and Threshold produces only NA's in the second column. At least this is what is happening when I do such multiplication with conventional matrices in R.

At the moment I figured out a workaround by defining thresholds separately for each variable and then concatenating them together. But I wonder if this can be avoided and if my code was wrong.

defAge1 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels="data.Zage_s_1", name="Age1")

defAge2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels="data.Zage_s_2", name="Age2")

defSex1 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels="data.sex_s_1", name="Sex1")

defSex2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels="data.sex_s_2", name="Sex2")

# Regression effects

B_Age <- mxMatrix( type="Full", nrow=1, ncol=nv, free=frMV, values=c(0,.1,.1,0), labels=betaLabs_age, name="bAge" )

B_Age_Thre <-mxMatrix( type="Full", nrow=nth, ncol=nvo, free=frTh, values=c(rep(.2,nth),0.2,NA,NA), labels=betaLabsThre_age, name="bAgeThre" )

B_Sex <- mxMatrix( type="Full", nrow=1, ncol=nv, free=frMV, values=c(0,.1,.1,0), labels=betaLabs_sex, name="bSex" )

B_Sex_Thre <-mxMatrix( type="Full", nrow=nth, ncol=nvo, free=frTh, values=c(rep(.2,nth),0.2,NA,NA), labels=betaLabsThre_sex, name="bSexThre" )

#setting up the regression

intercept <- mxMatrix( type="Full", nrow=1, ncol=nv, free=frMV, values=svMe, labels=meanLabs, name="intercept" )

expMean <- mxAlgebra( expression = cbind(intercept + (bAge%x%Age1) + (bSex%x%Sex1), intercept + (bAge%x%Age2) + (bSex%x%Sex2)) , name="expMean")

`inc <-mxMatrix( type="Lower",nrow=nth, ncol=nth, free=F, values=1, name="Low")`

threVar1 <-mxMatrix( type="Full", nrow=nth1, ncol=1, free=T, values=svTh, lbound=c(-3,rep(.001,ninc)), ubound=3, labels=c("healthp_s_thresh_1", "healthp_s_inc_1", "healthp_s_inc_2"), name="ThresholdVar1")

threVar2 <-mxMatrix( type="Full", nrow=nth1, ncol=1, free=c(rep(T,nth2),rep(F,nth1-nth2)), values=c(rep(svTh,nth2),rep(NA,nth1-nth2)), lbound=c(-3,rep(.001,ninc)), ubound=3, labels=c("ibs_gen_s_thresh_1", "ibs_gen_s_inc_1", "ibs_gen_s_inc_2"), name="ThresholdVar2")

threMatVar1 = mxAlgebra( expression= cbind(Low%*%ThresholdVar1), name='ThresholdMatvar1')

thresholds = mxAlgebra(cbind(ThresholdMatvar1, ThresholdVar2), name='ThresholdsInc')

expThre <-mxAlgebra( expression= cbind(ThresholdsInc + (bAgeThre%x%Age1) + (bSexThre%x%Sex1), ThresholdsInc + (bAgeThre%x%Age2) + (bSexThre%x%Sex2) ), name="expThreshold")

Thank you again!

Log in or register to post comments

In reply to Thank you for your reply! by Julia

## I am not sure what the

Does your workaround let you run your model, and do you get sensible-looking results?

Log in or register to post comments

In reply to I am not sure what the by AdminRobK

## It does run

Log in or register to post comments

## thresholds when not all vars have the same number of levels

You might well benefit from using

`umxThresholdMatrix()`

This handles all the details of building a thresholds matrix automatically, including sizing it, filling with reasonable defaults, and ensuring that levels are kept in order by using positive offsets instead of thresholds. It also handles twins properly. See

`?umx::umxThresholdMatrix`

for examples.Log in or register to post comments

In reply to thresholds when not all vars have the same number of levels by tbates

## leaving unused thresholds at zero or some other number

Log in or register to post comments