Hi everyone

I'm currently building an extended children-of-twins-and-siblings model that will have dichotomous data in the child generation and continuous data in the parent generation. The model includes both twins and their partners (co-parents), so that assortative mating can be taken into account. The number of free parameters (>10) make maximum likelihood unfeasible when using dichotomous data (i.e., the processing time is too long). I therefore use the WLS fit function instead.

I have quite a large sample size (> 10 000 sets of extended families), but the number of monozygotic twins with children is not that high (~250 female pairs and ~250 male pairs). I can double that number if I include twins where only one of them has children. For these pairs, the second twin's partner and children will be missing. The model seem to handle missing children (dichotomous data) without any problem, but runs into trouble when the second twin's partner is missing (continuous data). When examining the results (based on simulated data), it becomes apparent that the model underestimates the variance for the continuous variables with missing data, which leads to knock-on consequences as this error promulgates through the structure of equations.

I have now made a minimal example which reproduces this error. In this example, I have only simulated three variables (twin 1, twin 2, and child of twin 1). Here, it is the co-twin rather than the co-twin's partner that is missing in half the cases, but the error is the same. When doing this, I noticed that I get the exact same results if I set the missing values to 0 instead of NA. This leads me to suspect that this is what is happening under the hood with WLS. Is this a bug, or is it a mathematical consequence of how WLS works? What can I do about it?

Thanks in advance for all help!

Sincerely

Hans Fredrik

PS: (I have tried to include these twins as a separate group, and define separate expectation matrices for them, but then I simply return to the issue I had initially: That the subgroups are too small for WLS).

Here is the reproducible example:

