[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)

Frank Harrell f.harrell at vanderbilt.edu
Wed May 18 05:26:26 CEST 2011


I converted code using the Design package to be what rms expects but hadn't
bothered to re-run the code.  For the first note you sent, you can
attach(counties) (or better: use with(counties) inside the show.pts
function.  But I'm getting a warning I'll need to track down before it works
100%:

Warning messages:
1: In income - income.pt :
  longer object length is not a multiple of shorter object length
2: In income - income.pt :
  longer object length is not a multiple of shorter object length
3: In income - income.pt :
  longer object length is not a multiple of shorter object length
4: In income - income.pt :
  longer object length is not a multiple of shorter object length

For the nomogram call that needs to be re-worked.  The documentation to
nomogram in rms is pretty clear how to proceed.  Sorry about the problems. 
Note that you don't need windows().
Frank


John Smith-89 wrote:
> 
> Also I can not reproduce figure 7.11 by
> 
> f <- Newlabels(f, list(turnout='voter turnout (%)'))
> windows()
> nomogram(f, interact=list(income=incomes), turnout=seq(30,100,by=10),
>          lplabel='estimated % voting Democratic', cex.var=.8,
> cex.axis=.75)
> 
> 
> 
> 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)
>>
>>
>>
>>
>> Thanks
>>
>>
>>
> 
> 	[[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.
> 


-----
Frank Harrell
Department of Biostatistics, Vanderbilt University
--
View this message in context: http://r.789695.n4.nabble.com/can-not-use-plot-Predict-rms-reproduce-figure-7-8-from-Regression-Modeling-Strategies-http-biostat-m-tp3530561p3531488.html
Sent from the R help mailing list archive at Nabble.com.



More information about the R-help mailing list