Re: [R] Fitting a polynomial using lrm from the Design library

From: David Winsemius <dwinsemius_at_comcast.net>
Date: Fri, 18 Jun 2010 17:36:07 -0400

On Jun 18, 2010, at 5:16 PM, David Winsemius wrote:

>
> On Jun 18, 2010, at 5:13 PM, David Winsemius wrote:
>
>>
>> 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)
>
> Forgot to fix this:
> plot(mod.rcs3, PC1=NA)
>
>> # Perfectly sensible plot which includes the OR=0 line that would
>> be the theoretically ideal result.

That would be log(odds) = 0 or OR=1. I wonder how many other errors I committed?

-- 
DW

>>
>> # 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
>>
>> ______________________________________________
>> 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
>
David Winsemius, MD West Hartford, CT ______________________________________________ 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.
Received on Fri 18 Jun 2010 - 21:40:03 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Fri 18 Jun 2010 - 22:30:34 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive