# Bivariate Cholesky: error in do.call

Posted on

Daxide
Joined: 02/15/2016

Attachment | Size |
---|---|

cholesky (1).R | 6.2 KB |

newsimulatedData.csv | 166.97 KB |

Forums

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

# 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

## Weird. Try this.

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

In reply to Weird. Try this. by mhunter

## R version?

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

In reply to Weird. Try this. by mhunter

## Thank you very much! That

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

In reply to Thank you very much! That by Daxide

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

In reply to Try this code instead by AdminRobK

## Ok! Thanks so much, that code

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

Log in or register to post comments

In reply to Ok! Thanks so much, that code by Daxide

## correlation matrix

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)":1:23: unexpected ')'

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: solve(sqrt(I*A)) %*% A)

Log in or register to post comments

In reply to correlation matrix by Daxide

## Unbalanced parens

"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

In reply to Unbalanced parens by AdminRobK

## Yes you are right...I had not

Log in or register to post comments

In reply to Yes you are right...I had not by Daxide

## cov2cor

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

In reply to cov2cor by AdminRobK

## Thanks, I had spotted that

Log in or register to post comments

In reply to Thanks, I had spotted that by Daxide

## You're welcome.

Log in or register to post comments

In reply to You're welcome. by AdminRobK

## code working (kinda)

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

In reply to code working (kinda) by Daxide

## I think those results are right

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

In reply to I think those results are right by AdminRobK

## Ok thank you! Too bad I

Log in or register to post comments