no. observed statistics 0

8 replies [Last post]
lotgeels's picture
Offline
Joined: 12/22/2009

Hi all,

I'm trying to fit a saturated model to my data (1 binary variable, 5 zygosity groups), and want to estimate 10 thresholds (one for each twin in all zyg groups), and 5 twin correlations (1 per zyg group). The model runs, the right number of parameters are estimated and the estimations look good, but the number of observed statistics is 0, so I get a negative number of degrees of freedom... There are no error messages, and I've looked into different parts of the model, but can't find where it goes wrong. I hope you can help me with this, thanks in advance!

Lot Geels
Biological Psychology
VU University Amsterdam

neale's picture
Offline
Joined: 07/31/2009
Lot Without the script & data

Lot

Without the script & data posted it is difficult to debug this problem. If your own data are not suitable for posting to the website, perhaps use an example from the workshop?

Cheers
Mike

lotgeels's picture
Offline
Joined: 12/22/2009
Hi mike, here's my script: #

Hi mike,

here's my script:

# Original Program: UnivACE_MatrRawOrd_FR.R

setwd("D:/Documents and Settings/lm.geels/Desktop/adolalc")
getwd()

require(foreign)
require(OpenMx)
source("GenEpiHelperFunctions.R")

# Read data from '15_16jarigen.sav'
# Variabels: ntrid, zygdef alcooit1 alcooit2 alcfrq31 alcfrq32 alcfrq3_21 alcfrq3_22 alcpwk1 alcpwk2
# zyg: 1=mzm, 2=dzm, 3=mzf, 4=dzf, 5=dosmf
# Sex: 1=male, 2=female

allVars<- c('ntrid','zygdef','alcooit1','alcooit2','alcfrq31','alcfrq32','alcfrq3_21','alcfrq3_22','alcpwk1','alcpwk2')
groep1516=read.spss("15_16jarigen.sav", use.value.labels = F, to.data.frame =T)

head(groep1516)

summary(groep1516$alcooit1)
summary(groep1516$alcooit2)
table(groep1516$alcooit1,useNA="ifany")
table(groep1516$alcooit2,useNA="ifany")

nv <- 1
ntv <- nv*2
Vars <-('alcooit')
selVars <- c('alcooit1' , 'alcooit2')
mzmdata = subset(groep1516, zygdef==1, selVars)
dzmdata = subset(groep1516, zygdef==2, selVars)
mzfdata = subset(groep1516, zygdef==3, selVars)
dzfdata = subset(groep1516, zygdef==4, selVars)
dosmfdata = subset(groep1516, zygdef==5, selVars)

ls()

# Print Descriptive Statistics
# ---------------------------------------------------------------------
table(mzmdata$alcooit1,useNA="ifany")
table(mzmdata$alcooit2,useNA="ifany")

table(dzmdata$alcooit1,useNA="ifany")
table(dzmdata$alcooit2,useNA="ifany")

table(mzfdata$alcooit1,useNA="ifany")
table(mzfdata$alcooit2,useNA="ifany")

table(dzfdata$alcooit1,useNA="ifany")
table(dzfdata$alcooit2,useNA="ifany")

table(dosmfdata$alcooit1,useNA="ifany")
table(dosmfdata$alcooit2,useNA="ifany")

# Specify and Run Saturated Model (Tetrachoric correlations) with RawData and Matrix-style Input
# -----------------------------------------------------------------------
twinSatModel <- mxModel("twinSat",
mxModel("mzm",
# Matrix & Algebra for expected means vector (SND), Thresholds and correlation
mxMatrix( type="Zero", nrow=1, ncol=nv, name="M" ),
mxAlgebra( expression= cbind(M,M), name="expMean" ),
mxMatrix(type="Full", nrow=1, ncol=ntv, free=TRUE, values=-1.5, name="expThremzm", dimnames=list('th1',selVars) ),
mxMatrix(type="Stand", nrow=2, ncol=2, free=T, values=.5, lbound=-.99, ubound=.99, name="expCormzm"),
mxData(mzmdata, type="raw"),
mxFIMLObjective( covariance="expCormzm", means="expMean", dimnames=selVars, thresholds="expThremzm" )),

mxModel("dzm",
# Matrix & Algebra for expected means vector (SND), Thresholds and correlation
mxMatrix( type="Zero", nrow=1, ncol=nv, name="M" ),
mxAlgebra( expression= cbind(M,M), name="expMean" ),
mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=-1.5, name="expThredzm", dimnames=list('th1',selVars)),
mxMatrix(type="Stand", nrow=2, ncol=2, free=T, values=.5, lbound=-.99, ubound=.99, name="expCordzm"),
mxData(dzmdata, type="raw"),
mxFIMLObjective( covariance="expCordzm", means="expMean", dimnames=selVars, thresholds="expThredzm" )),

