R help - Jun 2010 - Fitting a polynomial using lrm from the Design library

Hi all,

I am looking to fit a logistic regression using the lrm function from the Design
library. I am interested in this function because I would like to obtain
"pseudo-R2" values (see
http://tolstoy.newcastle.edu.au/R/help/02b/1011.html).

Can anyone help me with the syntax? 

If I fit the model using the stats library, the code looks like this:
model <- glm(x$trait ~ x$PC1 + I((x$PC1)^2) + I((x$PC1)^3), family =
binomial)

What would be the equivalent syntax for the lrm function?

Thanks very much in advance,
-----------------------------------
Josh Banta, Ph.D
Center for Genomics and Systems Biology
New York University
100 Washington Square East
New York, NY 10003
Tel: (212) 998-8465
http://plantevolutionaryecology.org



      
	[[alternative HTML version deleted]]

On Jun 18, 2010, at 12:02 PM, Josh B wrote:
> Hi all,
>
> I am looking to fit a logistic regression using the lrm function  
> from the Design library. I am interested in this function because I  
> would like to obtain "pseudo-R2" values (see
http://tolstoy.newcastle.edu.au/R/help/02b/1011.html)
> .
>
> Can anyone help me with the syntax?
>
> If I fit the model using the stats library, the code looks like this:
> model <- glm(x$trait ~ x$PC1 + I((x$PC1)^2) + I((x$PC1)^3), family =  
> binomial)
>
> What would be the equivalent syntax for the lrm function?
Not sure if the code you gave above produces an orthogonal set, but  
perhaps this will be meaningful to some of r-help's readers (but not  
necessarily to me):

require(Design)
mod.poly3 <- lrm( trait ~ poly(PC1, 3), data=x)

This does report results, but I'm not sure how you would interpret.  
(See below for one attempt)

I think Harrell would probably recommend using restricted cubic  
splines, however.

mod.rcs3 <- lrm( trait ~ rcs(PC1, 3), data=x)

For plotting with Design/Hmisc functions, you will get better results  
with the datadist facilities.
 > ddx <- datadist(x)
 > options(datadist="ddx")
 > plot(mod3, PC1=NA)
# Perfectly sensible plot which includes the OR=0 line that would be  
the theoretically ideal result.

# Whereas plot.Design does not know how to plot the earlier result

 > plot(mod.poly3, PC1=NA)
Error in plot.Design(mod.poly3, PC1 = NA) :
   matrix or interaction factor may not be displayed

May still get meaningful results with predict:

plot(seq(-3, 2, by=.1), predict(mod.poly3, data.frame(PC1=seq(-3, 2,  
by=.1)) ) )

Bit it appears to be less satisfactory that the rcs fit, since it  
blows up at the extremes.

-- 
David.
>
> Thanks very much in advance,
> -----------------------------------
> Josh Banta, Ph.D
> Center for Genomics and Systems Biology
> New York University
> 100 Washington Square East
> New York, NY 10003
> Tel: (212) 998-8465
> http://plantevolutionaryecology.org
>
>
>
>
> 	[[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.
David Winsemius, MD
West Hartford, CT