.packageName <- "Zelig"
cmvglm<-function(formula, model, equal, zeros, ones, ndim){ 
  tt<-terms(formula)
  vars<-c("(Intercept)", attr(tt, "term.labels"))
  cm<-vector("list", length(vars))
  names(cm)<-vars
    for(i in 1:length(cm))
      cm[[i]]<-diag(1, ndim)
  if(!is.null(equal)){
    cnames<-names(equal)
    for(i in 1:length(equal)){
      ci<-sort(unique(equal[[i]]))
      indx<-pmatch(cnames[i], names(cm))
      for(j in 1:length(ci)){
        cm[[indx]][ci[j],ci[j]]<-0
        cm[[indx]][ci[j],ci[1]]<-1
      }
    }
  }
  if(!is.null(zeros)){
    cnames<-names(zeros)
    for(i in 1:length(zeros)){
      ci<-sort(unique(zeros[[i]]))
      indx<-pmatch(cnames[i], names(cm))
      for(j in 1:length(ci))
        cm[[indx]][ci[j],ci[j]]<-0
    }
  }
  if(!is.null(ones))
    for(i in 1:length(ones))
      for(j in 1:length(cm))
        if(names(cm)[j]!="(Intercept)")
          cm[[j]][ones[i],]<-matrix(0, ncol=ncol(cm[[j]]))
  for(i in 1:length(cm))
    cm[[i]]<-as.matrix(cm[[i]][,apply(cm[[i]], 2, sum)!=0])
  list("formula"=formula, "constraints"=cm)
}
gsource <- function(var.names= NULL, variables) {
  if (is.null(var.names)) {
    cat(variables, file = ".Rtmp.dat", sep = "\n")
    out <- read.table(".Rtmp.dat")
  }
  else {
    cat(var.names, variables, file = ".Rtmp.dat", sep = "\n")
    out <- read.table(".Rtmp.dat", header = TRUE)
  }
  unlink(".Rtmp.dat")
  return(as.data.frame(out))
}



help.zelig <- function(object) {
  # Elements of this function use functions from the Hmisc library for
  # R by Frank E Harrell Jr, distributed under the GNU GPL v.2.
  under.unix   <- !(version$os=='Microsoft Windows' ||
                    version$os=='Win32' || version$os=='mingw32')
  sys <- function (command, text = NULL) {
    cmd <- if (length(text)) 
      paste(command, text)
    else command
    if (under.unix) 
      system(cmd)
    else shell(cmd, wait = TRUE)
  }
  browser <- .Options$help.browser
  if(!length(browser)) browser <- .Options$browser
  if(!length(browser)) browser <- getOption("browser")
  url <- NULL
  if (missing(object)) 
    url <- c("http://gking.harvard.edu/zelig/docs/")
  else if (is.character(object)) { 
    z <- .libPaths()
    n <- urls <- array()
    for (i in 1:length(z))
      n[i] <- file.path(z[i], "Zelig", fsep = .Platform$file.sep)
    check <- file.exists(n)
    if (sum(check) > 1) {
      Zdir <- n[check][1]
      warning(paste("library Zelig found in two locations.  Using", 
Zdir))
    }
    else
      Zdir <- n[which(check)]
    Zdata <- file.path(Zdir, "data", fsep = .Platform$file.sep)
    files <- list.files(Zdata, pattern = "url")
    for (i in 1:length(files))
      urls[i] <- file.path(Zdata, files[i], fsep = .Platform$file.sep)
    data.path <- read.table(urls[1], header = FALSE)
    if (length(urls) > 1) {
      for (i in 2:length(urls)) { 
        tmp <- read.table(urls[i], header = FALSE)
        data.path <- rbind(data.path, tmp)
      }
    }
    url <- data.path[which(as.character(data.path[,1]) == object), 2]
  }
  else 
    stop("Please enclose the requested topic \n in quotes and try again.")
  if (is.null(url)) {
    cat("Warning: Requested topic not found in Zelig help.  \n If you are sure the topic exists, please check \n the full documentation at http://gking.harvard.edu/zelig.  \n Now searching R-help.\n\n")
    topic <- as.name(object)
    do.call("help", list(topic, htmlhelp = TRUE))
  }
  else {
    if (under.unix) {
      sys(paste(browser, as.character(url), '&'))
      invisible()
    }
    if (!under.unix) {
      browseURL(as.character(url), browser = browser)
      invisible("")
    }
  }
}  


".First.lib" <- function(lib, pkg) {
  cat("\nPlease refer to http://gking.harvard.edu/zelig for full documentation \n",
      "or help.zelig() for help with commands and models supported by Zelig.\n\n",
      sep='')
  
  if(!any(search()=="package:MASS"))
    require(MASS) 
  if(!any(search()=="package:boot"))
    require(boot) 
  
  options(digits = 4)
}

user.prompt <- function() 
  silent <- readline("\nPress <return> to continue: ")
"$.vglm"<-"$.summary.vglm"<-function(a, b){
  if(is.na(pmatch(b, slotNames(a)))){
    pmt<-pmatch(b, names.default(a))
    if(is.na(pmt))
      stop(paste("no such object exists in", a))
    else
      a[[pmt]]
  }
  else
    slot(a, b)
}

"$<-.vglm"<-"$<-summary.vglm"<-function(a, b, value){
  if(is.na(pmatch(b, slotNames(a)))){
    tmp<-class(a)
    class(a)<-NULL
    if(length(a)==0){
      a[[1]]<-value
      names(a)[1]<-b
    }
    else{
      pmt<-pmatch(b, names(a))
      if(is.na(pmt)){
        n<-length(names(a))
        a[[n+1]]<-value
        names(a)[n+1]<-b
      }
      else
        a[[pmt]]<-value
    }
    class(a)<-tmp
  }
  else
    slot(a, b)<-value
  a
}
model.end <- function(e, res) {
  variance <- -solve(res$hessian)
  lev <- eval(as.name("lev"), envir = e)
  Terms <- eval(as.name("Terms"), envir = e)
  call <- eval(as.name("m"), envir = e)   
  xlev <- eval(as.name("xlev"), envir = e)
  mf <- eval(as.name("mf"), envir = e)
  names(res$par) <- c("(Intercept)", attr(Terms, "term.labels"))
  colnames(variance) <- rownames(variance) <- names(res$par)
  fit <- list(coefficients = res$par, variance = variance, 
              lev = lev, terms = Terms, call = call, 
              convergence = res$convergence, xlevels = xlev) 
  attr(fit, "na.message") <- attr(mf, "na.message") 
  if (!is.null(attr(mf, "na.action"))) 
    fit$na.action <- attr(mf, "na.action") 
  fit
}
model.start <- function(e, m) {
  start <- m$start
  m$start <- m$... <- NULL
  m[[1]] <- as.name("model.frame")
  mf <- eval.parent(m)
  assign("mf", mf, envir = e)
  Terms <- attr(mf, "terms")
  assign("Terms", Terms, envir = e)
  x <- model.matrix(Terms, mf)
  assign("x", x, envir = e)
  yvar <- attr(Terms, "response")
  xvars <- as.character(attr(Terms, "response"))[-1]
  xlev <- if (length(xvars) > 0) {
    xlev <- lapply(mf[xvars], levels)
    xlev[!sapply(xlev, is.null)]
  }
  assign("xlev", xlev, envir = e)
  y <- model.response(mf)
  assign("y", y, envir = e)
  lev <- levels(y)
  assign("lev", lev, envir = e)
  if (is.null(start))
    start <- lm(y ~ x -1)$coef
  assign("start", start, envir = e)
}
names.relogit<-function(object,...){
  res<-list(default=names.default(object, ...),
            estimate = names.default(object$lower.estimate, ...))
  class(res)<-"names.relogit"
  res
}
names.summary.vglm <- function(object, ...){
  res <- list(default = slotNames(object, ...))
  class(res) <- "names.summary.vglm"
  res[[1]]
}
names.summary.zelig <- function(object,...) {
  res <- names.default(object, ...)[1:4]
  qi.names <- names.default(object$qi.stats, ...)
  for (i in 1:length(qi.names)) {
    qi.names[i] <- paste("qi.stats$", qi.names[i], sep = "")
  }
  res <- c(res, qi.names)
  res <- list(default = res)
  class(res)<-"names.zelig"
  res
}
names.summary.zelig.relogit <- function(object,...) {
  res <- names.default(object, ...)[1:5]
  qi.names <- names.default(object$qi.stats, ...)
  for (i in 1:length(qi.names)) {
    qi.names[i] <- paste("qi.stats$", qi.names[i], sep = "")
  }
  res <- c(res, qi.names)
  res <- list(default = res)
  class(res)<-"names.zelig"
  res
}
names.vglm <- function(object, ...)
  c(slotNames(object, ...), names.default(object, ...))
names.zelig <- function (object, ...) {
   res <- names.default(object, ...)[1:5]
   qi.names <- names.default(object$qi, ...)
   for (i in 1:length(qi.names)) {
     qi.names[i] <- paste("qi$", qi.names[i], sep = "")
  }
  res <- c(res, qi.names)
  res <- list(default = res)
  class(res) <- "names.zelig"
  res
}
names.zelig.relogit <- function (object, ...) {
   res <- names.default(object, ...)[1:5]
   qi.names <- names.default(object$qi, ...)
   for (i in 1:length(qi.names)) {
     qi.names[i] <- paste("qi$", qi.names[i], sep = "")
  }
  res <- c(res, qi.names)
  res <- list(default = res)
  class(res) <- "names.zelig"
  res
}
param<-function(object, num, ...)
  UseMethod("param")
param.default<-function(object, num, bootstrap = FALSE) {
  if (!bootstrap)
    res <- mvrnorm(num, mu=coef(object), Sigma=vcov(object))
  else
    res <- coef(object)
  res
}
param.glm <- function(object, num = NULL, bootstrap = FALSE) {
  if (!bootstrap) {
    coef <- mvrnorm(num, mu=coef(object), Sigma=vcov(object))
    if (object$zelig == "normal") {
      df <- object$df.residual
      sig2 <- summary(object)$dispersion
      alpha <- sqrt(df*sig2/rchisq(num, df=df))
      res <- cbind(coef, alpha)
    }
    else if (object$zelig == "gamma")  {
      rate <- gamma.shape(object) 
      alpha <- rnorm(num, mean = rate$alpha, sd = rate$SE)
      res <- cbind(coef, alpha)
    }
    else if (object$zelig == "negbin") {
      alpha <- object$theta
      res <- cbind(coef, c(alpha))
    }
    else
      res <- coef
  }
  else {
    coef <- coef(object)
    if (object$family$family == "gaussian") {
      alpha <- sum(object$residuals^2)/length(object$residuals)
      res <- c(coef, alpha)
    }
    else if (object$family$family == "Gamma") {
      alpha <- gamma.dispersion(object)
      res <- c(coef, alpha)
    }
    else if (object$family$family == "neg.bin") {
      alpha <- object$theta
      res <- c(coef, alpha)
    }
    else
      res <- coef
  }
  res
}




