You are here

NA values hit for thresholds

6 posts / 0 new
Last post
bhuibregtse's picture
Offline
Joined: 03/06/2014 - 22:15
NA values hit for thresholds

I am using a modified script for a bivariate twin analysis with two dichotomous variables.

I am stuck with the error message: "Error: In model 'Sat' I was expecting 1 thresholds in column 'nssi1' of matrix/algebra 'DZf.expThresDZf' but I hit NA values after only 0 thresholds. You need to increase the number of thresholds for 'nssi1' and give them values other than NA"

Another user had a similar issue- see previous discussion: (https://openmx.ssri.psu.edu/thread/2953). However, I do not believe my problem has to do with the mxFactor command, factors in general, or NAs in my data.

The error message implies there is something wrong with the start values I provided for the threshold ("error message... suggests that they are somehow started at NA"). I have tried a variety of different start values, but this does not seem to be the issue.

Both of my variables have incredibly low prevalence, although there are twins pairs in all groups (MZF, MZM, DZF, DZM) who are both affected (no missing cells in contingency table). Could low prevalence cause an issue with the thresholds?

Any help appreciated! Thanks.

tbates's picture
Offline
Joined: 07/31/2009 - 14:25
thresholds= expThresDZf

The error looks like the thresholds matrix you are passing in is the wrong shape, or has the wrong values (not NAs in the data).
i.e., where you call mxExpectationNormal(covariance, means, thresholds = expThresDZf), the thresholds matrix (or algebra?) "expThresDZf" needs to have as many rows as max thresholds, and as many columns as there are binary/ordinal variables in your data.

Attach your script so we can see?

bhuibregtse's picture
Offline
Joined: 03/06/2014 - 22:15
Thank you for the quick

Thank you for the quick reply.

Script is adapted from http://ibg.colorado.edu/cdrom2014/rijsdijk/Ordinal/Bivariate/BivThreshLi... in a few ways:
1)Has 2 binary vars rather than 1 binary and 1 ordinal
2)Uses more zyg groups to eventually test for sex differences
3)Switched to openMx 2.0 funtions

(read in dataset is mj_nssi)

#Make Factors
mj_nssi$mj1 <-mxFactor(mj_nssi$mj1, levels=c(0:1), exclude = NA)
mj_nssi$mj2 <-mxFactor(mj_nssi$mj2, levels=c(0:1), exclude = NA)
mj_nssi$nssi1 <-mxFactor(mj_nssi$nssi1, levels=c(0:1), exclude = NA)
mj_nssi$nssi2 <-mxFactor(mj_nssi$nssi2, levels=c(0:1), exclude = NA)

#Selects Vars for Analysis
vars <-c('nssi','mj')
selVars <-c('nssi1', 'mj1', 'nssi2', 'mj2' )
useVars <-c('nssi1', 'mj1', 'nssi2', 'mj2', 'age1', 'age2')

# Select Data for Analysis
mzfData <- subset(mj_nssi, zyg==1, useVars)
mzmData <- subset(mj_nssi, zyg==2, useVars)
dzfData <- subset(mj_nssi, zyg==3, useVars)
dzmData <- subset(mj_nssi, zyg==4, useVars)

nv <- 2 # number of variables per twin
ntv <- nv*2 # number of variables per pair
nth <- 1 # number of max thresholds

# 1) Fits a constrained Polychoric correlation model
# TH same across twins but different across zyg groups
# Age effect is different across variables, but same across thresholds within variables (if c>2) ###?
# There is one overall rPH between var1-2 and the x-trait x-twin correlations are symmetric
# ------------------------------------------------------------------------------------------------------------------------------

# CREATE LABELS & START VALUES as objects(to ease specification)
LabThMZf <-c('Tmz_11f','imz_11f') # THs for var 1 and 2 for a twin individual (mz)
LabThMZm <-c('Tmz_11m','imz_11m')
LabCorMZf <-c('r21f','rMZ1f','MZxtxtf','MZxtxtf','rMZ2f','r21f') #labels for correlation matrix
LabCorMZm <-c('r21m','rMZ1m','MZxtxtm','MZxtxtm','rMZ2m','r21m')

LabThDZf <-c('Tdz_11f','idz_11f')
LabThDZm <-c('Tdz_11m','idz_11m')
LabCorDZf <-c('r21f','rDZ1f','DZxtxtf','DZxtxtf','rDZ2f','r21f')
LabCorDZm <-c('r21m','rDZ1m','DZxtxtm','DZxtxtm','rDZ2m','r21m')

LabCovf <-c('BageThnnsif', 'BageThnnsif, 'BageThmjf', 'BageThmjf')
LabCovm <-c('BageThnnsim', 'BageThnnsim', 'BageTmjm', 'BageTmjm')

ThPat <-c(T,T)

StCorMZf <-c(.1, .6, .2, .2, .6, .1)
StCorMZm <-c(.1, .6, .2, .2, .6, .1)
StCorDZf <-c(.1, .2, .1, .1, .2, .1)
StCorDZm <-c(.1, .2, .1, .1, .2, .1)

StTH <-c(1.2816, 0.9741)

# Define definition variables to hold the Covariates
obsAge1 <- mxMatrix( type="Full", nrow=1, ncol=1, free=F, labels=c("data.age1"), name="Age1")
obsAge2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=F, labels=c("data.age2"), name="Age2")

# Matrix & Algebra for expected means (SND), Thresholds, effect of Age on Th and correlations
Mean <-mxMatrix( type="Zero", nrow=1, ncol=ntv, name="M" )

betaAf <-mxMatrix( type="Full", nrow=nth, ncol=nv, free=T, values=.2, labels=LabCovf, name="BageTHf" )
betaAm <-mxMatrix( type="Full", nrow=nth, ncol=nv, free=T, values=.2, labels=LabCovm, name="BageTHm" )

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

