tests/testthat/test-GREML_Error_Detection.R

In OpenMx: Extended Structural Equation Modelling

#
#   Copyright 2007-2020 by the individuals mentioned in the source code history
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
# 
#        http://www.apache.org/licenses/LICENSE-2.0
# 
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.


library(OpenMx)
library(testthat)
context("GREML Error Detection")

options(mxCondenseMatrixSlots=TRUE)
mxOption(NULL,"Analytic Gradients","Yes")
require(mvtnorm)

omxCheckError(mxExpectationGREML(V=1),
              "argument 'V' is not of type 'character' (the name of the expected covariance matrix)")

set.seed(1234)
dat <- cbind(rnorm(100),rep(1,100))
colnames(dat) <- c("y","x")
testmod <- mxModel(
  "GREMLtest",
  #mxData(observed=dat, type="raw", sort=F),
  mxMatrix(type = "Full", nrow = 1, ncol=1, free=T, values = 2, labels = "ve", lbound = 0.0001, name = "Ve"),
  mxMatrix("Iden",nrow=100,name="I",condenseSlots=T),
  mxAlgebra(I %x% Ve,name="V"),
  mxExpectationGREML(V="V",Xvars=list("x"),yvars="y",addOnes=F),
  mxFitFunctionGREML()
)
omxCheckError(mxRun(testmod),
              "the GREML expectation function does not have a dataset associated with it in model 'GREMLtest'")

testmod <- mxModel(
  testmod,
  mxData(observed=dat, type="raw", sort=F),
  mxMatrix("Full",1,1,F,labels="data.y",name="Z")
)
omxCheckError(mxRun(testmod),
              "definition variables are incompatible (and unnecessary) with GREML expectation")

testmod <- mxModel(
  "GREMLtest",
  mxData(observed=dat, type="raw", sort=F),
  mxMatrix(type = "Full", nrow = 1, ncol=1, free=T, values = 2, labels = "ve", lbound = 0.0001, name = "Ve"),
  mxMatrix("Full",nrow=100,ncol=99,free=F,values=diag(100)[,-100],name="V",condenseSlots=T),
  mxExpectationGREML(V="V",dataset.is.yX=T),
  mxFitFunctionGREML()
)
omxCheckError(mxRun(testmod),
              "'V' matrix is not square")

testmod$V <- mxMatrix("Full",nrow=99,ncol=99,free=F,values=diag(100)[-100,-100],name="V",condenseSlots=T)
omxCheckError(mxRun(testmod),
              "y and V matrices do not have equal numbers of rows")

testmod <- mxModel(
  "GREMLtest",
  mxData(observed=dat, type="raw", sort=F),
  mxMatrix(type = "Full", nrow = 1, ncol=1, free=T, values = 1, labels = "ve", lbound = 0.0001, name = "Ve"),
  mxMatrix("Iden",nrow=100,name="I",condenseSlots=T),
  mxAlgebra(I %x% Ve,name="V"),
  mxExpectationGREML(V="V",dataset.is.yX=TRUE,casesToDropFromV=101L),
  mxFitFunctionGREML()
)
omxCheckWarning(
  mxRun(testmod, suppressWarnings=TRUE),
  "casesToDrop vector in GREML expectation contains indices greater than the number of datapoints")

testmod <- mxModel(
  "GREMLtest",
  mxData(observed=dat, type="raw", sort=F),
  mxMatrix(type = "Full", nrow = 1, ncol=1, free=T, values = 2, labels = "ve", lbound = 0.0001, name = "Ve"),
  mxMatrix("Iden",nrow=100,name="I",condenseSlots=T),
  mxMatrix("Iden",nrow=99,name="J",condenseSlots=T),
  mxAlgebra(I %x% Ve,name="V"),
  mxExpectationGREML(V="V", dataset.is.yX = T),
  mxFitFunctionGREML(dV = c(ve="J"))
)
omxCheckError(mxRun(testmod),
              "all derivatives of V must have the same dimensions as V")

