MNP version 2.6-4 added

DESCRIPTION

@@ -0,0 +1,35 @@
+Package: MNP
+Version: 2.6-4
+Date: 2013-06-09
+Title: R Package for Fitting the Multinomial Probit Model
+Author: Kosuke Imai <kimai@princeton.edu>, David A. van Dyk
+        <dvd@uci.edu>.
+Maintainer: Kosuke Imai <kimai@princeton.edu>
+Depends: R (>= 2.1), MASS, utils
+Description: MNP is a publicly available R package that fits the
+        Bayesian multinomial probit model via Markov chain Monte Carlo.
+        The multinomial probit model is often used to analyze the
+        discrete choices made by individuals recorded in survey data.
+        Examples where the multinomial probit model may be useful
+        include the analysis of product choice by consumers in market
+        research and the analysis of candidate or party choice by
+        voters in electoral studies.  The MNP software can also fit the
+        model with different choice sets for each individual, and
+        complete or partial individual choice orderings of the
+        available alternatives from the choice set. The estimation is
+        based on the efficient marginal data augmentation algorithm
+        that is developed by Imai and van Dyk (2005). ``A Bayesian
+        Analysis of the Multinomial Probit Model Using the Data
+        Augmentation,'' Journal of Econometrics, Vol. 124, No. 2
+        (February), pp. 311-334. Detailed examples are given in Imai
+        and van Dyk (2005). ``MNP: R Package for Fitting the
+        Multinomial Probit Model.'' Journal of Statistical Software,
+        Vol. 14, No. 3 (May), pp. 1-32.
+LazyLoad: yes
+LazyData: yes
+License: GPL (>= 2)
+URL: http://imai.princeton.edu/software/MNP.html
+Packaged: 2013-06-09 17:45:34 UTC; kimai
+NeedsCompilation: yes
+Repository: CRAN
+Date/Publication: 2013-06-09 20:16:42

MD5

@@ -0,0 +1,30 @@
+7c923c417c32a5d1a9b62db575670762 *DESCRIPTION
+e5dbaa21eb5db71bef6d103cc66a5896 *MNP.pdf
+d30ea84d5abbc3f5323c8f384954ee69 *NAMESPACE
+97a7c7aa155788b4a62bdd3e3c02fab0 *R/coef.mnp.R
+127b786d97aa185eac38721c23aa2f90 *R/cov.mnp.R
+3f7ce384c42415e3e5fe81e7c74143aa *R/mnp.R
+5e40e20eb515081275508f808be52c3e *R/onAttach.R
+88edaa6559267b5e85bc8bba9b291b97 *R/predict.mnp.R
+e822d2e847f5daf2922a019fa2a4c491 *R/print.mnp.R
+a6f93ff28505296a3b6ed5c67a0286ab *R/print.summary.mnp.R
+a954cf57cf2e252b7837ac4739d87eab *R/summary.mnp.R
+cfdc76deb4c1e71bdb9422bf9f2cecd9 *R/xmatrix.mnp.R
+f594a4e2db80f456c9dee49d5f9efa6c *R/ymatrix.mnp.R
+c342e6c52a792ec58aac9efa0be703ae *data/detergent.txt.gz
+bb1181c0a078518436a2a7c9ee5594a9 *data/japan.txt.gz
+fd1fa970759414da81360ddd7a5fdec2 *man/coef.mnp.Rd
+4f44ccc2c56a236ec59793da04327863 *man/cov.mnp.Rd
+bb8f62b3e695266f35656484ff16c42d *man/detergent.Rd
+cc55ae92569bbafe1ac86019ee029382 *man/japan.Rd
+a2a40e00f7159946b4d573a8759d7ccf *man/mnp.Rd
+d349cc562db39915eafa1c6eaa8ecc21 *man/predict.mnp.Rd
+70c12884b2bd1fcb93b659af297dbdae *man/summary.mnp.Rd
+c0fbcfebfb2e1bdfe1a82ddcdfe22957 *src/MNP.c
+f009e46fcf131d28ea4ead122961b7bd *src/Makevars
+ca53d9c2191f31adfffde00b86d5ca59 *src/rand.c
+e79ea40e8f33e593ba24b0fff4cdbba7 *src/rand.h
+76dbc5c2f45a4f71cb3d279a7759ae67 *src/subroutines.c
+d7c8fe15ffb32a769555c331cfb6021b *src/subroutines.h
+b1260ae4ad8649146cb69704d32625e8 *src/vector.c
+87ab6c76b848c81968161e883914b2ce *src/vector.h

