R help - Jun 2013 - Nomogram (rms) for model with shrunk coefficients

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]

 ##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]]

You are using an informal shrinkage method.  It is much better to use 
penalized maximum likelihood estimation, built in to lrm.

If you really want to go the informal route, compute the linear 
predictor from your final estimates and use ols( ) to predict that from 
the component variables (you'll get an R^2 of 1.0 so specify sigma= to 
ols) and use nomogram on the result.

Frank


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

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