I am unable to figure out how to fix this. Check out this error:
Running Saturated DoC with 26 parameters Error in summary.MxModel(out, refModels = mxRefModels(out, run = T)) : Error : The job for model 'Saturated DoC' exited abnormally with the error message: fit is not finite (The continuous part of the model implied covariance (loc1) is not positive definite in data 'Saturated DZ.data' row 879. Detail: covariance = matrix(c( # 2x2 -4.20388567580781, 0.587702305205191 , 0.587702305205191, 0.955870642918316), byrow=TRUE, nrow=2, ncol=2) ) In addition: Warning message: In model 'Saturated DoC' Optimizer returned a non-zero status code 10. Starting values are not feasible. Consider mxTryHard() >
The expected cov for that particular dz model is full of zeros, which is what we want, no?
> mxGetExpected(satur$Saturated$`Saturated DZ`, "covariances") obese1 ht1 obese2 ht2 obese1 1 0 0 0 ht1 0 0 1 1 obese2 0 1 0 0 ht2 0 1 0 1 >
Finally, the eigen vectors become negative in one spot:
> mxGetExpected(satur$Saturated$`Saturated DZ`, "covariances") |> eigen() eigen() decomposition $values [1] 1.80194 1.00000 0.44504 -1.24698 $vectors [,1] [,2] [,3] [,4] [1,] 0.00000 1 0.00000 0.00000 [2,] -0.59101 0 -0.32799 0.73698 [3,] -0.32799 0 -0.73698 -0.59101 [4,] -0.73698 0 0.59101 -0.32799
Here is the MWE:
library(umx) library(dplyr) library(tidyr) options(mxByrow = TRUE) # IK sorry # Generating data data(twinData) glimpse(twinData) # Make "obese" variable with ~20% subjects categorised as overweight and 30% as obese obesityLevels = c('normal', 'overweight', 'obese', 'obese1') cutPoints = quantile(twinData[, "bmi1"], probs = c(.1, .2, .3), na.rm = TRUE) twinData$obese1 = cut(twinData$bmi1, breaks = c(-Inf, cutPoints, Inf), labels = obesityLevels) twinData$obese2 = cut(twinData$bmi2, breaks = c(-Inf, cutPoints, Inf), labels = obesityLevels) # Step 2: Make the ordinal variables into umxFactors (ordered, with the levels found in the data) selVars = c("obese1", "obese2") twinData <- mutate(twinData, obese1 = factor(obese1, ordered = T), obese2 = factor(obese2, ordered = T), ) mzData = subset(twinData, zygosity %in% c("MZFF", "MZMM")) dzData = subset(twinData, zygosity %in% c("DZFF", "DZMM")) pheno = c("obese", "ht") name = "DoC" sep = "" covar=NULL vnames = tvars(c(pheno), sep = sep) # Find ordinal variables colTypes = umx_is_ordered(xmu_extract_column(mzData, vnames), summaryObject= TRUE) if(colTypes$nOrdVars > 0){ ty = umxThresholdMatrix(rbind(mzData, dzData), fullVarNames = colTypes$ordVarNames, sep = sep, method="Mehta") ordinalPresent = TRUE } xmu_twin_check( selDVs = c(pheno), sep = sep, dzData = dzData, mzData = mzData, enforceSep = TRUE, nSib = 2, optimizer = "CSOLNP" ) top = mxModel("top", mxMatrix(name = "I", type = "Iden", nrow = 2, ncol = 2), umxMatrixFree("B", type = "Full", nrow = 2, ncol = 2, labels = c( NA, "g2", "g1", NA ) ), umxMatrix("psi_a", type = "Symm", ncol = 2, nrow = 2, byrow = TRUE, labels = c( "ax2", "ra", "ay2" ), free = c( T, T, T ), values = c( 0.1, 0.05, 0.1 ), ), umxMatrixFree("psi_c", type = "Symm", ncol = 2, nrow = 2, labels = c( "cx2", "rc", "cy2" ), ), umxMatrixFree("psi_e", type = "Symm", ncol = 2, nrow = 2, labels = c( "ex2", "re", "ey2" ), ), umxMatrix("lamba", type = "Diag", nrow = 2, ncol = 2, byrow = TRUE, free = c(F, F), labels = c("σ1", "σ2"), values = c(1, 1) ), mxAlgebra("A", expression = lamba %&% (solve(I - B) %&% psi_a)), mxAlgebra("C", expression = lamba %&% (solve(I - B) %&% psi_c)), mxAlgebra("E", expression = lamba %&% (solve(I - B) %&% psi_e)), mxAlgebra("CovMZ", expression = rbind( cbind(A + C + E, A + C), cbind(A + C, A + C + E) )), mxAlgebra(name = "CovDZ", expression = rbind( cbind(A + C + E, 0.5 %x% A + C), cbind(0.5 %x% A + C, A + C + E) )), mxAlgebra("VC", expression = cbind(A, C, E, A / (A + C + E), C / (A + C + E), E / (A + C + E)), dimnames = list( rep("VC", 2), rep(c("A", "C", "E", "SA", "SC", "SE"), each = 2) ) ) , # THE SECTION BELOW IS TO CALCULATE THE MEANS, THE OVERCOMPLICATED APPROACH IS LEFT OVER OF ADDING COVARS # BUT SHOULDNT AFFECT THE ISSUE # create algebra for expected mean & threshold matrices mxMatrix( type="Full", nrow=1, ncol=2, free= TRUE, labels=c("mean_Ph1","mean_Ph2"), name="meanT1" ), # in mz twins, prs_twin1 == prs_twin2 # so, twin2 in mz pairs does not have the prs variable mxMatrix( type="Full", nrow=1, ncol=2, free=TRUE, labels=c("mean_Ph1","mean_Ph2"), name="meanT2MZ"), # in dz twins, prs_twin1 != prs_twin2 mxMatrix( type="Full", nrow=1, ncol=2, free=TRUE, labels=c("mean_Ph1","mean_Ph2"), name="meanT2DZ" ) ) # This is leftovers of adding covars, but shouldn't affect the issue expMeanMZ <- mxAlgebra("expMeanMZ", expression = cbind(top.meanT1 , top.meanT2MZ )) expMeanDZ <- mxAlgebra("expMeanDZ", expression = cbind(top.meanT1 , top.meanT2DZ )) # Adding thresholds top = top+ty expMZ = mxExpectationNormal("top.CovMZ", means = "expMeanMZ", dimnames = vnames, thresholds="top.threshMat", threshnames = colTypes$ordVarNames ) expDZ =mxExpectationNormal("top.CovDZ", means = "expMeanDZ", dimnames = vnames, thresholds="top.threshMat", threshnames = colTypes$ordVarNames ) MZ = mxModel("MZ", expMeanMZ, expMZ, mxFitFunctionML()) DZ = mxModel("DZ", expMeanDZ, expDZ, mxFitFunctionML()) # Adding data MZ = MZ + mxData(mzData, type = "raw") DZ = DZ + mxData(dzData, type = "raw") model = mxModel(name, top, MZ, DZ, mxFitFunctionMultigroup(c("MZ","DZ") )) # Dropping g2 and re for identification DoC model = umxModify(model, update = c("g2", "re"), autoRun = F) model<- omxSetParameters(model, labels = c("obese_dev1","obese_dev2","obese_dev3"), free = c(F,F,T), value = c(0.1,0.1,0.2)) out <- mxRun(model) summary(out, refModels = mxRefModels(out, run = T)) satur <- mxRefModels(out) mxGetExpected(satur$Saturated$`Saturated DZ`, "covariances") |> eigen()
As I hit the submit buttom I got a suggestion from Mike:
Copying here:
I get to the same error using Prof. Maes example for a bi-phenotipic model, so I suppose the problem lies in the use of mxRefModels?
The error:
I can't edit the post above, the correct code from Prof Maes, which gives the error is:
To fix it, we can try harder with the optimization manually like this:
However, I note that it would be easier if mxRefModels() took an method= argument to use mxTryHard() or mxTryHardOrdinal() to fit the saturated and independence models in case of poor starting values.
Thanks.
This particular issue was resolved, sorry for not updating here. The culprit was the line:
options(mxByrow = TRUE)
This affects all code run after it. It should never be used this way. Such a noob error.