Ashley,
Looking at the code for the crq function it appears that with method="Portnoy"
the tau value is not passed to the fitting function (crq.fit.por) and that
function appears (I didn't delve into the code to make sure of this ...) to
calculate the grid of tau values given the nature of the data, formula, etc. ---
--- more telling though is the following snippet I copied from the
documentation of the crq function, specifically it states (in part) [I've split
apart the pertinent section first]:
"... Both the Portnoy and Peng-Huang estimators may be unable to compute
estimates of
the conditional quantile parameters in the upper tail of distribution.
... Like the Kaplan-Meier estimator, when censoring is heavy in the upper
tail the
estimated distribution is defective and quantiles are only estimable on a
sub-interval
of (0,1). The Peng and Huang estimator can be viewed as a generalization of
the Nelson
Aalen estimator of the cumulative hazard function, and can be formulated as
a variant
of the conventional quantile regression dual problem. See Koenker (2008) for
further
details. This paper is available from the package with vignette("crq"). ...
"
... so I'm guessing that using one of the other fitting methods may be get you
what you need --- HOWEVER I'm only looking at the underlying code and
documentation for same --- I'm completely ignorant of the statistical methods
you are trying to implement --- and whether or not one of the other fitting
methods is appropriate to your analyses.
Regardless it seems to me (IMHO) that the documentation for the crq function
should indicate that the "tau" or alternatively the "taus" supplied to the crq
function are only used with a subset of the method values.
Hope this helps-Allen
______________________________________
Allen Bingham
Bingham Statistical Consulting
aebingh...@gmail.com
This message and any attachments may contain confidential or privileged
information and are only for the use of the intended recipient of this message.
If you are not the intended recipient, please notify the sender by return
email, and delete or destroy this and all copies of this message and all
attachments. Any unauthorized disclosure, use, distribution, or reproduction of
this message or any attachments without permission from the sender is
prohibited and may be unlawful.
-----Original Message-----
From: Ashley Isaac Naimi, Mr [mailto:ashley.na...@mcgill.ca]
Sent: Tuesday, January 20, 2015 7:52 AM
To: r-help@r-project.org
Subject: [R] censored quantile regression
Hi there,
I've generated the following data:
library(quantreg)
library(survival)
set.seed(789)
N <- 2000
u <- runif(N)
x1 <- rbinom(N,1,.5)
x2 <- rbinom(N,1,.5)
x1x2<-x1*x2
lambda <- 1 + 1.5*x1 + 1.5*x2 + .5*x1x2
k <- 2
y <- lambda*((-log(1-u))^(1/k));max(y)
c <- runif(N,max=15)
event = as.numeric(y<=c)
mean(event);table(event)
cens <- 1-event
table(cens);mean(cens)
time <-as.matrix(ifelse(event==1,y,c))
St<-Surv(time,event,type="right")
To which I've fit the following censored quantile regression model:
q2 <- crq(St~x1 + x2 + x1x2,tau=.9,method="Portnoy")
summary(q2)
As one can see, I'm interested in the 0.9th quantile. But summary(q2) returns
only the 20th to 80th percentiles (by 20). How can I get only the 0.9th
quantile (aka the 90th percentile)?? All I actually need is the parameter
estimate for the interaction (x1x2) for the 90th percentile. I know how to
extract it if I were able to obtain estimates at tau=0.9, but no luck. Even
though I request the 90th percentile in crq (i.e., "tau=0.9"), the summary
function keeps returning the same set of (unwanted) percentiles (20th, 40th,
60th, 80th).
I’m using R version 3.1.2 (2014-10-31) -- "Pumpkin Helmet”, Platform:
x86_64-apple-darwin13.4.0 (64-bit).
Cheers,
Ashley I Naimi
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
