#
# 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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.