Nomogram (rms) for model with shrunk coefficients

classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

Nomogram (rms) for model with shrunk coefficients

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

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: Nomogram (rms) for model with shrunk coefficients

Frank Harrell
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

______________________________________________
[hidden email] 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 Harrell
Department of Biostatistics, Vanderbilt University
Reply | Threaded
Open this post in threaded view
|

Re: Nomogram (rms) for model with shrunk coefficients

David Winsemius
In reply to this post by Sander

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]]
>
> ______________________________________________
> [hidden email] 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

______________________________________________
[hidden email] 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.