Hi All,

I am using the new parametrization for a bivariate model, as suggested in "Type I Error Rates and Parameter Bias in Multivariate Behavioral Genetic Models".

The results (see below) include an rA correlation that is higher than one.

Am I doing something wrong or I have to interpret the rA correlation as 1?

The phenotypic correlation is .2. Is it be correct to say that is due to 1, -.96 and .2 for A, D and E respectively?

There are also negative values: there is no way to interpret these results, right?

Thank you so much in advance

$V mxAlgebra 'V' $formula: VA + VD + VE $result: [,1] [,2] [1,] 0.26015570 0.08651586 [2,] 0.08651586 0.64434885 dimnames: NULL $iSD mxAlgebra 'iSD' $formula: solve(sqrt(I * V)) $result: [,1] [,2] [1,] 1.9605744 0.0000000 [2,] 0.0000000 1.2457746 dimnames: NULL $rA mxAlgebra 'rA' $formula: solve(sqrt(I * VA)) %&% VA $result: [,1] [,2] [1,] 1.0000000 3.1877784 [2,] 3.1877784 1.0000000 dimnames: NULL $rD mxAlgebra 'rD' $formula: solve(sqrt(I * VD)) %&% VD $result: [,1] [,2] [1,] 1.00000000 -0.95732896 [2,] -0.95732896 1.00000000 dimnames: NULL $rE mxAlgebra 'rE' $formula: solve(sqrt(I * VE)) %&% VE $result: [,1] [,2] [1,] 1.00000000 0.20271077 [2,] 0.20271077 1.00000000 dimnames: NULL $US mxAlgebra 'US' $formula: cbind(VA, VD, VE, VA/V, VD/V, VE/V) $result: VA VA VD VD VE VE SA SA SD SD US 0.036015032 0.130102084 0.053615331 -0.096971671 0.170525336 0.053385447 0.13843645 1.503794611 0.2060894 -1.12085427 US 0.130102084 0.046249649 -0.096971671 0.191371989 0.053385447 0.406727210 1.50379461 0.071777344 -1.1208543 0.29700059 SE SE US 0.65547415 0.61705966 US 0.61705966 0.63122206

Hope that you are all well during this difficult situation

The simple answer is that, under the direct-symmetric parameterization, the

A,C, andEmatrices are not necessarily covariance matrices, and therefore when standardized, are not necessarily interpretable as correlation matrices.Thank you so much for your response. Super useful as always!

In that case, what would be your recommendation for interpreting the etiological correlations?

Thank you so much for your patience.

With all good wishes

Hi Juan

The large genetic correlation, rg greater than 1, is probably just some sampling error around an estimate of rg that is close to 1. Note that the genetic variances are small so rg may have confidence intervals that overlap 1.0 - perhaps by a lot. You could try a constraint to keep the correlation within the -1 to 1 bound and see how much difference it makes to the -2lnL.

To interpret it, I would inspect the correlations in the data. Likely what is happening is that the MZ cross-correlation between the traits is a bit bigger than expected, given how big the within-variable across MZ correlations are. So if looking at a correlation matrix that goes T1V1 T1V2 T2V1 T2V2 for columns and rows, it may be that rT1V1,T2V2 and/or rT1V2,T2V1 are a bit bigger than rT1V1,rT2V1 and rT1V2,rT2V2 (check their geometric mean).

Dear professor Neale,

Thank you so much for your response.

I have tried to constraint the rG <1 but I cannot make the model to work.

Here the lines that I have used and the error message.

MZconstraint <- mxConstraint(rA[2,1] < 1)

MZconstraint2 <- mxConstraint(rA[1,2] < 1)

modelADECONS <- mxModel( fitADE, name="ADECONS", MZconstraint,MZconstraint2)

fitADECONS <- mxTryHard( modelADECONS, intervals=F )

Warning message:

In mxTryHard(modelADECONS, intervals = F) :

Error : The following error occurred while evaluating the subexpression 'solve(sqrt(ADECONS.I * ADECONS.VA))' during the evaluation of 'rA' in model 'ADECONS' : Lapack routine dgesv: system is exactly singular: U[1,1] = 0

Regarding the correlations in the data, I think rT1V1,T2V2 are bigger than rT1V1,rT2V1 so probably that is the problem. I have calculated the geometric mean and there are differences.

Here the outputs in case useful.

$expCovMZ

mxAlgebra 'expCovMZ'

$formula: rbind(cbind(V, cMZ), cbind(t(cMZ), V))

$result:

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

[1,] 0.260155920 0.086515864 0.089630573 0.033130345

[2,] 0.086515864 0.644348796 0.033130345 0.237621696

[3,] 0.089630573 0.033130345 0.260155920 0.086515864

[4,] 0.033130345 0.237621696 0.086515864 0.644348796

dimnames: NULL

> mxEval(cov2cor(expCovMZ), fitADE$MZ, T)

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

[1,] 1.000000000 0.211309187 0.344526362 0.080918643

[2,] 0.211309187 1.000000000 0.080918643 0.368778055

[3,] 0.344526362 0.080918643 1.000000000 0.211309187

[4,] 0.080918643 0.368778055 0.211309187 1.000000000

Thank you so much in advance

Hi Juan

It is not necessary to add both constraints - they do the same thing because rA is symmetric. Just one is better (the two are perfectly correlated).

Second, VA for variable 1 is

set to, starting at, or fixed to zero, so when it tries to calculate rA, (remembering a correlation is a covariance divided b.y the square root of the product of the variances). In the present case with U[1,1]=0, it is trying to divide by zero, which isn't going to work.Third, I would make the constraint a bit more effective (prevent <-1 as well as >1) with

But in the end, I think there's something wrong with your VA parameters.

Dear Michael,

Thank you so much for your response. I have tried using the lines that you suggested but I got the same error message.

Error : The following error occurred while evaluating the subexpression 'solve(sqrt(ADECONS.I * ADECONS.VA))' during the evaluation of 'rA' in model 'ADECONS' : Lapack routine dgesv: system is exactly singular: U[1,1] = 0

I have also tried changing the starting values for VA but It does not work either.

Here you can see the values for VA from the ADE model without constraints:

$VA

SymmMatrix 'VA'

$labels

[,1] [,2]

[1,] "VA11" "VA21"

[2,] "VA21" "VA22"

$values

[,1] [,2]

[1,] 0.036015727 0.130102805

[2,] 0.130102805 0.046250591

Could it be because VA[2,1] is higher than VA[1,1] and VA[2,2]?

Thank you so much in advance for your help, I am totally stuck with this.

With all good wishes.

Hi again,

I have been trying to fix it but I cannot make it work.

I have checked the starting values for VA but I am not sure why VA is set to, starting at, or fixed to zero.

Here the lines for VA:

svPa <- .2 # start value for path coefficient

covA <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=valDiag(svPa,nv), label=labLower("VA",nv), name="VA" )

So, if I am not wrong, VA should start at .2. Am I doing something wrong?

I have also tried to find U[1,1] = 0 in the model but I cannot figure out where is that allocated.

Thank you so much for your help with this.

With all good wishes.