Interpret >1 and negative results from variance parametrisation of twin model
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,
Log in or register to post comments
In reply to The simple answer is that, by AdminRobK
Thank you so much for your
In that case, what would be your recommendation for interpreting the etiological correlations?
Thank you so much for your patience.
With all good wishes
Log in or register to post comments
Big genetic correlation > 1
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).
Log in or register to post comments
In reply to Big genetic correlation > 1 by AdminNeale
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
Log in or register to post comments
In reply to Dear professor Neale, by JuanJMV
Only one constraint needed
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
MZconstraint <- mxConstraint(abs(rA[2,1]) < 1)
But in the end, I think there's something wrong with your VA parameters.
Log in or register to post comments
In reply to Only one constraint needed by AdminNeale
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.
Log in or register to post comments
In reply to Dear Michael, by JuanJMV
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.
Log in or register to post comments