param.lm <-function(object, num, bootstrap = FALSE) {
  if (!bootstrap) {
    coef <- mvrnorm(num, mu=coef(object), Sigma=vcov.lm(object))
    df <- object$df.residual
    sig2 <- summary(object)$sigma^2
    alpha <- sqrt(df*sig2/rchisq(num, df=df))
    res <- cbind(coef, alpha)
  }
  else {
    coef <- coef(object)
    alpha <- summary(object)$sigma
    res <- c(coef, alpha)
  }
  res
}
param.multinom <- function(object, num, bootstrap = FALSE) {
  k <- length(object$lev)
  coef <- NULL
  tmp <- coef(object)
  for (i in 1:(k-1))
    coef <- c(coef, tmp[i,])
  if (!bootstrap) 
    sim.coef <- mvrnorm(num, mu=coef, Sigma=vcov(object))
  else
    sim.coef <- coef
  sim.coef
}
param.polr <- function(object, num, bootstrap=FALSE) {
  coef <- object$coefficients
  zeta <- object$zeta
  k <- length(coef)
  if (!bootstrap) {
    par <- mvrnorm(num, mu=c(coef,zeta), Sigma=vcov(object))
    coef <- par[,1:k]
    alpha <- par[,(k+1):ncol(par)]
    res <- as.matrix(cbind(coef, alpha))
  }
  else
    res <- c(coef, zeta)
  res
}
param.relogit <- function(object, num, x, bootstrap = FALSE, bootfn = NULL, ...) {
  tau <- eval(object$call$tau, sys.parent())
  if (length(tau) == 2) { # without full population information
    pping <- function(object, x, x1, num, bootstrap, bootfn,...) {
      if(!bootstrap) {
        par0 <- param.default(object$lower.estimate, num=num, bootstrap=bootstrap)
        par1 <- param.default(object$upper.estimate, num=num, bootstrap=bootstrap)
      }
      else {
        dta <- eval(object$data, sys.parent())
        dta <- dta[complete.cases(model.frame(object, dta)),]
        if (is.null(bootfn)) {
          bootfn <- function(data, i, obj) {
            d <- data[i,]
            obj$call$data <- d
            fit <- eval(obj$call, sys.parent())
            return(fit$coefficients)
          }
        }
        res0 <- boot(dta, bootfn, R=num, obj = object$lower.estimate, ...)
        res1 <- boot(dta, bootfn, R=num, obj = object$upper.estimate, ...)
        colnames(res0$t) <- names(res0$t0)
        colnames(res1$t) <- names(res1$t0)
        par0 <- res0$t
        par1 <- res1$t
      }
      sim00 <- qi.glm(object$lower.estimate, par0, x = x)
      P00 <- as.matrix(sim00$qi$ev)
      sim10 <- qi.glm(object$upper.estimate, par1, x = x)
      P10 <- as.matrix(sim10$qi$ev)
      test <- P00[,1] < P10[,1]
      par0 <- as.matrix(par0[test,])
      par1 <- as.matrix(par1[test,])
      list(par0 = par0, par1 = par1)
    }
    tmp <- pping(object, bootstrap=bootstrap, bootfn=bootfn, x=x, x1=x1, num=num, ...)
    par0 <- tmp$par0
    par1 <- tmp$par1
    while (nrow(par0) < num) {
      tmp <- pping(object, bootstrap=bootstrap, bootfn=bootfn, x=x, x1=x1, num=num,...)
      par0 <- rbind(par0, tmp$par0)
      par1 <- rbind(par1, tmp$par1)
    }
    if (nrow(par0) > num) {
      par0 <- par0[1:num,]
      par1 <- par1[1:num,]
    }
    par0 <- as.matrix(par0)
    par1 <- as.matrix(par1)
    rownames(par0) <- 1:nrow(par0)
    rownames(par1) <- 1:nrow(par1)
    return(list(par0 = par0, par1 = par1))
  }
  else { # with precise population info, or no population info
    if (!bootstrap) 
      simpar <- param.default(object, num=num, bootstrap = bootstrap)
    else {
      tt <- terms(object)
      dta <- eval(object$data, sys.parent())
      dta <- dta[complete.cases(model.frame(tt, dta)),]
      if(is.null(bootfn)) {
        bootfn <- function(data, i, object) {
          d <- data[i,]
          object$call$data <- d
          fit <- eval(object$call, sys.parent())
          return(param(fit, bootstrap = TRUE))
        }
      }
      res <- boot(dta, bootfn, R = num, object = object, ...)
      colnames(res$t) <- names(res$t0)
      simpar <- res$t
    }        
    return(simpar)
  }
}
























































param.survreg <- function(object, num, bootstrap = FALSE) {
  cov <- vcov(object)
  coef <- object$coefficients
  log.scale <- log(object$scale)
  k <- length(coef)
  if(!bootstrap) {
    if(ncol(cov)==k)
      res <- mvrnorm(num, mu=coef, Sigma=cov)
    else 
      res <- mvrnorm(num, mu=c(coef, log.scale), Sigma=cov)
  }
  else {
    if (ncol(cov) == k)
      res <- c(coef)
    else
      res <- c(coef, log.scale)
  }
  res
}
param.vglm <- function(object, num, bootstrap = FALSE) {
  cov <- vcov(object)
  res <- object@coefficients
  if (!bootstrap) 
    res <- mvrnorm(num, mu=res, Sigma=cov)
  res
}
plot.ci <- function(object, CI=95, qi = "ev", main = "",
                    ylab = NULL, xlab = NULL, xlim = NULL,
                    ylim = NULL, col = c("red", "blue"), ...) {
  if (class(object) != "zelig")
    stop(" plot.ci() works only for sim() output.")
  if (object$zelig.call$model %in% c("bprobit", "blogit", "ologit", "mlogit", "oprobit", "mloglm"))
    stop("\n    plot.ci() is valid only for non-categorical, univariate dependent variable models.")
  var <- var1 <- list()
  idx <- array()
  cip <- c((100-CI)/200, 1-(100-CI)/200)
  summarize <- function(z, cip){
    res <- NULL
    res <- cbind(res, apply(z, 2, quantile, prob=cip[1]))
    res <- cbind(res, apply(z, 2, quantile, prob=cip[2]))
    res
  }
  for (i in 1:ncol(object$x)) {
    var[[i]] <- unique(object$x[,i])
    if (length(var[[i]]) > 1)
      idx[i] <- i
    else
      idx[i] <- NA
  }
  if (!is.null(object$x1)) {
    idx1 <- NULL
    for (i in 1:ncol(object$x1)) {
      var1[[i]] <- unique(object$x1[,i])
      if (length(var1[[i]]) > 1) 
        idx1[i] <- i
      else
        idx1[i] <- NA
    }
    idx1 <- na.omit(idx1)
  }
  idx <- na.omit(idx)
  if (!is.null(idx1))
    if (!identical(idx, idx1))
      stop(" x and x1 in sim() object vary on different dimensions.")
  var <- var[[idx[1]]]
  q <- pmatch(qi, names.default(object$qi))
  qofi <- object$qi[[q]]
  sum.qi <- summarize(qofi, cip)
  if (!is.null(object$x1) && qi == "ev") {
    fd <- object$qi$fd
    ev1 <- fd + qofi
    sum.qi1 <- summarize(ev1, cip)
  }
  else
    sum.qi1 <- NULL
  if (is.null(ylab))
    ylab <- object$qi.name[[q]]
  if (is.null(xlab))
    xlab <- paste("Range of", colnames(object$x)[idx[1]])
  if (is.null(ylim)) {
    if (is.null(sum.qi1))
      ylim <- c(min(sum.qi), max(sum.qi))
    else
      ylim <- c(min(sum.qi, sum.qi1), max(sum.qi, sum.qi1))
  }
  if (is.null(xlim))
    xlim <- c(min(var), max(var))
  plot.default(var, type = "n", ylab = ylab, main = main, xlab = xlab, 
               xlim = xlim, ylim = ylim)
  for (i in 1:length(var)) {
    lines(c(var[i], var[i]), c(sum.qi[i,1], sum.qi[i,2]), col = col[1], ...)
    if (!is.null(object$x1) && qi == "ev")
      lines(c(var[i], var[i]), c(sum.qi1[i,1], sum.qi1[i,2]), col = col[2], ...)
  }
}










plot.zelig <- function(object, alt.col = "red", xlab = "", user.par = FALSE, ...) {
  model <- object$zelig.call$model
  object$qi$theta <- NULL
  k <- length(object$qi)
  op <- par(no.readonly = TRUE)
  if (!user.par) 
    par(mar = c(4,4,2,1), tcl = -0.25, mgp = c(2, 0.6, 0))
  if (dim(object$x)[1] > 1)
    stop("\n     default plot methods only supported for sim()\n     results from one-observation values of x.") 
  if (model %in% c("ls", "normal", "exp", "lognorm","weibull", "gamma")) {
    par(mfrow = c(k, 1))
    for (i in 1:k) {
      qi <- as.vector(object$qi[[i]])
      plot(density(qi), main = object$qi.name[[i]], xlab = xlab, ...)
    }
  }
  tau1 <- model == "relogit" && dim(object$qi[[1]])[2] == 1 
  if (model %in% c("logit", "probit") || tau1) {
    par(mfrow = c(k, 1))
    pr <- object$qi$pr
    y0 <- sum(pr == 0)
    y1 <- sum(pr == 1)
    xmax <- max(y0, y1)
    barplot(c(y0, y1), horiz = TRUE, col = alt.col, las = 1,
            names.arg = c("Y = 0", "Y = 1"), xlab = xlab,
            main = object$qi.name$pr, xlim = c(0, 1.15*xmax), ...)
    object$qi$pr <- object$qi.name$pr <- NULL
    for (i in 1:(k-1)) {
      qi <- as.vector(object$qi[[i]])
      plot(density(qi), main = object$qi.name[[i]], xlab = xlab, ...)
    }    
  }
  if (model == "relogit" && dim(object$qi[[1]])[2] == 2) {
    par(mfrow = c(k,1))
    for (i in 1:k) {
      qi <- object$qi[[i]]
      main <- as.character(object$qi.name[i])
      if (is.null(rownames(qi)))
        rownames(qi) <- 1:dim(qi)[1]
      idx <- as.integer(sample(rownames(qi), 100))
      tmp <- qi[idx,,1]
      xlim <- c(min(qi[,1,1]), max(qi[,2,1]))
      plot(xlim, type = "n", xlab = xlab, ylab = "Observations (n = 100)",
           main = main, ylim = c(0, 100), xlim = xlim, ...)
      for (j in 1:nrow(tmp))
        lines(c(tmp[j,1], tmp[j,2]), c(j,j), col = alt.col)
      abline(v = mean(qi[,1,1]))
      abline(v = mean(qi[,2,1]))
    }
  }
  if (model %in% c("blogit", "bprobit")) {
    if ("rr" %in% names(object$qi)) {
      object$qi$rr <- object$qi.name$rr <- NULL
      k <- k - 1
    }
    main <- object$qi.name
    par(mfrow = c(k, 1))
    for (i in 1:k) {
      qi <- as.matrix(object$qi[[i]])
      if (is.character(qi[1])) {
        y00 <- sum(as.integer(qi[,1]))
        y01 <- sum(as.integer(qi[,2]))
        y10 <- sum(as.integer(qi[,3]))
        y11 <- sum(as.integer(qi[,4]))
        xmax <- max(y00, y01, y10, y11)
        labels <- c("(0,0)", "(0,1)","(1,0)", "(1,1)")
        barplot(c(y00, y01, y10, y11), horiz = TRUE, col = alt.col,
                names.arg = labels, xpd = TRUE, main = main[[i]],
                las = 1, xlim = c(0, 1.15*xmax), ...)
      }
      else if (is.numeric(qi[1])) {
        y1 <- qi[, 3] + qi[, 4]
        y2 <- qi[, 2] + qi[, 4]
        contour(kde2d(y1, y2), xlab = "Pr(Y1 = 1)", 
                ylab = "Pr(Y2 = 1)", main = main[[i]], ...)
      }
    }
  }
  if(model %in% c("negbin", "poisson")) {
    par(mfrow = c(k,1))
    truehist(object$qi$pr, main = object$qi.name$pr, x0 = -0.25,
             xlab = xlab, ylab = "Probability", ...)
    plot(density(object$qi$ev), main = object$qi.name$ev, xlab = xlab, ...) 
    if (k == 3)
      plot(density(object$qi$fd), main = object$qi$fd, xlab = xlab, ...)
  }
  if (model %in% c("oprobit", "ologit", "mlogit")) {
    if (!is.null(object$qi$rr)) {
      k <- k - 1
      object$qi$rr <- object$qi.name$rr <- NULL
    }
    par(mfrow = c(k,1))
    pr <- object$qi$pr
    spr <- array()
    lev <- sort(unique(pr))
    K <- length(lev)
    for (i in 1:K)
      spr[i] <- sum(pr == lev[i])
    xmax <- max(spr)
    labels <- paste("Y=", lev, sep = "")
    barplot(spr, horiz = TRUE, col = rainbow(K), names.arg = labels,
            las = 1, main = object$qi.name$pr, xlim = c(0, 1.15*xmax), ...)
    object$qi$pr <- object$qi.name$pr <- NULL
    main <- object$qi.name
    for (i in 1:(k-1)) {
      qi <- object$qi[[i]]
      if (length(dim(qi)) == 3 && dim(qi)[3] == 1)
        qi <- qi[,,1]
      dens <- list()
      xmax <- ymax <- array()
      for (j in 1:ncol(qi)) {
        dens[[j]] <- density(qi[,j])
        xmax[j] <- max(dens[[j]]$x)
        ymax[j] <- max(dens[[j]]$y)
      }
      plot(dens[[1]], col = rainbow(K)[1],
           xlim = c(min(min(qi), 0), max(xmax)),
           xlab = xlab, main = "",
           ylim = c(0, max(ymax)), ...)
      for (j in 2:ncol(qi)) 
        lines(dens[[j]], col = rainbow(K)[j])
      title(main = object$qi.name[[i]][1])
    }
  }
  par(op)
}

















print.names.relogit<-function(obj, ...){
  if (length(obj$tau) == 2) {
    print(obj$default, ...)
    cat(paste("Additional objects available in lower.estimate and upper.estimate: \n"))
    print(obj$estimate)
  }
  else
    print(obj$default, ...)
}