testmod <- mxModel(
  "GREMLtest",
  mxData(observed=dat, type="raw", sort=F),
  mxMatrix(type = "Full", nrow = 1, ncol=1, free=T, values = 2, labels = "ve", lbound = 0.0001, name = "Ve"),
  mxMatrix("Iden",nrow=100,name="I",condenseSlots=T),
  mxAlgebra(I %x% Ve,name="V"),
  mxExpectationGREML(V="V", dataset.is.yX=T),
  mxFitFunctionGREML(dV = c(ve="J"))
)
omxCheckError(mxRun(testmod),
              "The reference 'J' does not exist.  It is used by named reference 'GREMLtest.fitfunction' .")


testmod <- mxModel(
  "GREMLtest",
  mxData(observed = matrix(dat[,1],1,100,dimnames=list(NULL,paste("y",1:100,sep=""))), type="raw"),
  mxMatrix(type = "Full", nrow = 1, ncol=1, free=T, values = 2, labels = "ve", lbound = 0.0001, name = "Ve"),
  mxMatrix("Iden",nrow=100,name="I",condenseSlots=T),
  mxMatrix("Zero",1,100,name="Zm"),
  mxAlgebra(I %x% Ve,name="V"),
  mxExpectationNormal(covariance="V",means="Zm", dimnames=paste("y",1:100,sep="")),
  mxFitFunctionGREML()
)
omxCheckError(mxRun(testmod),
              "GREML fitfunction is currently only compatible with GREML expectation")


testmod <- mxModel(
	"GREMLtest",
	mxData(observed=dat, type="raw", sort=F),
	mxMatrix(type = "Full", nrow = 1, ncol=1, free=T, values = 2, labels = "ve", lbound = 0.0001, name = "Ve"),
	mxMatrix("Iden",nrow=100,name="I",condenseSlots=T),
	mxAlgebra(I %x% Ve,name="V"),
	mxExpectationGREML(V="V",Xvars=list("x"),yvars="y",addOnes=F),
	mxFitFunctionGREML()
)
omxCheckError(mxRefModels(testmod),
							"Reference models for GREML expectation are not implemented")


omxCheckError(mxModel(
	"GREMLtest",
	mxData(observed=dat, type="raw", sort=F),
	mxMatrix(type = "Full", nrow = 1, ncol=1, free=T, values = 2, labels = "ve", lbound = 0.0001, name = "Ve"),
	mxMatrix("Iden",nrow=100,name="I",condenseSlots=T),
	mxAlgebra(I %x% Ve,name="V"),
	mxExpectationGREML(V="V",Xvars=list("x"),yvars="y",addOnes=F),
	mxFitFunctionGREML(autoDerivType="muneric")
),"'muneric' should be one of 'semiAnalyt' and 'numeric'")


testmod <- mxModel(
	"GREMLtest",
	mxData(observed=dat, type="raw", sort=F),
	mxMatrix(type = "Unit", nrow = 100, ncol=100, name = "V", condenseSlots = T),
	mxExpectationGREML(V="V",Xvars=list("x"),yvars="y",addOnes=F),
	mxFitFunctionGREML()
)
omxCheckError(mxRun(testmod),
							"Expected covariance matrix is non-positive-definite at initial values")



testmod <- mxModel(
	"GREMLtest",
	mxData(observed=dat, type="raw", sort=F),
	mxMatrix(type = "Unit", nrow = 100, ncol=100, name = "V", condenseSlots = T),
	mxExpectationGREML(V="V",Xvars=list("x"),yvars="y",addOnes=F),
	mxFitFunctionML()
)
omxCheckError(mxRun(testmod),
							"Expected covariance matrix is non-positive-definite at initial values")



