skip_on_cran()
if (require(future)) {
N <- 2000
library(semtree)
library(future)
future::plan(multisession, workers=5)
observed <- MASS::mvrnorm(n=N, mu=c(0,0), Sigma=matrix(c(1,0.2,0.2,1),nrow=2))
pred1 <- rnorm(N)
pred2 <- rnorm(N)
noise <- rnorm(N)
observed[,1] <- observed[,1]+ifelse(pred1>0,2,0)
observed[,2] <- observed[,2]+ifelse(pred2>0,2,0)
datf <- data.frame(x1=observed[,1], x2=observed[,2],pred1, pred2, noise)
#
manifests<-c("x1","x2")
model.biv <- mxModel("Bivariate_Model",
type="RAM",
manifestVars = manifests,
latentVars = c(),
mxPath(from="x1",to=c("x1","x2"),
free=c(TRUE,TRUE), value=c(1.0,.2) ,
arrows=2, label=c("VAR_x1","COV_x1_x2") ),
mxPath(from="x2",to=c("x2"), free=c(TRUE),
value=c(1.0) , arrows=2, label=c("VAR_x2") ),
mxPath(from="one",to=c("x1","x2"), label=c("mu1","mu2"),
free=TRUE, value=0, arrows=1),
mxData(datf, type = "raw")
);
result <- mxRun(model.biv)
summary(result)
# ---
forest <- semforest(model=model.biv, data=datf,
constraints = semtree.constraints(focus.parameters="mu2"),
control =
semforest.control(num.trees = 30,
control=semtree.control(min.N = 50, min.bucket = 100,alpha=1,
method="score")))
vim_naive <- varimp(forest)
vim_focus <- varimp(forest, method="permutationFocus")
vimdat <- data.frame( vim=rep(c("naive","focus"),each=3),
param=rep(c("pred1","pred2","noise"),2),
vals=c( semtree:::aggregateVarimp(vim_naive),semtree:::aggregateVarimp(vim_focus)))
library(ggplot2)
library(tidyr)
vimdat %>% ggplot(aes(x=vim,y=vals,group=param,fill=param))+geom_col(position="dodge")
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.