[R] Probit regression with limited parameter space

[Ben Bolker]

Sally Luo <shali623 <at> gmail.com> writes:

> 
> Dear R helpers,
> 
> I need to estimate a probit model with box constraints placed on several of
> the model parameters.  I have the following two questions:
> 
> 1) How are the standard errors calclulated in glm
> (family=binomial(link="probit")?  I ran a typical probit model using the
> glm probit link and the nlminb function with my own coding of the
> loglikehood, separately. As nlminb does not produce the hessian matrix, I
> used hessian (numDeriv) to calculate it.  However, the standard errors
> calculated using hessian function are quite different from the ones
> generated by the glm function, although the parameter estimates are very
> close.  I was wondering what makes this difference in the estmation of
> standard errors and how this computation is carried out in glm.
> 
> 2) Does any one know how to estimate a constrained probit model in R (to be
> specific, I need to restrain the range of three parameters to [-1,1])?
> Among the optimation functions, so far nlminb and spg work for my problem,
> but neither produces a hessian matrix.  As I mentioned above, if I use
> hessian funciton and calculate standard errors manually, the standard
> errors seem not right.
> 

   I'm a little biased, but I think the bbmle package is the
easiest way to get this done -- it provides convenient wrappers
for a range of optimizers including nlminb.
   I would warn however that you should be very careful interpreting
the meaning of the Hessian matrix if some of your parameters lie
on the boundary of the feasible space ...

set.seed(101)
x <- runif(100)
p <- pnorm(1+3*x)
y <- rbinom(100,p,size=1)
d <- data.frame(x,y)
glmfit <- glm(y~x,family=binomial(link="probit"),data=d)
coef(summary(glmfit))

library(bbmle)
mle2fit <- mle2(y~dbinom(pnorm(pprobit),size=1),
             parameters=list(pprobit~x),
             start=list(pprobit=0),
             data=d)
coef(summary(mle2fit))

## now add constraints
mle2fit_constr <- mle2(y~dbinom(pnorm(pprobit),size=1),
             parameters=list(pprobit~x),
             start=list(pprobit=0),
             optimizer="nlminb",
             lower=c(2,0.15),
             data=d)

coef(summary(mle2fit_constr))