Login

News

no. observed statistics 0

Posted on
Picture of user. lotgeels Joined: 12/22/2009
Forums
Behavioral Genetics Models
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

Replied on Mon, 03/15/2010 - 10:10
Picture of user. neale Joined: 07/31/2009

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

Log in or register to post comments
Replied on Mon, 03/15/2010 - 10:38
Picture of user. lotgeels Joined: 12/22/2009

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

Log in or register to post comments
Replied on Mon, 03/15/2010 - 11:09
Picture of user. neale Joined: 07/31/2009

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?
Log in or register to post comments
Replied on Mon, 03/15/2010 - 11:55
Picture of user. lotgeels Joined: 12/22/2009

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

Log in or register to post comments
  Login

News