print.names.zelig<-function(obj, ...){
  print(obj$default, ...)
}
print.relogit <- function(object, digits = getOption("digits"), ...) {
  if (names(object)$default[1] == "lower.estimate") {
    cat(paste("\nEstimate for tau =",
              round(object$lower.estimate$call$tau, digits = 5), "\n"))
    print.glm(object$lower.estimate, digits = digits, ...)
    cat(paste("\n\nEstimate for tau =",
              round(object$upper.estimate$call$tau, digits = 5), "\n"))    
    print.glm(object$upper.estimate, digits = digits, ...)
  }
  else
    print.glm(object, digits  = digits, ...)
}
print.setx<-function(object, digits=getOption("digits"),
                     ...){
  attributes(object)$class<-NULL
  attributes(object)$assign<-NULL
  attributes(object)$contrasts<-NULL
  print.matrix(object, digits=digits, ...)
}
print.summary.MI <- function(object, subset = NULL, ...){
  m <- length(object$all)
  if (m == 0)
    m <- 1
  if (any(subset > max(m)))
    stop("the subset selected lies outside the range of available \n        observations in the MI regression output.")
  cat("\n  Model:", object$zelig)
  cat("\n  Number of multiply imputed data sets:", m, "\n")
  if (is.null(subset)) {
    cat("\nCombined results:\n\n")
    cat("Call:\n")
    print(object$call)
    cat("\nCoefficients:\n")
    print(object$coefficients)
    cat("\nFor combined results from datasets i to j, use summary(object, subset = i:j).\nFor separate results, use print(summary(object), subset = i:j).\n\n")
  }
  else {
    if (is.function(subset))
      M <- 1:m
    if (is.numeric(subset))
      M <- subset
    for(i in M){
      cat(paste("\nResult with dataset", i, "\n"))
      print(object$all[[i]], ...)
    }
  }
}

print.summary.multiple <- function(object, subset = NULL, ...){
  if (is.null(subset))
    m <- length(object$M)
  else if (any(subset > max(m)))
    stop("the subset selected lies outside the range of available \n        sets of regression output.")
  else
    m <- subset
  cat("\n  Model:", object$zelig)
  cat("\n  Number of subsets evaluated:", m, "\n")
  for (i in 1:m) {
    cat(paste("\nResults for", object$call$by[1], "=", object$lev[i], "\n"))
    print(object[[i]])
    cat("\n")
  }
}

print.summary.relogit <- function(object, digits = getOption("digits")) {
  if (length(object) == 2)
    call <- object[[1]]$call
  else
    call <- object$call
  tau <- eval(call$tau, sys.parent())
  if (length(tau) == 2) {
    printing <- function(obj, tau, digits) {
      print.glm(obj, digits = digits)
      cat(paste("Prior correction is made with tau = ",
                round(tau, digits = max(digits, 5)),"\n"))
      if (call$bias.correct)
        cat("Rare events bias correction was made.")
      cat("\n\n")
    }
    printing(object$lower.summary, tau[1], digits)
    printing(object$upper.summary, tau[2], digits)
  }
  if (length(tau) == 1 || length(tau) == 0){
    print.glm(object, digits)
    if(!is.null(tau))
      cat(paste("Prior correction is made with tau = ",
                round(tau, digits = max(digits, 5)),"\n"))
    if(call$bias.correct)
      cat("Rare events bias correction was made.\n\n")
  }
}
print.summary.zelig<-function(object, digits=getOption("digits"),
                              print.x=FALSE, ...){
  cat("\n  Model:", object$zelig, "\n")
  if (!is.null(object$num))
      cat("  Number of simulations:", object$num, "\n")
  if (!is.null(object$x)) {
    if (print.x || nrow(object$x) == 1 || is.null(dim(object$x))) {
      if (class(object$x) == "setx.cond")
        cat("\nObserved Data \n")
      else
        cat("\nValues of X \n")
      print.setx(object$x, digits=digits, ...)
      if(!is.null(object$x1)){
        cat("\nValues of X1 \n")
        print.setx(object$x1, digits=digits, ...)
      }
    }
    else {
      if (class(object$x) == "setx.cond")
        cat("\nMean Values of Observed Data (n = ", nrow(object$x), ") \n", sep = "")
      else
        cat("\nMean Values of X (n = ", nrow(object$x), ") \n", sep = "")
      print.setx(apply(object$x, 2, mean), digits=digits, ...)
      if (!is.null(object$x1)) {
        cat("\nMean Values of X1 (n = ", nrow(object$x1), ") \n", sep = "")
        print.setx(apply(object$x1, 2, mean), digits=digits, ...) 
      }
    }
  }
  for (i in 1:length(object$qi.name)){
    indx <- pmatch(names(object$qi.name[i]), names(object$qi.stats))
    tmp <- object$qi.stats[[indx]]
    if (names(object$qi.name)[indx] == "pr" && colnames(tmp)[1] != "mean")
      lab <- paste(object$qi.name[[i]], "(percentage of simulations)", sep = " ")
    else
      lab <- object$qi.name[[i]]
    cat("\n", lab, "\n", sep = "")
    if (length(dim(tmp)) == 3) {
        for (j in 1:dim(tmp)[3]){
          cat("\n  Observation", dimnames(tmp)[[3]][j], "\n")
          if (is.null(rownames(tmp[,,j])))
            rownames(tmp[,,j]) <- 1:nrow(tmp[,,j])
          if (!is.null(names(tmp[,,j])))
            names(tmp[,,j]) <- NULL
          print.matrix(tmp[,,j], digits=digits, ...)
        }
      }
    else {
      if (is.null(rownames(tmp)))
        rownames(tmp) <- 1:nrow(tmp)
      if (!is.null(names(tmp)))
        names(tmp) <- NULL
      print.matrix(tmp, digits=digits, ...)
    }
  }
}












print.summary.zelig.multiple <- function(object, subset = NULL, ...){
  if (is.null(subset))
    m <- length(object)
  else if (any(subset > max(m)))
    stop("the subset selected lies outside the range of available \n        sets of regression output.")
  else
    m <- subset
  for (i in 1:m) {
    cat(paste("\nResults for", names(object)[i], "\n"))
    print(object[[i]])
    cat("\n")
  }
}

print.summary.zmatchit <- function(object, digits=getOption("digits"), ...){
  M<-length(object)
  for(i in 1:(M/2)){
    if(M>2)
      cat("\nSubclass", i, "\n")
    cat("\nTreatment Group:\n")
    print(object[[(i-1)*2+1]], ...)
    cat("\nControl Group:\n")
    print(object[[i*2]], ...)
  }
}
print.zelig <- function (object, digits = max(3, getOption("digits") - 3), ...) 
{
    cat("\nModel:", object$zelig.call$model, "\n", sep = " ")
    cat("Number of simulations:", object$call$num, "\n\n", sep = " ")
    idx <- unlist(object$qi.name)
    cat("Available Quantities of Interest: \n")
    for (i in 1:length(object$qi.name)){
      cat(paste("  qi$", names(object$qi)[i], " = ", idx[[i]], "\n", sep = ""))
    }
    cat("\nPlease use summary() to obtain more information. \n")
    invisible(object)
}
print.zmatchit <- function(object, digits=getOption("digits"), ...){
  M<-length(object)-1
  for(i in 1:(M/2)){
    if(M>2)
      cat("\nSubclass", i)
    cat("\n\nTreatment Group:\n")
    print(object[[(i-1)*2+1]], ...)
    cat("\nControl Group:\n")
    print(object[[i*2]], ...)
  }
}
qi<-function(object, ...)
  UseMethod("qi")
qi.glm <- function(object, par, x, x1 = NULL, y = NULL) {
  check <- FALSE
  model <- object$zelig
  k <- length(object$coef)
  coef <- par[,1:k]
  if (k < ncol(par)) 
    alpha <- par[,(k+1):ncol(par)]
  eta <- coef %*% t(x)
  theta <- matrix(object$family$linkinv(eta), nrow = nrow(coef))
  pr <- ev <- matrix(NA, nrow = nrow(theta), ncol = ncol(theta))
  dimnames(pr) <- dimnames(ev) <- dimnames(theta)
  if (model %in% c("logit", "probit", "relogit")) {
    check <- TRUE
    ev <- theta
    for (i in 1:ncol(theta)) 
      pr[,i] <- as.character(rbinom(length(ev[,i]), 1, ev[,i]))
  }
  else if (model == "normal") {
    ev <- theta
    for (i in 1:nrow(ev)) 
      pr[i,] <- rnorm(length(ev[i,]), mean = ev[i,], sd = alpha[i])
  }
  else if (model == "gamma") {
    ev <- theta * 1/alpha
    for (i in 1:nrow(ev))  
      pr[i,] <- rgamma(length(ev[i,]), shape = theta[i,], scale = 1/alpha[i])
  }
  else if (model == "poisson") {
    ev <- theta
    for (i in 1:ncol(ev))
      pr[,i] <- rpois(length(ev[,i]), lambda = ev[,i])
  }
  else if (model == "negbin") {  
    ev <- theta
    for (i in 1:nrow(ev)) 
      pr[i,] <- rnegbin(length(ev[i,]), mu = ev[i,], theta = alpha[i])
  }
  qi <- list(ev = ev, pr = pr)
  qi.name <- list(ev = "Expected Values: E(Y|X)",
                  pr = "Predicted Values: Y|X")
  if (!is.null(x1)){
    theta1 <- matrix(object$family$linkinv(coef %*% t(as.matrix(x1))),
                     nrow = nrow(coef))
    if (model == "gamma")
      ev1 <- theta1 * 1/alpha
    else
      ev1 <- theta1
    qi$fd <- ev1-ev
    qi.name$fd <- "First Differences in Expected Values: E(Y|X1)-E(Y|X)"
    if (model %in% c("logit", "probit", "relogit")) {
      qi$rr <- ev1/ev
      qi.name$rr <- "Risk Ratios: P(Y=1|X1)/P(Y=1|X)"
    }
  }
  if (!is.null(y)) {
    yvar <- matrix(rep(y, nrow(par)), nrow = nrow(par), byrow = TRUE)
    tmp.ev <- yvar - qi$ev
    if (check)
      tmp.pr <- yvar - as.integer(qi$pr)
    else
      tmp.pr <- yvar - qi$pr
    qi$ate.ev <- matrix(apply(tmp.ev, 1, mean), nrow = nrow(par))
    qi$ate.pr <- matrix(apply(tmp.pr, 1, mean), nrow = nrow(par))
    qi.name$ate.ev <- "Average Treatment Effect: Y - EV"
    qi.name$ate.pr <- "Average Treatment Effect: Y - PR"
  }
  list(qi=qi, qi.name=qi.name)
}






qi.lm <- function(object, par, x, x1 = NULL, y = NULL) {
  k <- length(object$coef)
  coef <- par[,1:k]
  alpha <- par[,(k+1):ncol(par)]
  ev <- coef %*% t(x)
  pr <- matrix(NA, nrow = nrow(ev), ncol = ncol(ev))
  dimnames(pr) <- dimnames(ev) 
  for (i in 1:nrow(pr)) 
    pr[i,] <- rnorm(length(ev[i,]), mean = ev[i,], sd = alpha[i])
  qi <- list(ev=ev, pr=pr)
  qi.name <- list(ev="Expected Values: E(Y|X)",
                  pr = "Predicted Values: Y|X")
  if(!is.null(x1)){
    ev1 <- coef %*% t(x1)
    qi$fd <- ev1-ev
    qi.name$fd <-
      "First Differences in Expected Values: E(Y|X1)-E(Y|X)"
  }
  if (!is.null(y)) {
    yvar <- matrix(rep(y, nrow(par)), nrow = nrow(par), byrow = TRUE)
    tmp.ev <- yvar - qi$ev
    tmp.pr <- yvar - qi$pr
    qi$ate.ev <- matrix(apply(tmp.ev, 1, mean), nrow = nrow(par))
    qi$ate.pr <- matrix(apply(tmp.pr, 1, mean), nrow = nrow(par))
    qi.name$ate.ev <- "Average Treatment Effect: Y - EV"
    qi.name$ate.pr <- "Average Treatment Effect: Y - PR"
  }
  list(qi=qi, qi.name=qi.name)
}



