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)

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 &lt;zmring at gmail.com&gt;
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]]