Welcome to the Documentation forum

Yes, I know that I can make up my own, or copy and paste from the documentation to get the equivalent, but having a working example to run always helps.

Bill

	likelihood <- model@output$Minus2LogLikelihood
	saturated <- model@output$SaturatedLikelihood
	chi <- likelihood - saturated
	DoF <- retval$degreesOfFreedom
	retval[['AIC']] <- likelihood - 2 * DoF
	retval[['BIC']] <- 0.5 * (likelihood - DoF * log(data@numObs))
	rmseaSquared <- (chi / DoF - 1) / data@numObs
	if (rmseaSquared < 0) {
		retval[['RMSEA']] <- 0
	} else {
		retval[['RMSEA']] <- sqrt(rmseaSquared)
	}

A few questions:
- Is chi [likelihood - saturated] or should it be [saturated - likelihood] ?
- Degrees of freedom does not currently take into account non-linear constraints. To compute D.o.F. I should be subtracting the number of non-linear constraints?

Congratulations on getting underway!!

I started a bit late because when I upgraded R to 2.9.1 from 2.8 and updated my packages, it caused SweaveListingUtils, a package I use, to stop functioning. It is near the start of the term and I am constructing slides with LaTeX,SWeave, R, Beamer, and SWeaveListingUtils. Suddenly nothing would run.
It was a crisis. It took a lot of troubleshooting to figure out what was happening.

The temporary resolution was to simply kill the "updated" SWeaveListingUtils and go back to the version I was using with 2.8. It took me a while to figure out what was going on, and by then I was due to leave on vacation, from which I have just returned.

I loaded up the system, and started in on the examples.

Although there is some excellent potential there, my sense is that these tutorials will not work as currently set up, because they lack the appropriate conceptual setup, AND they don't work in the way that most people are used to having introductory R examples work.

Let me tell you what I did. I installed (worked great!) powered up R and said

>library(OpenMx)
Loading required package: snowfall
Loading required package: snow
snowfall 1.70 initialized: sequential execution, one CPU.

Notice, I'm showing the output. In my opinion, examples should ALWAYS show the output. If you are not developing your docs in SWeave and LaTeX, you should be, because that way, you can immediately embed output in your examples *and* see if any of them are erroneous.

I then started the examples. I did what every R user is accustomed to --- I copied in the example by cutting and pasting:

> mxMatrix(
+ type="Full",
+ nrow=3,
+ ncol=1,
+ values=c(1,2,3),
+ name='A'
+ )
FullMatrix 'A'

Labels matrix: No labels assigned.

Values matrix:
[,1]
[1,] 1
[2,] 2
[3,] 3

Free matrix: No free parameters.

Lbound matrix: No lower bounds assigned.

Ubound matrix: No upper bounds assigned.

> mxMatrix(
+ type="Full",
+ nrow=3,
+ ncol=1,
+ values=c(1,2,3),
+ name='A'
+ )
FullMatrix 'A'

Labels matrix: No labels assigned.

Values matrix:
[,1]
[1,] 1
[2,] 2
[3,] 3

Free matrix: No free parameters.

Lbound matrix: No lower bounds assigned.

Ubound matrix: No upper bounds assigned.

So far, so good --- apparently. It was not clear from the documentation what was supposed to happen, because no output was shown.

I had no idea what some of the output meant, but I continued on.

> mxAlgebra(
+ A %*% t(A),
+ name='q1'
+ )
mxAlgebra 'q1'
formula: A %*% t(A)
result:
<0 x 0 matrix>
dimnames: NULL
> mxAlgebra(
+ t(A) %*% A,
+ name='q2'
+ )
mxAlgebra 'q2'
formula: t(A) %*% A
result:
<0 x 0 matrix>
dimnames: NULL
> mxAlgebra(
+ A * A,
+ name='q3'
+ )
mxAlgebra 'q3'
formula: A * A
result:
<0 x 0 matrix>
dimnames: NULL
> mxAlgebra(
+ A + B,
+ name='q4'
+ )
mxAlgebra 'q4'
formula: A + B
result:
<0 x 0 matrix>
dimnames: NULL
>

Note that it tells me that the result is a 0x0 matrix!

This was very disconcerting, and I went back over what I had done several times trying to detect and error. I could find none. I spent a lot of time.

In my opinion, this example puts the cart WAYYY before the horse.

What is a structural model? It is a statement that the data (or at least moments of the data) is a function of a vector of parameters.

Sigma = Sigma.Implied

where Sigma.Implied = M(theta)

Sigma can be large. If theta is small and M is a function with some interpretability, then we may have learned something by fitting M() to Sigma.

Often M is a matrix product, involving matrices with free parameters in them. So we need a system to tell us (a)
what the matrices are, (b) where the free parameters are in the matrices, and (c) how to combine them functionally to generate M().

In other words, you need a lot more in the introduction for it to make sense to the typical statistician. And it should be more conceptual.

Finally I realized that you actually didn't intend me to run this code at all -- you were just throwing some code snippets in to illustrate some ideas. Imagine how that made me feel (dumb) after I wasted quite a bit of time double checking. (I do things sequentially...)

So I restarted R and moved on to the next example. I copied in the code, and here's what happened:

