Login

News

Bivariate Cholesky: error in do.call

Posted on
No user picture. Daxide Joined: 02/15/2016
Forums
Behavioral Genetics Models
Hi! I am a novice, trying to move from univariate to multivariate analysis. So I downloaded the script cholesky.R from Sarah Medland's webpage (https://www.genepi.qimr.edu.au/staff/sarahMe/Workshop13.html), along with the data file. I tried to adapt the script to my data but I ran into an error. So I ran the script as I downloaded it, on the supplied simulated data, without modifying it (see eattachment) and I still encounter the same error, at this point. I do not know what this error means and I do not know why it mentions "what" since it does not exist in the script. I would greatly appreciate your help as I am currently stuck and I do not have any colleagues that use OpenMx.
# To create Labels for Lower Triangular Matrices
> aLabs <- paste("a", do.call(c, sapply(seq(1, nv), function(x){ paste(x:nv, x,sep="") })), sep="")
Error in do.call(c, sapply(seq(1, nv), function(x) { :
'what' must be a character string or a function
Replied on Wed, 02/17/2016 - 11:25
Picture of user. mhunter Joined: 07/31/2009

Thanks for providing all the information we need to replicate the behavior! Sorry that this is causing problems.

Here's some explanation first; the solution comes later. The line


aLabs <- paste("a", do.call(c, sapply(seq(1, nv), function(x){ paste(x:nv, x,sep="") })), sep="")

creates labels for the "a" matrix later. It does this through nested function calls. The outer-most call is paste("a", stuff) where "stuff" is do.call(c, moreStuff). The error is thrown by "stuff", or more precisely the do.call() function. The first argument of the do.call() function is "what". The error message is saying that it doesn't think 'c' is a valid 'what' argument. I think the error must be a Mac vs Windows issue because the line worked on my machine. Hopefully that helps a bit.

The thing produced by this line is


aLabs
[1] "a11" "a21" "a31" "a41" "a51" "a61" "a22" "a32" "a42" "a52" "a62" "a33"
[13] "a43" "a53" "a63" "a44" "a54" "a64" "a55" "a65" "a66"

If the line doesn't work for you, let's just use different code to produce the same output. Here's an example.


aLabs <- paste("a", outer(1:nv, 1:nv, paste, sep="")[lower.tri(matrix(NA, nv, nv), TRUE)], sep="")

A last thing you could try is to replace that line with the same line, but put quotes around "c".


aLabs <- paste("a", do.call("c", sapply(seq(1, nv), function(x){ paste(x:nv, x,sep="") })), sep="")

Log in or register to post comments
Replied on Wed, 02/17/2016 - 11:57
Picture of user. AdminRobK Joined: 01/24/2014

FWIW, I'm not able to reproduce the error either. Under 32-bit Windows, the error doesn't occur in R version 3.1.1 or 2.15.3. It doesn't occur either under x86-64 Linux with R version 3.2.3.

Anyhow, yet another way of getting that line to run would be

aLabs <- paste("a", do.call(base::c, sapply(seq(1, nv), function(x){ paste(x:nv, x,sep="") })), sep="")

Log in or register to post comments
Replied on Thu, 02/18/2016 - 04:33
No user picture. Daxide Joined: 02/15/2016

Thank you very much! That worked (putting the quotes around "c")!
I ran the rest of the script and encountered other problems. I had to replace deprecated functions (I guess the script is old) and that fixed it. However, I encountered one more issue towards the end of the script, when I get this error message:

CholAceFit <- mxRun(CholAceModel)
Error: The references 'MZ.fitfunction' and 'DZ.fitfunction' are unknown in the algebra named 'CholACE.m2LL'
In addition: Warning messages:
1: In is.na(x) : is.na() applied to non-(list or vector) of type 'symbol'
2: In is.na(x) : is.na() applied to non-(list or vector) of type 'symbol'

I guess I am getting this because earlier in the script I had to replace deprecated function "mxAlgebraObjective" (as suggested by an error message) like so:
obj <- mxAlgebraObjective( "m2LL" )
replaced with
obj <- mxFitFunctionAlgebra( "m2LL" )

I suspect I have to modify the following line which appears after the previous line, to include references "MZ.fitfunction" and "DZ.fitfunction", but I have no clue how to do that:
CholAceModel <- mxModel( "CholACE", pars, modelMZ, modelDZ, minus2ll, obj )

To summarize, the concerned lines of code are these (culimating in error message at the last line):

obj <- mxFitFunctionAlgebra("m2LL" )
CholAceModel <- mxModel( "CholACE", pars, modelMZ, modelDZ, minus2ll, obj )
CholAceFit <- mxRun(CholAceModel)

Log in or register to post comments
Replied on Thu, 02/18/2016 - 10:19
Picture of user. AdminRobK Joined: 01/24/2014

Try this code instead:

# Objective objects for Multiple Groups
xpecMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selVars )
xpecDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars )
twinfit <- mxFitFunctionML()

# Combine Groups
pars <- list( pathA, pathC, pathE, covA, covC, covE, covP, matI, invSD, invSDa, invSDe, estVars ) #invSDc,
modelMZ <- mxModel( pars, inclusions, covMZ, dataMZ, xpecMZ, twinfit, name="MZ" )
modelDZ <- mxModel( pars, inclusions, covDZ, dataDZ, xpecDZ, twinfit, name="DZ" )
minus2ll <- mxFitFunctionMultigroup(c("MZ.fitfunction","DZ.fitfunction"))
CholAceModel <- mxModel( "CholACE", pars, modelMZ, modelDZ, minus2ll)

Edit: BTW, I get status code 6 when I run this model on my system. Try re-running it a few times, with CholAceFit <- mxRun(CholAceFit).

Log in or register to post comments
Replied on Thu, 02/18/2016 - 16:30
No user picture. Daxide Joined: 02/15/2016

I was able to get the genetic and environmental influences as stCov. However, I get the following error when trying to compute the correlation matrix, and I do not know what it's caused by?I used the following code:


ACEcorMatrices<-c(
"solve(sqrt(I*A)) %*% A)",
"solve(sqrt(I*C)) %*% C)",
"solve(sqrt(I*E)) %*% E)",
"solve(sqrt(I*V)) %*% V)")
ACEcorLabels=c("CorA", "CorC", "CorE", "CorP")
formatOutputMatrices(CholAceFit ,ACEcorMatrices,ACEcorLabels, Vars, 4)

and this is the error I get:

[1] "Matrix solve(sqrt(I*A)) %*% A)"
Error in print(formatOutputMatrix(evalQuote(matricesList[[k]], fittedModel), :
error in evaluating the argument 'x' in selecting a method for function 'print': Error in parse(text = expstring) : :1:23: unexpected ')'
1: solve(sqrt(I*A)) %*% A)

Log in or register to post comments
Replied on Fri, 02/19/2016 - 11:04
Picture of user. AdminRobK Joined: 01/24/2014

I'm going to guess the error occurs because, in the following,

"solve(sqrt(I*A)) %*% A)",
"solve(sqrt(I*C)) %*% C)",
"solve(sqrt(I*E)) %*% E)",
"solve(sqrt(I*V)) %*% V)"

each line has two open-parentheses but three close-parentheses. Try deleting one of the close-parentheses...probably the last one in each line (I'm not sure what this code is supposed to do)?
Log in or register to post comments
Replied on Fri, 02/19/2016 - 11:25
Picture of user. AdminRobK Joined: 01/24/2014

This code is supposed to convert a covariance into a correlation matrix.

Then shouldn't it be

"solve(sqrt(I*A)) %&% A"

etc., i.e. with the quadratic product operator %&% instead of %*%?

If all you want to do is turn a covariance matrix into a correlation matrix, consider cov2cor().

Log in or register to post comments
Replied on Fri, 02/19/2016 - 15:20
No user picture. Daxide Joined: 02/15/2016

Sorry to bother again! I will use the novice excuse one last time.
The code now seems to work, but only for a variable pair and not for the other. I attach both the code and the dataset.
The first bivariate analysis (Po-VarX) yields meaningful results , but the other (Zic-VarX) yields a strange matrix. E.g.:
CorA1 CorA2
VarX 1.0000 1.0000
Zic 1.0000 1.0000
CorC1 CorC2
VarX 1.0000 -1.0000
Zic -1.0000 1.0000

CorE1 CorE2
VarX 1.0000 -0.2285
Zic -0.2285 1.0000

I have no idea what this is due to as I used the same code and just changed variable names. I also ran univariate on "Zic" and it produced consistent results. I don't even get any error messages.

Log in or register to post comments
Replied on Sat, 02/20/2016 - 17:57
Picture of user. AdminRobK Joined: 01/24/2014

You only have 135 rows of data, and therefore, you're going to be getting pretty noisy estimates of those correlations. To obtain confidence intervals for those correlations, I modified part of your script as follows:

CholAceModel <- mxModel(
"CholACE", pars, modelMZ, modelDZ, minus2ll,
mxAlgebra(cov2cor(A),name="CorA"),
mxAlgebra(cov2cor(C),name="CorC"),
mxAlgebra(cov2cor(E),name="CorE"),
mxAlgebra(rbind(CorA[2,1],CorC[2,1],CorE[2,1]),name="Cors"),
mxCI(c("Cors"))
)
CholAceFit <- mxRun(CholAceModel,intervals=T)
CholAceSumm <- summary(CholAceFit,verbose=T)
CholAceSumm

I get (displaying only the first 4 columns of 'CI details'):

confidence intervals:
lbound estimate ubound note
CholACE.Cors[1,1] 0.9998999 1.0000000 NA !!!
CholACE.Cors[2,1] NA -1.0000000 1.00000000 !!!
CholACE.Cors[3,1] -0.4112075 -0.2284786 -0.01972023

CI details:
parameter value side fit
1 CholACE.Cors[1,1] 1.00000000 upper 1650.391
2 CholACE.Cors[1,1] 0.99989994 lower 1654.233
3 CholACE.Cors[2,1] 1.00000000 upper 1654.233
4 CholACE.Cors[2,1] -1.00000000 lower 1650.391
5 CholACE.Cors[3,1] -0.01972023 upper 1654.233
6 CholACE.Cors[3,1] -0.41120753 lower 1654.233

The upper confidence limit for the A correlation is 1.0, the same as the point estimate; the lower confidence limit is very close to the point estimate, but the change in fit is approximately correct for a 95% profile-likelihood confidence limit (under most circumstances--I won't get into the complication that the point estimate is on the boundary of the parameter space). On the other hand, your estimate for the C correlation has no statistical precision: the confidence interval spans the entire parameter space of [-1,1].

Log in or register to post comments
  Login

News