R help - Jun 2004 - ordered probit or logit / recursive regression

On Thu, 10 Jun 2004, VUILLEUMIER Mathieu wrote:
> I make a study in health econometrics and have a categorical dependent
> variable (take value 1-5). I would like to fit an ordered probit or
> ordered logit but i didn't find a command or package who make that.
Does
> anyone know if it's exists ?
Does polr in package MASS or lrm in package Design do what you want?

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595

> I make a study in health econometrics and have a categorical
> dependent variable (take value 1-5). I would like to fit an ordered
> probit or ordered logit but i didn't find a command or package who
> make that. Does anyone know if it's exists ?
R is very fancy. You won't get mundane things like ordered probit off
the shelf. (I will be very happy if someone will show how to use glm()
to do a vanilla probit!) But you will get markov chain monte carlo for
the ordered probit model! :-) It's cool.
> help.search("ordered probit")
MCMCoprobit(mcmcpack-0.4-8)
                        Markov chain Monte Carlo for Ordered Probit
                        Regression
> library(MCMCpack)
##
## Markov chain Monte Carlo Package (MCMCpack)
## Copyright (C) 2003 Andrew D. Martin and Kevin M. Quinn
##
Loading required package: coda 
Loading required package: MASS 
> ?MCMCoprobit
The example in this function suggests:

        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;
        posterior <- MCMCoprobit(y ~ x1 + x2, tune=0.3, mcmc=20000)
        plot(posterior)
        summary(posterior)

In our standard econometrics notation for the ordered probit model, we
use beta for the coefficients and tau for the vector of cutoffs, I
think they are using beta = (1,.1,-.5), T=100, y* = beta'x + u,
sigma(u) = 1, tau = (0, 1, 1.5). (They use 'z' for our y*).

When I run this, I get:
> summary(posterior)
Iterations = 1:19996
Thinning interval = 5 
Number of chains = 1 
Sample size per chain = 4000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

                Mean     SD Naive SE Time-series SE
(Intercept)  1.10010 0.1580 0.002497       0.003518
x1           0.09309 0.1238 0.001958       0.001953
x2          -0.53336 0.1194 0.001887       0.001943
gamma2       1.08807 0.1565 0.002474       0.004187
gamma3       1.54971 0.1818 0.002875       0.004952

2. Quantiles for each variable:

               2.5%      25%      50%     75%   97.5%
(Intercept)  0.8037  0.99397  1.09560  1.2053  1.4253
x1          -0.1503  0.01075  0.09285  0.1784  0.3393
x2          -0.7671 -0.61606 -0.53358 -0.4534 -0.2977
gamma2       0.7989  0.98142  1.08194  1.1901  1.4077
gamma3       1.2020  1.42604  1.54146  1.6618  1.9291


-- 
Ajay Shah                                                   Consultant
ajayshah at mayin.org                      Department of Economic Affairs
http://www.mayin.org/ajayshah           Ministry of Finance, New Delhi

On Sat, 12 Jun 2004, Ajay Shah wrote:
>
> R is very fancy. You won't get mundane things like ordered probit off
> the shelf. (I will be very happy if someone will show how to use glm()
> to do a vanilla probit!)

glm(y~x+z, family=binomial(probit))


Be happy,

	-thomas