mxModel("mzf",
# Matrix & Algebra for expected means vector (SND), Thresholds and correlation
mxMatrix( type="Zero", nrow=1, ncol=nv, name="M" ),
mxAlgebra( expression= cbind(M,M), name="expMean" ),
mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=-1.5, name="expThremzf", dimnames=list('th1',selVars)),
mxMatrix(type="Stand", nrow=2, ncol=2, free=T, values=.7, lbound=-.99, ubound=.99, name="expCormzf"),
mxData(mzfdata, type="raw"),
mxFIMLObjective( covariance="expCormzf", means="expMean", dimnames=selVars, thresholds="expThremzf" )),

mxModel("dzf",
# Matrix & Algebra for expected means vector (SND), Thresholds and correlation
mxMatrix( type="Zero", nrow=1, ncol=nv, name="M" ),
mxAlgebra( expression= cbind(M,M), name="expMean" ),
mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=-1.7, name="expThredzf", dimnames=list('th1',selVars)),
mxMatrix(type="Stand", nrow=2, ncol=2, free=T, values=.5, lbound=-.99, ubound=.99, name="expCordzf"),
mxData(dzfdata, type="raw"),
mxFIMLObjective( covariance="expCordzf", means="expMean", dimnames=selVars, thresholds="expThredzf" )),

mxModel("dosmf",
# Matrix & Algebra for expected means vector (SND), Thresholds and correlation
mxMatrix( type="Zero", nrow=1, ncol=nv, name="M" ),
mxAlgebra( expression= cbind(M,M), name="expMean" ),
mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=-1.7, name="expThredosmf", dimnames=list('th1',selVars)),
mxMatrix(type="Stand", nrow=2, ncol=2, free=T, values=.5, lbound=-.99, ubound=.99, name="expCordosmf"),
mxData(dosmfdata, type="raw"),
mxFIMLObjective( covariance="expCordosmf", means="expMean", dimnames=selVars, thresholds="expThredosmf" )),
mxAlgebra( mzm.objective + dzm.objective + mzf.objective + dzf.objective + dosmf.objective, name="-2sumll" ),
mxAlgebraObjective("-2sumll")

)
twinSatFit <- mxRun(twinSatModel)
twinSatSumm <- summary(twinSatFit)
twinSatSumm

I don't think my supervisors will appreciate it if I put the data online, but maybe I can email them to you if you need them?

Lot

neale's picture
Offline
Joined: 07/31/2009
Yes, by all means email me

Yes, by all means email me the data. However, I suspect that what you are seeing is a bug in which the number of statistics is not calculated correctly with an mxAlgebraObjective. Can you post the results of summary(twinSatFit) to illustrate?

lotgeels's picture
Offline
Joined: 12/22/2009
of course! here it

of course! here it is:

$parameters
name matrix row col Estimate Std.Error
1 mzm.expThremzm th1 alcooit1 -1.1448220 0.06942189
2 mzm.expThremzm th1 alcooit2 -1.1977616 0.07237408
3 mzm.expCormzm alcooit1 alcooit2 0.9223936 0.02712722
4 dzm.expThredzm th1 alcooit1 -1.3544697 0.08880781
5 dzm.expThredzm th1 alcooit2 -1.3676424 0.08930811
6 dzm.expCordzm alcooit1 alcooit2 0.6877297 0.09314534
7 mzf.expThremzf th1 alcooit1 -1.2179928 0.06120047
8 mzf.expThremzf th1 alcooit2 -1.1834995 0.06136057
9 mzf.expCormzf alcooit1 alcooit2 0.9496579 0.01782034
10 dzf.expThredzf th1 alcooit1 -1.5202752 0.08934453
11 dzf.expThredzf th1 alcooit2 -1.3797419 0.07745697
12 dzf.expCordzf alcooit1 alcooit2 0.6528512 0.09502078
13 dosmf.expThredosmf th1 alcooit1 -1.2377992 0.06080836
14 dosmf.expThredosmf th1 alcooit2 -1.3489826 0.06004347
15 dosmf.expCordosmf alcooit1 alcooit2 0.7187523 0.05144278

$estimatedParameters
[1] 15

$observedStatistics
[1] 0

$degreesOfFreedom
[1] -15

$Minus2LogLikelihood
[1] 1630.903

I just discussed this with a colleague (Sanja Franic) and she told me that when she simulated some data and fit a multivariate model to it she had the same problem (good estimates but 0 observed statistics).

thanks!

Lot

mspiegel's picture
Offline
Joined: 07/31/2009
Could you run mxVersion() and

Could you run mxVersion() and make sure you're running 0.2.10? There was an overhaul in the summary function recently.

lotgeels's picture
Offline
Joined: 12/22/2009
this is what I

this is what I get:

mxVersion()
[1] "0.2.3-1006"

neale's picture
Offline
Joined: 07/31/2009
Looks like you need to

Looks like you need to upgrade... At the R prompt

source('http://openmx.psyc.virginia.edu/getOpenMx.R')

lotgeels's picture
Offline
Joined: 12/22/2009
that works, thanks very

that works, thanks very much!

Lot