## Generate some fake data
require(mvtnorm)
set.seed(12345)
N <- 2000
Df <- cbind(1, rmvnorm(N, mean=rep(0, 5)))
## Set coefficients
beta.a <- c(-2.00, 0.33, 0.66, 1.00)
beta.b <- c(0.00, 1.50, -0.25)
## Generate path probabilities following a normal model.
y.a <- as.vector(pnorm(tcrossprod(beta.a, Df[, 1:4])))
y.b <- as.vector(pnorm(tcrossprod(beta.b, Df[, c(1, 5, 6)])))
## AND and OR-model
or <- function(x, y) { x + y - x * y }
and <- function(x, y) { x * y }
y.star.OR <- or(y.a, y.b)
y.star.AND <- and(y.a, y.b)
## Observed responses
y.OR <- rbinom(N, 1, y.star.OR)
y.AND <- rbinom(N, 1, y.star.AND)
## Set up data.frame for estimation
Df <- cbind(1, Df)
Df <- as.data.frame(Df)
Df[,1] <- y.OR
Df[,2] <- y.AND
names(Df) <- c("y.OR", "y.AND", "x1", "x2", "x3", "x4", "x5")
## Before estimating, boolean models need to be specified using the
## boolprep function.
## OR model, specified to use a probit link for each causal path. This
## matches the data generating process above.
mod.OR <- boolprep(y.OR ~ (a | b), a ~ x1 + x2 + x3, b ~ x4 + x5,
data = Df, family=list(binomial("probit")))
## Fit a model using the nlminb optimizer (the default). Verbose output is
## requested by specifying trace=1 as a control parameter.
(fit <- boolean(mod.OR, method="nlminb", control=list(trace=1)))
## Multiple optimizers can be specified in a single call to boolean. Here
## we fit with the nlm and nlminb optimizers.
(fit1 <- boolean(mod.OR, method=c("nlm", "nlminb")))
## The summary function will report the detailed results for each
## optimization method.
summary(fit1)
Run the code above in your browser using DataLab