Take the 2-minute tour ×
Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It's 100% free, no registration required.

I have an outcome with right censoring like this:

y<-c(rep(2.83,3),
     rep(3.17,4),
     rep(3.83,4),
     rep(4.17,5),
     rep(4.83,8),
     rep(5.5,3),
     rep(7.17,5),
     rep(8.17,7),
     rep(8.83,12),
     rep(9.5, 12),
     rep(9.83,17),
     rep(10.17,30),
     rep(10.50,100)) 

where y=10.5 are right censoring values. Then, I would try to use quantreg::crq to fit a censored quantile regression model and start with a binary intervention variable:

set.seed(123)
require(quantreg)
yc<-rep(10.5, length(y))
treat<-rbinom(length(y), 1, 0.5)
age<-as.integer(rnorm(length(y), 50, 2))

model1

fit1<-crq(Curv(y, yc, "right")~treat, taus=(1:4)/5, , method="Powell")
Error in solve.default(x[h, ]) : 
  Lapack routine dgesv: system is exactly singular: U[2,2] = 0
Error in crq.fit.pow(X, y, cen, tau = taus[i], weights, left = left, ...) : 
  Singular basic solution generated by 'start'

model2

fit2<-crq(Curv(y, yc, "right")~treat+age, taus=(1:4)/5)
Error in solve.default(x[h, ]) : 
  Lapack routine dgesv: system is exactly singular: U[2,2] = 0
Error in crq.fit.pow(X, y, cen, tau = taus[i], weights, left = left, ...) : 
  Singular basic solution generated by 'start'

May somebody here know what's wrong of the models? Is it because of the ties in y? Is there a solution?

share|improve this question

1 Answer 1

up vote 13 down vote accepted

In such artificial data problems the default starting values for the Powell method are not very propitious. Here is what is happening: crq.fit.pow naively starts by trying to find an rq solution by ignoring the censoring. In your case, since your covariates are independent of the response and one of the covariates is binary, this is likely to yield a solution with a hard zero treatment coefficient. Then the algorithm tries to start at this solution and finds that this basic solution (the pair of observations that characterize the initial fit) both have treatment indicator 0, (or 1), and at that point, trying to solve for the starting value yields a singular linear system and you get your error.

So the problem arises from a rather nasty conspiracy of problems that have to do with your replicated data the lack of a model signal, and, frankly, a rather naïve choice of a protocol for choosing a starting value. If you really want to force R to produce an answer you can use start = "global" -- and (at least for small problems like this) crq will produce a globally optimal solution. But I suspect that the better path is to change the model somewhat.

share|improve this answer
8  
Thank you very much for joining us to answer this question! We are grateful to have such authoritative advice. –  whuber Jun 27 at 16:26

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.