mcmcoprobit.rd

使用R语言的马尔科夫链蒙特卡洛模拟（MCMC）源代码程序。

\name{MCMCoprobit}
\alias{MCMCoprobit}
\title{Markov Chain Monte Carlo for Ordered Probit Regression}
\description{
  This function generates a sample from the posterior distribution of
  an ordered probit regression model using the data augmentation 
  approach of Cowles (1996). The user supplies data and priors,
  and a sample from the posterior distribution is returned as an mcmc
  object, which can be subsequently analyzed with functions 
  provided in the coda package.
  }
  
\usage{
MCMCoprobit(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
   thin=1, tune = NA, tdf = 1, verbose = 0, seed = NA, beta.start = NA,
   b0 = 0, B0 = 0, a0 = 0, A0 = 0, mcmc.method = c("Cowles", "AC"), ...) }

\arguments{
    \item{formula}{Model formula.}

    \item{data}{Data frame.}

    \item{burnin}{The number of burn-in iterations for the sampler.}

    \item{mcmc}{The number of MCMC iterations for the sampler.}

    \item{thin}{The thinning interval used in the simulation.  The number of
    Gibbs iterations must be divisible by this value.}

    \item{tune}{The tuning parameter for the Metropolis-Hastings
      step. Default of NA corresponds to a choice of 0.05 divided by the
      number of categories in the response variable.}

    \item{tdf}{Degrees of freedom for the multivariate-t proposal 
     distribution when \code{mcmc.method} is set to "IndMH". Must be
     positive. }

    \item{verbose}{A switch which determines whether or not the progress of
    the sampler is printed to the screen.  If \code{verbose} is greater
    than 0 the iteration
    number, the beta vector, and the Metropolis-Hastings acceptance rate
    are printed to the screen every \code{verbose}th iteration.}
  
    \item{seed}{The seed for the random number generator.  If NA, the Mersenne
      Twister generator is used with default seed 12345; if an integer is 
      passed it is used to seed the Mersenne twister.  The user can also
      pass a list of length two to use the L'Ecuyer random number generator,
      which is suitable for parallel computation.  The first element of the
      list is the L'Ecuyer seed, which is a vector of length six or NA (if NA 
      a default seed of \code{rep(12345,6)} is used).  The second element of 
      list is a positive substream number. See the MCMCpack 
      specification for more details.}
  
    \item{beta.start}{The starting value for the \eqn{\beta}{beta} vector.   
    This can either be a scalar or a column vector with dimension equal
    to the number of betas. If this takes a scalar value, then that
    value will serve as the starting value for all of the betas. The
    default value of NA will use rescaled estimates from an ordered
    logit model.}

    \item{b0}{The prior mean of \eqn{\beta}{beta}.  This can either be a 
    scalar or a column vector with dimension equal to the number of
    betas. If this takes a scalar value, then that value will serve as
    the prior mean for all of the betas.}
    
    \item{B0}{The prior precision of \eqn{\beta}{beta}.  This can either be a 
    scalar or a square matrix with dimensions equal to the number of
    betas.  If this takes a scalar value, then that value times an
    identity matrix serves as the prior precision of \eqn{\beta}{beta}.
    Default value of 0 is equivalent to  an improper uniform prior on
    \eqn{\beta}{beta}. } 

    \item{a0}{The prior mean of \eqn{\alpha}{alpha}.  This can either be a 
    scalar or a column vector with dimension equal to the number of
    betas. If this takes a scalar value, then that value will serve as
    the prior mean for all of the betas.}
    
    \item{A0}{The prior precision of \eqn{\alpha}{alpha}.  This can either be a 
    scalar or a square matrix with dimensions equal to the number of
    betas.  If this takes a scalar value, then that value times an
    identity matrix serves as the prior precision of \eqn{\beta}{beta}.
    Default value of 0 is equivalent to  an improper uniform prior on
    \eqn{\beta}{beta}. } 

   \item{mcmc.method}{Can be set to either "Cowles" (default) or "AC" to perform
posterior sampling of cutpoints based on Cowles (1996) or Albert and Chib (2001) respectively.}
    
    \item{...}{further arguments to be passed}       
}

\value{
   An mcmc object that contains the posterior sample.  This 
   object can be summarized by functions provided by the coda package.
}

\details{
\code{MCMCoprobit} simulates from the posterior distribution of a
ordered probit regression model using data augmentation. The simulation
proper is done in compiled C++ code to maximize efficiency.  Please
consult the coda documentation for a comprehensive list of functions
that can be used to analyze the posterior sample. 
  
  The observed variable \eqn{y_i}{y_i} is ordinal with a total of \eqn{C}{C} 
  categories, with distribution
  governed by a latent variable:
  \deqn{z_i = x_i'\beta + \varepsilon_i}{z_i = x_i'beta + epsilon_i}
  The errors are assumed to be from a standard Normal distribution.  The 
  probabilities of observing each outcome is governed by this latent
  variable and \eqn{C-1}{C-1} estimable cutpoints, which are denoted
  \eqn{\gamma_c}{gamma_c}.  The probability that individual \eqn{i}{i}
  is in category \eqn{c}{c} is computed by:
  \deqn{
   \pi_{ic} = \Phi(\gamma_c - x_i'\beta) - \Phi(\gamma_{c-1} - x_i'\beta)
   }{
   pi_ic = Phi(gamma_c - x_i'beta) - Phi(gamma_(c-1) - x_i'beta)
   }
   These probabilities are used to form the multinomial distribution
   that defines the likelihoods.
   
   \code{MCMCoprobit} provides two ways to sample the cutpoints. Cowles (1996) proposes a sampling scheme that groups sampling of a latent variable with cutpoints.  In this case, for identification the first element  \eqn{\gamma_1}{gamma_1} is normalized to zero. Albert and Chib (2001) show that we can sample cutpoints indirectly without constraints by transforming cutpoints into real-valued parameters (\eqn{\alpha}{alpha}).	
}
  
\references{
  Albert, James and Siddhartha Chib. 2001. ``Sequential Ordinal Modeling with Applications to Survival Data." \emph{Biometrics.} 57: 829-836. 
  
  M. K. Cowles. 1996. ``Accelerating Monte Carlo Markov Chain Convergence for
  Cumulative-link Generalized Linear Models." \emph{Statistics and Computing.}
  6: 101-110.
     
  Valen E. Johnson and James H. Albert. 1999. \emph{Ordinal Data Modeling}. 
  Springer: New York.
      
   Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin.  2007.  
   \emph{Scythe Statistical Library 1.0.} \url{http://scythe.wustl.edu}.
   
  Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2002.
   \emph{Output Analysis and Diagnostics for MCMC (CODA)}.
   \url{http://www-fis.iarc.fr/coda/}
}


\examples{
   \dontrun{
   x1 <- rnorm(100); x2 <- rnorm(100);
   z <- 1.0 + x1*0.1 - x2*0.5 + rnorm(100);
   y <- z; y[z < 0] <- 0; y[z >= 0 & z < 1] <- 1;
   y[z >= 1 & z < 1.5] <- 2; y[z >= 1.5] <- 3;
   out1 <- MCMCoprobit(y ~ x1 + x2, tune=0.3)
   out2 <- MCMCoprobit(y ~ x1 + x2, tune=0.3, tdf=3, verbose=1000, mcmc.method="AC")
   summary(out1)
   summary(out2)
   plot(out1)
   plot(out2)
   }
}

\keyword{models}

\seealso{\code{\link[coda]{plot.mcmc}},\code{\link[coda]{summary.mcmc}}}