library(OpenMx) mxOption(key="Default optimizer", value="NPSOL") # ---------------------------------------------------------------------------------------------------------------------- # PREPARE DATA (SIMULATION) # Set parameters n <- 10000 # n of families per zygosity VA1 <- 0.4 # VE1 = 1-VA1 cut <- 1 # threshold # Simulates genetic components A1_mz <- MASS::mvrnorm(n,c(0,0,0), matrix(c( 1, 1,.5, 1, 1,.5, .5,.5, 1 ),3,3,T)) A1_dz <- MASS::mvrnorm(n,c(0,0,0), matrix(c( 1,.5,.5, .5, 1,.25, .5,.25, 1),3,3,T)) # Simulates environmental components E1_mz <- MASS::mvrnorm(n,c(0,0,0), matrix(c( 1, 0, 0, 0, 1, 0, 0, 0, 1 ),3,3,T)) E1_dz <- MASS::mvrnorm(n,c(0,0,0), matrix(c( 1, 0, 0, 0, 1, 0, 0, 0, 1 ),3,3,T)) # Weighs and combines genetic and environmental components. dt_mz <- as.data.frame( sqrt(VA1)*A1_mz + sqrt(1-VA1)*E1_mz ) dt_dz <- as.data.frame( sqrt(VA1)*A1_dz + sqrt(1-VA1)*E1_dz ) # Sets variable names colnames(dt_mz) <- c("proband","cotwin","child") colnames(dt_dz) <- c("proband","cotwin","child") # Randomly removes half of cotwin observations dt_mz[sample(nrow(dt_mz),nrow(dt_mz)/2),"cotwin"] <- NA dt_dz[sample(nrow(dt_dz),nrow(dt_dz)/2),"cotwin"] <- NA # This dichotomizes the child phenotype and turns it into a factor dt_mz[,"child"] <- mxFactor(ifelse(dt_mz[,"child"] > cut, 1, 0), levels = c(0,1)) dt_dz[,"child"] <- mxFactor(ifelse(dt_dz[,"child"] > cut, 1, 0), levels = c(0,1)) # ---------------------------------------------------------------------------------------------------------------------- # PREPARE MODEL # Select Variables for Analysis selVars <- c("proband","cotwin","child") # Create Data Objects dataMZ <- mxData( observed=dt_mz, type="raw") dataDZ <- mxData( observed=dt_dz, type="raw" ) ##################################################### #### Creates Free Parameters #### ##################################################### varA1 <- mxMatrix(name="VA1", type="Symm", nrow=1, ncol=1, free=TRUE, values=.5, label="est_VA1" ) varE1 <- mxMatrix(name="VE1", type="Symm", nrow=1, ncol=1, free=TRUE, values=.5, label="est_VE1" ) varE2 <- mxMatrix(name="VE2", type="Symm", nrow=1, ncol=1, free=TRUE, values=.5, label="est_VE2" ) ####################################################### #### Creates Algebra for Expected (Co-)Variances #### ####################################################### # Parent variance var1 <- mxAlgebra(expression = VA1+VE1, name="var1" ) #Child variance var2 <- mxAlgebra(expression = VA1+VE2, name="var2" ) cons <- mxConstraint(expression = var2 == 1, name="var2_cons") # Twin covariance covTW_MZ <- mxAlgebra(expression = VA1, name="covTW_MZ") covTW_DZ <- mxAlgebra(expression = 0.5*VA1, name="covTW_DZ") # Parent-offspring covariance covPO <- mxAlgebra(expression = 0.5*VA1, name="covPO") # Avuncular covariance covAV_MZ <- mxAlgebra(expression = 0.5*VA1, name="covAV_MZ") covAV_DZ <- mxAlgebra(expression =0.25*VA1, name="covAV_DZ") ####################################################### #### Creates Expectation Objects #### ####################################################### expCovMZ <- mxAlgebra( expression= rbind( cbind( var1, covTW_MZ, covPO), cbind( covTW_MZ, var1, covAV_MZ), cbind( covPO, covAV_MZ, var2)), name="expCov_MZ" ) expCovDZ <- mxAlgebra( expression= rbind( cbind( var1, covTW_DZ, covPO), cbind( covTW_DZ, var1, covAV_DZ), cbind( covPO, covAV_DZ, var2)), name="expCov_DZ" ) meanG <- mxMatrix( type="Full", nrow=1, ncol=3, free=c(T,T,F), values=0, labels=c("mean1","mean1","mean2"), name="meanG" ) threG <- mxMatrix( type="Full", nrow=1, ncol=1, free=T, values=0, labels="threshold", name="threG" ) expMZ <- mxExpectationNormal( covariance = "expCov_MZ", means = "meanG", thresholds = "threG", dimnames = selVars, threshnames = "child") expDZ <- mxExpectationNormal( covariance = "expCov_DZ", means = "meanG", thresholds = "threG", dimnames = selVars, threshnames = "child") ####################################################### #### Creates Model #### ####################################################### fitFun <- mxFitFunctionWLS("WLS") # (Using ML on the toy example yields the correct results) # fitFun <- mxFitFunctionML() # Create Model Objects modelMZ <- mxModel(varA1, varE1, varE2, var1, var2, cons, covTW_MZ, covPO, covAV_MZ, meanG, threG, expCovMZ, dataMZ, expMZ, fitFun, name="modelMZ" ) modelDZ <- mxModel(varA1, varE1, varE2, var1, var2, cons, covTW_DZ, covPO, covAV_DZ, meanG, threG, expCovDZ, dataDZ, expDZ, fitFun, name="modelDZ" ) multi <- mxFitFunctionMultigroup( c("modelMZ", "modelDZ") ) model1 <- mxModel(varA1, varE1, varE2, modelMZ,modelDZ,multi, name = "Model1") ####################################################### #### Run Model #### ####################################################### fit1 <- mxRun(model1) # Check estimates (should be VA1 = .40, VE1/2 = .60) round(fit1$output$estimate,2) # check calculated covariance # Notice how .[2,2] is ~.5, even though it should be ~1 fit1$output$data[[1]]$cov # Check estimated variance (Should be var1 = 1) fit1$modelMZ$var1

I'm not able to run that syntax as-is. I get an error at around line 54: "Error: The name 'co-twin' is illegal because it contains the '-' character in mxData(observed = dt_mz, type = "raw")". What is your

`mxVersion()`

output?Anyhow, if I globally change "co-twin" to "coTwin", I can run the script. I reproduce the underestimation of the variance for co-twins. I don't know enough about the internals of WLS in OpenMx to diagnose what is going on, though.

Ah, thanks for pointing that out. I got the same error. I must have changed the variable names into something more sensible last-minute, and forgot to double check that it actually ran. I have edited the supplied code now.

Good catch! That is, indeed, a bug. WLS is treating NAs as zeros for continuous variables. We'll be fixing this shortly and making a new release. Thanks for the great detective work and thorough example!

An update on this: I found the bug. It was not treating missing values as zeros, but it was dividing by the total number of rows instead of the total number of non-missing rows. In R-style code, it's

I'm testing the fix now. The OpenMx dev team will make a release shortly thereafter.

Fantastic news!

Thank you for finding the bug and fixing it quickly. It really helps my PhD project staying on track.