z <- matrix(-1,100,2)
colnames(z) <- c("z1","z2")
dat2 <- cbind(dat,z)
testmod <- mxModel(
	"GREMLtest",
	mxData(observed=dat2, type="raw", sort=F),
	mxMatrix(type = "Full", nrow = 1, ncol=1, free=T, values = 2, labels = "ve", lbound = 0.0001, name = "Ve"),
	mxMatrix("Iden",nrow=100,name="I",condenseSlots=T),
	mxAlgebra(I %x% Ve,name="V"),
	mxExpectationGREML(V="V",Xvars=list(c("x","z1","z2")),yvars="y",addOnes=F),
	mxFitFunctionGREML()
)
omxCheckError(mxRun(testmod),
	"Cholesky factorization failed at initial values; possibly, the matrix of covariates is rank-deficient")


testmod <- mxModel(
	"GREMLtest",
	mxData(observed=dat2, type="raw", sort=F),
	mxMatrix(type = "Full", nrow = 1, ncol=1, free=T, values = 2, labels = "ve", lbound = 0.0001, name = "Ve"),
	mxMatrix("Iden",nrow=100,name="I",condenseSlots=T),
	mxAlgebra(I %x% Ve,name="V"),
	mxExpectationGREML(V="V",Xvars=list(c("x","z1","z2")),yvars="y",addOnes=F),
	mxFitFunctionML()
)
omxCheckError(mxRun(testmod),
							"Cholesky factorization failed at initial values; possibly, the matrix of covariates is rank-deficient")


set.seed(476)
A1 <- matrix(0,100,100)  
A1[lower.tri(A1)] <- runif(4950, -0.025, 0.025)
A1 <- A1 + t(A1)
diag(A1) <- runif(100,0.95,1.05)
A2 <- matrix(0,100,100)  
A2[lower.tri(A2)] <- runif(4950, -0.025, 0.025)
A2 <- A2 + t(A2)
diag(A2) <- runif(100,0.95,1.05)
y <- t(rmvnorm(1,sigma=A1*0.25)+rmvnorm(1,sigma=A2*0.25))  
y <- y + rnorm(100,sd=sqrt(0.5))
x <- rnorm(100) 
dat3 <- cbind(y,x)
rm(x,y)
colnames(dat3) <- c("y","x")
testmod <- mxModel(
	"GREMLtest",
	mxMatrix(type = "Full", nrow = 1, ncol=1, free=T, values =0.5, labels = "ve", lbound = 0.0001, 
					 name = "Ve"),
	mxMatrix(type = "Full", nrow = 1, ncol=1, free=T, values = 0.25, labels = "va1", name = "Va1"),
	mxMatrix(type = "Full", nrow = 1, ncol=1, free=T, values = 0.25, labels = "va2", name = "Va2"),
	mxData(observed = dat3, type="raw", sort=FALSE),
	mxExpectationGREML(V="V",yvars="y", Xvars="x", addOnes=T),
	mxMatrix("Iden",nrow=100,name="I"),
	mxMatrix("Symm",nrow=100,free=F,values=A1,name="A1"),
	mxMatrix("Symm",nrow=100,free=F,values=A2,name="A2"),
	mxAlgebra((A1%x%Va1) + (A2%x%Va2) + (I%x%Ve), name="V"),
	mxMatrix(type="Full",nrow=1,ncol=1,free=F,values=0.64,name="aug"),
	mxMatrix(type="Zero",nrow=1,ncol=1,name="Zilch"),
	mxFitFunctionGREML(dV=c(ve="I",va1="A1",va2="A2"),aug="aug",augHess="Zilch")
)
omxCheckError(
	mxRun(testmod),
	"if argument 'augHess' has nonzero length, then argument 'augGrad' must as well")


testmod$fitfunction <- mxFitFunctionGREML(dV=c(ve="I",va1="A1",va2="A2"),aug="aug")
omxCheckError(
	mxRun(testmod),
	"if arguments 'dV' and 'aug' have nonzero length, then 'augGrad' must as well")


testmod$fitfunction <- mxFitFunctionGREML(dV=c(ve="I",va1="A1",va2="A2",va3="I"))
omxCheckError(
	mxRun(testmod),
	"length of argument 'dV' is greater than the number of explicit free parameters")



