R Under development (unstable) (2023-11-25 r85635 ucrt) -- "Unsuffered Consequences" Copyright (C) 2023 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > # > # testing invariance > # > > # simulation: > # factor structure holds over SES but not over age > # both SES and age predict a mean difference in cognitive outcome > > set.seed(123) > > require("semtree") Loading required package: semtree Loading required package: OpenMx > require("lavaan") Loading required package: lavaan This is lavaan 0.6-16 lavaan is FREE software! Please report any bugs. > # > lambda <- list() > age <- c(0,0,1,1) # 0=young, 1 =old > ses <- c(0,1,0,1) # 0=low, 1=high > lambda[[1]] <- c(1,0.9,0.8,0.8) > lambda[[2]] <- c(1,0.9,0.8,0.8) > lambda[[3]] <- c(1,0.4,0.2,0.9) > lambda[[4]] <- c(1,0.4,0.2,0.9) > cogmean <- 80-age*30 + ses*20 > cogsd <- 1 > errsd <- 1 > > Nsub <- 200 # persons per agexSES group > > cbind(age,ses,cogmean) age ses cogmean [1,] 0 0 80 [2,] 0 1 100 [3,] 1 0 50 [4,] 1 1 70 > > # simulate data from 4 groups > data <- c() > for (i in 1:4) { + cogsim <- rnorm(n = Nsub,mean = cogmean[i],cogsd) + scores <- as.matrix(t(outer(lambda[[i]],cogsim))) + rnorm(Nsub*4,0,errsd) + data <- rbind(data,scores) + } > > # rescale data > data <- scale(data) > > # put data into dataframe and label observed variables > data <- data.frame(data) > names(data) <- paste0("x",1:4) > > # add predictors to create full data set > fulldata <- cbind(data, age=factor(rep(age,each=Nsub)),ses=factor(rep(ses,each=Nsub))) > > # specify model > model<-" + ! regressions + F=~1.0*x1 + F=~l2*x2 + F=~l3*x3 + F=~l4*x4 + ! residuals, variances and covariances + x1 ~~ VAR_x1*x1 + x2 ~~ VAR_x2*x2 + x3 ~~ VAR_x3*x3 + x4 ~~ VAR_x4*x4 + F ~~ 1.0*F + ! means + F~1 + x1~0*1; + x2~0*1; + x3~0*1; + x4~0*1; + "; > result<-lavaan(model, data=data, fixed.x=FALSE, missing="FIML"); Warning message: In lav_object_post_check(object) : lavaan WARNING: some estimated ov variances are negative > > manifests<-c("x1","x2","x3","x4") > latents<-c("F") > model <- mxModel("Unnamed_Model", + type="RAM", + manifestVars = manifests, + latentVars = latents, + mxPath(from="F",to=c("x1","x2","x3","x4"), free=c(FALSE,TRUE,TRUE,TRUE), + value=c(1.0,1.0,1.0,1.0) , arrows=1, label=c("F__x1","l2","l3","l4") ), + mxPath(from="one",to=c("F"), free=c(TRUE), value=c(1.0) , arrows=1, label=c("const__F") ), + mxPath(from="one",to=c("x2","x3","x4"), free=c(TRUE,TRUE,TRUE), value=c(1.0,1.0,1.0) , arrows=1, label=c("const__x2","const__x3","const__x4") ), + mxPath(from="x1",to=c("x1"), free=c(TRUE), value=c(1.0) , arrows=2, label=c("VAR_x1") ), + mxPath(from="x2",to=c("x2"), free=c(TRUE), value=c(1.0) , arrows=2, label=c("VAR_x2") ), + mxPath(from="x3",to=c("x3"), free=c(TRUE), value=c(1.0) , arrows=2, label=c("VAR_x3") ), + mxPath(from="x4",to=c("x4"), free=c(TRUE), value=c(1.0) , arrows=2, label=c("VAR_x4") ), + mxPath(from="F",to=c("F"), free=c(FALSE), value=c(1.0) , arrows=2, label=c("VAR_F") ), + mxPath(from="one",to=c("x1"), free=F, value=0, arrows=1), + mxData(data[1:50,], type = "raw") + ); > > result <- mxRun(model) Running Unnamed_Model with 11 parameters Warning message: In model 'Unnamed_Model' Optimizer returned a non-zero status code 6. The model does not satisfy the first-order optimality conditions to the required accuracy, and no improved point for the merit function could be found during the final linesearch (Mx status RED) > summary(result) Summary of Unnamed_Model The model does not satisfy the first-order optimality conditions to the required accuracy, and no improved point for the merit function could be found during the final linesearch (Mx status RED) free parameters: name matrix row col Estimate Std.Error A 1 l2 A x2 F 0.110681832 0.0700818456 2 l3 A x3 F 0.258542310 0.0784017931 3 l4 A x4 F 0.577944014 0.1992752561 4 VAR_x1 S x1 x1 0.004467686 0.0035868467 5 VAR_x2 S x2 x2 0.001213089 0.0002540675 6 VAR_x3 S x3 x3 0.001042346 0.0003111305 7 VAR_x4 S x4 x4 0.006968799 0.0018057267 8 const__x2 M 1 x2 0.629667879 0.0201839921 9 const__x3 M 1 x3 0.641788444 0.0224711009 10 const__x4 M 1 x4 -0.081926307 0.0570797980 11 const__F M 1 F 0.278899151 0.1417367901 Model Statistics: | Parameters | Degrees of Freedom | Fit (-2lnL units) Model: 11 189 -387.4425 Saturated: 14 186 NA Independence: 8 192 NA Number of observations/statistics: 50/200 Information Criteria: | df Penalty | Parameters Penalty | Sample-Size Adjusted AIC: -765.4425 -365.4425 -358.4951 BIC: -1126.8149 -344.4103 -378.9374 To get additional fit indices, see help(mxRefModels) timestamp: 2023-11-26 15:11:32 Wall clock time: 0.1311638 secs optimizer: SLSQP OpenMx version number: 2.21.10 Need help? See help(mxSummary) > > subset <- data[1:50, ] > > > ctr <- semtree.control(verbose=TRUE) > ctr$exclude.heywood <- FALSE > # naive tree should find both effects, age & ses, with age having the stronger effect > tree <- semtree(model = result, data=fulldata, control=ctr) Detected OpenMx model. OpenMx model estimation selected! ❯ Growing level 0 n = 800 Beginning initial fit attempt Fit attempt 0, fit=2843.02725567401, new current best! (was 1281526.09932206) ❯ Testing Predictor: age ❯ Testing Predictor: ses ✔ Best split is age with statistic = 8512.9 ❯ Growing level 1 n = 400 Beginning initial fit attempt Fit attempt 0, fit=-2928.40985127395, new current best! (was 1727.38609524129) Beginning fit attempt 1 of at maximum 10 extra tries Fit attempt 1, fit=-2928.40985182956, new current best! (was -2928.40985127395) Beginning fit attempt 2 of at maximum 10 extra tries Fit attempt 2, fit=-2928.40985182957, new current best! (was -2928.40985182956) Beginning fit attempt 3 of at maximum 10 extra tries Beginning fit attempt 4 of at maximum 10 extra tries Beginning fit attempt 5 of at maximum 10 extra tries Beginning fit attempt 6 of at maximum 10 extra tries Beginning fit attempt 7 of at maximum 10 extra tries Beginning fit attempt 8 of at maximum 10 extra tries Beginning fit attempt 9 of at maximum 10 extra tries Beginning fit attempt 10 of at maximum 10 extra tries ❯ Testing Predictor: age ❯ Testing Predictor: ses ✔ Best split is ses with statistic = 219.28 ❯ Growing level 2 n = 200 Beginning initial fit attempt Fit attempt 0, fit=-1577.85039131669, new current best! (was -1479.85059188343) Beginning fit attempt 1 of at maximum 10 extra tries Fit attempt 1, fit=-1577.85039192277, new current best! (was -1577.85039131669) Beginning fit attempt 2 of at maximum 10 extra tries Beginning fit attempt 3 of at maximum 10 extra tries Beginning fit attempt 4 of at maximum 10 extra tries Beginning fit attempt 5 of at maximum 10 extra tries Beginning fit attempt 6 of at maximum 10 extra tries Beginning fit attempt 7 of at maximum 10 extra tries Beginning fit attempt 8 of at maximum 10 extra tries Beginning fit attempt 9 of at maximum 10 extra tries Beginning fit attempt 10 of at maximum 10 extra tries ❯ Testing Predictor: age ❯ Testing Predictor: ses ❯ Best Likelihood Ratio was NULL. Stop splitting ❯ Growing level 2 n = 200 Beginning initial fit attempt Fit attempt 0, fit=-1569.84006709087, new current best! (was -1448.5592599462) Beginning fit attempt 1 of at maximum 10 extra tries Fit attempt 1, fit=-1569.84006709889, new current best! (was -1569.84006709087) Beginning fit attempt 2 of at maximum 10 extra tries Fit attempt 2, fit=-1569.84006768963, new current best! (was -1569.84006709889) Beginning fit attempt 3 of at maximum 10 extra tries Beginning fit attempt 4 of at maximum 10 extra tries Beginning fit attempt 5 of at maximum 10 extra tries Beginning fit attempt 6 of at maximum 10 extra tries Beginning fit attempt 7 of at maximum 10 extra tries Beginning fit attempt 8 of at maximum 10 extra tries Beginning fit attempt 9 of at maximum 10 extra tries Beginning fit attempt 10 of at maximum 10 extra tries ❯ Testing Predictor: age ❯ Testing Predictor: ses ❯ Best Likelihood Ratio was NULL. Stop splitting ❯ Growing level 1 n = 400 Beginning initial fit attempt Fit attempt 0, fit=-2741.46632160248, new current best! (was 1115.64116043285) Beginning fit attempt 1 of at maximum 10 extra tries Beginning fit attempt 2 of at maximum 10 extra tries Beginning fit attempt 3 of at maximum 10 extra tries Beginning fit attempt 4 of at maximum 10 extra tries Beginning fit attempt 5 of at maximum 10 extra tries Beginning fit attempt 6 of at maximum 10 extra tries Beginning fit attempt 7 of at maximum 10 extra tries Beginning fit attempt 8 of at maximum 10 extra tries Beginning fit attempt 9 of at maximum 10 extra tries Beginning fit attempt 10 of at maximum 10 extra tries ❯ Testing Predictor: age ❯ Testing Predictor: ses ✔ Best split is ses with statistic = 289.57 ❯ Growing level 2 n = 200 Beginning initial fit attempt Fit attempt 0, fit=-1496.78521317663, new current best! (was -1362.46535446136) Beginning fit attempt 1 of at maximum 10 extra tries Fit attempt 1, fit=-1496.7852132728, new current best! (was -1496.78521317663) Beginning fit attempt 2 of at maximum 10 extra tries Fit attempt 2, fit=-1496.78521340084, new current best! (was -1496.7852132728) Beginning fit attempt 3 of at maximum 10 extra tries Beginning fit attempt 4 of at maximum 10 extra tries Beginning fit attempt 5 of at maximum 10 extra tries Beginning fit attempt 6 of at maximum 10 extra tries Beginning fit attempt 7 of at maximum 10 extra tries Beginning fit attempt 8 of at maximum 10 extra tries Beginning fit attempt 9 of at maximum 10 extra tries Beginning fit attempt 10 of at maximum 10 extra tries ❯ Testing Predictor: age ❯ Testing Predictor: ses ❯ Best Likelihood Ratio was NULL. Stop splitting ❯ Growing level 2 n = 200 Beginning initial fit attempt Fit attempt 0, fit=-1534.25522765984, new current best! (was -1379.00096714098) Beginning fit attempt 1 of at maximum 10 extra tries Fit attempt 1, fit=-1534.25522768047, new current best! (was -1534.25522765984) Beginning fit attempt 2 of at maximum 10 extra tries Beginning fit attempt 3 of at maximum 10 extra tries Beginning fit attempt 4 of at maximum 10 extra tries Beginning fit attempt 5 of at maximum 10 extra tries Beginning fit attempt 6 of at maximum 10 extra tries Beginning fit attempt 7 of at maximum 10 extra tries Beginning fit attempt 8 of at maximum 10 extra tries Beginning fit attempt 9 of at maximum 10 extra tries Beginning fit attempt 10 of at maximum 10 extra tries ❯ Testing Predictor: age ❯ Testing Predictor: ses ❯ Best Likelihood Ratio was NULL. Stop splitting ✔ Tree construction finished [took 4s]. > #plot(tree) > > # invariance tree should exclude splits wrt age and only splir wrt to ses > ctr$report.level <- 99 > ctr$alpha.invariance <- 0.01 > > cnst <- semtree.constraints(local.invariance = c("l2","l3","l4")) > #cnst <- semtree.constraints(local.invariance=) > tree2 <- semtree(model = result, data=fulldata, control=ctr, constraints = cnst ) Detected OpenMx model. OpenMx model estimation selected! ❯ Growing level 0 n = 800 Growing tree level 0 Beginning initial fit attempt Fit attempt 0, fit=2843.02725567401, new current best! (was 1281526.09932206) |--Estimating baseline likelihood: 2843.02725567297 ❯ Testing Predictor: age |--Estimating baseline likelihood: 2843.02725567297 ❯ Testing Predictor: ses ❯ Filtering possible splits by invariance LL.within: 8512.9034285494,6670.01382649603 LL.max: NA within.split: 1,1 split max NA ❯ Best Likelihood Ratio was NULL. Stop splitting ✔ Tree construction finished [took 1s]. Warning messages: 1: In model 'sharedModel' Optimizer returned a non-zero status code 6. The model does not satisfy the first-order optimality conditions to the required accuracy, and no improved point for the merit function could be found during the final linesearch (Mx status RED) 2: In model 'sharedModel' Optimizer returned a non-zero status code 6. The model does not satisfy the first-order optimality conditions to the required accuracy, and no improved point for the merit function could be found during the final linesearch (Mx status RED) > > #plot(tree2) > > proc.time() user system elapsed 7.35 0.28 7.64