NAMESPACE

@@ -0,0 +1,10 @@
+useDynLib(MNP)
+
+importFrom(MASS, mvrnorm)
+export(mnp, coef.mnp, cov.mnp, predict.mnp, summary.mnp, print.summary.mnp)
+
+S3method(coef, mnp)
+S3method(summary, mnp)
+S3method(predict, mnp)
+S3method(print, mnp)
+S3method(print, summary.mnp)

R/coef.mnp.R

@@ -0,0 +1,12 @@
+coef.mnp <- function(object, subset = NULL, ...) {
+  param <- object$param
+  p <- object$n.alt
+  n.cov <- ncol(param) - p*(p-1)/2
+  n <- nrow(param)
+  if (is.null(subset))
+    return(param[,1:n.cov])
+  else if (subset > n)
+    stop(paste("invalid input for `subset.' only", nrow(param), "draws are stored.")) 
+  else
+    return(param[subset, 1:n.cov])                 
+}

R/cov.mnp.R

@@ -0,0 +1,30 @@
+cov.mnp <- function(object, subset = NULL, ...) {
+  if (is.null(subset))
+    subset <- 1:nrow(object$param)
+  else if (max(subset) > nrow(object$param))
+    stop(paste("invalid input for `subset.' only", nrow(param), "draws are stored."))
+
+  p <- object$n.alt
+  n <- length(subset)
+  Sigma <- array(0, c(p-1, p-1, n))
+  param <- object$param
+  n.cov <- ncol(param) - p*(p-1)/2
+  cov <- param[,(n.cov+1):ncol(param)]
+  for (i in 1:n) {
+    count <- 1
+    for (j in 1:(p-1)) {
+      Sigma[j,j:(p-1),i] <- cov[subset[i],count:(count+p-j-1)]
+      count <- count + p - j
+    }
+    diag(Sigma[,,i]) <- diag(Sigma[,,i]/2)
+    Sigma[,,i] <- Sigma[,,i] + t(Sigma[,,i])
+  }
+  tmp <- list()
+  tmp[[1]] <- tmp[[2]] <- object$alt[-pmatch(object$base, object$alt)]
+  tmp[[3]] <- as.character(subset)
+  dimnames(Sigma) <- tmp
+  if (n > 1)
+    return(Sigma)
+  else
+    return(Sigma[,,1])
+}

R/mnp.R

