# Create and fit a WLS model using RAM, and then using matrices.
library(OpenMx)
# Simulate some data where y = .5x + error
x = rnorm(1000, mean = 0, sd = 1)
y = 0.5*x + rnorm(1000, mean = 0, sd = 1)
tmpFrame = data.frame(x, y)
varNames = names(tmpFrame)
# =======================
# = A RAM model example =
# =======================
m1 = mxModel("my_first_WLS", type = "RAM",
manifestVars = c("x", "y"),
mxPath(c("x", "y"), arrows = 2, values = 1, labels = c("xVar", "yVar")),
mxPath("x", to = "y", labels = "x_to_y"),
mxFitFunctionWLS(),
mxData(tmpFrame, 'raw')
)
m1 = mxRun(m1)
summary(m1)$parameters
# Here are the cov, acov and Weight matrices:
print(m1$data$observedStats)
# Use a different weight matrix
m2 = m1
os <- m1$data$observedStats
os$asymCov <- solve(rWishart(n=1, df= nrow(tmpFrame), Sigma= diag(3))[,,1])
os$useWeight <- solve(os$asymCov * nrow(tmpFrame))
m2$data$observedStats <- os
# Set verbose to check if our new weights are used
m2$data$verbose <- 1L
# Run model
m2 <- mxRun(m2)
# SE indeed changed due to new weights
print(m2$output$standardErrors - m1$output$standardErrors)
# ==========================
# = A matrix-based example =
# ==========================
# Define matrices for Symmetric (S) and Asymmetric (A) paths and an Identity matrix.
S <- mxMatrix(type = "Full", nrow = 2, ncol = 2, values=c(1,0,0,1),
free=c(TRUE,FALSE,FALSE,TRUE), labels=c("Vx", NA, NA, "Vy"), name = "S")
A <- mxMatrix(type = "Full", nrow = 2, ncol = 2, values=c(0,1,0,0),
free=c(FALSE,TRUE,FALSE,FALSE), labels=c(NA, "b", NA, NA), name = "A")
I <- mxMatrix(type="Iden", nrow=2, ncol=2, name="I")
# Build the model
tmpModel <- mxModel(model="exampleModel",
# Add the S, A, and I matrices constructed above
S, A, I,
# Define the expectation
mxAlgebra(name="expCov", solve(I-A) %*% S %*% t(solve(I-A))),
# Choose a normal expectation and WLS as the fit function
mxExpectationNormal(covariance= "expCov", dimnames= varNames),
mxFitFunctionWLS(),
# Add the data
mxData(tmpFrame, 'raw')
)
# Fit the model and print a summary
tmpModel <- mxRun(tmpModel)
summary(tmpModel)
Run the code above in your browser using DataLab