demo/focusParameter.R

In semtree: Recursive Partitioning for Structural Equation Models

#
# focusParameter.R
# ----------------
# demo script
#
# we simulate three subgroups with bivariate observations. Observations are equally correlated (r=0.5) in all three 
# groups and have equal variance. Groups only differ w.r.t their means.
#
#
require("semtree")
require("future")
plan(multisession)

N.sub <- 200

# simulate multivariate data

require("MASS")
smu <- 3
sub1 <- mvrnorm(n=N.sub,mu = c(0,0), Sigma=matrix(data = c(1,0.5,
                                                   0.5,1),nrow=2,byrow = TRUE))
sub2 <- mvrnorm(n=N.sub,mu = c(1*smu,1*smu), Sigma=matrix(data = c(1,0.5,
                                                           0.5,1),nrow=2,byrow = TRUE))
sub3 <- mvrnorm(n=N.sub,mu = c(-1*smu,-1*smu), Sigma=matrix(data = c(1,0.5,
                                                           0.5,1),nrow=2,byrow = TRUE))

groups <- rep(1:3, each=N.sub)

pred1 <- factor(rep(c(1,1,0), each=N.sub))
pred2 <- factor(rep(c(1,0,1), each=N.sub))
pred3 <- sample(c(0,1), N.sub*3, replace = TRUE)

data <- data.frame(rbind(sub1,sub2,sub3), pred1, pred2, pred3)
names(data) <- c("x1","x2","p1","p2","p3")

opar <- par(no.readonly = TRUE)

par(mfrow=c(2,1))

plot(data)
plot(data,col=c("orange","black","red")[groups],pch=20,cex=1)

par(opar)



#
# This model specification was automatically generated by Onyx
# 
require("OpenMx");
manifests<-names(data)[1:2]
latents<-c()
model <- mxModel("Unnamed_Model", 
                 type="RAM",
                 manifestVars = manifests,
                 latentVars = latents,
                 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","x1"), free=c(TRUE,TRUE), value=c(1.0,0.1) , arrows=2, label=c("VAR_x2","COV_x2_x1") ),
                 mxPath(from="one",to=c("x1","x2"), free=TRUE, value=0, arrows=1,
                        label=c("mux1","mux2")),
                 mxData(data, type = "raw")
);

result <- mxRun(model)
summary(result)


ctrl <- semforest.control(num.trees = 210)

constraints1 <- semtree.constraints(focus.parameters = "mux1")
constraints2 <- semtree.constraints(focus.parameters = "VAR_x1")


forest <- semforest(model = model, data = data, 
                    constraints=NULL, control = ctrl)

forest.constrained1 <- semforest(model = model, data = data, 
                                constraints=constraints1, control = ctrl)
forest.constrained2 <- semforest(model = model, data = data, 
                                 constraints=constraints2, control = ctrl)


vim.constrained1 <- varimp(forest.constrained1)
vim.constrained2 <- varimp(forest.constrained2)

vim <- varimp(forest)

opar <- par(no.readonly = TRUE)
par(mfrow=c(3,1))
plot(vim)
plot(vim.constrained1)
plot(vim.constrained2)
par(opar)