@@ -0,0 +1,187 @@
+mnp <- function(formula, data = parent.frame(), choiceX = NULL,
+                cXnames = NULL, base = NULL, latent = FALSE,
+                invcdf = FALSE, trace = TRUE, n.draws = 5000, p.var = "Inf", 
+                p.df = n.dim+1, p.scale = 1, coef.start = 0,
+                cov.start = 1, burnin = 0, thin = 0, verbose = FALSE) {   
+  call <- match.call()
+  mf <- match.call(expand.dots = FALSE)
+  mf$choiceX <- mf$cXnames <- mf$base <- mf$n.draws <- mf$latent <-
+    mf$p.var <- mf$p.df <- mf$p.scale <- mf$coef.start <- mf$invcdf <-
+      mf$trace <- mf$cov.start <- mf$verbose <- mf$burnin <- mf$thin <- NULL   
+  mf[[1]] <- as.name("model.frame")
+  mf$na.action <- 'na.pass'
+  mf <- eval.parent(mf)
+
+  ## fix this parameter
+  p.alpha0 <- 1
+
+  ## obtaining Y
+  tmp <- ymatrix.mnp(mf, base=base, extra=TRUE, verbose=verbose)
+  Y <- tmp$Y
+  MoP <- tmp$MoP
+  lev <- tmp$lev
+  base <- tmp$base
+  p <- tmp$p
+  n.dim <- p - 1
+  if(verbose)
+    cat("\nThe base category is `", base, "'.\n\n", sep="") 
+  if (p < 3)
+    stop("The number of alternatives should be at least 3.")
+  if(verbose) 
+    cat("The total number of alternatives is ", p, ".\n\n", sep="") 
+  if(verbose) {
+    if (trace)
+      cat("The trace restriction is used instead of the diagonal restriction.\n\n")
+    else
+      cat("The diagonal restriction is used instead of the trace restriction.\n\n")
+  }
+
+  ### obtaining X
+  tmp <- xmatrix.mnp(formula, data=eval.parent(data),
+                     choiceX=call$choiceX, cXnames=cXnames, 
+                     base=base, n.dim=n.dim, lev=lev, MoP=MoP,
+                     verbose=verbose, extra=TRUE)
+  X <- tmp$X
+  coefnames <- tmp$coefnames
+  n.cov <- ncol(X) / n.dim
+
+  ## listwise deletion for X
+  na.ind <- apply(is.na(X), 1, sum)
+  if (ncol(Y) == 1)
+    na.ind <- na.ind + is.na(Y)
+  Y <- Y[na.ind==0,]
+  X <- X[na.ind==0,]
+  n.obs <- nrow(X)
+
+  if (verbose) {
+    cat("The dimension of beta is ", n.cov, ".\n\n", sep="")
+    cat("The number of observations is ", n.obs, ".\n\n", sep="")
+    if (sum(na.ind>0)>0) {
+      if (sum(na.ind>0)==1)
+        cat("The observation ", (1:length(na.ind))[na.ind>0], " is dropped due to missing values.\n\n", sep="")
+      else {
+        cat("The following ", sum(na.ind>0), " observations are dropped due to missing values:\n", sep="")
+        cat((1:length(na.ind))[na.ind>0], "\n\n")
+      }
+    }
+  } 
+
+  ## checking the prior for beta
+  p.imp <- FALSE 
+  if (p.var == Inf) {
+    p.imp <- TRUE
+    p.prec <- diag(0, n.cov)
+    if (verbose)
+      cat("Improper prior will be used for beta.\n\n")
+  }
+  else if (is.matrix(p.var)) {
+    if (ncol(p.var) != n.cov || nrow(p.var) != n.cov)
+      stop("The dimension of `p.var' should be ", n.cov, " x ", n.cov, sep="")
+    if (sum(sign(eigen(p.var)$values) < 1) > 0)
+      stop("`p.var' must be positive definite.")
+    p.prec <- solve(p.var)
+  }
+  else {
+    p.var <- diag(p.var, n.cov)
+    p.prec <- solve(p.var)
+  }
+  p.mean <- rep(0, n.cov)
+
+  ## checking prior for Sigma
+  p.df <- eval(p.df)
+  if (length(p.df) > 1)
+    stop("`p.df' must be a positive integer.")
+  if (p.df < n.dim)
+    stop(paste("`p.df' must be at least ", n.dim, ".", sep=""))
+  if (abs(as.integer(p.df) - p.df) > 0)
+    stop("`p.df' must be a positive integer.")
+  if (!is.matrix(p.scale))  
+    p.scale <- diag(p.scale, n.dim)
+  if (ncol(p.scale) != n.dim || nrow(p.scale) != n.dim)
+    stop("`p.scale' must be ", n.dim, " x ", n.dim, sep="")
+  if (sum(sign(eigen(p.scale)$values) < 1) > 0)
+    stop("`p.scale' must be positive definite.")
+  else if ((trace == FALSE) & (p.scale[1,1] != 1)) {
+    p.scale[1,1] <- 1
+    warning("p.scale[1,1] will be set to 1.")
+  }
+  Signames <- NULL
+  for(j in 1:n.dim)
+    for(k in 1:n.dim)
+      if (j<=k)
+        Signames <- c(Signames, paste(if(MoP) lev[j] else lev[j+1],
+                                      ":", if(MoP) lev[k] else lev[k+1], sep="")) 
+
+  ## checking starting values
+  if (length(coef.start) == 1)
+    coef.start <- rep(coef.start, n.cov)
+  else if (length(coef.start) != n.cov)
+    stop(paste("The dimenstion of `coef.start' must be  ",
+               n.cov, ".", sep=""))
+  if (!is.matrix(cov.start)) {
+    cov.start <- diag(n.dim)*cov.start
+    if (!trace)
+      cov.start[1,1] <- 1
+  }
+  else if (ncol(cov.start) != n.dim || nrow(cov.start) != n.dim)
+    stop("The dimension of `cov.start' must be ", n.dim, " x ", n.dim, sep="")
+  else if (sum(sign(eigen(cov.start)$values) < 1) > 0)
+    stop("`cov.start' must be a positive definite matrix.")
+  else if ((trace == FALSE) & (cov.start[1,1] != 1)) {
+    cov.start[1,1] <- 1
+    warning("cov.start[1,1] will be set to 1.")
+  }
+
+  ## checking thinnig and burnin intervals
+  if (burnin < 0)
+    stop("`burnin' should be a non-negative integer.") 
+  if (thin < 0)
+    stop("`thin' should be a non-negative integer.")
+  keep <- thin + 1
+
+  ## running the algorithm
+  if (latent)
+    n.par <- n.cov + n.dim*(n.dim+1)/2 + n.dim*n.obs
+  else
+    n.par <- n.cov + n.dim*(n.dim+1)/2
+  if(verbose)
+    cat("Starting Gibbs sampler...\n")
+  # recoding NA into -1
+  Y[is.na(Y)] <- -1
+
+  param <- .C("cMNPgibbs", as.integer(n.dim),
+              as.integer(n.cov), as.integer(n.obs), as.integer(n.draws),
+              as.double(p.mean), as.double(p.prec), as.integer(p.df),
+              as.double(p.scale*p.alpha0), as.double(X), as.integer(Y), 
+              as.double(coef.start), as.double(cov.start), 
+              as.integer(p.imp), as.integer(invcdf),
+              as.integer(burnin), as.integer(keep), as.integer(trace),
+              as.integer(verbose), as.integer(MoP), as.integer(latent),
+              pdStore = double(n.par*floor((n.draws-burnin)/keep)),
+              PACKAGE="MNP")$pdStore 
+  param <- matrix(param, ncol = n.par,
+                  nrow = floor((n.draws-burnin)/keep), byrow=TRUE)
+  if (latent) {
+    W <- array(as.vector(t(param[,(n.par-n.dim*n.obs+1):n.par])),
+               dim = c(n.dim, n.obs, floor((n.draws-burnin)/keep)),
+               dimnames = list(lev[!(lev %in% base)], rownames(Y), NULL))
+    param <- param[,1:(n.par-n.dim*n.obs)]
+    }
+  else
+    W <- NULL
+  colnames(param) <- c(coefnames, Signames)
+
+  ##recoding -1 back into NA
+  Y[Y==-1] <- NA
+
+  ## returning the object
+  res <- list(param = param, x = X, y = Y, w = W, call = call, alt = lev,
+              n.alt = p, base = base, invcdf = invcdf, trace = trace,
+              p.mean = if(p.imp) NULL else p.mean, p.var = p.var, 
+              p.df = p.df, p.scale = p.scale, burnin = burnin, thin = thin) 
+  class(res) <- "mnp"
+  return(res)
+}
+
+
+

R/onAttach.R

@@ -0,0 +1,7 @@
+".onAttach" <- function(lib, pkg) {
+  mylib <- dirname(system.file(package = pkg))
+  title <- packageDescription(pkg, lib.loc = mylib)$Title
+  ver <- packageDescription(pkg, lib.loc = mylib)$Version
+  author <- packageDescription(pkg, lib.loc = mylib)$Author
+  packageStartupMessage(pkg, ": ", title, "\nVersion: ", ver, "\nAuthors: ", author, "\n")
+}