polytomousEffects <- function(mod, newdata, ci=c("logits", "probabilities"),
    level=.95){
    
#   last modified 1 Jan 05 by J. Fox
#
# mod: a model of class "multinom" or "polr"
# newdata: a data frame with rows at which the effects are to be estimated
# ci: compute confidence intervals for the fitted probabilities using the
#     the standard errors of the logits or of the probabilities
# level: confidence level
#
# Returns a data frame with newdata plus the fitted probabilities and logits,
#   their standard errors, and the lower and upper bounds of the confidence
#   intervals for the response-category probabilities

    #   define some local functions:

    eff <- function(x0, mod, ...){
        UseMethod("eff", mod)
        }
        
    eff.multinom <- function(x0, mod, ...){
        d <- array(0, c(m, m - 1, p))
        exp.x0.B <- as.vector(exp(x0 %*% B))
        sum.exp.x0.B <- sum(exp.x0.B)
        for (j in 1:(m-1)){
           d[m, j,] <- - exp.x0.B[j]*x0
           for (jj in 1:(m-1)){
               d[j, jj,] <- if (jj != j)
                   - exp(x0 %*% (B[,jj] + B[,j]))*x0
                    else exp.x0.B[j]*(1 + sum.exp.x0.B - exp.x0.B[j])*x0
               }
           }
        d <- d/(1 + sum.exp.x0.B)^2
        V.mu <- rep(0, m)
        for (j in 1:m){
           dd <- as.vector(t(d[j,,]))
           for (s in 1:r){
               for (t in 1:r){
                   V.mu[j] <- V.mu[j] + V[s,t]*dd[s]*dd[t]
                   }
               }
           }
        mu <- exp(x0 %*% B)
        mu <- mu/(1 + sum(mu))
        mu[m] <- 1 - sum(mu)
        logits <- log(mu/(1 - mu))
        V.logits <- V.mu/(mu^2 * (1 - mu)^2)
        list(p=mu, std.err.p=sqrt(V.mu), logits=logits,
            std.error.logits=sqrt(V.logits))
        }
        
    eff.polr <- function(x0, mod, ...){
        eta0 <- x0 %*% b
        mu <- rep(0, m)
        mu[1] <- 1/(1 + exp(alpha[1] + eta0))
        for (j in 2:(m-1)){
            mu[j] <- exp(eta0)*(exp(alpha[j - 1]) - exp(alpha[j]))/
                ((1 + exp(alpha[j - 1] + eta0))*(1 + exp(alpha[j] + eta0)))
            }
        mu[m] <- 1 - sum(mu)
        d <- matrix(0, m, r)
        d[1, 1] <- - exp(alpha[1] + eta0)/(1 + exp(alpha[1] + eta0))^2
        d[1, m:r] <- - exp(alpha[1] + eta0)*x0/(1 + exp(alpha[1] + eta0))^2
        for (j in 2:(m-1)){
            d[j, j-1] <- exp(alpha[j-1] + eta0)/(1 + exp(alpha[j-1] + eta0))^2
            d[j, j]   <- - exp(alpha[j] + eta0)/(1 + exp(alpha[j] + eta0))^2
            d[j, m:r] <- exp(eta0)*(exp(alpha[j]) - exp(alpha[j-1]))*
                (exp(alpha[j-1] + alpha[j] + 2*eta0) - 1) * x0 /
                    (((1 + exp(alpha[j-1] + eta0))^2)*
                    ((1 + exp(alpha[j] + eta0))^2))
            }
        d[m, m-1] <- exp(alpha[m-1] + eta0)/(1 + exp(alpha[m-1] + eta0))^2
        d[m, m:r] <- exp(alpha[m-1] + eta0)*x0/(1 + exp(alpha[m-1] + eta0))^2
        V.mu <- rep(0, m)
        for (j in 1:m){
           dd <- d[j,]
           for (s in 1:r){
               for (t in 1:r){
                   V.mu[j] <- V.mu[j] + V[s,t]*dd[s]*dd[t]
                   }
               }
           }
        logits <- log(mu/(1 - mu))
        V.logits <- V.mu/(mu^2 * (1 - mu)^2)
        list(p=mu, std.err.p=sqrt(V.mu), logits=logits,
            std.error.logits=sqrt(V.logits))
        }
        
    logit2p <- function(logit) 1/(1 + exp(-logit))
    
    # refit model to produce 'safe' predictions when the model matrix includes
    #   terms -- e.g., poly(), bs() -- whose basis depends upon the data
    
    formula.rhs <- formula(mod)[c(1,3)]
    new <- newdata
    newdata[[as.character(formula(mod)[2])]] <- rep(mod$lev[1], nrow(newdata))
    extras <- setdiff(all.vars(formula(mod)), names(model.frame(mod)))
    X <- if (length(extras) == 0) model.frame(mod)
        else {
            if (is.null(mod$call$data))
                mod$call$data <- environment(formula(mod))
            expand.model.frame(mod, extras)
            }
    nrow.X <- nrow(X)
    data <- rbind(X[,names(newdata),drop=FALSE], newdata)
    data$wt <- rep(0, nrow(data))
    data$wt[1:nrow.X] <- 1
    mod.matrix.all <- model.matrix(formula.rhs, data=data)
    X0 <- mod.matrix.all[-(1:nrow.X),]
    resp.names <- make.names(mod$lev, unique=TRUE)
    if (inherits(mod, "multinom")){
        resp.names <- c(resp.names[-1], resp.names[1]) # make the last level
                                                       # the reference level
        mod <- multinom(formula(mod), data=data, Hess=TRUE, weights=wt)
        B <- t(coef(mod))
        V <- vcov(mod)
        m <- ncol(B) + 1
        p <- nrow(B)
        r <- p*(m - 1)
        }
    else {
        mod <- polr(formula(mod), data=data, Hess=TRUE, weights=wt)
        X0 <- X0[,-1]
        b <- coef(mod)
        p <- length(b)  # corresponds to p - 1 in the text
        alpha <- - mod$zeta  # intercepts are negatives of thresholds
        m <- length(alpha) + 1
        r <- m + p - 1
        indices <- c((p+1):r, 1:p)
        V <- vcov(mod)[indices, indices]
        for (j in 1:(m-1)){  # fix up the signs of the covariances
            V[j,] <- -V[j,]  #  for the intercepts
            V[,j] <- -V[,j]
            }
        }
        
    n <- nrow(X0)
    z <- qnorm(1 - (1 - level)/2)
    ci <- match.arg(ci)
    Lower <- Upper <- P <- Logit <- SE.P <- SE.Logit <- matrix(0, n, m)
    colnames(Lower) <-  paste("L.", resp.names, sep="")
    colnames(Upper) <-  paste("U.", resp.names, sep="")
    colnames(P) <-  paste("p.", resp.names, sep="")
    colnames(Logit) <-  paste("logit.", resp.names, sep="")
    colnames(SE.P) <-  paste("se.p.", resp.names, sep="")
    colnames(SE.Logit) <-  paste("se.logit.", resp.names, sep="")
    for (i in 1:n){
        res <- eff(X0[i,], mod) # compute effects
        P[i,] <- prob <- res$p # fitted probabilities
        SE.P[i,] <- se.p <- res$std.err.p # std. errors of fitted probs
        Logit[i,] <- logit <- res$logits # fitted logits
        SE.Logit[i,] <- se.logit <- res$std.error.logits # std. errors of logits
        if (ci == "probabilities"){ # confidence intervals
            Lower[i,] <- prob - z*se.p
            Upper[i,] <- prob + z*se.p
            }
        else{
            Lower[i,] <- logit2p(logit - z*se.logit)
            Upper[i,] <- logit2p(logit + z*se.logit)
            }
        }
    cbind(new, P, Logit, SE.P, SE.Logit, Lower, Upper)
    }