qi.polr <- function(object, par, x = x, x1 = NULL) {
  m <- ncol(x)
  lev <- object$lev
  k <- length(lev)
  ev.multinom<-function(num, simpar, x, lev, k, m){
    eta <- NULL
    for (i in 1:(k-1))
      eta <- cbind(eta, simpar[,(m*(i-1)+1):(m*i)]%*%t(x))
    eta <- exp(cbind(rep(0, num), eta))
    ev <- eta/apply(eta, 1, sum)
    colnames(ev) <- lev
    ev
  }
  ev <- ev.multinom(num, simpar, x, lev, k, m)
  Ipr <- sim.cut <- matrix(NA, nrow = nrow(ev), ncol = ncol(ev))
  colnames(Ipr) <- colnames(sim.cut) <- lev
  sim.cut[,1] <- ev[,1]
  for (j in 2:k) 
    sim.cut[,j] <- sim.cut[,(j-1)] + ev[,j]
  tmp <- runif(nrow(ev), 0, 1)
  for (l in 1:k) 
    Ipr[,l] <- tmp > sim.cut[,l]
  pr <- NULL
  for (m in 1:nrow(Ipr))
    pr[m] <- 1 + sum(Ipr[m,])
  res$qi <- list(ev=ev, pr=pr)
  res$qi.name <- list(ev = "Predicted Probabilities: P(Y=j|X)",
                      pr = "Predicted Values: Y|X")
  if(!is.null(x1)){
    ev1 <- ev.multinom(num, simpar, x1, lev, k, m)
    res$qi$fd <- ev1-ev
    res$qi.name$fd <- "First Differences: P(Y=j|X1)-P(Y=j|X)"
  }
  list(qi=qi, qi.name=qi.name)
}
qi.polr <- function(object, par, x, x1 = NULL, y = NULL) {
  m <- length(coef(object))
  sim.coef <- par[,1:m]
  sim.zeta <- par[,(m+1):ncol(par)]
  k <- length(object$zeta) + 1
  lev <- object$lev
  eta <- array()
  if (ncol(x) == 1)
    eta <- 0
  else 
    eta <- x[,-1] %*% t(sim.coef)
  ev.polr <- function(ev, eta, sim.zeta, lev) {
    for (j in 1:ncol(sim.zeta)) 
      ev[,j,] <- 1 / (1 + exp(-sim.zeta[,j] + eta))
    for (j in 0:(ncol(sim.zeta)-1))
      for (i in 1:dim(ev)[3])
        ev[,ncol(ev)-j,i] <- ev[,ncol(ev)-j,i]-ev[,(ncol(ev)-j-1),i]
    ev
  }
  ev <- array(1, dim = c(nrow(sim.coef), length(lev), nrow(x)),
              dimnames = list(NULL, lev, rownames(x)))
  n <- nrow(sim.coef)
  z <- nrow(x)
  ev <- ev.polr(ev, eta, sim.zeta, lev)
  Ipr <- array(NA, dim(ev), dimnames(ev))
  sim.cut <- array(NA, dim(ev)[1:2], dimnames(ev)[1:2])
  tmp <- array()
  pr <- matrix(NA, nrow = n, ncol = z)
  pr.polr <- function(ev, Ipr, lev) {
    sim.cut[,1] <- ev[,1]
    for (j in 2:length(lev)) 
      sim.cut[,j] <- sim.cut[,(j-1)] + ev[,j]
    tmp <- runif(dim(ev)[1], 0, 1)
    for (j in 1:length(lev)) 
      Ipr[,j] <- tmp > sim.cut[,j]
    for (l in 1:nrow(Ipr))
      tmp[l] <- 1 + sum(Ipr[l,])
    factor(tmp, levels = sort(unique(tmp)), labels = lev)
  }
  for (i in 1:z)
    pr[,i] <- as.character(pr.polr(ev[,,i], Ipr[,,i], lev))
  colnames(pr) <- rownames(x)
  qi <- list(ev=ev, pr=pr)
  qi.name <- list(ev="Expected Values: P(Y=j|X)",
                  pr="Predicted Values: Y|X")
  if(!is.null(x1)){
    ev1 <- array(1, dim = c(nrow(sim.coef), length(lev), nrow(x)),
                 dimnames = list(NULL, lev, rownames(x)))
    eta1 <- x1[,-1] %*% t(sim.coef)
    ev1 <- ev.polr(ev1, eta1, sim.zeta, lev)
    qi$fd <- ev1-ev
    qi$rr <- ev1/ev
    qi.name$fd <- "First Differences: P(Y=j|X1)-P(Y=j|X)"
    qi.name$rr <- "Risk Ratio: P(Y=j|X1)-P(Y=j|X)"
  }
  if (!is.null(y)) {
    yvar <- matrix(NA, nrow = length(y), ncol = length(lev))
    tmp.ev <- tmp.pr <- array(NA, dim = dim(qi$ev))
    pr.idx <- array(NA, dim = c(nrow(pr), length(lev), nrow(x)))
    qi$ate.ev <- qi$ate.pr <- matrix(NA, dim(qi$ev)[1], dim(qi$ev)[2])
    for (i in 1:length(lev)) {
      yvar[,i] <- as.integer(y == lev[i])
      pr.idx[,i,] <- as.integer(pr[,i] == lev[i])
    }
    colnames(yvar) <- lev 
    for (j in 1:ncol(yvar)) {
      for (i in 1:n) {
        tmp.ev[i,j,] <- yvar[,j] - qi$ev[i,j,]
        tmp.pr[i,j,] <- yvar[,j] - pr.idx[i,j,]
      }
      qi$ate.ev[,j] <- apply(tmp.ev[,j,], 1, mean)
      qi$ate.pr[,j] <- apply(tmp.pr[,j,], 1, mean)
    }
    colnames(qi$ate.ev) <- colnames(qi$ate.pr) <- lev
    qi.name$ate.ev <- "Average Treatment Effect: Y - EV"
    qi.name$ate.pr <- "Average Treatment Effect: Y - PR"
  }
  list(qi=qi, qi.name=qi.name)
}











qi.relogit <- function(object, par, x, x1 = NULL, y = NULL) {
  tau <- eval(object$call$tau, sys.parent())
  num <- nrow(par$par0)
  if (length(tau) == 2){
    low <- qi.glm(object$lower.estimate, par$par0, x, x1)
    up <- qi.glm(object$upper.estimate, par$par1, x, x1)
    PP <- PR <- array(NA, dim = c(num, 2, nrow(x)),
                      dimnames = list(NULL, c("Lower Bound", "Upper Bound"),
                        rownames(x)))
    PP[,1,] <- P00 <- low$qi$ev
    PP[,2,] <- P10 <- up$qi$ev
    qi <- list(ev = PP)
    qi.name <- list(ev = "Expected Values: E(Y|X)")
    if (!is.null(x1)) {
      FD <- RR <- array(NA, dim = c(num, 2, nrow(x)),
                        dimnames = list(NULL, d2 = c("Lower Bound", "Upper Bound"), 
                          rownames(x)))
      sim01 <- qi.glm(object$lower.estimate, par0, x = x1, x1 = NULL)
      sim11 <- qi.glm(object$upper.estimate, par1, x = x1, x1 = NULL)
      tau0 <- object$lower.estimate$call$tau
      tau1 <- object$upper.estimate$call$tau
      P01 <- as.matrix(sim01$qi$ev)
      P11 <- as.matrix(sim11$qi$ev)
      OR <- (P10/(1-P10)) / (P00/(1-P00))
      RR[,1,] <- pmin(as.matrix(P01/P00), as.matrix(P11/P10))
      RR[,2,] <- pmax(as.matrix(P01/P00), as.matrix(P11/P10))
      RD0 <- as.matrix(P01-P00)
      RD1 <- as.matrix(P11-P10)
      RD <- as.matrix((sqrt(OR)-1) / (sqrt(OR)+1))
      ## checking monotonicity
      y.bar <- mean(object$lower.estimate$y)
      beta0.e <- object$lower.estimate$coefficients
      beta1.e <- object$upper.estimate$coefficients
      ## evaluating RD at tau0 and tau1
      RD0.p <- 1/(1+exp(-t(beta0.e) %*% t(x1))) - 1/(1+exp(-t(beta0.e) %*% t(x)))
      RD1.p <- 1/(1+exp(-t(beta1.e) %*% t(x1))) - 1/(1+exp(-t(beta1.e) %*% t(x)))
      ## evaluating RD at tau0+e and tau1+e
      e <- 0.001
      beta0.e["(Intercept)"] <- beta0.e["(Intercept)"]+log(1-tau0)-log(tau0) -
        log(1-tau0-0.001)+log(tau0+0.001)
      beta1.e["(Intercept)"] <- beta1.e["(Intercept)"]+log(1-tau1)-log(tau1) -
        log(1-tau1-e)+log(tau1+e)
      RD0.e <- 1/(1+exp(-t(beta0.e) %*% t(x1))) - 1/(1+exp(-t(beta0.e) %*% t(x)))
      RD1.e <- 1/(1+exp(-t(beta1.e) %*% t(x1))) - 1/(1+exp(-t(beta1.e) %*% t(x)))
      ## checking the sign and computing the bounds
      check <- sum((RD1.e-RD1.p) * (RD0.e-RD0.p))
      if (check > 0) {
        FD[,1,] <- pmin(RD0, RD1)
        FD[,2,] <- pmax(RD0, RD1)
      }
      else {
        FD[,1,] <- pmin(RD0, RD1, RD)
        FD[,2,] <- pmax(RD0, RD1, RD)
      }
      qi$fd <- FD
      qi$rr <- RR
      qi.name$fd <- "First Differences: P(Y=1|X1) - P(Y=1|X)"
      qi.name$rr <- "Risk Ratios: P(Y=1|X1) / P(Y=1|X)"
    }
    if (!is.null(y)) {
      yvar <- matrix(rep(y, num), nrow = num, byrow = TRUE)
      tmp.ev <- yvar - qi$ev
      tmp.pr <- yvar - as.integer(qi$pr)
      qi$ate.ev <- matrix(apply(tmp.ev, 1, mean), nrow = num)
      qi$ate.pr <- matrix(apply(tmp.pr, 1, mean), nrow = num)
      qi.name$ate.ev <- "Average Treatment Effect: Y - EV"
      qi.name$ate.pr <- "Average Treatment Effect: Y - PR"
    }
    out <- list(qi = qi, qi.name = qi.name)
  }
  else 
    out <- qi.glm(object = object, par = par, x = x, x1 = x1, y = y)
  out
}
qi.survreg <- function(object, par, x, x1 = NULL, y = NULL, cond.data = NULL) {
  model <- object$zelig
  k <- length(object$coef)
  sim.coef <- as.matrix(par[,1:k])
  if (model %in% c("weibull", "Weibull", "lognorm")) #{
    sim.scale <- par[,(k+1):ncol(par)]
  if (!is.null(y)) {
    status <- x[,1]
    x <- as.data.frame(x[,2:ncol(x)])
  }
  link <- survreg.distributions[[object$dist]]$itrans
  ev.pr.surv <- function(sim.coef, sim.scale, x, link) {
    eta <- sim.coef %*% t(x)
    theta <- as.matrix(apply(eta, 2, link))
    ev <- pr <- matrix(NA, ncol=ncol(theta), nrow=nrow(theta))
    dimnames(pr) <- dimnames(ev) <- dimnames(theta) 
    if (model == "exp") {
      ev <- theta
      for (i in 1:nrow(pr))
        pr[i,] <- rexp(length(ev[i,]), rate = 1/ev[i,])
    }
    else if (model == "weibull" || model == "Weibull") {
      ev <- theta * gamma(1 + exp(sim.scale))
      for (i in 1:nrow(pr))
        pr[i,] <- rweibull(length(ev[i,]), shape=1/exp(sim.scale[i]),
                           scale=theta[i,])
    }
    else if (model == "lognorm") {
      ev <- exp(log(theta) + 0.5*(1/exp(sim.scale))^2)
      for (i in 1:nrow(pr)) 
        pr[i,] <- rlnorm(length(ev[i,]), meanlog = log(theta[i,]),
                         sdlog = 1/exp(sim.scale[i]))
    }
    list(ev=ev, pr=pr)
  } 
  qi <- ev.pr.surv(sim.coef, sim.scale, x, link)
  if (!is.null(y)) {
    if (is.null(cond.data))
      stop("`cond.data' is required for the exponential, Weibull, and lognormal models.")
    call <- object$call
    call$by <- NULL
    call$data <- cond.data
    est0 <- eval(call, sys.frame())
    par0 <- matrix(param(est0, bootstrap = TRUE), nrow = 1)
    coef0 <- matrix(par0[,1:k], nrow = 1)
    if (model %in% c("weibull", "Weibull", "lognorm"))
      scale0 <- rep(par0[,(k+1):ncol(par0)], nrow(x))
    qi0 <- ev.pr.surv(coef0, scale0, x, link)
    y1 <- y
    tmp <- status
    y[status == 0] <- qi0$pr[(status == 0), 1]
    status[status == 0] <- as.integer(y1[status == 0] < y[status == 0])
    while (sum(status) < length(status)) {
      qi0 <- ev.pr.surv(coef0, scale0, x, link)
      y[status == 0] <- qi0$pr[status == 0, 1]
      status[status == 0] <- as.integer(y1[status == 0] < y[status == 0])
    }
  }
  qi.name <- list(ev="Expected Values: E(Y|X)",
                  pr="Predicted Values: Y|X")
  if (!is.null(x1)) {
    eta1.sim <- sim.coef %*% t(x1)
    ev1 <- apply(eta1.sim, 2, link)
    qi$fd <- ev1 - qi$ev
    qi.name$fd <- "First Differences: E(Y|X1)-E(Y|X)"
  }
  if (!is.null(y)) {
    yvar <- matrix(rep(y, nrow(par)), nrow = nrow(par), byrow = TRUE)
    tmp.ev <- yvar - qi$ev
    tmp.pr <- yvar - qi$pr
    qi$ate.ev <- matrix(apply(tmp.ev, 1, mean), nrow = nrow(par))
    qi$ate.pr <- matrix(apply(tmp.pr, 1, mean), nrow = nrow(par))
    qi.name$ate.ev <- "Average Treatment Effect: Y - EV"
    qi.name$ate.pr <- "Average Treatment Effect: Y - PR"
  }
  list(qi=qi, qi.name=qi.name)
}  















