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]]
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
aebingham2 at 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.naimi at mcgill.ca]
Sent: Tuesday, January 20, 2015 7:52 AM
To: r-help at 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]]
Hi Allen,
Thanks for your response. Regarding the issue that Portnoy and PH estimators may
not compute parameters in the tail, I believe this is dependent on the amount of
censoring (e.g., in the extreme case, if all observations are censored in the
upper tail, one cannot obtain estimates in the absence of parametric
assumptions). Furthermore, it could not explain the problem I encountered, since
other fitting methods did not provide what I needed, nor was able to obtain
parameter estimates at the desired quantile (0.9th) in the absence of censoring.
Dr Koenker offered a simple solution, which was to replace the last line with:
summary(q2, taus = 1:9/10)
instead of requesting the desired taus in the crq() function.
Regards,
Ashley
Ashley I Naimi, PhD
Assistant Professor
Department Obstetrics and Gynecology
McGill University
???????????
web: www.ashleyisaacnaimi.com<http://www.ashleyisaacnaimi.com>
email: ashley.naimi at mcgill.ca<mailto:ashley.naimi at mcgill.ca>
phone: 514.574.8946
fax: 514.834.1678
-----------------------
On Feb 6, 2015, at 2:51 PM, Allen Bingham <aebingham2 at
gmail.com<mailto:aebingham2 at gmail.com>> wrote:
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
aebingham2 at gmail.com<mailto:aebingham2 at 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.naimi at mcgill.ca]
Sent: Tuesday, January 20, 2015 7:52 AM
To: r-help at 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]]
[[alternative HTML version deleted]]