[R] can not use plot.Predict {rms} reproduce figure 7.8 from Regression Modeling Strategies (http://biostat.mc.vanderbilt.edu/wiki/pub/Main/RmS/course2.pdf)
David Winsemius
dwinsemius at comcast.net
Wed Mar 14 22:52:47 CET 2012
On Mar 14, 2012, at 4:09 PM, John Smith wrote:
> With most current version of R and RMS, the 4 curves are drew in
> 4 separate panels. Can anyone show me how can I get the figure
> exactly like
> the figure 7.8 in *Regression Modeling Strategies* (
> http://biostat.mc.vanderbilt.edu/wiki/pub/Main/RmS/course2.pdf)
In case anyone else is scratching their head wondering what connection
there might be between figure 7.8 in that document (which happens to
be a nomogram), they should realize that the questioner is most
probably asking about 7.8 in the RMS *book* and that there is no such
Figure in chapter 7 of the pdf document . (Nor, for that matter, does
a search for "democrat" bring up any hits in the document so I don't
think the code comes from there either.)
I suspect the data is here:
http://biostat.mc.vanderbilt.edu/wiki/pub/Main/DataSets/counties.sav
--
David.
>
>
>
>
> On Tue, May 17, 2011 at 4:04 PM, John Smith <zmring at gmail.com> wrote:
>
>> Dear R-users,
>>
>> I am using R 2.13.0 and rms 3.3-0 , but can not reproduce figure
>> 7.8 of
>> the handouts *Regression Modeling Strategies* (
>> http://biostat.mc.vanderbilt.edu/wiki/pub/Main/RmS/course2.pdf) by
>> the
>> following code. Could any one help me figure out how to solve this?
>>
>>
>> setwd('C:/Rharrell')
>> require(rms)
>> load('data/counties.sav')
>>
>> older <- counties$age6574 + counties$age75
>> label(older) <- '% age >= 65, 1990'
>> pdensity <- logb(counties$pop.density+1, 10)
>> label(pdensity) <- 'log 10 of 1992 pop per 1990 miles^2'
>> counties <- cbind(counties, older, pdensity) # add 2 vars. not in
>> data
>> frame
>> dd <- datadist(counties)
>> options(datadist='dd')
>>
>> f <- ols(democrat ~ rcs(pdensity,4) + rcs(pop.change,3) +
>> rcs(older,3) + crime + rcs(college,5) + rcs(income,4) +
>> rcs(college,5) %ia% rcs(income,4) +
>> rcs(farm,3) + rcs(white,5) + rcs(turnout,3), data=counties)
>>
>> incomes <- seq(22900, 32800, length=4)
>> show.pts <- function(college.pts, income.pt)
>> {
>> s <- abs(income - income.pt) < 1650
>> # Compute 10th smallest and
>> 10th
>> largest % college
>> # educated in counties with
>> median
>> family income within
>> # $1650 of the target income
>> x <- college[s]
>> x <- sort(x[!is.na(x)])
>> n <- length(x)
>> low <- x[10]; high <- x[n-9]
>> college.pts >= low & college.pts <= high
>> }
>> windows()
>> plot(Predict(f, college, income=incomes, conf.int=FALSE),
>> xlim=c(0,35), ylim=c(30,55), lty=1, lwd=c(.25,1.5,3.5,6),
>> col=c(1,1,2,2), perim=show.pts)
>>
>>
>
David Winsemius, MD
West Hartford, CT
More information about the R-help
mailing list