[R] Nomogram (rms) for model with shrunk coefficients

[David Winsemius]

On Jun 24, 2013, at 12:00 PM, Sander van Kuijk wrote:

> Dear R-users,
> 
> I have used the nomogram function from the rms package for a logistic
> regresison model made with lrm(). Everything works perfectly (r version
> 2.15.1 on a mac). My question is this: if my final model is not the one
> created by lrm, but I internally validated the model and 'shrunk' the
> regression coefficients and computed a new intercept, how can I build a
> nomogram using that object? Please see the simplified code for details:
> 
> library(rms)
> 
> x1<-rnorm(100,1,1)
> x2<-rnorm(100,1,1)
> y<-rbinom(100,1,0.5)
> 
> d<-data.frame(x1,x2,y)
> 
> attach(d)
> 
> ddist<-datadist(d)
> options(datadist='ddist')
> 
> model<-lrm(y~x1+x2, x=TRUE, y=TRUE, data=d)
> plot(nomogram(model))
> ##Nomogram is printed, as expected
> 
> ##Now the model is internally validated, and regression coefficients are
> penalized
> bootstrap<-validate(model, bw=FALSE, B=100)
> shrinkage<-round(bootstrap[4,5],2)
> final<-round(model$coef*shrinkage, 3)
> final.lp<-cbind(model$x)%*%final[-1]
> final["Intercept"]<-round(lrm.fit(y=d$y, offset=final.lp)$coef,3)
> final.lp<-final[1]+model$x%*%final[-1]
> 

I cannot speak to this with authority...just another interested user. My seat-of-the-pants notion is that the Intercept estimate should move toward the unadjusted odds for the outcome in the population while the coefficients should move toward the Null value. It's not clear to me that your code is accomplshing that goal for the intercept. 

When I replace model$coefficients with final.lp (which has the same overall structure) the change in the plotted version of nomgram shifts too radically to be supportive of the notion that you have properly calculated these values:

> str(final)
 Named num [1:3] 0.015 0.016 0.009
 - attr(*, "names")= chr [1:3] "Intercept" "x1" "x2"

> str(model$coefficients)
 Named num [1:3] -0.2045 0.1603 0.0889
 - attr(*, "names")= chr [1:3] "Intercept" "x1" "x2"

It appears to me that the "Intercept" value has been shifted too far. It should have stayed near zero. Nonetheless, the nomogram function does not throw an error with this code:

> mod2 <- model
> mod2$coefficients <- final
> plot(nomogram(mod2))

But the linear predictor scale is radically shifted. And the "Total Points" scale was compressed inappropriately.
-- 
David.


> ##The object 'final' now contains all model parameters, yet in a different
> fashion than the lrm-created model. Does anyone have an idea how to make
> this compatible again with nomogram()?
> 
> Thanks in advance for any thoughts on it.
> 
> Best wishes,
> Sander
> 
> 	[[alternative HTML version deleted]]
> 
> ______________________________________________
> 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.

David Winsemius
Alameda, CA, USA