Nothing
#
# Copyright 2007-2019 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.
#------------------------------------------------------------------------------
# Frontpage model in LISREL form
# with factor scores
library(OpenMx)
library(testthat)
context("FactorScores")
suppressWarnings(RNGversion("3.5"))
skip_on_cran()
data(demoOneFactor)
manifests <- names(demoOneFactor)
latents <- c("G")
factorModelL <- mxModel("OneFactor",
type="LISREL",
manifestVars=list(exo=manifests),
latentVars=list(exo=latents),
mxPath(from=latents, to=manifests),
mxPath(from=manifests, arrows=2),
mxPath(from=latents, arrows=2, free=FALSE, values=1.0),
mxPath(from='one', to=manifests),
mxData(observed=demoOneFactor, type="raw"))
summary(factorRunL <- mxRun(factorModelL))
r1 <- mxFactorScores(factorRunL, 'ML')
r2 <- mxFactorScores(factorRunL, 'Regression')
r3 <- mxFactorScores(factorRunL, 'WeightedML')
omxCheckCloseEnough(cor(r1[,,1], r2[,,1]), 1)
#pdf('plotFactorScores.pdf')
#plot(r2[,,1], r1[,,1], main='Factor Scoring Methods', xlab='Regression Score', ylab='Likelihood Score')
#legend('bottomright', legend=c('ML', 'Weighted ML'), pch=1, col=c('black', 'blue'), lty=1)
#points(r2[,,1], r3[,,1], col='blue')
#abline(v=0, h=0)
lmS <- lm(r3[,,1] ~ r2[,,1])
lmR <- lm(r1[,,1] ~ r2[,,1])
#abline(lmS, col='blue')
#abline(lmR, col='black')
#dev.off()
summary(lmS)
summary(lmR)
omxCheckCloseEnough(coef(lmS), c(0, .5), 0.01)
omxCheckCloseEnough(coef(lmR), c(0, 1), 0.03)
if (mxOption(key="Default optimizer") != 'NPSOL') {
# NPSOL can add a warning about code green
#Test warning about no SEs:
#factorRun <- mxOption(factorRun,"Standard Errors","No")
mxOption(NULL,"Standard Errors","No")
omxCheckWarning(
mxFactorScores(factorRunL,"ML"),
"factor-score standard errors not available from MxModel 'OneFactor' because calculating SEs is turned off for that model (possibly due to one or more MxConstraints)")
omxCheckWarning(
mxFactorScores(factorRunL,"WeightedML"),
"factor-score standard errors not available from MxModel 'OneFactor' because calculating SEs is turned off for that model (possibly due to one or more MxConstraints)")
mxOption(NULL,"Standard Errors","Yes")
}
#TODO compare standard errors
#cbind(r2[,,2], r1[,,2], r3[,,2])
#------------------------------------------------------------------------------
# Ordinal example similar to thresholdModel1Factor3Variate.R
# Step 1: load libraries
require(OpenMx)
#
# Step 2: set up simulation parameters
# Note: nVariables>=3, nThresholds>=1, nSubjects>=nVariables*nThresholds (maybe more)
# and model should be identified
#
nVariables<-3
nFactors<-1
nThresholds<-3
nSubjects<-500
isIdentified<-function(nVariables,nFactors) as.logical(1+sign((nVariables*(nVariables-1)/2) - nVariables*nFactors + nFactors*(nFactors-1)/2))
# if this function returns FALSE then model is not identified, otherwise it is.
isIdentified(nVariables,nFactors)
# Step 3: simulate multivariate normal data
set.seed(1234)
quants <- quantile(rnorm(1001), probs = c((1:nThresholds)/(nThresholds+1)))
invL <- matrix(0, nThresholds, nThresholds)
invL[lower.tri(invL, diag=TRUE)] <- 1
thresholdStart <- solve(invL) %*% quants
fruitynames<-paste("banana",1:nVariables,sep="")
lis <- mxModel("thresholdModel",
mxMatrix("Full", nVariables, nFactors, values=0.7, free=TRUE, lbound=-.99, ubound=.99, name="L", dimnames=list(fruitynames, 'F')),
mxMatrix("Unit", nVariables, 1, name="vectorofOnes"),
mxMatrix("Zero", nVariables, 1, name="M", dimnames=list(fruitynames, 'F')),
mxAlgebra(vec2diag(vectorofOnes - (diag2vec(L %*% t(L)))) , name="E"),
mxMatrix('Iden', 1, 1, name='P', dimnames=list('F', 'F')),
mxMatrix('Zero', 1, 1, name='K', dimnames=list('F', c())),
mxMatrix("Full",
name="thresholdDeviations", nrow=nThresholds, ncol=nVariables,
values=thresholdStart,
free = TRUE,
lbound = rep( c(-Inf,rep(.01,(nThresholds-1))) , nVariables),
dimnames = list(c(), fruitynames)),
mxMatrix("Lower",nThresholds,nThresholds,values=1,free=F,name="unitLower"),
mxAlgebra(unitLower %*% thresholdDeviations, name="thresholdMatrix"),
mxFitFunctionML(),mxExpectationLISREL(LX='L', TX='M', PH='P', TD='E', KA='K', thresholds="thresholdMatrix")
)
lis$thresholdDeviations$ubound <- 1
# Generate data based on model
ordinalData <- mxGenerateData(lis, nSubjects)
# Put the simulated data in the model and run
lis <- mxModel(lis, mxData(observed=ordinalData, type='raw'))
lisr <- mxRun(lis)
omxCheckCloseEnough(lisr$output$fit, 3858.052, .01)
# Compute factor scores for the model
lism <- mxFactorScores(lisr)
lisw <- mxFactorScores(lisr, 'WeightedML')
omxCheckError(lisreg <- mxFactorScores(lisr, 'Regression'), "Regression factor scores cannot be computed when there are thresholds (ordinal data).")
mask <- abs(lism[,,1]) < 5
omxCheckCloseEnough(cor(lism[mask,,1], lisw[mask,,1]), 1, 0.01)
#pdf('plotOrdinalFactorScores.pdf')
#plot(lism[,,1], main='Ordinal Factor Scoring Methods', xlab='Sorted Data Row', ylab='Factor Score')
#points(lisw[,,1], col='blue')
#legend('bottomright', legend=c('ML', 'Weighted ML', paste('r =', round(cor(lism[,,1], lisw[,,1]), 2))), pch=c(1, 1, NA), col=c('black', 'blue', 'black'))
#dev.off()
#------------------------------------------------------------------------------
# Frontpage model in RAM form
require(OpenMx)
data(demoOneFactor)
manifests <- names(demoOneFactor)
latents <- c("G")
factorModel <- mxModel("OneFactor",
type="RAM",
manifestVars = manifests,
latentVars = latents,
mxPath(from=latents, to=manifests),
mxPath(from=manifests, arrows=2),
mxPath(from=latents, arrows=2,
free=FALSE, values=1.0),
mxPath(from='one', to=manifests),
mxData(observed=demoOneFactor, type="raw"))
summary(factorRamRun <- mxRun(factorModel))
rr1 <- mxFactorScores(factorRamRun, 'ML')
rr2 <- mxFactorScores(factorRamRun, 'Regression')
rr3 <- mxFactorScores(factorRamRun, 'WeightedML')
# Compare RAM factor scores to LISREL
omxCheckCloseEnough(cor(rr1[,,1], r1[,,1]), 1, 1e-8)
omxCheckCloseEnough(cor(rr2[,,1], r2[,,1]), 1, 1e-8)
omxCheckCloseEnough(cor(rr3[,,1], r3[,,1]), 1, 1e-8)
rms <- function(x, y){sqrt(mean((x-y)^2))}
omxCheckCloseEnough(rms(rr1[,,1], r1[,,1]), 0, .1)
omxCheckCloseEnough(rms(rr2[,,1], r2[,,1]), 0, .1)
omxCheckCloseEnough(rms(rr3[,,1], r3[,,1]), 0, .1)
# Compare RAM standard errors to LISREL
omxCheckCloseEnough(rms(rr1[,,2], r1[,,2]), 0, .001)
omxCheckCloseEnough(rms(rr2[,,2], r2[,,2]), 0, .001)
omxCheckCloseEnough(rms(rr3[,,2], r3[,,2]), 0, .001)
#------------------------------------------------------------------------------
# Multiple group example
factorModel2 <- mxModel(factorModel, name='ConstrainedResidual',
mxPath(from=manifests, arrows=2, labels='resid')
)
twoGroup <- mxModel(model='multipleGroup', factorModel, factorModel2,
mxFitFunctionMultigroup(c('OneFactor', 'ConstrainedResidual'))
)
twoGroupRun <- mxRun(twoGroup)
tg1 <- mxFactorScores(twoGroupRun, 'ML')
tg2 <- mxFactorScores(model=twoGroupRun,"Regression")
omxCheckCloseEnough(rms(tg1[[1]], rr1), 0, .001)
omxCheckCloseEnough(rms(tg1[[2]], rr1), 0, .01)
omxCheckCloseEnough(rms(tg2[[1]], rr2), 0, .001)
omxCheckCloseEnough(rms(tg2[[2]], rr2), 0, .01)
#-------------------------------------------------------------------
#Ensure regression RAM scoring does not fail in the presence of missing data:
demoOneFactor.miss <- as.matrix(demoOneFactor)
demoOneFactor.miss[sample(1:2500,size=100,replace=F)] <- NA
demoOneFactor.miss[100,] <- NA
factorModel.miss <- mxModel("OneFactor",
type="RAM",
manifestVars = manifests,
latentVars = latents,
mxPath(from=latents, to=manifests),
mxPath(from=manifests, arrows=2),
mxPath(from=latents, arrows=2,
free=FALSE, values=1.0),
mxPath(from='one', to=manifests),
mxData(observed=demoOneFactor.miss, type="raw"))
factorRamRun.miss <- mxRun(factorModel.miss)
for (type in c("Regression", "WeightedML", "ML")) {
omxCheckError(mxFactorScores(model=factorRamRun.miss,type),
"mxFactorScores: row 8 has missing data. Hence, you must specify minManifests")
}
regs <- mxFactorScores(model=factorRamRun.miss,"Regression", minManifests=3)
omxCheckTrue(is.na(regs[100,1,1]))
omxCheckTrue(is.na(regs[100,1,2]))
omxCheckTrue(cor(regs[,,1], rr2[,,1], use="complete.obs") > 0.95)
#Ensure regression LISREL scoring does not fail in the presence of missing data:
factorModelL.miss <- mxModel(name="OneFactor", factorModelL,
mxData(observed=demoOneFactor.miss, type="raw"))
factorRunL.miss <- mxRun(factorModelL.miss)
for (type in c("Regression", "WeightedML", "ML")) {
omxCheckError(mxFactorScores(model=factorRunL.miss,type),
"mxFactorScores: row 8 has missing data. Hence, you must specify minManifests")
}
regsl <- mxFactorScores(model=factorRunL.miss,"Regression", minManifests=3)
omxCheckTrue(is.na(regsl[100,1,1]))
omxCheckTrue(is.na(regsl[100,1,2]))
omxCheckTrue(cor(regsl[,,1], r2[,,1], use="complete.obs") > 0.95)
#------------------------------------------------------------------------------
# Check some data generation for LISREL too
nvar <- 6
varnames <- paste("x",1:nvar,sep="")
# specify model with out means and try to generate data
ssModelLisrel <- mxModel(model="lisrel",
mxMatrix("Full",6,6,TRUE,.2,name="BE"),
mxMatrix("Full",6,6,TRUE,.5,name="LY",dimnames=list(varnames,varnames)),
mxMatrix("Diag",6,6,FALSE,1,name="PS"),
mxMatrix("Diag",6,6,FALSE,1,name="TE"),
mxExpectationLISREL(BE="BE",LY="LY",PS="PS",TE="TE"),
mxFitFunctionML()
)
expect_warning(mxGetExpected(ssModelLisrel, 'means'),
"Means requested, but model has no means")
# Add means and it works fine
ssModelLisrel <- mxModel(ssModelLisrel, name="lisrel",
mxMatrix("Full",6,1,FALSE,2,name="AL"),
mxMatrix("Full",6,1,FALSE,2,name="TY"),
mxExpectationLISREL(BE="BE",LY="LY",PS="PS",TE="TE",AL="AL",TY="TY")
)
# Check expectations
omxCheckCloseEnough(mxGetExpected(ssModelLisrel, 'means'), matrix(-28, 6), 1e-10)
cc <- matrix(37.5, 6, 6)
diag(cc) <- 38.5
omxCheckCloseEnough(mxGetExpected(ssModelLisrel, 'covariance'), cc, 1e-10)
# Check data generation
set.seed(72)
ssDataLisrel <- mxGenerateData(ssModelLisrel, 200)
omxCheckEquals(dim(ssDataLisrel), c(200, 6))
rms <- function(x, y){sqrt(mean((x-y)^2))}
omxCheckTrue(rms(colMeans(ssDataLisrel), rep(-28, 6)) < 1)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.