qi.vglm <- function (object, sim.coef, x, x1=NULL, y = NULL) {
  model <- object$zelig
  cm <- object@constraints
  if (model=="mlogit" || model=="oprobit")
    ndim <- (ncol(object@y)-1)
  else if (model=="blogit" || model=="bprobit")
    ndim <- 3
  else if (model=="sur")
    ndim <- ncol(object@y)
  v <- rep(list(NULL), ndim)
  for(i in 1:length(cm)) {
    if(ncol(cm[[i]])==1){
      for(j in 1:ndim)
        if(sum(cm[[i]][j,])==1)
          v[[j]] <- c(v[[j]], names(cm)[i])
    }
    else {
      for (j in 1:ndim)
        if (sum(cm[[i]][j,])==1)
          v[[j]] <- c(v[[j]], paste(names(cm)[i], ":", j, sep=""))
    }
  }
  all.coef <- NULL
  for(i in 1:ndim)
    all.coef <- c(all.coef, list(sim.coef[,v[[i]]]))
  if (model=="mlogit" || model=="oprobit"){
    if(is.null(colnames(object@y)))
      cnames <- ynames <- seq(1,ndim+1,1)
    else
      cnames <- ynames <- colnames(object@y)
    for(i in 1:(ndim+1))
      cnames[i] <- paste("Pr(Y=", ynames[i],")",sep="")
  }
  else if(model=="blogit" || model=="bprobit") 
    cnames <- c("Pr(Y1=0, Y2=0)", "Pr(Y1=0, Y2=1)",
                "Pr(Y1=1, Y2=0)", "Pr(Y1=1, Y2=1)")
  else if(model=="sur"){
    if(is.null(colnames(object@y)))
      cnames <- ynames <- seq(1,ndim,1)
    else
      ynames <- cnames <- colnames(object@y)
    for(i in 1:ndim)
      cnames[i] <- paste("E(Y_", ynames[i],"|X)",sep="")
  }
  else
    stop(paste(model, "is not supported"))
  pp.vglm <- function(object, cm, all.coef, x, ndim, cnames){
    xm <- rep(list(NULL), ndim)
    for (i in 1:length(cm))
      for (j in 1:length(xm))
        if (sum(cm[[i]][j,]) == 1)
          xm[[j]] <- c(xm[[j]], x[,names(cm)[i]])
    sim.eta <- NULL
    for (i in 1:ndim)
      sim.eta <- cbind(sim.eta, all.coef[[i]] %*% as.matrix(xm[[i]]))
    ev <- object@family@inverse(sim.eta)
    colnames(ev) <- cnames
    ev
  }
  pr.vglm <- function(object, ev, ynames) { # To assign predicted values.
    model <- object$zelig
    if (model == "mlogit" || model == "oprobit") {
      k <- ncol(ev)
      Ipr <- sim.cut <- matrix(NA, nrow = nrow(ev), ncol = ncol(ev))
      colnames(Ipr) <- colnames(sim.cut) <- colnames(ev)
      sim.cut[,1] <- ev[,1]
      for (j in 2:k) 
        sim.cut[,j] <- sim.cut[,(j-1)] + ev[,j]
      tmp <- runif(nrow(ev), 0, 1)
      for (l in 1:k) 
        Ipr[,l] <- tmp > sim.cut[,l]
      pr <- NULL
      for (m in 1:nrow(Ipr))
        pr[m] <- 1 + sum(Ipr[m,])
      pr <- factor(pr, levels = sort(unique(pr)), labels = ynames)
      if (model == "mlogit")
        pr <- factor(pr, ordered = FALSE)
      pr <- matrix(pr, nrow = dim(ev)[1])
    }
    else if (model == "blogit" || model == "bprobit") {
      mpr <- cbind((ev[,3] + ev[,4]), (ev[,2] + ev[,4]))
      index <- matrix(NA, ncol = 2, nrow = nrow(mpr))
      for (i in 1:2)
        index[,i] <- rbinom(length(mpr[,i]), 1, mpr[,i])
      pr <- matrix(NA, nrow(sim.coef), 4)
      pr[,1] <- as.character(as.integer(index[,1] == 0 & index[,2] == 0))
      pr[,2] <- as.character(as.integer(index[,1] == 0 & index[,2] == 1))
      pr[,3] <- as.character(as.integer(index[,1] == 1 & index[,2] == 0))
      pr[,4] <- as.character(as.integer(index[,1] == 1 & index[,2] == 1))
      colnames(pr) <- c("(Y1=0, Y2=0)", "(Y1=0, Y2=1)", "(Y1=1, Y2=0)",
                        "(Y1=1, Y2=1)")
    }
#    else if (model == "sur") { # Need to test/add pr.  
#      pr[i,] <- mvrnorm(length(ev[i,], mu = ev[i,], vcov = vcov(object)))
#    }
    pr
  }
  if (nrow(x) == 1) {
    ev <- pp.vglm(object, cm, all.coef, x, ndim, cnames)
    pr <- pr.vglm(object, ev, ynames)
  }
  else {
    ev <- array(dim = c(nrow(sim.coef), if(model=="sur") ndim else (ndim+1), nrow(x)))
    if (model == "mlogit" || model == "oprobit")
      pr <- matrix(nrow=nrow(sim.coef), ncol=nrow(x))
    else if (model == "blogit" || model == "bprobit")
      pr <- array(dim = c(nrow(sim.coef), 4, nrow(x)))
    else if (model == "sur")   #### Need to change this
      pr <- NULL  
    for (i in 1:nrow(x)){
      tmp <- matrix(x[i,], nrow=1)
      colnames(tmp) <- colnames(x)
      rownames(tmp) <- rownames(x)[i]
      tmp.ev <- pp.vglm(object, cm, all.coef, tmp, ndim, cnames)
      tmp.pr <- pr.vglm(object, tmp.ev, ynames)
      ev[,,i] <- tmp.ev
      if (model == "blogit" || model == "bprobit")
        pr[,,i] <- tmp.pr
      else if (model == "mlogit" || model == "oprobit")
        pr[,i] <- tmp.pr
      else if (model == "sur") ##### Placeholder NEED TO FIX!!!!
        pr <- NULL 
    }
    dimnames(ev) <- list(rownames(tmp.ev), colnames(tmp.ev), NULL)
    if (model == "blogit" || model == "bprobit")
      dimnames(pr) <- list(rownames(tmp.pr), colnames(tmp.pr), NULL)
    else if (model == "mlogit" || model == "oprobit")
      dimnames(pr) <- list(c(1:nrow(sim.coef)), NULL)
    else if (model == "sur")  ### Need to change.
      dimnames(pr) <- NULL
  }
  if (model=="mlogit" || model=="oprobit") {
    qi <- list(ev=ev, pr=pr)
    qi.name <- list(ev="Predicted Probabilities: Pr(Y=k|X)",
                    pr="Predicted Values: Y=k|X")
  }
  else if (model=="blogit" || model=="bprobit") {
    qi <- list(ev=ev, pr=pr)
    qi.name <- list(ev="Predicted Probabilities: Pr(Y1=k,Y2=l|X)",
                    pr="Predicted Values: (Y1,Y2)|X")
  }
  else if (model=="sur") {
    qi <- list(ev=ev, pr=pr)
    qi.name <- list(ev="Expected Values: E(Yj|X)",
                    pr="Predicted Values: Yj|X")
  }
  if (!is.null(x1)) {
    if (nrow(x1) == 1)
      ev1 <- pp.vglm(object, cm, all.coef, x1, ndim, cnames)
    else {
      ev1 <- array(dim = c(nrow(sim.coef), if(model=="sur") ndim else (ndim+1), nrow(x1)))
      for (i in 1:nrow(x1)) {
        tmp <- matrix(x1[i,], nrow=1)
        colnames(tmp) <- colnames(x1)
        rownames(tmp) <- rownames(x1)[i]
        tmp <- pp.vglm(object, cm, all.coef, tmp, ndim, cnames)
        ev1[,,i] <- tmp
      }
      dimnames(ev1) <- list(rownames(tmp), colnames(tmp), NULL)
    }
    qi$fd <- ev1 - ev
    if (model=="mlogit" || model=="oprobit") {
      qi$rr <- ev1 / ev
      qi.name$fd <- "First Differences: Pr(Y=k|X1) - Pr(Y=k|X)"
      qi.name$rr <- "Risk Ratios: Pr(Y=k|X1) / Pr(Y=k|X)"
    }
    else if (model=="blogit" || model=="bprobit") {
      qi$rr <- ev1/ev
      qi.name$fd <- "First Differences: Pr(Y1=k,Y2=l|X1) - Pr(Y1=k,Y2=l|X)"
      qi.name$rr <- "Risk Ratios: Pr(Y1=k,Y2=l|X1) / Pr(Y1=k,Y2=l|X)"
    }
    else if(model=="sur")
      qi.name$fd <- "First Differences: E(Yj|X1j) - E(Yj|Xj)"
  }
  if (!is.null(y)) {
    tmp.ev <- tmp.pr <- array(NA, dim = dim(qi$ev))
    qi$ate.ev <- qi$ate.pr <- matrix(NA, dim(qi$ev)[1], dim(qi$ev)[2])
    if (model=="mlogit" || model=="oprobit") {
      yvar <- matrix(NA, nrow = length(y), ncol = length(ynames))
      pr.idx <- array(NA, dim = c(nrow(pr), length(ynames), nrow(x)))
      for (i in 1:length(ynames)) {
        yvar[,i] <- as.integer(y == ynames[i])
        pr.idx[,i,] <- as.integer(pr[,i] == ynames[i])
      }
      colnames(qi$ate.ev) <- colnames(qi$ate.pr) <- ynames
    }
    else if (model=="blogit" || model=="bprobit") {
      yvar <- matrix(NA, nrow = nrow(y), ncol = 4)
      yvar[,1] <- as.integer(y[,1] == 0 & y[,2] == 0)
      yvar[,2] <- as.integer(y[,1] == 0 & y[,2] == 1)
      yvar[,3] <- as.integer(y[,1] == 1 & y[,2] == 0)
      yvar[,4] <- as.integer(y[,1] == 1 & y[,2] == 1)
      pr.idx <- array(NA, dim = c(nrow(pr), 4, nrow(x)))
      for (i in 1:4)
        pr.idx[,i,] <- as.integer(pr[,i,])
      colnames(qi$ate.ev) <- colnames(qi$ate.pr) <-
        c("(Y1=0, Y2=0)", "(Y1=0, Y2=1)",
          "(Y1=1, Y2=0)", "(Y1=1, Y2=1)")
    }
    for (j in 1:ncol(yvar)) {
      for (i in 1:nrow(sim.coef)) {
        tmp.ev[i,j,] <- yvar[,j] - qi$ev[i,j,]
        tmp.pr[i,j,] <- yvar[,j] - pr.idx[i,j,]
      }
      qi$ate.ev[,j] <- apply(tmp.ev[,j,], 1, mean)
      qi$ate.pr[,j] <- apply(tmp.pr[,j,], 1, mean)
    }
    qi.name$ate.ev <- "Average Treatment Effect: Y - EV"
    qi.name$ate.pr <- "Average Treatment Effect: Y - PR"
  }
  list(qi=qi, qi.name=qi.name)
}
















relogit <- function(formula, data=sys.parent(), tau=NULL,
                    bias.correct=TRUE, ...){
  mf <- match.call()
  mf$tau <- mf$bias.correct <- NULL
  if (!is.null(tau))
    tau <- unique(tau)
  relog <- function(formula=formula, data=data, tau=tau,
                    bias.correct=bias.correct, ...) {
    mf[[1]] <- as.name("glm")
    mf$family <- as.name("binomial")
    res <- eval(as.call(mf))
    ## prior correction 
    y <- res$y  
    ybar <- mean(y)
    #w1 <- tau/ybar
    #w0 <- (1-tau)/(1-ybar)
    #wi <- w1*y + w0*(1-y)
    if (!is.null(tau)){      
      if (tau <= 0 || tau >= 1) 
        stop("\ntau needs to be between 0 and 1.\n") 
      res$coefficients["(Intercept)"] <- res$coefficients["(Intercept)"] - 
        log(((1-tau)/tau) * (ybar/(1-ybar)))
    }
    X <- model.matrix(res)
    if (bias.correct){
      pihat <- fitted(res)
     # if (is.null(tau))
     #   W <- pihat * (1 - pihat) * wi
     # else 
        W <- pihat * (1 - pihat)
      Qdiag <- lm.influence(glm(y ~ X-1, weights=W))$hat
      #if (is.null(tau)) 
      #  xi <- 0.5 * Qdiag * ((1+w0)*pihat-w0)
      #else 
        xi <- 0.5 * Qdiag * (2*pihat - 1)
      res$coefficients <- res$coefficients -
        glm(xi ~ X - 1, weights=W)$coefficients 
    }
    res$linear.predictors <- t(res$coef) %*% t(X) 
    res$fitted.values <- 1/(1+exp(-res$linear.predictors))
    res
  }
  if (length(tau) > 2)
    stop("tau must be a vector of length less than or equal to 2")
  else if (is.null(tau) || length(tau) == 1) 
    res <- relog(mf, data, tau, bias.correct)
  else if(length(tau) == 2){ # without full population information
    lower <- relog(mf, data, min(tau), bias.correct)
    upper <- relog(mf, data, max(tau), bias.correct)
    lower$zelig <- upper$zelig <- "relogit"
    res <- list(lower.estimate=lower, upper.estimate=upper,
                terms = lower$terms)
  }
  res$zelig <- "relogit"
  class(res) <- "relogit"
  res
}

