[R] Help - lm, glm, aov results inconsistent with other stati stical package

Christoph Buser buser at stat.math.ethz.ch
Wed Mar 1 10:08:14 CET 2006


Dear Ben

Berwin is correct in his answer about different
parameterizations.

After changing the contrasts in R from treatment to sum

options(contrasts = c("contr.sum", "contr.poly" ))
test.lm<-lm(y~A*x)
summary(test.lm)

I got similar results to JMP. Be careful in doing a correct
interpretation of your coefficients, especially when you have an
interaction in your model.

Regards,

Christoph Buser

--------------------------------------------------------------
Christoph Buser <buser at stat.math.ethz.ch>
Seminar fuer Statistik, LEO C13
ETH (Federal Inst. Technology)	8092 Zurich	 SWITZERLAND
phone: x-41-44-632-4673		fax: 632-1228
http://stat.ethz.ch/~buser/
--------------------------------------------------------------




Ben Ridenhour writes:
 > Okay, I took the data to SAS and it gives me the same answer that R does.  I don't why JMP is giving me an incorrect answer, but it seems to be.  (Either that or I have made the same mistake in SAS and R.) Any ideas what JMP might be doing?
 >  
 >  Ben
 >  
 > -------------------------------
 > Benjamin Ridenhour
 > School of Biological Sciences
 > Washigton State University
 > P.O. Box 644236
 > Pullman, WA 99164-4236
 > Phone (509)335-7218
 > -------------------------------- 
 > "Nothing in biology makes sense except in the light of evolution."
 > -T. Dobzhansky
 > 
 > 
 > ----- Original Message ----
 > 
 > To: "Liaw, Andy" <andy_liaw at merck.com>; r-help at stat.math.ethz.ch
 > Sent: Tuesday, February 28, 2006 10:50:18 PM
 > Subject: Re: [R] Help - lm, glm, aov results inconsistent with other stati stical package
 > 
 > Alright, I'll try to give some sample code.
 >  
 >  # create A with 2 levels - 2 and 4
 >  > A<-c(rep(2,times=30),rep(4,times=42))     
 >  
 >  # make A a factor
 >  > A<-as.factor(A)                                          
 >  
 >    #generate 72 random x points
 >  > x<-rnorm(72)                                               
 >  
 >   #create different slopes & intercepts for A=2 and A=4
 >  #add a random error term 
 >  > y<-c(x[(1:30)]+1,2*x[(31:72)]-2)+rnorm(72,mean=0,sd=2)           
 >  
 >  #use model y~A*x for lm (or glm)
 >  > test.lm<-lm(y~A*x)
 >  
 >  #check the output
 >  > summary(test.lm)
 >  
 >  This essentially creates something like my data set and uses the same model.  In response to (1), I was just using 0/1 because these are codings for the 2 levels, correct? (i.e., when A=2 the dummy variable=0, when A=4 the dummy variable=1?).  In response to (2), yes, I do want different slopes for the two categories (that is what I am interested in testing).
 >  
 >  If I export the data created above to JMP and run what I believe to be the same model, I get a different answer for my equations :(
 >  
 >  
 > -------------------------------
 > Benjamin Ridenhour
 > School of Biological Sciences
 > Washigton State University
 > P.O. Box 644236
 > Pullman, WA 99164-4236
 > Phone (509)335-7218
 > -------------------------------- 
 > "Nothing in biology makes sense except in the light of evolution."
 > -T. Dobzhansky
 > 
 > 
 > ----- Original Message ----
 > From: "Liaw, Andy" <andy_liaw at merck.com>
 > 
 > Sent: Tuesday, February 28, 2006 5:14:57 PM
 > Subject: RE: [R] Help - lm, glm, aov results inconsistent with other stati stical package
 > 
 > 1. You have levels(A) as "2" and "4", yet you showed equations for A=0 and
 > A=1?
 > 
 > 2. y = A + X + A*X means you're allowing the different groups of A to have
 > different slopes.  Probably not what you intended.
 > 
 > 3. It's probably best to provide a small sample of the data (and R code) so
 > we know how you got what you got.
 > 
 > Andy
 > 
 > From: Ben Ridenhour
 > > 
 > > Hello,
 > >  
 > >  I 'm sure there must a be a simple explanation for what I'm 
 > > doing wrong but I am stumped.  I am a novice R user and this 
 > > has shaken my confidence in what I'm doing! I am trying to 
 > > run a simple ANCOVA using the model y~A*x, where y and x are 
 > > continuous and A has two levels.  Everything seems fine until 
 > > I compare the output with what I get from another statistical 
 > > package (JMP).  JMP has the model y=A+x+A*x (this should be 
 > > the same as what I specified to R, correct?).  In comparison 
 > > I get the line equations
 > >  
 > >  y = 7072.09-1024.94 x (for A=0) and
 > >  y = 7875.58-970.088 x (for A=1)
 > >  
 > >  from JMP.  And from R I get
 > >  
 > >  y = 6276.7-1259.8 x (for A=0) and
 > >  y = 7867.5-1150.1 x (for A=1).
 > >  
 > >  Obviously, these aren't even close to the same answer.  I've 
 > > tried this using glm(), lm(), and aov() and as expected they 
 > > all give the same answer.  If I do
 > >  
 > >  >levels(A)
 > >  [1] "2" "4"
 > >  
 > >  which are the two levels of A.  Why can't I get the same 
 > > answer from JMP as in R?  This is very disturbing to me!
 > >  
 > >  Thanks,
 > >  Ben
 > >  
 > >  
 > > -------------------------------
 > > Benjamin Ridenhour
 > > School of Biological Sciences
 > > Washigton State University
 > > P.O. Box 644236
 > > Pullman, WA 99164-4236
 > > Phone (509)335-7218
 > > -------------------------------- 
 > > "Nothing in biology makes sense except in the light of evolution."
 > > -T. Dobzhansky
 > > 
 > > 
 > > 
 > >     [[alternative HTML version deleted]]
 > > 
 > > ______________________________________________
 > > R-help at stat.math.ethz.ch mailing list
 > > https://stat.ethz.ch/mailman/listinfo/r-help
 > > PLEASE do read the posting guide! 
 > > http://www.R-project.org/posting-guide.html
 > > 
 > > 
 > 
 > 
 > ------------------------------------------------------------------------------
 > 
 > ------------------------------------------------------------------------------
 > 
 > 
 > 
 > 
 >     [[alternative HTML version deleted]]
 > 
 > ______________________________________________
 > R-help at stat.math.ethz.ch mailing list
 > https://stat.ethz.ch/mailman/listinfo/r-help
 > PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
 > 
 > 
 > 
 > 
 > 	[[alternative HTML version deleted]]
 > 
 > ______________________________________________
 > R-help at stat.math.ethz.ch mailing list
 > https://stat.ethz.ch/mailman/listinfo/r-help
 > PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html




More information about the R-help mailing list