# Problems with variance estimation with missing continuous data: WLS seem to replace NA with 0

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

## mxVersion()

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.

Log in or register to post comments

In reply to mxVersion() by AdminRobK

## Thanks, edited the code

Log in or register to post comments

## That's a bug

Log in or register to post comments

In reply to That's a bug by AdminHunter

## Bug Found

# Sort of treats missing values as zeros

sum(x, na.rm=TRUE)/length(x)

# versus

# Same sum, but remove missing count from demoninator

sum(x, na.rm=TRUE)/(length(x) - sum(is.na(x)))

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

Log in or register to post comments

## THANKS!

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

Log in or register to post comments