repl <- function(object, ...)
  UseMethod("repl")
repl.default <- function(object, data=NULL, ...) {
  if (!is.null(data))
    object$call$data <- data
  eval(object$call, sys.parent())
}
repl.zelig <- function(object, prev=NULL, data=NULL, x=NULL,
  x1=NULL, bootfn=NULL, ...) {
  rep.zelig <- eval(object$zelig.call, sys.parent())
  object$call$object <- rep.zelig
  if (is.null(x))
    object$call$x <- object$x
  else
    object$call$x <- x
  if (is.null(x1))
    object$call$x1 <- object$x1
  else
    object$call$x1
  if (!is.null(prev))
    object$call$prev <- prev
  eval(object$call, sys.parent())
}




rocplot <- function(y1, y2, fitted1, fitted2,
                    cutoff = seq(from=0, to=1, length=100), lty1="solid",
                    lty2="dashed", lwd1=par("lwd"), lwd2=par("lwd"),
                    col1=par("col"), col2=par("col"), main="ROC Curve",
                    xlab = "Proportion of 1's Correctly Predicted",
                    ylab="Proportion of 0's Correctly Predicted", ...) {
  op <- par(no.readonly = TRUE)
  par(mar = c(4,4,2,1), tcl = -0.25, mgp = c(2, 0.6, 0))
  ones <- zeros <- matrix(NA, nrow = length(cutoff), ncol = 2)
  for (i in 1:length(cutoff)) {
    ones[i,1] <- mean(fitted1[y1==1] >= cutoff[i]) 
    ones[i,2] <- mean(fitted2[y2==1] >= cutoff[i])
    zeros[i,1] <- mean(fitted1[y1==0] < cutoff[i])
    zeros[i,2] <- mean(fitted2[y2==0] < cutoff[i])
  }
  plot(0:1, 0:1, type = "n", xaxs = "i", yaxs = "i",
       main=main, xlab=xlab, ylab=ylab, ...)
  lines(cbind(ones[,1], zeros[,1]), lty = lty1, lwd = lwd1, col=col1)
  lines(cbind(ones[,2], zeros[,2]), lty = lty2, lwd = lwd2, col=col2)
  abline(1, -1, lty = "dotted")
  par(op)
}
















setx.MI <- function(object, ...) {
  M <- length(object)
  data <- NULL
  ## combine data
  for (i in 1:M) {
    tmp <- as.data.frame(eval(object[[i]]$call$data, sys.parent()))
    data <- rbind(data, tmp)
  }
  model <- object[[1]]
  model$call$data <- as.data.frame(data)
  ## call setx()
  mf <- match.call()
  mf$object <- model
  mf[[1]] <- as.name("setx")
  eval(as.call(mf), sys.frame(sys.parent()))
}

setx<-function(object, ...)
  UseMethod("setx")

setx.default <- function(object, fn = list(numeric = mean, ordered =
                                   median, other = mode),
                         data = NULL, cond = FALSE, ...){
  mf <- match.call()
  mode <- function(x){
    tb <- tapply(x, x, length)
    if(is.factor(x))
      value <- factor(unlist(labels(tb[seq(along=tb)[tb==max(tb)]])),
                      levels=levels(x))
    else if (is.logical(x))
      value <- as.logical(unlist(labels(tb[seq(along=tb)[tb==max(tb)]])))
    else if (is.character(x))
      value <- as.character(unlist(labels(tb[seq(along=tb)[tb==max(tb)]])))
    else
      stop(paste(vars[i], "is not a supported variable type."))
    if (length(value)>1) {
      warning("There is more than one mode. The first level is selected.")
      value <- sort(value)[1]
    }
    return(value)
  }
  median.default <- median
  median <- function(x) {
    if(is.numeric(x))
      value <- median.default(x)
    else if (is.ordered(x))
      value <- factor(levels(x)[median.default(as.integer(x))],
                      levels=levels(x)) 
    else
      stop("median cannot be calculated for this data type")
    return(value)
  }
  max.default <- max
  max <- function(x, na.rm=FALSE) {
    if(is.numeric(x))
      value <- max.default(x, na.rm=na.rm)
    else if (is.ordered(x)) 
      value <- factor(levels(x)[length(levels(x))], levels=levels(x))
    else
      stop("max cannot be calculated for this data type")
    return(value)
  }
  min.default <- min
  min <- function(x, na.rm=FALSE) {
    if(is.numeric(x))
      value <- min.default(x, na.rm = na.rm)
    else if (is.ordered(x))
      value <- factor(levels(x)[1], levels=levels(x))
    else
      stop("min cannot be calculated for this data type")
    return(value)
  }
  if (any(class(object)=="vglm"))
    tt <- object@terms$terms
  else
    tt <- terms(object)
  if (is.null(data)) 
    dta <- eval(object$call$data, sys.parent())
  else
    dta <- data
  data <- dta <- dta[complete.cases(model.frame(tt, dta)), ]
  vars <- names(dta)
  if (cond) {
    if (is.null(data)) 
      stop(paste("if cond = TRUE, you must specify the data frame."))
    if (!is.null(fn)) {
      warning(paste("when cond = TRUE, fn is coerced to NULL"))
      fn <- NULL
    }
    if (class(object)[1] == "vglm")
      Y <- model.response(model.frame(tt, data=dta), "any")
    else
      Y <- model.extract(model.frame(tt, data=dta), "response")
  }
  else if (!is.null(fn)) {
    if (is.null(fn$numeric) || !is.function(fn$numeric)) {
      warning("fn$numeric coerced to mean().")
      fn$numeric <- mean
    }
    if (is.null(fn$ordered) || !is.function(fn$ordered) || 
        identical(mean, fn$ordered)) {
      warning("fn$ordered coreced to median().")
      fn$ordered <- median
    }
    else if (identical(min.default, fn$ordered)) 
      fn$ordered <- min
    else if (identical(max.default, fn$ordered)) 
      fn$ordered <- max
    else if (identical(median.default, fn$ordered)) 
      fn$ordered <- median
    if (is.null(fn$other) || !is.function(fn$other)) { 
      warning("the only available fn for other is mode.")
      fn$other <- mode
    }
    for (i in 1:ncol(dta)) {
      v <- na.omit(dta[,i])
      if (any(vars[i] == names(model.frame(tt, dta)))) {
        if (is.numeric(v))
          value <- lapply(list(v), fn$numeric)[[1]]
        else if (is.ordered(v)) 
          value <- lapply(list(v), fn$ordered)[[1]]
        else 
          value <- lapply(list(v), fn$other)[[1]]
        data[,i] <- value
      }
    }
    opt <- vars[na.omit(pmatch(names(mf), vars))]
    maxl <- 1
    if (length(opt) > 0)
      for (i in 1:length(opt)) {
        value <- eval(mf[[opt[i]]], sys.parent())
        lv <- length(value)
        if (lv>1)
          if (maxl==1 || maxl==lv) {
            maxl <- lv
            data <- data[1:lv,]
          }
          else
            stop("vector inputs should have the same length.")
        if (is.factor(data[,opt[i]]))
          data[,opt[i]] <- list(as.factor(value))
        else if (is.numeric(data[,opt[i]]))
          data[,opt[i]] <- list(as.numeric(value))
        else if (is.logical(data[,opt[i]]))
          data[,opt[i]] <- list(as.logical(value))
        else
          data[,opt[i]] <- list(value)
      }
    data <- data[1:maxl,]
  }
  X <- model.matrix(delete.response(tt), data = data, contrasts =
                    if (length(object$contrasts)>0) object$contrasts
                    else NULL, xlev = object$xlevels, ...)
  if (cond) {
    X <- cbind(Y, X)
    class(X) <- "setx.cond"
  }
  else 
    class(X) <- "setx"
  row.names(X) <- 1:nrow(X)
  return(as.matrix(X))
}










setx.multiple <- function(object, data = NULL, cond = FALSE, ...){
  dta <- eval(object[[1]]$call$data, sys.parent())
  d <- x <- list()
  by <- object[[1]]$call$by
  idx <- pmatch(by, names(dta))[1]
  lev <- sort(unique(eval(dta[[idx]], sys.parent())))
  if (cond) {
    for (i in 1:length(lev)) {
      d[[i]] <- as.data.frame(dta[dta[[idx]] == lev[i],])
      x[[i]] <- setx.default(object[[i]], data = d[[i]], cond = cond, ...)
      class(x[[i]]) <- "setx.cond"
    }
  }
  else {
    x <- setx.default(object[[1]], data = data, cond = cond, ...)
    class(x) <- "setx"
  }
  x 
}




sim.MI <- function(object, x, x1 = NULL, num = c(1000, 100), prev =
                   NULL, bootstrap = FALSE, bootfn = NULL, ...){
  M <- length(object)
  simpar <- simqi <- res <- list()
  if (is.null(prev)) {
    if (!bootstrap) {
      if (length(num) == 2)
        num <- num[1]
      numM <- round(num/M)
      for (i in 1:M) 
        simpar[[i]] <- param(object[[i]], num=numM, bootstrap=bootstrap)
    }
    else {
      if (length(num) == 2)
        num <- num[2]
      numM <- round(num/M)
      if (any(class(object[[1]]) == "vglm")) 
        tt <- object[[1]]@terms$terms
      else
        tt <- terms(object[[1]])
      if (missing(bootfn)) {
        bootfn <- function(data, i, object) {
          d <- data[i,]
          object$call$data <- d
          fit <- eval(object$call, sys.parent())
          l1 <- length(param(fit, bootstrap = TRUE))
          l <- length(param(object, bootstrap = TRUE))
          while (l > l1) {
            d1 <- data[sample(nrow(data), replace = TRUE),]
            fit <- update(object, data = d1)
            l1 <- length(param(fit, bootstrap = TRUE))
          }
          return(param(fit, bootstrap = TRUE))
        }
      }
      for (i in 1:M) {
        dta <- eval(object[[i]]$data, sys.parent())
        res <- boot(dta, bootfn, R = num, object = object[[i]],  ...)
        colnames(res$t) <- names(res$t0)
        simpar[[i]] <- res$t
      }
    }
    par <- as.matrix(simpar[[1]])
    for (j in 2:M)
      par <- rbind(par, as.matrix(simpar[[j]]))
  }
  else {
    if (bootstrap)
      stop("Error: Choosing 'bootstrap = TRUE' generates new parameters.  \nIf you wish to use previously generated parameters, \nplease specify only 'prev'.")  
    else
      simpar <- prev
  }
  simqi <- qi(object[[1]], par, x=x, x1=x1)
  c <- match.call()
  c$num <- num
  res <- list(x=x, x1=x1, call = c, zelig.call = object[[1]]$call, 
              par = simpar, qi=simqi$qi, qi.name = simqi$qi.name)
  class(res) <- "zelig"
  res
}







sim <- function(object, ...)
  UseMethod("sim")
