# Bivariate Cholesky: error in do.call

15 posts / 0 new
Offline
Joined: 02/15/2016 - 14:08
Bivariate Cholesky: error in do.call
AttachmentSize
6.2 KB
166.97 KB

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

Offline
Joined: 07/31/2009 - 15:26
Weird. Try this.

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="")
Offline
Joined: 01/24/2014 - 12:15
R version?

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="")
Offline
Joined: 02/15/2016 - 14:08
Thank you very much! That

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'
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)
Offline
Joined: 01/24/2014 - 12:15

# 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).

Offline
Joined: 02/15/2016 - 14:08
Ok! Thanks so much, that code

Ok! Thanks so much, that code did the trick! I also had to add a path to GenEpiHelperFunctions.R:

source("http://www.vipbg.vcu.edu/~vipbg/Tc24/GenEpiHelperFunctions.R")

Offline
Joined: 02/15/2016 - 14:08
correlation matrix

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(IA)) %% 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(IA)) %% A)

Offline
Joined: 01/24/2014 - 12:15
Unbalanced parens

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)?

Offline
Joined: 02/15/2016 - 14:08
Yes you are right...I had not

Yes you are right...I had not noticed that. Now it works! This code is supposed to convert a covariance into a correlation matrix.

Offline
Joined: 01/24/2014 - 12:15
cov2cor
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().

Offline
Joined: 02/15/2016 - 14:08

Thanks, I had spotted that error and was gonna add that. Thanks for pointing it out though. Your help has been invaluable.

Offline
Joined: 01/24/2014 - 12:15
You're welcome.

You're welcome.

Offline
Joined: 02/15/2016 - 14:08
code working (kinda)

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.

Offline
Joined: 01/24/2014 - 12:15
I think those results are right

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].

Offline
Joined: 02/15/2016 - 14:08
Ok thank you! Too bad I

Ok thank you! Too bad I cannot collect more data!