From: Duncan Murdoch <murdoch_at_stats.uwo.ca>

Date: Mon, 24 Mar 2008 09:11:16 -0400

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 Mon 24 Mar 2008 - 13:22:57 GMT

Date: Mon, 24 Mar 2008 09:11:16 -0400

On 3/24/2008 9:03 AM, zhijie zhang wrote:

> Dear Murdoch,

*> Thanks very much for your rapid response. It helps me greatly.
**> If i want to write the model formula according to the estimets from
**> R, Is it correct for the below formula? I'm not very sure about it, and
**> i also hope that you can recommend a good book on this area. I want to
**> learn it systematically. Thanks a million.
**> Logit/p/=12.10+5.815*x+6.654*elevation-6.755*elevation^2
**> -1.291*distance_1 -10.348*distance2-3.53*distance3
**> _ -6.828*_y1 _ -4.426*_y2 _ -11.216*_y3 _
*

I would use predict() instead. What you have there doesn't look as though it uses the B-spline basis.

The reference given in the ?bs help page is a reasonable starting point, but just about any book that covers splines should handle the B-spline basis and the linear case.

Duncan Murdoch

> #R code and the estimate results--Knots for distance are 16.13 and 24,

*> respectively, and Knots for y are -0.4357 and -0.3202
**> m.glm<-glm(mark~x+poly(elevation,2)+bs(distance,degree=1,knots=c(16.13,24))
**>
**> +bs(y,degree=1,knots=c(-0.4357,-0.3202)),family=binomial(logit),data=point)
**> summary(m.glm)
**>
**> Estimate Std.
**> Error z value Pr(>|z|)
**> (Intercept) 12.104
**> 3.138 3.857 0.000115 ***
**> x 5.815
**> 1.987 2.926 0.003433 **
**> poly(elevation, 2)1 6.654
**> 4.457 1.493 0.135444
**> poly(elevation, 2)2 -6.755
**> 3.441 -1.963 0.049645 *
**> bs(distance, degree = 1, knots = c(16.13, 24))1 -1.291
**> 1.139 -1.133 0.257220
**> bs(distance, degree = 1, knots = c(16.13, 24))2 -10.348
**> 2.025 -5.110 3.22e-07 ***
**> bs(distance, degree = 1, knots = c(16.13, 24))3 -3.530
**> 3.752 -0.941 0.346850
**> bs(y, degree = 1, knots = c(-0.4357, -0.3202))1 -6.828
**> 1.836 -3.719 0.000200 ***
**> bs(y, degree = 1, knots = c(-0.4357, -0.3202))2 -4.426
**> 1.614 -2.742 0.006105 **
**> bs(y, degree = 1, knots = c(-0.4357, -0.3202))3 -11.216
**> 2.861 -3.920 8.86e-05 ***
**>
**> Thanks again.
**>
**>
**>
**> On Mon, Mar 24, 2008 at 8:08 PM, Duncan Murdoch <murdoch_at_stats.uwo.ca
**> <mailto:murdoch_at_stats.uwo.ca>> wrote:
**>
**> On 24/03/2008 7:06 AM, zhijie zhang wrote:
**> > Dear Murdoch,
**> > "Compare the predictions, not the coefficients.", this is the
**> key point.
**> > I have checked their predicted probability and their results
**> are the
**> > same.
**> > What do you mean by "You're using different bases"? Could you
**> please give
**> > me a little more info on it, so that i can go to find some materials?
**>
**> There are many ways to represent a piecewise linear function. R and
**> your SAS code have both chosen linear combinations of basis functions,
**> but you have used different basis functions. R uses the B-spline basis.
**> You have what is sometimes called the truncated power basis, but maybe
**> not commonly in the linear context. I don't know if it has a name here.
**>
**> Duncan Murdoch
**>
**> > Thanks a lot.
**> >
**> >
**> > On 3/24/08, Duncan Murdoch <murdoch_at_stats.uwo.ca
**> <mailto:murdoch_at_stats.uwo.ca>> wrote:
**> >> On 24/03/2008 5:23 AM, zhijie zhang wrote:
**> >>> Dear Rusers,
**> >>> I am now using R and SAS to fit the piecewise linear
**> functions, and
**> >> what
**> >>> surprised me is that they have a great differrent result. See
**> below.
**> >> You're using different bases. Compare the predictions, not the
**> >> coefficients.
**> >>
**> >> To see the bases that R uses, do this:
**> >>
**> >> matplot(distance, bs(distance,degree=1,knots=c(16.13,24)))
**> >> matplot(y, bs(y,degree=1,knots=c(-0.4357,-0.3202)))
**> >>
**> >> Duncan Murdoch
**> >>
**> >>> #R code--Knots for distance are 16.13 and 24, respectively, and
**> Knots
**> >> for y
**> >>> are -0.4357 and -0.3202
**> >>>
**> m.glm<-glm(mark~x+poly(elevation,2)+bs(distance,degree=1,knots=c(16.13
**> >> ,24))
**> >>> +bs(y,degree=1,knots=c(-0.4357,-0.3202
**> >>> )),family=binomial(logit),data=point)
**> >>> summary(m.glm)
**> >>>
**> >>> Coefficients:
**> >>> Estimate Std.
**> >> Error z
**> >>> value Pr(>|z|)
**> >>> (Intercept) 12.104
**> >> 3.138
**> >>> 3.857 0.000115 ***
**> >>> x 5.815
**> 1.987
**> >>> 2.926 0.003433 **
**> >>> poly(elevation, 2)1 6.654
**> 4.457
**> >>> 1.493 0.135444
**> >>> poly(elevation, 2)2 -6.755
**> >> 3.441 -
**> >>> 1.963 0.049645 *
**> >>> bs(distance, degree = 1, knots = c(16.13, 24))1 -1.291
**> >> 1.139 -
**> >>> 1.133 0.257220
**> >>> bs(distance, degree = 1, knots = c(16.13, 24))2 -10.348
**> >> 2.025 -
**> >>> 5.110 3.22e-07 ***
**> >>> bs(distance, degree = 1, knots = c(16.13, 24))3 -3.530
**> 3.752
**> >> -
**> >>> 0.941 0.346850
**> >>> bs(y, degree = 1, knots = c(-0.4357, -0.3202))1 -6.828
**> 1.836
**> >> -
**> >>> 3.719 0.000200 ***
**> >>> bs(y, degree = 1, knots = c(-0.4357, -0.3202))2 -4.426
**> 1.614
**> >> -
**> >>> 2.742 0.006105 **
**> >>> bs(y, degree = 1, knots = c(-0.4357, -0.3202))3 -11.216
**> >> 2.861 -
**> >>> 3.920 8.86e-05 ***
**> >>>
**> >>> #SAS codes
**> >>> data b;
**> >>> set a;
**> >>> if distance > 16.13 then d1=1; else d1=0;
**> >>> distance2=d1*(distance - 16.13);
**> >>> if distance > 24 then d2=1; else d2=0;
**> >>> distance3=d2*(distance - 24);
**> >>> if y>-0.4357 then dd1=1; else dd1=0;
**> >>> y2=dd1*(y+0.4357);
**> >>> if y>-0.3202 then dd2=1; else dd2=0;
**> >>> y3=dd2*(y+0.3202);
**> >>> run;
**> >>>
**> >>> proc logistic descending data=b;
**> >>> model mark =x elevation elevation*elevation distance distance2
**> >> distance3 y
**> >>> y2 y3;
**> >>> run;
**> >>>
**> >>>
**> >>> The LOGISTIC Procedure Analysis of Maximum
**> Likelihood
**> >>> Estimates
**> >>>
**> >>> Standard
**> Wald
**> >>> Parameter DF Estimate Error
**> >>> Chi-Square Pr > ChiSq
**> >>>
**> >>> Intercept 1 -2.6148 2.1445
**> >>> 1.4867
**> >>> 0.2227
**> >>> x 1 5.8146 1.9872
**> >>> 8.5615
**> >>> 0.0034
**> >>> elevation 1 0.4457 0.1506
**> >>> 8.7545
**> >>> 0.0031
**> >>> elevation*elevation 1 -0.0279 0.0142
**> >>> 3.8533
**> >>> 0.0496
**> >>> distance 1 -0.1091 0.0963
**> >>> 1.2836
**> >>> 0.2572
**> >>> distance2 1 -1.0418 0.2668
**> >>> 15.2424
**> >>> <.0001
**> >>> distance3 1 2.8633 0.7555
**> >>> 14.3625
**> >>> 0.0002
**> >>> y 1 -16.2032 4.3568
**> >>> 13.8314
**> >>> 0.0002
**> >>> y2 1 36.9974 10.3219
**> >>> 12.8476
**> >>> 0.0003
**> >>> y3 1 -58.4938 14.0279
**> >>> 17.3875
**> >>> <.0001
**> >>> Q: What is the problem? which one is correct for the piecewise
**> linear
**> >>> function?
**> >>> Thanks very much.
**> >>>
**> >>>
**> >>
**> >
**> >
**>
**>
**>
**>
**> --
**> With Kind Regards,
**>
**> oooO:::::::::
**> (..):::::::::
**> :\.(:::Oooo::
**> ::\_)::(..)::
**> :::::::)./:::
**> ::::::(_/::::
**> :::::::::::::
**> [***********************************************************************]
**> Zhi Jie,Zhang ,PHD
**> Tel:+86-21-54237149
**> Dept. of Epidemiology,School of Public Health,Fudan University
**> Address:No. 138 Yi Xue Yuan Road,Shanghai,China
**> Postcode:200032
**> Email:epistat_at_gmail.com <mailto:Email:epistat_at_gmail.com>
**> Website: www.statABC.com <http://www.statABC.com>
**> [***********************************************************************]
**> oooO:::::::::
**> (..):::::::::
**> :\.(:::Oooo::
**> ::\_)::(..)::
**> :::::::)./:::
**> ::::::(_/::::
**> :::::::::::::
*

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 Mon 24 Mar 2008 - 13:22:57 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 Mon 24 Mar 2008 - 18:30:24 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.
*