Errors in fitting a censored quantile regression model

Question

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?

Answer 1

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.