[R] glm

[Frank E Harrell Jr]

On 06/23/2010 06:30 AM, Samuel Okoye wrote:
> Thank you ver much.
>
> Is there is a function in R which is doing penalized cubic
> regression, say spl.plr(), that if I have weeks = 1:9 I can use
> somthing like pp<- spl.plr(weeks,c(1,3,5,7)) and for 8 and 9 will be
> linear? Is rcs() library(Design)  doing this?

rcs in Design (which you should replace with rms before long) does 
unpenalized cubic splines that are linear in the tails.  You can also 
penalize such fits using the lrm and ols functions in rms/Design.  Also 
see the pspline function in the survival package, the splines package, 
and others.

Frank


>
> Many thanks, Samuel
>
> --- On Tue, 22/6/10, Joris Meys<jorismeys at gmail.com>  wrote:
>
> From: Joris Meys<jorismeys at gmail.com> Subject: Re: [R] glm To:
> "Samuel Okoye"<samuoko at yahoo.com> Cc: r-help at r-project.org Date:
> Tuesday, 22 June, 2010, 9:50
>
> On Tue, Jun 22, 2010 at 1:00 AM, Samuel Okoye<samuoko at yahoo.com>
> wrote:
>> Hi,
>>
>> I have the following data
>>
>> data1<- data.frame(count = c(0,1,1,2,4,5,13,16,14), weeks = 1:9,
>> treat=c(rep("1mg",3),rep("5mg",3),rep("10mg",3))) and I am using
>>
>> library(splines)
>>
>> to fit
>>
>> glm.m<-
>> glm(count~weeks)+as.factor(treat),family=poisson,data=data1)
>>
>> and I am interested in predicting the count variale for the weeks
>> 10, 11 and 12 with treat 10mg and 15mg.
>
> bad luck for you.
>
> newdat<-data.frame( weeks=rep(10:12,each=2),
> treat=rep(c("5mg","10mg"),times=3) )
>
> preds<- predict(glm.m,type="response",newdata=newdat,se.fit=T)
> cbind(newdat,preds)
>
> gives as expected : Warning message: In bs(weeks, degree = 3L, knots
> = numeric(0), Boundary.knots = c(1L,  : some 'x' values beyond
> boundary knots may cause ill-conditioned bases
>
> weeks treat       fit    se.fit residual.scale 1    10   5mg
> 5.934881  5.205426              1 2    10  10mg 12.041639  9.514347
> 1 3    11   5mg  4.345165  6.924663              1 4    11  10mg
> 8.816168 15.805171              1 5    12   5mg  2.781063  8.123436
> 1 6    12  10mg  5.642667 18.221007              1
>
>
> Watch the standard errors on the predicted values. No, you shouldn't
> predict outside your data space, especially when using splines. And
> when interested in 15mg, well, you shouldn't put treatment as a
> factor to start with.
>
> Cheers Joris
>
>
>
>
> ______________________________________________ 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.


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