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

Frank E Harrell Jr f.harrell at Vanderbilt.Edu
Sat Jun 19 00:15:25 CEST 2010


On 06/18/2010 04:36 PM, David Winsemius wrote:
>
> 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?

Best to convert to the rms package - see 
http://biostat.mc.vanderbilt.edu/Rrms for differences with Design.

If using ordinary polynomials, use e.g.

mod.poly3 <- lrm(trait ~ pol(PC1, 3), data=x) # not pol not poly

Frank


-- 
Frank E Harrell Jr   Professor and Chairman        School of Medicine
                      Department of Biostatistics   Vanderbilt University



More information about the R-help mailing list