library(OpenMx)

# cohort
cohort <- sample(0:5, size=1000, replace=TRUE)

# x 
x <- rnorm(1000)

# y 
y <- .1 * cohort * x + .2 * x + rnorm(1000, 0, .3)

# try the regression
summary(lm(y~x*cohort-cohort))

# let's make openmx do this
defData <- data.frame(x, y, cohort)

# let's make a model
defModel <- mxModel("for Pablo",
  # data and RAM
  mxData(defData, type="raw"),
  type="RAM",
  manifestVars=c("x", "y"),
  # simple reg paths
  mxPath("x", arrows=2, values=1, labels="varX"),
  mxPath("y", arrows=2, values=1, labels="varYresid"),
  mxPath("one", c("x", "y"), labels=c("meanX", "intercept")),
  # matrices
  mxMatrix("Full", 1, 2, values=.2, free=TRUE, labels=c("beatX", "betaInt"), name="B"),
  mxMatrix("Full", 2, 1, values=1, labels=c("unity", "data.cohort"), name="D"),
  # algebra
  mxAlgebra(B%*%D, name="regResult"),
  # last path
  mxPath("x", "y", free=FALSE, labels="regResult[1,1]")
  )
# run it
defRes <- mxRun(defModel)


# try the regression
summary(lm(y~x*cohort-cohort))
summary(defRes)