Nothing
## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
## ----setup, warning=FALSE, message=FALSE,echo=FALSE---------------------------
library(semtree)
library(MASS)
## ----eval=FALSE, message=FALSE, warning=FALSE, results="hide"-----------------
# library(semtree)
#
# cnst <- semtree.constraints(local.invariance=NULL,
# global.invariance=NULL,
# focus.parameters=NULL)
#
# semtree(model.x, data=df, constraints=cnst)
## ----difsim-------------------------------------------------------------------
N <- 2000
p1 <- sample(size = N,
x = c(0, 1),
replace = TRUE)
p2 <- sample(size = N,
x = c(0, 1),
replace = TRUE)
lat <- rnorm(N, mean = 0 + p1)
loadings <- c(.5, .8, .7, .9)
observed <- lat %*% t(loadings) + rnorm(N * length(loadings), sd = .1)
observed[, 3] <- observed[, 3] + p2 * 0.5 * lat
cfa.sim <- data.frame(observed, p1 = factor(p1), p2 = factor(p2))
names(cfa.sim)[1:4] <- paste0("x", 1:4)
## ----cfadefinition, echo=TRUE-------------------------------------------------
require("OpenMx");
manifests<-c("x1","x2","x3","x4")
latents<-c("F")
model.cfa <- mxModel("CFA", type="RAM", manifestVars = manifests,
latentVars = latents,
mxPath(from="F",to=c("x1","x2","x3","x4"),
free=c(TRUE,TRUE,TRUE,TRUE), value=c(1.0,1.0,1.0,1.0) ,
arrows=1, label=c("F__x1","F__x2","F__x3","F__x4") ),
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="one",to=c("F"), free=c(TRUE),
value=c(1.0) , arrows=1, label=c("const__F") ),
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(cfa.sim, type = "raw")
);
## ----message=FALSE, echo=TRUE, warning=FALSE----------------------------------
tree.gc <- semtree(model.cfa, data=cfa.sim, constraints=
semtree.constraints(global.invariance =
c("F__x1","F__x2","F__x3","F__x4")))
## ----fig.height=3-------------------------------------------------------------
plot(tree.gc)
## ----warning=FALSE, message=FALSE, error=FALSE, echo=TRUE---------------------
tree.lc <- semtree(model.cfa, data=cfa.sim, constraints=
semtree.constraints(
local.invariance= c("F__x1","F__x2","F__x3","F__x4")))
## -----------------------------------------------------------------------------
plot(tree.lc)
## -----------------------------------------------------------------------------
set.seed(123)
N <- 1000
grp1 <- sample(x = c(0,1), size=N, replace=TRUE)
grp2 <- sample(x = c(0,1), size=N, replace=TRUE)
Sigma <- matrix(byrow=TRUE,
nrow=2,c(2,0.2,
0.2,1))
obs <- MASS::mvrnorm(N,mu=c(0,0),
Sigma=Sigma)
obs[,1] <- obs[,1] + ifelse(grp1,3,0)
obs[,2] <- obs[,2] + ifelse(grp2,3,0)
df.biv <- data.frame(obs, grp1=factor(grp1), grp2=factor(grp2))
names(df.biv)[1:2] <- paste0("x",1:2)
## -----------------------------------------------------------------------------
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(df.biv, type = "raw")
);
## -----------------------------------------------------------------------------
result <- mxRun(model.biv)
summary(result)
## ----growbivtree, warning=FALSE, message=FALSE--------------------------------
tree.biv <- semtree(model.biv, data=df.biv)
## -----------------------------------------------------------------------------
# default white color for all nodes
cols <- rep("black", semtree:::getNumNodes(tree.biv))
cols[as.numeric(row.names(semtree:::getTerminalNodes(tree.biv)))] <- viridis:::viridis_pal()(4)
plot(tree.biv, border.col=cols)
## ----warning=FALSE, message=FALSE---------------------------------------------
require("ggplot2")
ggplot(data = df.biv, aes(x=x1, y=x2))+
geom_density_2d()+
theme_classic()
## -----------------------------------------------------------------------------
df.biv.pred <- data.frame(df.biv,
leaf=factor(getLeafs(tree=tree.biv, data = df.biv)))
ggplot(data = df.biv.pred, aes(x=x1, y=x2))+
geom_density_2d(aes(colour=leaf))+
viridis::scale_color_viridis(discrete=TRUE)+
theme_classic()
## ----bv2, warning=FALSE, message=FALSE----------------------------------------
tree.biv2 <- semtree(model.biv, df.biv, constraints=
semtree.constraints(focus.parameters = "mu1"))
## ----plbv2--------------------------------------------------------------------
plot(tree.biv2)
## ----message=FALSE, warning=FALSE---------------------------------------------
tree.biv3 <- semtree(model.biv, df.biv, constraints=
semtree.constraints(focus.parameters = "mu2"))
## -----------------------------------------------------------------------------
plot(tree.biv3)
## -----------------------------------------------------------------------------
tree.biv4 <- semtree(model.biv, df.biv, constraints=
semtree.constraints(focus.parameters = "VAR_x2"))
plot(tree.biv4)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.