######################################################################## ######################################################################## ######################################################################## ###################################### ### ### Optimization problems in OpenMx ### ###################################### ### Libraries used library( OpenMx ) library( polycor ) # For independent calculation of tetrachoric correlation mxOption(NULL,"Default optimizer","SLSQP") set.seed(161103) ### Generate data datS <- mxFactor( data.frame( x1=c(rep(0,1600000),rep(1,800),rep(0,810),rep(1,195),rep(0,810),rep(1,3),rep(0,9),rep(1,2),rep(0,780),rep(1,6),rep(0,6),rep(1,2),rep(0,200),rep(1,2),rep(1,1)) , x2=c(rep(0,1600000),rep(0,800),rep(1,810),rep(1,195),rep(0,810),rep(0,3),rep(1,9),rep(1,2),rep(0,780),rep(0,6),rep(1,6),rep(1,2),rep(0,200),rep(0,2),rep(1,1)) , x3=c(rep(0,1600000),rep(0,800),rep(0,810),rep(0,195),rep(1,810),rep(1,3),rep(1,9),rep(1,2),rep(0,780),rep(0,6),rep(0,6),rep(0,2),rep(1,200),rep(1,2),rep(1,1)) , x4=c(rep(0,1600000),rep(0,800),rep(0,810),rep(0,195),rep(0,810),rep(0,3),rep(0,9),rep(0,2),rep(1,780),rep(1,6),rep(1,6),rep(1,2),rep(1,200),rep(1,2),rep(1,1)) ) , levels=c(0,1) ) varNames <- names(datS) ### Observed tetrachoric correlations obsMat <- matrix(1,4,4) for(i in 1:3){ for(j in (i+1):4){ print( table(datS[,c(i,j)] ) ) # Print each table obsMat[i,j] <- obsMat[j,i] <- polychor(table(datS[,c(i,j)]) , ML=T ) } } round( obsMat , 3 ) ### Observed prevalences colMeans( 1*(datS==1) ) ### Problem inverting observed correlation matrix? round( solve( obsMat ) , 4 ) ### Very simple model simpMod <- mxModel( 'SimpMod', ### Expected covariance matrix (variance==1 due to binary variables) mxMatrix(type='Stand',nrow=4,ncol=4,free=T,values=obsMat[lower.tri(obsMat)],name='eC',lbound=0,ubound=0.99), ### Expected Means mxMatrix(type='Full',nrow=1,ncol=4,free=F,values=0,name='eM'), ### Expected thresholds mxMatrix(type='Full',nrow=1,ncol=4,free=T,values=qnorm(colMeans( 1*(datS==1) ),lower.tail=F),name='eT'), ### Data and likelihood mxModel( 'Datmod', mxData(datS , type='raw'), mxExpectationNormal(means='SimpMod.eM',covariance='SimpMod.eC',thresholds='SimpMod.eT',threshnames=varNames,dimnames=varNames), mxFitFunctionML() ), ### Add the likelihoods (only one) mxFitFunctionMultigroup( c('Datmod.fitfunction') ) ) ### Fit model simpModFit <- mxTryHardOrdinal( simpMod, scale=0.1 ) summary(simpModFit) fitfunc <- function(model,state){ b <- model$Betas$values out <- 0 #0000: out <- out - 1600000*2*( dbinom(0,1,pnorm(b[1,1],lower.tail=F),log=T) + dbinom(0,1,pnorm(b[2,2],lower.tail=F),log=T) + dbinom(0,1,pnorm(b[3,3],lower.tail=F),log=T) + dbinom(0,1,pnorm(b[4,4],lower.tail=F),log=T) ) #0001: out <- out - 780*2*( dbinom(0,1,pnorm(b[1,1],lower.tail=F),log=T) + dbinom(0,1,pnorm(b[2,2],lower.tail=F),log=T) + dbinom(0,1,pnorm(b[3,3],lower.tail=F),log=T) + dbinom(1,1,pnorm(b[4,4],lower.tail=F),log=T) ) #0010: out <- out - 810*2*( dbinom(0,1,pnorm(b[1,1],lower.tail=F),log=T) + dbinom(0,1,pnorm(b[2,2],lower.tail=F),log=T) + dbinom(1,1,pnorm(b[3,3],lower.tail=F),log=T) + dbinom(0,1,pnorm(b[4,4]+b[4,3],lower.tail=F),log=T) ) #0011: out <- out - 200*2*( dbinom(0,1,pnorm(b[1,1],lower.tail=F),log=T) + dbinom(0,1,pnorm(b[2,2],lower.tail=F),log=T) + dbinom(1,1,pnorm(b[3,3],lower.tail=F),log=T) + dbinom(1,1,pnorm(b[4,4]+b[4,3],lower.tail=F),log=T) ) #0100: out <- out - 810*2*( dbinom(0,1,pnorm(b[1,1],lower.tail=F),log=T) + dbinom(1,1,pnorm(b[2,2],lower.tail=F),log=T) + dbinom(0,1,pnorm(b[3,3]+b[3,2],lower.tail=F),log=T) + dbinom(0,1,pnorm(b[4,4]+b[4,2],lower.tail=F),log=T) ) #0101: out <- out - 6*2*( dbinom(0,1,pnorm(b[1,1],lower.tail=F),log=T) + dbinom(1,1,pnorm(b[2,2],lower.tail=F),log=T) + dbinom(0,1,pnorm(b[3,3]+b[3,2],lower.tail=F),log=T) + dbinom(1,1,pnorm(b[4,4]+b[4,2],lower.tail=F),log=T) ) #0110: out <- out - 9*2*( dbinom(0,1,pnorm(b[1,1],lower.tail=F),log=T) + dbinom(1,1,pnorm(b[2,2],lower.tail=F),log=T) + dbinom(1,1,pnorm(b[3,3]+b[3,2],lower.tail=F),log=T) + dbinom(0,1,pnorm(b[4,4]+b[4,3]+b[4,2],lower.tail=F),log=T) ) #1000: out <- out - 800*2*( dbinom(1,1,pnorm(b[1,1],lower.tail=F),log=T) + dbinom(0,1,pnorm(b[2,2]+b[2,1],lower.tail=F),log=T) + dbinom(0,1,pnorm(b[3,3]+b[3,1],lower.tail=F),log=T) + dbinom(0,1,pnorm(b[4,4]+b[4,1],lower.tail=F),log=T) ) #1001: out <- out - 6*2*( dbinom(1,1,pnorm(b[1,1],lower.tail=F),log=T) + dbinom(0,1,pnorm(b[2,2]+b[2,1],lower.tail=F),log=T) + dbinom(0,1,pnorm(b[3,3]+b[3,1],lower.tail=F),log=T) + dbinom(1,1,pnorm(b[4,4]+b[4,1],lower.tail=F),log=T) ) #1010: out <- out - 3*2*( dbinom(1,1,pnorm(b[1,1],lower.tail=F),log=T) + dbinom(0,1,pnorm(b[2,2]+b[2,1],lower.tail=F),log=T) + dbinom(1,1,pnorm(b[3,3]+b[3,1],lower.tail=F),log=T) + dbinom(0,1,pnorm(b[4,4]+b[4,3]+b[4,1],lower.tail=F),log=T) ) #1011: out <- out - 2*2*( dbinom(1,1,pnorm(b[1,1],lower.tail=F),log=T) + dbinom(0,1,pnorm(b[2,2]+b[2,1],lower.tail=F),log=T) + dbinom(1,1,pnorm(b[3,3]+b[3,1],lower.tail=F),log=T) + dbinom(1,1,pnorm(b[4,4]+b[4,3]+b[4,1],lower.tail=F),log=T) ) #1100: out <- out - 195*2*( dbinom(1,1,pnorm(b[1,1],lower.tail=F),log=T) + dbinom(1,1,pnorm(b[2,2]+b[2,1],lower.tail=F),log=T) + dbinom(0,1,pnorm(b[3,3]+b[3,2]+b[3,1],lower.tail=F),log=T) + dbinom(0,1,pnorm(b[4,4]+b[4,2]+b[4,1],lower.tail=F),log=T) ) #1101: out <- out - 2*2*( dbinom(1,1,pnorm(b[1,1],lower.tail=F),log=T) + dbinom(1,1,pnorm(b[2,2]+b[2,1],lower.tail=F),log=T) + dbinom(0,1,pnorm(b[3,3]+b[3,2]+b[3,1],lower.tail=F),log=T) + dbinom(1,1,pnorm(b[4,4]+b[4,2]+b[4,1],lower.tail=F),log=T) ) #1110: out <- out - 2*2*( dbinom(1,1,pnorm(b[1,1],lower.tail=F),log=T) + dbinom(1,1,pnorm(b[2,2]+b[2,1],lower.tail=F),log=T) + dbinom(1,1,pnorm(b[3,3]+b[3,2]+b[3,1],lower.tail=F),log=T) + dbinom(0,1,pnorm(b[4,4]+b[4,3]+b[4,2]+b[4,1],lower.tail=F),log=T) ) #1111: out <- out - 2*( dbinom(1,1,pnorm(b[1,1],lower.tail=F),log=T) + dbinom(1,1,pnorm(b[2,2]+b[2,1],lower.tail=F),log=T) + dbinom(1,1,pnorm(b[3,3]+b[3,2]+b[3,1],lower.tail=F),log=T) + dbinom(1,1,pnorm(b[4,4]+b[4,3]+b[4,2]+b[4,1],lower.tail=F),log=T) ) return(out) } simpMod2 <- mxModel( "SimpMod2", mxMatrix(type="Lower",nrow=4,free=T,byrow=T,name="Betas", values=c( 3, 0,3, 0,0,3, 0,0,0,3), labels=c( "beta11", "beta21","beta22", "beta31","beta32","beta33", "beta41","beta42","beta43","beta44" )), mxFitFunctionR(fitfunc) ) fitfunc(simpMod2) ### Fit model set.seed(161103) simpModFit2 <- mxTryHard(simpMod2) summary(simpModFit2)