R help - Dec 2008 - Logistic regression with rcs() and inequality constraints?

Dear guRus,

I am doing a logistic regression using restricted cubic splines via 
rcs(). However, the fitted probabilities should be nondecreasing with 
increasing predictor. Example:

predictor <- seq(1,20)
y <- c(rep(0,9),rep(1,10),0)
model <- glm(y~rcs(predictor,n.knots=3),family="binomial")
print(1/(1+exp(-predict(model))))

The last expression should be a nondecreasing sequence, as fitted 
probabilities make no sense to decrease in this case.

Is there any elegant way to constrain the model fit?

Right now, I am trying to do this by defining a "large" set of linear 
inequalities and using constrOptim() on the log-likelihood. However, 
this seems to be very sensitive to the starting values, and I am not 
really happy with the results.

BTW: yes, working with lrm() would probably be more natural, but lrm() 
breaks down on my data because of a rank-deficient information matrix.

Thank you for your time!
Stephan Kolassa

Stephan Kolassa wrote:
> Dear guRus,
> 
> I am doing a logistic regression using restricted cubic splines via 
> rcs(). However, the fitted probabilities should be nondecreasing with 
> increasing predictor. Example:
> 
> predictor <- seq(1,20)
> y <- c(rep(0,9),rep(1,10),0)
> model <- glm(y~rcs(predictor,n.knots=3),family="binomial")
> print(1/(1+exp(-predict(model))))
> 
> The last expression should be a nondecreasing sequence, as fitted 
> probabilities make no sense to decrease in this case.
You may be mixing the concepts of population parameters and sample 
estimates.  Also, the confidence intervals for the fit will include 
almost anything because of the limited sample size.
> 
> Is there any elegant way to constrain the model fit?
Not with these methods.  Bootstrap bumping can be used though.
> 
> Right now, I am trying to do this by defining a "large" set of
linear
> inequalities and using constrOptim() on the log-likelihood. However, 
> this seems to be very sensitive to the starting values, and I am not 
> really happy with the results.
> 
> BTW: yes, working with lrm() would probably be more natural, but lrm() 
> breaks down on my data because of a rank-deficient information matrix.
You can try different tolerances for singularities with lrm.

Frank
> 
> Thank you for your time!
> Stephan Kolassa
> 
> ______________________________________________
> R-help at r-project.org mailing list
> 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.
> 

-- 
Frank E Harrell Jr   Professor and Chair           School of Medicine
                      Department of Biostatistics   Vanderbilt University