demo/partialDependence2.R

In brandmaier/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)

## - single tree -
tree<-semtree(model, data)

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

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

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

#semforest(model = model, data = data, constraints=constraints)

vim <- varimp(forest)

print(vim, aggregate="median")
print(vim, aggregate="mean")

pd1.mu <- partialDependence(forest, "p1", "mux1")
pd1.var <- partialDependence(forest, "p1", "VAR_x1")
pd3.mu <- partialDependence(forest, "p3", "mux1")
pd3.var <- partialDependence(forest, "p3", "VAR_x1")
plot(pd)

opar <- par(no.readonly = TRUE)

par(mfrow=c(2,2))
plot(pd1.mu)
title("Mu pred1")
plot(pd1.var)
title("Var pred1")
plot(pd3.mu)
title("Mu pred3")
plot(pd3.var)
title("Var pred3")

par(opar)