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

David Winsemius dwinsemius at comcast.net
Fri Jun 18 23:36:07 CEST 2010

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