Problems with variance estimation with missing continuous data: WLS seem to replace NA with 0
hansfredriksunde
Joined: 12/10/2021
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