[R] Simple censored quantile regression question

Ajay Shah ajayshah at mayin.org
Mon Sep 15 18:59:39 CEST 2008


I start by doing a simple gaussian tobit by MLE:

x1 <- runif(1000)           # E() = 0.5
x2 <- runif(1000)*2         # E() = 1
x3 <- runif(1000)*4         # E() = 2
ystar <- -7  + 4*x1 + 5*x2 + rnorm(1000)           # is mean 0
y <- ystar
censored <- ystar <= 0
y[censored] <- 0
library(AER)
m <- tobit(y ~ x1 + x2, left=0, data=D)
summary(m)

Which gives:

Call:
tobit(formula = y ~ x1 + x2, left = 0, data = D)

Observations:
         Total  Left-censored     Uncensored Right-censored 
          1000            500            500              0 

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -7.19392    0.20772 -34.632   <2e-16 ***
x1           4.17418    0.15678  26.625   <2e-16 ***
x2           5.10423    0.11441  44.612   <2e-16 ***
Log(scale)   0.02672    0.03149   0.849    0.396    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Scale: 1.027 

And we're doing fine; compared with the true paraemters of (-7,4,5,1)
we've recovered (-7.2, 4.2, 5.1, 1.03).

I now try to replicate this using crq():

yc <- rep(0,1000)
m <- crq(Curv(y, yc) ~ x1 + x2, tau=0.5, method="Pow", data=D)

This gives:

Coefficients:
[1] -7.154124  4.278266  5.051516

These numerical values look okay but:

> summary(m)

tau: [1] 0.5

Coefficients:
            Value    Std. Error t value  Pr(>|t|)
(Intercept) -7.15412  0.00000       -Inf  0.00000
x1           4.27827  0.00000        Inf  0.00000
x2           5.05152  0.00000        Inf  0.00000

The standard errors are all 0.

What am I missing?


While I have your attention, could someone explain the `Surv()' object
to me? How would I use method="Por" and method="Pen" in doing a
censored quantile regression?

-- 
Ajay Shah                                      http://www.mayin.org/ajayshah  
ajayshah at mayin.org                             http://ajayshahblog.blogspot.com
<*(:-? - wizard who doesn't know the answer.



More information about the R-help mailing list