Tmzf <-mxMatrix( type="Full", nrow=nth, ncol=nv, free=ThPat, values=StTH, lbound= -3, ubound=3, labels=LabThMZf, name="ThMZf")
Tmzm <-mxMatrix( type="Full", nrow=nth, ncol=nv, free=ThPat, values=StTH, lbound= -3, ubound=3, labels=LabThMZm, name="ThMZm") #BMH now give only one lower limit
ThresMZf <-mxAlgebra( expression= cbind(Low%*%ThMZf + BageTHf%x%Age1, Low%*%ThMZf + BageTHf%x%Age2), name="expThresMZf")
ThresMZm <-mxAlgebra( expression= cbind(Low%*%ThMZm + BageTHm%x%Age1, Low%*%ThMZm + BageTHm%x%Age2), name="expThresMZm")
CorMZf <-mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=StCorMZf, labels=LabCorMZf, lbound=-.99, ubound=.99, name="expCorMZf")
CorMZm <-mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=StCorMZm, labels=LabCorMZm, lbound=-.99, ubound=.99, name="expCorMZm")

Tdzf <-mxMatrix( type="Full", nrow=nth, ncol=nv, free=ThPat, values=StTH, lbound= -3, ubound=3, labels=LabThDZf, name="ThDZf") #BMH now give only one lower limit
Tdzm <-mxMatrix( type="Full", nrow=nth, ncol=nv, free=ThPat, values=StTH, lbound= -3, ubound=3, labels=LabThDZm, name="ThDZm")
ThresDZf <-mxAlgebra( expression= cbind(Low%*%ThDZf + BageTHf%x%Age1, Low%*%ThDZf + BageTHf%x%Age2), name="expThresDZf")
ThresDZm <-mxAlgebra( expression= cbind(Low%*%ThDZm + BageTHm%x%Age1, Low%*%ThDZm + BageTHm%x%Age2), name="expThresDZm")
CorDZf <-mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=StCorDZf, labels=LabCorDZf, lbound=-.99, ubound=.99, name="expCorDZf")
CorDZm <-mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=StCorDZm, labels=LabCorDZm, lbound=-.99, ubound=.99, name="expCorDZm")

# Data objects for Multiple Groups
dataMZf <- mxData( observed=mzfData, type="raw" )
dataMZm <- mxData( observed=mzmData, type="raw" )
dataDZf <- mxData( observed=dzfData, type="raw" )
dataDZm <- mxData( observed=dzmData, type="raw" )

expMZf <- mxExpectationNormal( covariance="expCorMZf", means="M", dimnames=selVars, thresholds="expThresMZf" )
expMZm <- mxExpectationNormal( covariance="expCorMZm", means="M", dimnames=selVars, thresholds="expThresMZm" )
expDZf <- mxExpectationNormal( covariance="expCorDZf", means="M", dimnames=selVars, thresholds="expThresDZf" )
expDZm <- mxExpectationNormal( covariance="expCorDZm", means="M", dimnames=selVars, thresholds="expThresDZm" )

#FitFunctions
fun <- mxFitFunctionML()
fitFunction <- mxFitFunctionMultigroup(c("MZf.fitfunction", "MZm.fitfunction", "DZf.fitfunction", "DZm.fitfunction"))

# Combine Groups
modelMZf <- mxModel( obsAge1, obsAge2, Mean, betaAf, Tmzf, inc, ThresMZf, CorMZf, dataMZf, expMZf, fun, name="MZf" )
modelMZm <- mxModel( obsAge1, obsAge2, Mean, betaAm, Tmzm, inc, ThresMZm, CorMZm, dataMZm, expMZm, fun, name="MZm" )
modelDZf <- mxModel( obsAge1, obsAge2, Mean, betaAf, Tdzf, inc, ThresDZf, CorDZf, dataDZf, expDZf, fun, name="DZf" )
modelDZm <- mxModel( obsAge1, obsAge2, Mean, betaAm, Tdzm, inc, ThresDZm, CorDZm, dataDZm, expDZm, fun, name="DZm" )

SatModel <- mxModel( "Sat", modelMZf, modelMZm, modelDZf, modelDZm, fitFunction)

# -------------------------------------------------------------------------------------------------------------------------------
# 1) RUN Saturated Model

SatFit <- mxRun(SatModel, intervals=F)
(SatSumm <- summary(SatFit))

bhuibregtse's picture
Offline
Joined: 03/06/2014 - 22:15
Two changes I am not

Two changes I am not completely confident in since nv=2 instead of nv=1

Changing start values for thresholds to 2 values from 4
(StTh<-c(1.2816, 0.9741) from StTH <-c(1.5, 0, -.8, 1.4 ))

and giving only one lower bound:

(Tmz<-mxMatrix( type="Full", nrow=nth, ncol=nv, free=ThPat, values=1, lbound= -3, ubound=3, labels=LabThMZf, name="ThMZf") from
Tmz <-mxMatrix( type="Full", nrow=nth, ncol=nv, free=ThPat, values=StTH, lbound=c(-3,.001), ubound=3, labels=LabThMZ, name="ThMZ"))

Any ideas are appreciated! Thanks.

tbates's picture
Offline
Joined: 07/31/2009 - 14:25
Replace NA with any number

Your threshold algebra is being made here:

ThresMZf = mxAlgebra(name= "expThresMZf",
    cbind(Low %*% ThMZf + BageTHf %x% Age1, Low %*% ThMZf + BageTHf %x% Age2)
)

Replace the NAs in ThMZf with any legal number and let us know it that solves the issue?

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
thread?

Did you mean to post that comment in another thread?