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

[David Winsemius]

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