Murray Jorgensen
2002-Sep-30 04:51 UTC
[R] Decompose numerical factor into orthog. poly parts
Consider the following analysis of a class experiment done as a Latin Square:> spinner <- gl(4,4,16,label=c("Murray","Angela","Shasha","Stephen")) > order <- gl(4,1,16) > treat <- scan()1: 1 2 4 3 5: 4 3 1 2 9: 3 4 2 1 13: 2 1 3 4 17: Read 16 items> coin <- factor(treat,label=c("5c","10c","20c","50c")) > time <- scan()1: 12.0 13.3 13.0 8.9 5: 9.7 7.2 7.3 8.6 9: 9.2 9.0 7.6 7.4 13: 10.1 5.7 7.5 8.9 17: Read 16 items> diam <- c(19,24,28,33) > treat[1] 1 2 4 3 4 3 1 2 3 4 2 1 2 1 3 4> diam[treat][1] 19 24 33 28 33 28 19 24 28 33 24 19 24 19 28 33> d <- diam[treat] > > lsmod <- lm(time ~ spinner + order + poly(d,3)) > anova(lsmod)Analysis of Variance Table Response: time Df Sum Sq Mean Sq F value Pr(>F) spinner 3 39.367 13.122 9.9476 0.009594 ** order 3 7.587 2.529 1.9172 0.228001 poly(d, 3) 3 14.208 4.736 3.5900 0.085602 . Residuals 6 7.915 1.319 --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 I'm a bit disappointed: I thought that the poly(d,3) sum of squares would be decomposed into its three parts. How would I do that? (Usual apologies for mental laziness.) Dr Murray Jorgensen http://www.stats.waikato.ac.nz/Staff/maj.html Department of Statistics, University of Waikato, Hamilton, New Zealand Email: maj at waikato.ac.nz Fax +64-7 838 4155 Phone +64-7 838 4773 wk +64 7 849 6486 home Mobile 021 395 862 -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help 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-help-request at stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
John Maindonald
2002-Sep-30 05:58 UTC
[R] Decompose numerical factor into orthog. poly parts
The following perhaps ought to work: d <- ordered(diam[treat]) lsmod <- lm(time ~ spinner + order + d) anova(lsmod) Unfortunately, it gives exactly the same result as before. The following, which is inelegant, does split up the orthogonal terms: d <- ordered(diam[treat]) p3 <- model.matrix(~ -1 + d) lsmod <- lm(time ~ spinner + order + p3[,1]+p3[,2]+p3[,3]) anova(lsmod) Note that orthogonal polynomial contrasts are the default for ordered factors. John Maindonald>Consider the following analysis of a class experiment done as a Latin Square: > > > spinner <- gl(4,4,16,label=c("Murray","Angela","Shasha","Stephen")) >> order <- gl(4,1,16) >> treat <- scan() >1: 1 2 4 3 >5: 4 3 1 2 >9: 3 4 2 1 >13: 2 1 3 4 >17: >Read 16 items >> coin <- factor(treat,label=c("5c","10c","20c","50c")) >> time <- scan() >1: 12.0 13.3 13.0 8.9 >5: 9.7 7.2 7.3 8.6 >9: 9.2 9.0 7.6 7.4 >13: 10.1 5.7 7.5 8.9 >17: >Read 16 items >> diam <- c(19,24,28,33) >> treat > [1] 1 2 4 3 4 3 1 2 3 4 2 1 2 1 3 4 >> diam[treat] > [1] 19 24 33 28 33 28 19 24 28 33 24 19 24 19 28 33 >> d <- diam[treat] >> >> lsmod <- lm(time ~ spinner + order + poly(d,3)) > > anova(lsmod) >Analysis of Variance Table > >Response: time > Df Sum Sq Mean Sq F value Pr(>F) >spinner 3 39.367 13.122 9.9476 0.009594 ** >order 3 7.587 2.529 1.9172 0.228001 >poly(d, 3) 3 14.208 4.736 3.5900 0.085602 . >Residuals 6 7.915 1.319 >--- >Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 > >I'm a bit disappointed: I thought that the poly(d,3) sum of squares would >be decomposed into its three parts. How would I do that? > >(Usual apologies for mental laziness.) > > > >Dr Murray Jorgensen http://www.stats.waikato.ac.nz/Staff/maj.html >Department of Statistics, University of Waikato, Hamilton, New Zealand >Email: maj at waikato.ac.nz Fax +64-7 838 4155 >Phone +64-7 838 4773 wk +64 7 849 6486 home Mobile 021 395 862 > >-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- >r-help 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-help-request at stat.math.ethz.ch >_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._-- John Maindonald email: john.maindonald at anu.edu.au phone : +61 2 (6125)3473 fax : +61 2(6125)5549 Centre for Bioinformation Science, Room 1194, John Dedman Mathematical Sciences Building (Building 27) Australian National University, Canberra ACT 0200. -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help 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-help-request at stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Bill.Venables@cmis.csiro.au
2002-Sep-30 07:26 UTC
[R] Decompose numerical factor into orthog. poly parts
Multiple degrees of freedom terms can be split up in the anova table using lsmod <- aov(time ~ spinner + order + d) # using aov( ) makes things a bit easier. summary(lsmod, split = list(.....)) ie the splitting is done using summary.aov and the split= argument allows the split to be controlled in an orderly way. The details are just too long to fit into this margin, though... Bill Venables.> -----Original Message----- > From: John Maindonald [mailto:john.maindonald at anu.edu.au] > Sent: Monday, September 30, 2002 3:58 PM > To: Murray Jorgensen > Cc: r-help at stat.math.ethz.ch > Subject: Re: [R] Decompose numerical factor into orthog. poly parts > > The following perhaps ought to work: > > d <- ordered(diam[treat]) > lsmod <- lm(time ~ spinner + order + d) > anova(lsmod) > > Unfortunately, it gives exactly the same result as before. > The following, which is inelegant, does split up the orthogonal > terms: > > d <- ordered(diam[treat]) > p3 <- model.matrix(~ -1 + d) > lsmod <- lm(time ~ spinner + order + p3[,1]+p3[,2]+p3[,3]) > anova(lsmod) > > Note that orthogonal polynomial contrasts are the default for > ordered factors. > > John Maindonald > > >Consider the following analysis of a class experiment done as a Latin > Square: > > > > > spinner <- gl(4,4,16,label=c("Murray","Angela","Shasha","Stephen")) > >> order <- gl(4,1,16) > >> treat <- scan() > >1: 1 2 4 3 > >5: 4 3 1 2 > >9: 3 4 2 1 > >13: 2 1 3 4 > >17: > >Read 16 items > >> coin <- factor(treat,label=c("5c","10c","20c","50c")) > >> time <- scan() > >1: 12.0 13.3 13.0 8.9 > >5: 9.7 7.2 7.3 8.6 > >9: 9.2 9.0 7.6 7.4 > >13: 10.1 5.7 7.5 8.9 > >17: > >Read 16 items > >> diam <- c(19,24,28,33) > >> treat > > [1] 1 2 4 3 4 3 1 2 3 4 2 1 2 1 3 4 > >> diam[treat] > > [1] 19 24 33 28 33 28 19 24 28 33 24 19 24 19 28 33 > >> d <- diam[treat] > >> > >> lsmod <- lm(time ~ spinner + order + poly(d,3)) > > > anova(lsmod) > >Analysis of Variance Table > > > >Response: time > > Df Sum Sq Mean Sq F value Pr(>F) > >spinner 3 39.367 13.122 9.9476 0.009594 ** > >order 3 7.587 2.529 1.9172 0.228001 > >poly(d, 3) 3 14.208 4.736 3.5900 0.085602 . > >Residuals 6 7.915 1.319 > >--- > >Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 > > > >I'm a bit disappointed: I thought that the poly(d,3) sum of squares would > >be decomposed into its three parts. How would I do that? > > > >(Usual apologies for mental laziness.) > > > > > > > >Dr Murray Jorgensen http://www.stats.waikato.ac.nz/Staff/maj.html > >Department of Statistics, University of Waikato, Hamilton, New Zealand > >Email: maj at waikato.ac.nz Fax +64-7 838 4155 > >Phone +64-7 838 4773 wk +64 7 849 6486 home Mobile 021 395 862 > > > >-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- > .-.-.- > >r-help 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-help-request at stat.math.ethz.ch > >_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._ > ._._._ > > -- > John Maindonald email: john.maindonald at anu.edu.au > phone : +61 2 (6125)3473 fax : +61 2(6125)5549 > Centre for Bioinformation Science, Room 1194, > John Dedman Mathematical Sciences Building (Building 27) > Australian National University, Canberra ACT 0200. > -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-. > -.-.- > r-help 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-help-request at stat.math.ethz.ch > _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._. > _._._-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help 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-help-request at stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._