sim.default <- function(object, x, x1=NULL, num=c(1000, 100), prev =
      NULL, bootstrap = FALSE, bootfn=NULL, cond.data = NULL, ...) {
  if (is.null(dim(x)))
    x <- matrix(x, nrow = 1)
  if (class(x) == "setx.cond") {
    if (!is.null(x1)) {
      warning("First Differences are not calculated in conditional prediction models.")
      x1 <- NULL
    }
    if (object$call$model %in% c("bprobit", "blogit")) {
      yvar <- x[,1:2]
      x <- x[,3:ncol(x)]
    }
    else {
      yvar <- x[,1]
      x <- x[,2:ncol(x)]
    }
    class(x) <- "setx.cond"
  }
  if (any(class(object) == "MCMCZelig"))
    num <- nrow(object$coefficients)
  else if (!bootstrap)
    num <- num[1]
  else
    num <- num[2]
  if (is.null(prev)) {
    if (!bootstrap & any(class(object) != "relogit"))
      simpar <- param(object, num=num, bootstrap=bootstrap)
    else if (class(object) == "relogit") 
      simpar <- param.relogit(object, num=num, x=x, bootstrap=bootstrap, bootfn=bootfn, ...) 
    else {
      if (any(class(object)=="vglm")) 
        tt <- object@terms$terms
      else
        tt <- terms(object)
      dta <- eval(object$data, sys.parent())
      dta <- dta[complete.cases(model.frame(tt, dta)),]
      if (is.null(bootfn)) {
        bootfn <- function(data, i, object) {
          d <- data[i,]
          object$call$data <- d
          fit <- eval(object$call, sys.parent())
          l <- length(param(object, bootstrap = TRUE))
          l1 <- length(param(fit, bootstrap = TRUE))
          while (l > l1) {
            object$call$data <- data[sample(nrow(data), replace=TRUE),]
            fit <- eval(object$call, sys.parent())
            l1 <- length(param(fit, bootstrap = TRUE))
          }
          return(param(fit, bootstrap = TRUE))
        }
      }
      res <- boot(dta, bootfn, R = num, object = object, ...)
      colnames(res$t) <- names(res$t0)
      simpar <- res$t
    }
  }
  else {
    if (bootstrap)
      stop("Error: Choosing 'bootstrap = TRUE' generates new parameters.  \nIf you wish to use previously generated parameters, \nplease specify only 'prev'.")
    else
      simpar <- prev
  }
  if (class(x) == "setx.cond") {
    if (class(object)[1] == "survreg") 
      simqi <- qi(object, simpar, x = x, x1 = x1, y = yvar, cond.data = cond.data)
    else
      simqi <- qi(object, simpar, x = x, x1 = x1, y = yvar)
    class(x) <- "setx.cond"
  }
  else 
    simqi <- qi(object, simpar, x = x, x1 = x1, y = NULL)
  c <- match.call()
  c$num <- num
  res <- list(x=x, x1=x1, call = c, zelig.call = object$call,
              par = simpar, qi=simqi$qi, qi.name=simqi$qi.name)
  class(res) <- "zelig"
  res
}
















sim.multiple <- function(object, x, x1 = NULL, num = 100, ...){
  m <- length(object)
  dta <- eval(object[[1]]$data, sys.parent())
  res <- simpar <- list()
  by <- object[[1]]$call$by
  idx <- pmatch(by, names(dta))
  lev <- sort(unique(eval(dta[[idx]], sys.parent())))
  for (i in 1:length(lev)) {
    if (class(x) == "list") {
      object[[i]]$data <- eval(dta[[idx]] == lev[i])
      res[[i]] <- sim(object[[i]], x = x[[i]], x1 = x1[[i]], num = num, ...)
      res[[i]]$nx <- nrow(x[[i]])/nrow(dta)
    }
    else if (class(x) == "setx") {
      res[[i]] <- sim(object[[i]], x = x, x1 = x1, num = num, ...)
      res[[i]]$nx <- nrow(dta[dta[[idx]] == lev[i],])/nrow(dta)
    }
  }
  class(res) <- "zelig.multiple"
  names(res) <- names(object)
  res
}







summary.MI <- function(object, subset = NULL, ...){
  res <- list()
  if(is.null(subset))
    M <- 1:length(object)
  else
    M <- c(subset)
  for (i in M)
    res[[i]] <- summary(object[[i]])
  ans <- res[[1]]
  ans$zelig <- object[[1]]$zelig
  ans$call <- object[[1]]$call
  if (length(M) > 1) {
    ans$all <- res
    coef <- se <- NULL
    for (i in M){
      if (any(class(object[[1]]) == "vglm")) {
        coef <- cbind(coef, res[[i]]@coef3[,1])
        se <- cbind(se, res[[i]]@coef3[,2])
      }
      else {
        coef <- cbind(coef, coef(res[[i]])[,1])
        se <- cbind(se, coef(res[[i]])[,2])
      }
    }
    Q <- apply(coef, 1, mean)
    U <- apply(se^2, 1, mean)
    B <- apply((coef-Q)^2, 1, sum)/(length(M)-1)
    var <- U+(1+1/length(M))*B
    nu <- (length(M)-1)*(1+U/((1+1/length(M))*B))^2
    if (any(class(object[[1]]) == "vglm")) {
      coef.table <- matrix(NA, nrow(ans@coef3), 4)
      dimnames(coef.table) <- list(c(dimnames(ans@coef3)[[1]]),
                                   c("Value", "Std. Error", "t-stat", "p-value"))
    }
    else {
      coef.table <- matrix(NA, nrow = nrow(coef(ans)), ncol = 4)
      dimnames(coef.table) <- list(rownames(coef(ans)),
                                   c("Value", "Std. Error", "t-stat", "p-value"))
    }
    coef.table[,1] <- Q
    coef.table[,2] <- sqrt(var)
    coef.table[,3] <- Q/sqrt(var)
    coef.table[,4] <- pt(abs(Q/sqrt(var)), df=nu, lower.tail=F)*2
    ans$coefficients <- coef.table
    ans$cov.scaled <- ans$cov.unscaled <- NULL
    if (!any(class(object[[1]]) == "vglm")) {
      for (i in 1:length(ans)) {
        if (is.numeric(ans[[i]]) && !names(ans)[i] %in% c("coefficients")){
          tmp <- NULL
          for (j in M) 
            tmp <- cbind(tmp, res[[j]][[pmatch(names(ans)[i], names(res[[j]]))]])
          ans[[i]] <- apply(tmp, 1, mean) 
        }
      }
    }
    class(ans) <- "summary.MI"
  }
  else
    ans <- summary(object[[M]])
  ans
}

summary.multiple <- function(object, subset = NULL, ...) {
  res <- list()
  if(is.null(subset))
    M <- c(1:length(object))
  else
    M <- c(subset)
  for (i in 1:length(M))
    res[[i]] <- summary(object[[i]])
  names(res) <- paste("summary.", object[[1]]$call$by, names(object)[M], sep = "")
  res$zelig <- object[[1]]$zelig
  res$call <- object[[1]]$call
  res$lev <- names(object)[M]
  res$M <- M
  class(res) <- "summary.multiple"
  res
}



summary.relogit <- function(object, ...) {
  if (!is.numeric(object$call$tau))
    tau <- eval(object$call$tau, sys.frame(sys.parent()))
  else
    tau <- object$call$tau
  dta <- eval(object$call$data, sys.parent())
  tt <- terms(object)
  dta <- dta[complete.cases(model.frame(tt, dta)),]
  n <- nrow(dta)
  k <- ncol(dta)
  bias.cor <- function(res.obj, bc, n, k){
    if (bc) {
      res.obj$cov.unscaled <- res.obj$cov.unscaled * (n/(n+k))^2
      res.obj$cov.scaled <- res.obj$cov.unscaled * res.obj$dispersion
      res.obj$coef[,2] <- sqrt(diag(res.obj$cov.scaled))
      res.obj$coef[,3] <- res.obj$coef[,1] / res.obj$coef[,2]
      res.obj$coef[,4 ] <- 2*pt(-abs(res.obj$coef[,3]), res.obj$df.residual)
    }
    res.obj
  }
  if(length(tau) == 1 || length(tau) == 0){
    res <- summary.glm(object, ...)
    res <- bias.cor(res, object$call$bias.correct, n, k)
    res$call <- object$call
  }
  else if (length(tau) == 2) {
    res <- list()
    res0 <- summary.glm(object$lower.estimate, ...)
    res0$call <- object$call
    res$lower.summary <- bias.cor(res0, object$call$bias.correct, n, k)
    res1 <- summary.glm(object$upper.estimate, ...)
    res1$call <- object$call
    res$upper.summary <- bias.cor(res1, object$call$bias.correct, n, k)
  }
  class(res)<-c("summary.relogit", class(object))
  res
}












summary.setx <- function(x) {
  res <- matrix(NA, nrow(x), ncol(x))
  for (i in 1:ncol(x))
    res[,i] <- x[,i]
  colnames(res) <- colnames(x)
  res <- as.data.frame(res)
  summary(res)
}
summary.setx.cond <- function(x) {
  res <- matrix(NA, nrow(x), ncol(x))
  for (i in 1:ncol(x))
    res[,i] <- x[,i]
  colnames(res) <- colnames(x)
  res <- as.data.frame(res)
  summary(res)
}
summary.zelig<-function(object, subset = NULL, CI=95, 
                        stats=c("mean", "sd", "min", "max"), ...){
  cip <- c((100-CI)/200, 1-(100-CI)/200)
  model <- object$zelig.call$model
  summarize <- function(x, cip, stats, model){
    res <- NULL
    if (is.numeric(x)) {
      if (!is.null(stats))
        for (i in 1:length(stats)) 
          res <- cbind(res, apply(x, 2, stats[i]))
      res <- cbind(res, apply(x, 2, quantile, prob=cip[1]))
      res <- cbind(res, apply(x, 2, quantile, prob=cip[2]))
      colnames(res) <- c(stats, paste(cip[1]*100, "%", sep=""),
                         paste(cip[2]*100, "%", sep=""))
    }
    else {
      res <- t(apply(if(is.matrix(x)) x else as.matrix(x), 2, table))
      if (length(unique(x))==1) {
        colnames(res) <- unique(x)
        if(model=="relogit")
          if(colnames(res)==1) {
            res <- cbind(0, res)
            colnames(res)[1] <-0
          }
          else {
            res <- cbind(res, 0)
            colnames(res)[2] <-1
          }
      }
      res <- res/nrow(x) * 100
    }
    return(res)
  }
  qi.stats<-list()
  X <- object$x
  X1 <- object$x1
  if (is.null(dim(X))) {
    X <- matrix(X, ncol = 1)
    colnames(X) <- "(Intercept)"
    if (!is.null(X1)) {
      X1 <- matrix(X1, ncol = 1)
      colnames(X1) <- "(Intercept)"
    }
  }
  for (i in 1:length(object$qi)) {
    qi <- object$qi[[i]]
    if (nrow(X) > 1) {
      if (length(dim(qi)) == 3) {
        if (is.function(subset)){  # subset = all; all is class "function"
          tmp <- summarize(qi[,,1], cip, stats, model)
          tmp <- array(NA, dim=c(nrow(tmp), ncol(tmp), dim(qi)[3]),
                       dimnames=list(dimnames(qi)[[2]], dimnames(tmp)[[2]],
                       rownames(object$x)))
          for (j in 1:dim(qi)[3])
            tmp[,,j] <- summarize(qi[,,j], cip, stats, model)
          qi.stats[[i]] <- tmp
        }
        if(is.null(subset)){# subset = NULL; summarizes all obs at once. 
          tmp <-NULL
          for (j in 1:dim(qi)[3])
            tmp <- rbind(tmp, qi[,,j])
          qi.stats[[i]] <- summarize(tmp, cip, stats, model)
          object$qi.name[[i]] <- paste("Pooled", object$qi.name[[i]])
        }
        if (is.numeric(subset)) { # subset=integer, summarizes identified obs
          k <- dim(qi)[[3]]
          if (subset < k || subset == k) {
            tmp <- summarize(qi[,,1], cip, stats, model)
            tmp <- array(NA, dim=c(nrow(tmp), ncol(tmp), length(subset)),
                         dimnames=list(dimnames(qi)[[2]],
                           dimnames(tmp)[[2]], subset)) 
            for (l in 1:length(subset)){
              tmp[,,l] <- summarize(qi[,,subset[l]], cip, stats, model)
              qi.stats[[i]] <- tmp
              X <- object$x[subset,]
              if (!is.null(object$x1))
              X1 <- object$x1[subset,]
            }
          }
          else 
            print("Error: Selected observation to subset is not in available range of observations.")
        }
      }
      if (length(dim(qi)) == 2) {
        if (is.function(subset)) 
          qi.stats[[i]] <- summarize(qi, cip, stats, model)
        if (is.null(subset)) {
          tmp <- matrix(c(qi), ncol=1)
          tmp <- summarize(tmp, cip, stats, model)
          qi.stats[[i]] <- as.matrix(tmp)
          names(qi.stats[[i]]) <- colnames(tmp)
          object$qi.name[[i]] <- paste("Pooled", object$qi.name[[i]])
        }
        if (is.numeric(subset)) {
          k <- dim(qi)[[2]]
          if (subset < k || subset == k){
            qi.stats[[i]] <- summarize(as.matrix(qi[,subset]), cip,
                                       stats, model)
            X <- X[subset,]
            if(!is.null(object$x1))
              X1 <- X1[subset,]
          }
          else
            print("Error: Selected observation to subset is not in available range of observations.")
        }
      }
    }
    else
      qi.stats[[i]] <- summarize(qi, cip, stats, model)
  }
  names(qi.stats) <- names(object$qi)
  res <- list(zelig=model, num=object$call$num, x=X, x1=X1,
              qi.stats=qi.stats, qi.name=object$qi.name) 
  class(res)<-"summary.zelig"
  res
}




