> algebraExercises <- mxModel(
+ mxMatrix(type="Full", values=c(1,2,3), nrow=3, ncol=1, name='A'),
+ mxMatrix(type="Full", values=c(1,2,3), nrow=3, ncol=1, name='B'),
+ mxAlgebra(A%*%t(A), name='q1'),
+ mxAlgebra(t(A)%*%A, name='q2'),
+ mxAlgebra(A*A, name='q3'),
+ mxAlgebra(A+B, name='q4'))
>
> answers <- mxRun(algebraExercises)
Running untitled1
> answers@algebras
$q1
mxAlgebra 'q1'
formula: A %*% t(A)
result:
[,1] [,2] [,3]
[1,] 1 2 3
[2,] 2 4 6
[3,] 3 6 9
dimnames: NULL

$q2
mxAlgebra 'q2'
formula: t(A) %*% A
result:
[,1]
[1,] 14
dimnames: NULL

$q3
mxAlgebra 'q3'
formula: A * A
result:
[,1]
[1,] 1
[2,] 4
[3,] 9
dimnames: NULL

$q4
mxAlgebra 'q4'
formula: A + B
result:
[,1]
[1,] 2
[2,] 4
[3,] 6
dimnames: NULL

> result
[[1]]
[,1] [,2] [,3]
[1,] 1 2 3
[2,] 2 4 6
[3,] 3 6 9

[[2]]
[,1]
[1,] 14

[[3]]
[,1]
[1,] 1
[2,] 4
[3,] 9

[[4]]
[,1]
[1,] 2
[2,] 4
[3,] 6

Great.

But you missed a significant opportunity in that little example. I had a fundamental question that the tutorial example *might* have answered. Are Mx matrices entered columnwise or rowwise by default, as in R? I loaded up the referenc manual and found that indeed they are, and that the mxMatrix command has a byrow option just like the R matrix command. Great. However, if you had included a slightly more ambitious example, you could have pointed that out in passing and smoothed the transition, while saving the reader some time.

I move on to the Optimization Script example. I copied in the code, and here is what happened:

> bivCorModel <- mxModel("bivCor",
+ mxMatrix(
+ type="Full",
+ nrow=1,
+ ncol=2,
+ free=True,
+ values=c(0,0),
+ dimnames=list(NULL, selVars),
+ name="expMean"
+ ),
+ mxMatrix(
+ type="Full",
+ nrow=2,
+ ncol=2,
+ free=c(T,T,F,T),
+ values=c(1,.2,0,1),
+ dimnames=list(selVars, selVars),
+ name="Chol"
+ ),
+ mxAlgebra(
+ expression=Chol %*% t(Chol),
+ name="expCov",
+ dimnames=list(selVars, selVars)
+ ),
+ mxData(
+ observed=testData,
+ type="raw"
+ ),
+ mxFIMLObjective(
+ covariance="expCov",
+ means="expMean")
+ )
Error in matrixCheckErrors(type, values, free, labels, lbound, ubound, :
object 'True' not found

Note there is an error. This could not have happened had you developed your docs with Sweave and LaTeX, could it?

Anyway, we are off and running Steve, congratulations!!

--Jim

p13. optimization script
EC <- bivCorFit[[’expCov’]]@values
should have @results instead of @values

p 15, script has data initially defined as MZ and DZ, but then changes names mid-script

p 16, script for MZ model has missing comma
mxMatrix(
type="Lower",
nrow=2,
ncol=2,
free=TRUE,

Section on # Generate Saturated Model Comparison Output
this line:
LL_Sub2 <- mxEval(objective, twinSatFitSub1)
should be:
LL_Sub2 <- mxEval(objective, twinSatFitSub2)

btw, is there any way of adding line numbering to R scripts, so one can easily refer to a line of code?

It states there that to equate variances across twin order and zygosity, we can use a statement such as:
mxMatrix("Lower", 2, 2, T, .5, labels= c("var","MZcov","var"), name="CholMZ"),
mxModel("DZ",
mxMatrix("Full", 1, 2, T, 0, "mean", name="expMeanDZ"),
mxMatrix("Lower", 2, 2, T, .5, labels= c("var","DZcov","var"), name="CholDZ"))

However, this seems to assume that the 3 entries on the Cholesky matrix correspond to variances and covariance, but they don't, because the ultimate covariance matrix from Cholesky multiplication of lower matrix [a b c] is:
[ a^2 ab
ab b^2 + c^2]

I'd be curious to know how you can equate covariances - but it's clear that the script above is wrong - as I discovered when I ran it and obtained an enormously significant different from the saturated model.

1) In mxModelMZ:
mxFIMLObjective("twinACE.expCovMZ", "twinACE.expMean, selVars")),
should be:
mxFIMLObjective("twinACE.expCovMZ", "twinACE.expMean", selVars)),

2) Script doesn't run for me unless the dimnames (ie selVars) are also provided for the parallel mxFIMLObjective for DZs

3) mxAlgebra statement:
mxAlgebra(rbind(cbind(A+C+E, .5%x%A+C), cbind(.5%x%A+C , A+C+E)), name="expCovDZ"),

should have * rather than x