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)

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

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)

John Smith-89
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]]

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

John Smith-89
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 <[hidden email]> 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]]

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

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

John Smith-89
In reply to this post by John Smith-89
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)

Thanks




On Tue, May 17, 2011 at 4:04 PM, John Smith <[hidden email]> 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]]

______________________________________________
[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: 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

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

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

John Smith-89
Sorry. I am still using the 9-11 March 2011 version of course2.pdf.



On Wed, Mar 14, 2012 at 5:52 PM, David Winsemius <[hidden email]>wrote:

>
> 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<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<http://biostat.mc.vanderbilt.edu/wiki/pub/Main/DataSets/counties.sav>
>
>
> --
>
> David.
>
>>
>>
>>
>>
>> On Tue, May 17, 2011 at 4:04 PM, John Smith <[hidden email]> 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<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
>
>

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