summary.zelig.multiple <-function(object, subset = NULL, CI=95, 
                        stats=c("mean", "sd", "min", "max"), ...){
  m <- length(object)
  if(is.null(subset)) { # summary for all strata together
    w <- NULL
    for (j in 1:m) 
      w <- c(w, object[[j]]$nx)
    w <- round(w*object[[1]]$call$num)
    qi <- object[[1]]$qi
    for (i in 1:length(qi)){
      qi1i <- object[[1]]$qi[[i]]
      if(class(object[[1]]$x)=="setx.cond") # conditional prediction
        if(length(dim(qi1i))==3){
          tdim <- dim(qi1i)
          for (j in 2:m)
            tdim[2:3] <- tdim[2:3] + dim(object[[j]]$qi[[i]])[2:3]
          tmp <- array(NA, dim=tdim)
          tmp[,1:dim(qi1i)[2],1:dim(qi1i)[3]] <- qi1i
          tdim <- dim(qi1i)
        }
        else
          tmp <- as.matrix(qi1i)
      else # unconditional prediction
        if(length(dim(qi1i))==3){
          tmp <- array(NA, dim=c(sum(w), dim(qi1i))[2:3])
          tmp[1:w[1],,] <- qi1i[1:w[1],,] 
        }
        else
          tmp <- as.matrix(qi1i[1:w[1],])
      for (j in 2:m) {
        qiji <- object[[j]]$qi[[i]]
        if(class(object[[j]]$x)=="setx.cond") # conditional prediction
          if(length(dim(qi1i))==3) {
            tmp[,(tdim[2]+1):(tdim[2]+dim(qiji[2])),(tdim[3]+1):(tdim[3]+dim(qiji[3]))] <- qiji
            tdim <- tdim + dim(qiji)
          }
          else
            tmp <- cbind(tmp, qiji)
        else{ # unconditional prediction
          if(length(dim(qi1i))==3)
            tmp[(sum(w[1:(j-1)])+1):sum(w[1:j]),,] <- qiji[1:w[j],,]
          else
            tmp <- rbind(tmp, as.matrix(qiji[1:w[j],]))
        }
      }
      qi[[i]] <- tmp
    }
    c <- match.call()
    c$num <- object[[1]]$call$num
    res <- list(qi=qi, qi.name=object[[1]]$qi.name,
                x=object[[1]]$x, x1=NULL, call=c,
                zelig.call=object[[1]]$zelig)
  }
  else { # summary for each strata
    res <- list()
    if(is.function(subset)){
      m <- length(object)
      subset <- 1:m
    }
    else
      m <- length(subset)
    for (i in 1:m) {
      res[[i]] <- summary.zelig(object[[subset[i]]], subset=NULL, ...)
      names(res)[i] <- names(object)[i]
    }
    class(res) <- "summary.zelig.multiple"
    return(res)
  }
  return(summary.zelig(res, ...))  
}



















ternarypoints <- function(object, pch = 19, col = "blue", ...){
    s <- rowSums(object)
    if (any(s <= 0))
        stop("each row of the input `object' must have a positive sum")
    object <- object/s
    top <- sqrt(3)/2
    xp <- object[, 2] + object[, 3]/2
    yp <- object[, 3] * top
    points(xp, yp, pch = pch, col = col, ...)
}
vcov.lm <- function(object) {
  so <- summary.lm(object, corr = FALSE)
  so$sigma^2 * so$cov.unscaled
}
vcov.relogit <- function(object) {
  dta <- eval(object$call$data, sys.parent())
  tt <- terms(object)
  dta <- dta[complete.cases(model.frame(tt, dta)),]
  n <- nrow(dta)
  k <- ncol(dta)
  bias.cor <- function(res.obj, bc, n, k){
    if (bc) {
      res.obj$cov.unscaled <- res.obj$cov.unscaled * (n/(n+k))^2
      res.obj$cov.scaled <- res.obj$cov.unscaled * res.obj$dispersion
      res.obj$coef[,2] <- sqrt(diag(res.obj$cov.scaled))
      res.obj$coef[,3] <- res.obj$coef[,1] / res.obj$coef[,2]
      res.obj$coef[,4 ] <- 2*pt(-abs(res.obj$coef[,3]), res.obj$df.residual)
    }
    res.obj
  }
  res <- summary.glm(object)
  res <- bias.cor(res, object$call$bias.correct, n, k)
  res$cov.scaled
}
zelig <- function(formula, model, data, by = NULL, ...) {
  fn <- paste("zelig2", model, sep = "")
  if (!exists(fn))
    stop("Error:", model, "not supported. Type help.zelig(\"models\") to list supported models.")
  mf <- match.call()
  if (missing(by))
    by <- NULL
  if (!is.null(by)) {
    if (!is.data.frame(data))
      stop("Error: subsetting using 'by' is currently not supported for MI datasets.")  
    temp <- list()
    if (as.character(by) %in% c(formula[[2]], formula[[3]]))
      stop("Error: the variable selected for subsetting cannot be called in the formula.")
    if (as.character(by) %in% names(data)) {
      mf$by <- as.name(by)
      by <- eval(mf$by, data)
      mf$by <- NULL
      lev <- sort(unique(by))
      for (i in 1:length(lev))
        temp[[i]] <- data[by == lev[i],]
      names(temp) <- lev
      data <- temp
    }
  }
  object <- list()
  M <- 1
  if (!is.data.frame(data)) {
    M <- length(data)
    dta <- data
  }
  else
    d <- data
  mf <- do.call(fn, list(formula, model, data, M, ...))
  for (i in 1:M){
    if (M > 1) 
      d <- eval(dta[[i]], sys.parent())
    if (is.data.frame(d)){
      d <- d[complete.cases(model.frame(as.formula(formula), d)),]
      mf$data <- d
      res <- eval(as.call(mf))
      res$call <- match.call()
      res$data <- res$call$data
      res$zelig <- model
      if (M > 1) 
        object[[i]] <- res
      else
        object <- res
    }
    else
      stop(paste("Error:", class(data), "is not a supported data format."))
  }
  if (model == "relogit") {
    if (is.null(object$call$bias.correct))
      object$call$bias.correct <- TRUE
    else
      object$call$bias.correct <- FALSE
    tau <- eval(object$call$tau, sys.parent())
    if (length(tau) == 2) {
      object$lower.estimate$call <- object$upper.estimate$call <- object$call
      object$lower.estimate$call$tau <- round(min(tau),
                                              digits = getOption("digits"))
      object$upper.estimate$call$tau <- round(max(tau),
                                              digits = getOption("digits"))
    }
  }
  if (M > 1)
    if (is.null(by)) 
      class(object) <- "MI"
    else {
      names(object) <- lev
      class(object) <- "multiple"
    }
  object
}
zelig2blogit <- function(formula, model, data, M, ...) {
  require(VGAM)
  mf <- match.call(expand.dots = TRUE)
  if(is.null(mf$ones))
    mf$ones <- 3
  mf[[1]] <- as.name("vglm")
  mf$family <- as.name("blogit")
  tmp <- cmvglm(mf$formula, mf$model, mf$equal, mf$zeros, mf$ones, 3)
  mf$formula <- tmp$formula 
  mf$constraints <- tmp$constraints
  blogit <<- function() binom2.or(zero=NULL)
  mf$model <- mf$equal <- mf$ones <- mf$zeros <- mf$M <- NULL
  as.call(mf)
}
zelig2bprobit <- function(formula, model, data, M, ...) {
  require(VGAM)
  mf <- match.call(expand.dots = TRUE)
  if(is.null(mf$ones))
    mf$ones <- 3
  mf[[1]] <- as.name("vglm")
  mf$family <- as.name("bprobit")
  tmp <- cmvglm(mf$formula, mf$model, mf$equal, mf$zeros, mf$ones, 3)
  mf$formula <- tmp$formula 
  mf$constraints <- tmp$constraints
  bprobit <<- function() binom2.rho(zero=NULL)
  mf$model <- mf$equal <- mf$ones <- mf$zeros <- mf$M <- NULL
  as.call(mf)
}
zelig2exp <- function(formula, model, data, M, ...) {
  mf <- match.call(expand.dots = TRUE)
  mf$model <- mf$M <- NULL
  require(survival)
  mf[[1]] <- as.name("survreg")
  mf$dist <- as.character("exponential")
  as.call(mf)
}
zelig2gamma <- function(formula, model, data, M, ...) {
  mf <-  match.call(expand.dots = TRUE)
  mf$model <- mf$M <- NULL
  mf[[1]] <- as.name("glm")
  mf$family <- as.name("Gamma")
  as.call(mf)
}
zelig2logit <- function(formula, model, data, M, ...) {
  mf <- match.call(expand.dots = TRUE)
  mf$model <- mf$M <- NULL
  mf[[1]] <- as.name("glm")
  mf$family <- as.name("binomial")
  as.call(mf)
}
zelig2lognorm <- function(formula, model, data, M, ...) {
  mf <- match.call(expand.dots = TRUE)
  mf$model <- mf$M <- NULL
  require(survival)
  mf[[1]] <- as.name("survreg")
  mf$dist <- as.character("lognormal")
  as.call(mf)
}
zelig2ls <- function(formula, model, data, M, ...) {
  mf <- match.call(expand.dots = TRUE)
  mf$model <- mf$M <- NULL
  mf[[1]] <- as.name("lm")
  as.call(mf)
}
zelig2mlogit <- function(formula, model, data, M, ...) {
  require(VGAM)
  mf <- match.call(expand.dots = TRUE)
  mf[[1]] <- as.name("vglm") 
  mf$family <- as.name("multinomial")
  if (mf$M == 1)
    ndim <- length(unique(na.omit((eval(mf$formula[[2]], mf$data))))) - 1
  if (mf$M > 1) {
   ntmp <- array()
   for (i in 1:mf$M)  
     ntmp[i] <- length(unique(na.omit((eval(mf$formula[[2]], mf$data[[i]]))))) $
   ndim <- max(ntmp)
  }
  tmp <- cmvglm(mf$formula, mf$model, mf$equal, mf$zeros, mf$ones, ndim)
  mf$formula <- tmp$formula  
  mf$constraints <- tmp$constraints
  mf$model <- mf$equal <- mf$ones <- mf$zeros <- mf$M <- NULL
  as.call(mf)
}
zelig2mloglm <- function(formula, model, data, M, ...) {
  mf$model <- mf$M <- mf$... <- NULL
  require(nnet)
  mf[[1]] <- as.name("multinom")
  mf$Hess <- TRUE
  as.call(mf)
}
zelig2negbin <- function(formula, model, data, M, ...) {
  mf <- match.call(expand.dots = TRUE)
  mf$model <- mf$M <- NULL
  mf[[1]] <- as.name("glm.nb")
  as.call(mf)
}
zelig2normal <- function(formula, model, data, M, ...) {
  mf <- match.call(expand.dots = TRUE)
  mf$model <- mf$M <- NULL
  mf[[1]] <- as.name("glm")
  mf$family <- as.name("gaussian")
  as.call(mf)
}
zelig2ologit <- function(formula, model, data, M, ...) {
  mf <- match.call(expand.dots = TRUE)
  mf$model <- mf$M <- NULL
  mf[[1]] <- as.name("polr")
  mf$Hess <- TRUE
  as.call(mf)
}
zelig2oprobit <- function(formula, model, data, M, ...) {
  require(VGAM)
  mf <- match.call(expand.dots = TRUE)
  if(is.null(mf$equal)) 
    ordinal.probit <<- function() cumulative(link=probit, parallel=T)
  else
    ordinal.probit <<- function() cumulative(link=probit)
  mf[[1]] <- as.name("vglm")
  mf$family <- as.name("ordinal.probit")
  if (mf$M == 1)
    ndim <- length(unique(na.omit((eval(mf$formula[[2]], mf$data))))) - 1
  if (mf$M > 1) {
    ntmp <- array()
    for (i in 1:mf$M) 
      ntmp[i] <- length(unique(na.omit((eval(mf$formula[[2]], mf$data[[i]]))))) - 1
    ndim <- max(ntmp)
  }
  tmp <- cmvglm(mf$formula, mf$model, mf$equal, mf$zeros, mf$ones, ndim)
  mf$formula <- tmp$formula
  mf$constraints <- tmp$constraints
  mf$model <- mf$equal <- mf$ones <- mf$zeros <- mf$M <- NULL
  as.call(mf)
}
zelig2poisson <- function(formula, model, data, M, ...) {
  mf <- match.call(expand.dots = TRUE)
  mf$model <- mf$M <- NULL
  mf[[1]] <- as.name("glm")
  mf$family <- as.name("poisson")
  as.call(mf)
}
zelig2probit <- function(formula, model, data, M, ...) {
  mf <- match.call(expand.dots = TRUE)
  mf$model <- mf$M <- NULL
  mf[[1]] <- as.name("glm")
  mf$family <- as.name("binomial.probit")
  binomial.probit <<- function() binomial(link="probit")
  as.call(mf)
}
zelig2relogit <- function(formula, model, data, M, ...) {
  mf <- match.call(expand.dots = TRUE)
  mf$model <- mf$M <- NULL
  if (is.null(mf$bias.correct))
    mf$bias.correct <- TRUE
  mf[[1]] <- as.name("relogit")
  as.call(mf)
}
zelig2weibull <- function(formula, model, data, M, ...) {
  mf <- match.call(expand.dots = TRUE)
  mf$model <- mf$M <- NULL
  require(survival)
  mf[[1]] <- as.name("survreg")
  mf$dist <- as.character("weibull")
  as.call(mf)
}
