How to get estimates of latent random slopes and random intercept
Here is code to simulate HLM data
library(dplyr)
library(tidyr)
library(OpenMx)
# Simulate Data -------------
J <- 100 #number of clusters
n_j <- 30#number of observations
gamma_00 <- 5
gamma_10 <- 3
mu_x <- 4
sim.dat <- data.frame(J_ID = 1:J) |>
mutate(U_0j = rnorm(J, mean = 0 , sd = 1),
U_1j = rnorm(J, mean = 0, sd = 1),
gamma_00 = gamma_00,
gamma_10 = gamma_10,
n_j = n_j) |>
rowwise() |>
mutate(X = list(rnorm(n_j, mu_x, sd = 1)),
R = list(rnorm(n_j, 0, sd = 1))) |>
unnest(cols = c(X, R)) |>
mutate(Y = gamma_00 + U_0j + (gamma_10 + U_1j) * X + R)
sim.dat.obs <- sim.dat|>
select(J_ID, X, Y) |>
mutate(J_ID = as.factor(J_ID))
Here is code for openmx
## J level (Teacher) --------
J_obs <- sim.dat.obs |>
select(J_ID) |>
distinct()
J_cluster_name <- "J_ID"
J_latent <- c("U_0", "U_1")
J_data <- mxData(J_obs, type = "raw",
primaryKey = J_cluster_name)
J_var <- mxPath(from = J_latent,
arrows = 2, values = 1,
connect = "single")
J_cov <- mxPath(from = J_latent,
arrows = 2, values = 0,
connect = "unique.bivariate",
free = FALSE)
J_model <- mxModel("Level_J", type = "RAM",
latentVars = J_latent,
J_data,
J_var,
J_cov)
## I level (Student) --------
I_obs <- mxData(sim.dat.obs, type = "raw")
I_manifest <- c("X", "Y")
I_means <- mxPath(from = "one",
to = c("X"),
labels = c("mu_x"),
value = 1)
I_var <- mxPath(from = c("X", "Y"),
labels = c("sd_x", "R_ysd"),
arrow = 2,
value = 1)
I_fixed <- mxPath(from = c("one","X"),
to = "Y",
value = 1,
labels = c("gamma_00", "gamma_10"))
I_randomslope <- mxPath(from = c("Level_J.U_1"),
to = c("Y"),
labels = c("data.X"),
free = FALSE,
joinKey = J_cluster_name)
I_randomint <- mxPath(from = c("Level_J.U_0"),
to = c("Y"),
value = 1,
free = FALSE,
joinKey = J_cluster_name)
I_model<- mxModel("Level_I", type = "RAM",
manifestVars = I_manifest,
J_model,
I_obs,
I_means,
I_var,
I_fixed,
I_randomslope,
I_randomint)
I.fit <- mxRun(I_model)
#> Running Level_I with 7 parameters
summary(I.fit)
#> Summary of Level_I
#>
#> free parameters:
#> name matrix row col Estimate Std.Error A
#> 1 gamma_10 A Y X 3.0702452 0.09300885
#> 2 sd_x S X X 0.9924763 0.02562636
#> 3 R_ysd S Y Y 1.0691329 0.02857744
#> 4 mu_x M 1 X 4.0054015 0.01818861
#> 5 gamma_00 M 1 Y 5.0011756 0.13307218
#> 6 Level_J.S[1,1] Level_J.S U_0 U_0 1.1217698 0.24712431 !
#> 7 Level_J.S[2,2] Level_J.S U_1 U_1 0.8267078 0.12055658
#>
#> Model Statistics:
#> | Parameters | Degrees of Freedom | Fit (-2lnL units)
#> Model: 7 5993 17907.01
#> Saturated: NA NA NA
#> Independence: NA NA NA
#> Number of observations/statistics: 3100/6000
#>
#> Information Criteria:
#> | df Penalty | Parameters Penalty | Sample-Size Adjusted
#> AIC: 5921.013 17921.01 17921.05
#> BIC: -30271.658 17963.29 17941.04
#> CFI: NA
#> TLI: 1 (also known as NNFI)
#> RMSEA: 0 [95% CI (NA, NA)]
#> Prob(RMSEA <= 0.05): NA
#> To get additional fit indices, see help(mxRefModels)
#> timestamp: 2024-05-03 20:11:30
#> Wall clock time: 9.318122 secs
#> optimizer: SLSQP
#> OpenMx version number: 2.21.8
#> Need help? See help(mxSummary)
mxFactorScores(I.fit, "WeightedML")
#> Running Container with 0 parameters
#> Error in `[.data.frame`(got, , paste0(names(coef(fit)), "SE"), drop = FALSE): undefined columns selected
mxFactorScores(I.fit@submodels$Level_J, "WeightedML")
#> Error in mxEvalByName(model$expectation$M, model, compute = TRUE): 'name' argument must be a character argument
One work around I tried was to specify dummy latent variables which are the Level_J.U1 times one at the level 1 but the results were not coherent. Observations within a cluster had different values of U_1.
I also asked on Stackoverflow
Log in or register to post comments
Not Implemented
Unfortunately, no such feature is currently implemented. In your specific case, it wouldn't be terribly hard for OpenMx to figure out the factor scores, but I don't think the general case (arbitrarily many levels with multiple groups and different factor structures everywhere) is solved. I don't see this feature being implemented anytime soon.
I've created an Issue on GitHub to give a better error message though: [https://github.com/OpenMx/OpenMx/issues/391](https://github.com/OpenMx/OpenMx/issues/391)
Log in or register to post comments