Skip to content

Possible bug (PR#633)

2 messages · agusalon@mcr.ucm.es, Peter Dalgaard

#
This is a multi-part message in MIME format.

------=_NextPart_000_0007_01C00569.70069B40
Content-Type: multipart/alternative;
	boundary="----=_NextPart_001_0008_01C00569.70069B40"


------=_NextPart_001_0008_01C00569.70069B40
Content-Type: text/plain;
	charset="iso-8859-1"
Content-Transfer-Encoding: quoted-printable

Dear Sirs:

I am reading the book Mixed-Effects Models in S and S-Plus, by Pinheiro =
and Bates,  (Springer-Verlag New York, 2000) but using the nlme library =
for R (version:3.1-7, dated: 2000-07-07), and for the Linux (Mandrake =
7.0) operating system.

Working the first example of the book (pp. 4 - 12), I reproduce the =
reported results but with one exception: the model on page 6, called =
fm2Rail.lm. The estimated coefficients that I get are not those reported =
on the book.

In the attached file (report), I am sending you the short program (which =
reproduces the one in the book) and the results for estimating the =
fm2Rail.lm model. The estimated coefficients which I get are:

    Rail^0  Rail^1       Rail^2   Rail^3   Rail^4       Rail^5
    66.5    54.303    -4.692    -2.658    -0.567    11.192
which do not seem correct.

Those reported in the book are:

    Rail2    Rail5    Rail1    Rail6       Rail3    Rail4
    31.667    50      54        82.667    84.667    96

I have tried several times and other examples but the function lm seems =
not to work correctly, and perhaps the problem is in the adaptation of =
nlme from S to R.

Would you, please, let me know what is the reason for the estimated =
coefficients which I get?

Thanking you for your attention, I am

Sincerely yours


Agustin Alonso-Rodriguez
Prof. of Econometrics
Real C. Univ. "Escorial-Maria Cristina"
28200 San Lorenzo del Escorial(Madrid)
Spain
aalonso@mcr.ucm.es



------=_NextPart_001_0008_01C00569.70069B40
Content-Type: text/html;
	charset="iso-8859-1"
Content-Transfer-Encoding: quoted-printable

<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN">
<HTML><HEAD>
<META content=3D"text/html; charset=3Diso-8859-1" =
http-equiv=3DContent-Type>
<META content=3D"MSHTML 5.00.2014.210" name=3DGENERATOR>
<STYLE></STYLE>
</HEAD>
<BODY bgColor=3D#ffffff>
<DIV><FONT size=3D2>Dear Sirs:</FONT></DIV>
<DIV>&nbsp;</DIV>
<DIV><FONT size=3D2>I am reading the book <U>Mixed-Effects Models in S =
and=20
S-Plus</U>, by Pinheiro and Bates,&nbsp; (Springer-Verlag New York, =
2000) but=20
using the nlme library for R (version:3.1-7, dated: 2000-07-07), and for =
the=20
Linux (Mandrake 7.0) operating system.</FONT></DIV>
<DIV>&nbsp;</DIV>
<DIV><FONT size=3D2>Working the first example of the book (pp. 4 - 12), =
I=20
reproduce the reported results but with one exception: the model on page =
6,=20
called fm2Rail.lm. The estimated coefficients that I get are not those =
reported=20
on the book.</FONT></DIV>
<DIV>&nbsp;</DIV>
<DIV><FONT size=3D2>In the attached file (report), I am sending you the =
short=20
program (which reproduces the one in the book) and the results for =
estimating=20
the fm2Rail.lm model. The estimated coefficients which I get =
are:</FONT></DIV>
<DIV>&nbsp;</DIV>
<DIV><FONT size=3D2>&nbsp;&nbsp;&nbsp; Rail^0&nbsp; Rail^1&nbsp;&nbsp;=20
&nbsp;&nbsp;&nbsp; Rail^2&nbsp;&nbsp; Rail^3&nbsp;&nbsp; =
Rail^4&nbsp;&nbsp;=20
&nbsp;&nbsp;&nbsp; Rail^5</FONT></DIV>
<DIV><FONT size=3D2>&nbsp;&nbsp;&nbsp; 66.5&nbsp;&nbsp;&nbsp;=20
54.303&nbsp;&nbsp;&nbsp; -4.692&nbsp;&nbsp;&nbsp; =
-2.658&nbsp;&nbsp;&nbsp;=20
-0.567&nbsp;&nbsp;&nbsp; 11.192</FONT></DIV>
<DIV><FONT size=3D2>which do not seem correct.</FONT></DIV>
<DIV>&nbsp;</DIV>
<DIV><FONT size=3D2>Those reported in the book are:</FONT></DIV>
<DIV>&nbsp;</DIV>
<DIV><FONT size=3D2>&nbsp;&nbsp;&nbsp; Rail2&nbsp;&nbsp;&nbsp;=20
Rail5&nbsp;&nbsp;&nbsp; Rail1&nbsp;&nbsp;&nbsp;=20
Rail6&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; Rail3&nbsp;&nbsp;&nbsp;=20
Rail4</FONT></DIV>
<DIV><FONT size=3D2>&nbsp;&nbsp;&nbsp; 31.667&nbsp;&nbsp;&nbsp;=20
50&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 54&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; =

82.667&nbsp;&nbsp;&nbsp; 84.667&nbsp;&nbsp;&nbsp; 96</FONT></DIV>
<DIV>&nbsp;</DIV>
<DIV><FONT size=3D2>I have tried several times and other examples but =
the function=20
<STRONG>lm </STRONG>seems not to work correctly, and perhaps the problem =
is in=20
the adaptation of nlme from S to R.</FONT></DIV>
<DIV>&nbsp;</DIV>
<DIV><FONT size=3D2>Would you, please, let me know what is the reason =
for the=20
estimated coefficients which I get?</FONT></DIV>
<DIV>&nbsp;</DIV>
<DIV><FONT size=3D2>Thanking you for your attention, I am</FONT></DIV>
<DIV>&nbsp;</DIV>
<DIV><FONT size=3D2>Sincerely yours</FONT></DIV>
<DIV>&nbsp;</DIV>
<DIV>&nbsp;</DIV>
<DIV><FONT size=3D2>Agustin Alonso-Rodriguez</FONT></DIV>
<DIV><FONT size=3D2>Prof. of Econometrics</FONT></DIV>
<DIV><FONT size=3D2>Real C. Univ. "Escorial-Maria Cristina"</FONT></DIV>
<DIV><FONT size=3D2>28200 San Lorenzo del Escorial(Madrid)</FONT></DIV>
<DIV><FONT size=3D2>Spain</FONT></DIV>
<DIV><FONT size=3D2><A=20
href=3D"mailto:aalonso@mcr.ucm.es">aalonso@mcr.ucm.es</A></FONT></DIV>
<DIV>&nbsp;</DIV>
<DIV>&nbsp;</DIV></BODY></HTML>

------=_NextPart_001_0008_01C00569.70069B40--

------=_NextPart_000_0007_01C00569.70069B40
Content-Type: application/octet-stream;
	name="report"
Content-Transfer-Encoding: quoted-printable
Content-Disposition: attachment;
	filename="report"

1.- Program=0A=
=0A=
library(nlme)=0A=
data(Rail)=0A=
Rail=0A=
=0A=
# Fixed-effects model=0A=
fm2Rail.lm <- lm(travel ~ Rail - 1, data =3D Rail)=0A=
fm2Rail.lm=0A=
summary(fm2Rail.lm)=0A=
=0A=
=0A=
2.- Results=0A=
=0A=
=0A=
R : Copyright 2000, The R Development Core Team=0A=
Version 1.1.0  (June 15, 2000)=0A=
=0A=
R is free software and comes with ABSOLUTELY NO WARRANTY.=0A=
You are welcome to redistribute it under certain conditions.=0A=
Type	"?license" or "?licence" for distribution details.=0A=
=0A=
R is a collaborative project with many contributors.=0A=
Type	"?contributors" for a list.=0A=
=0A=
Type	"demo()" for some demos, "help()" for on-line help, or=0A=
    	"help.start()" for a HTML browser interface to help.=0A=
Type	"q()" to quit R.=0A=
=0A=
Loading required package: nls =0A=
Grouped Data: travel ~ 1 | Rail=0A=
   Rail travel=0A=
1     1     55=0A=
2     1     53=0A=
3     1     54=0A=
4     2     26=0A=
5     2     37=0A=
6     2     32=0A=
7     3     78=0A=
8     3     91=0A=
9     3     85=0A=
10    4     92=0A=
11    4    100=0A=
12    4     96=0A=
13    5     49=0A=
14    5     51=0A=
15    5     50=0A=
16    6     80=0A=
17    6     85=0A=
18    6     83=0A=
=0A=
Call:=0A=
lm(formula =3D travel ~ Rail - 1, data =3D Rail)=0A=
=0A=
Coefficients:=0A=
Rail^0  Rail^1  Rail^2  Rail^3  Rail^4  Rail^5  =0A=
66.500  54.303  -4.692  -2.658  -0.567  11.192  =0A=
=0A=
=0A=
Call:=0A=
lm(formula =3D travel ~ Rail - 1, data =3D Rail)=0A=
=0A=
Residuals:=0A=
    Min      1Q  Median      3Q     Max =0A=
-6.6667 -1.0000  0.1667  1.0000  6.3333 =0A=
=0A=
Coefficients:=0A=
       Estimate Std. Error t value Pr(>|t|)    =0A=
Rail^0  66.5000     0.9477  70.169  < 2e-16 ***=0A=
Rail^1  54.3032     2.3214  23.392 2.22e-11 ***=0A=
Rail^2  -4.6917     2.3214  -2.021 0.066161 .  =0A=
Rail^3  -2.6584     2.3214  -1.145 0.274458    =0A=
Rail^4  -0.5669     2.3214  -0.244 0.811181    =0A=
Rail^5  11.1919     2.3214   4.821 0.000418 ***=0A=
---=0A=
Signif. codes:  0  `***'  0.001  `**'  0.01  `*'  0.05  `.'  0.1  ` '  1 =0A=
=0A=
Residual standard error: 4.021 on 12 degrees of freedom=0A=
Multiple R-Squared: 0.9978,	Adjusted R-squared: 0.9967 =0A=
F-statistic: 916.6 on 6 and 12 degrees of freedom,	p-value: 2.998e-15 =0A=
=0A=
------=_NextPart_000_0007_01C00569.70069B40--


-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
#
agusalon@mcr.ucm.es writes:
This is a matter of differing conventions in handling

lm(travel ~ Rail - 1, data = Rail)

when Rail is an ordered factor. R uses polynomial contrasts in the
no-intercept case as it does when the intercept is present, whereas
Splus treats the ordered factor as if it were an unordered factor. If
you do something like the below, you get your coefficients.

I suspect that this could indeed be argued to be a bug in R, but
others may know better.
Call:
lm(formula = travel ~ factor(Rail, ordered = F) - 1, data = Rail)

Coefficients:
factor(Rail, ordered = F)2  factor(Rail, ordered = F)5  
                     31.67                       50.00  
factor(Rail, ordered = F)1  factor(Rail, ordered = F)6  
                     54.00                       82.67  
factor(Rail, ordered = F)3  factor(Rail, ordered = F)4  
                     84.67                       96.00