testmod <- mxModel(
	"GREMLtest",
	mxMatrix(type = "Full", nrow = 1, ncol=1, free=T, values =0.5, labels = "ve", lbound = 0.0001, 
					 name = "Ve"),
	mxMatrix(type = "Full", nrow = 1, ncol=1, free=T, values = 0.25, labels = "va1", name = "Va1"),
	mxMatrix(type = "Full", nrow = 1, ncol=1, free=T, values = 0.25, labels = "va2", name = "Va2"),
	mxData(observed = dat3, type="raw", sort=FALSE),
	mxExpectationGREML(V="V",yvars="y", Xvars="x", addOnes=T),
	mxMatrix("Iden",nrow=100,name="I"),
	mxMatrix("Symm",nrow=100,free=F,values=A1,name="A1"),
	mxMatrix("Symm",nrow=100,free=F,values=A2,name="A2"),
	mxAlgebra((A1%x%Va1) + (A2%x%Va2) + (I%x%Ve), name="V"),
	mxComputeSequence(steps=list(
		mxComputeNewtonRaphson(fitfunction="fitfunction"),
		mxComputeOnce('fitfunction', c('fit','gradient','hessian','ihessian')),
		mxComputeStandardError(),
		mxComputeReportDeriv(),
		mxComputeReportExpectation()
	)),
	mxFitFunctionGREML(dV=c(ve="I",va1="A1",va2="A2",va3="I"))
)
omxCheckError(
	mxRun(testmod),
	"length of argument 'dV' is greater than the number of explicit free parameters")



testmod <- mxModel(
	"GREMLtest",
	mxMatrix(type = "Full", nrow = 1, ncol=1, free=T, values =0.5, labels = "ve", lbound = 0.0001, 
					 name = "Ve"),
	mxMatrix(type = "Full", nrow = 1, ncol=1, free=T, values = 0.25, labels = "va1", name = "Va1"),
	mxMatrix(type = "Full", nrow = 1, ncol=1, free=T, values = 0.25, labels = "va2", name = "Va2"),
	mxData(observed = dat3, type="raw", sort=FALSE),
	mxExpectationGREML(V="V",yvars="y", Xvars="x", addOnes=T),
	mxMatrix("Iden",nrow=100,name="I"),
	mxMatrix("Symm",nrow=100,free=F,values=A1,name="A1"),
	mxMatrix("Symm",nrow=100,free=F,values=A2,name="A2"),
	mxAlgebra((A1%x%Va1) + (A2%x%Va2) + (I%x%Ve), name="V"),
	mxComputeSequence(steps=list(
		mxComputeNewtonRaphson(fitfunction="fitfunction"),
		mxComputeOnce('fitfunction', c('fit','gradient','hessian')),
		mxComputeStandardError(),
		mxComputeReportDeriv(),
		mxComputeReportExpectation()
	)),
	mxMatrix(type="Full",nrow=1,ncol=1,free=F,values=2.1,name="aug"),
	mxMatrix(type="Zero",nrow=2,ncol=1,free=F,name="ag"),
	mxMatrix(type="Zero",nrow=3,ncol=3,free=F,name="ah"),
	mxFitFunctionGREML(dV=c(ve="I",va1="A1",va2="A2"),aug="aug",augGrad="ag",augHess="ah")
)
omxCheckError(
	mxRun(testmod),
	"matrix referenced by 'augGrad' must have as many elements as there are explicit free parameters")



testmod$ag <- mxMatrix(type="Zero",nrow=3,ncol=1,free=F,name="ag")
testmod$ah <- mxMatrix(type="Zero",nrow=2,ncol=3,free=F,name="ah")
omxCheckError(
	mxRun(testmod),
	"matrix referenced by 'augHess' must be square (instead of 2x3)")



testmod$ah <- mxMatrix(type="Zero",nrow=2,ncol=2,free=F,name="ah")
omxCheckError(
	mxRun(testmod),
	"Augmentation derivatives non-conformable (gradient is size 3 and Hessian is 2x2)")

mxOption(reset=TRUE)