Pooled results work differently than in Zelig 3.5. This needs both a tweak and a documentation change
library(Zelig)
obs <- 1000
dat <- data.frame(cbind(rnorm(obs), runif(obs)))
names(dat) <- c("myX1", "myX2")
dat$myY <- dat$myX1 + 2*dat$myX2 - dat$myX1*dat$myX2 + rnorm(obs)
m1 <- zelig(myY ~ myX1 + myX2 + myX1:myX2, data = dat, model = "ls")
x.out <- setx(m1, myX1 = 2:3, myX2 = .4)
x.out2 <- setx(m1, myX1 = 2:3, myX2 = .6)
s.out <- sim(m1, x = x.out, x1 = x.out2)
s.out$qi$fd