Re: [R] BBH2 and FrF2 packages
Hi Dennis, Thank you for your answer. Are you sure that the following commands are working? library(FrF2) plan - FrF2(8, 3, factor.names=c(T,C,K), default.level=c(-,+), randomize = FALSE) y - c(60,72,54,68,52,83,45,80) plan - add.response(plan, y) MEPlot(plan) IAPlot(plan) I get the following errors: library(FrF2) Loading required package: BsMD Loading required package: scatterplot3d Loading required package: igraph Loading required package: sfsmisc Loading required package: DoE.base Loading required package: relimp Loading Tcl/Tk interface ... done Loading required package: tcltk Loading required package: vcd Loading required package: MASS Loading required package: grid Loading required package: colorspace Attaching package: 'DoE.base' The following object(s) are masked from 'package:relimp': showData The following object(s) are masked from 'package:stats': lm The following object(s) are masked from 'package:graphics': plot.design The following object(s) are masked from 'package:utils': fix Attaching package: 'FrF2' The following object(s) are masked from 'package:BsMD': DanielPlot Warning message: In fun(...) : no valid postscript previewer found; consider setting options(eps_view= ) yourself plan - FrF2(8, 3, factor.names=c(T,C,K), default.level=c(-,+), randomize = FALSE) creating full factorial with 8 runs ... y - c(60,72,54,68,52,83,45,80) plan - add.response(plan, y) MEPlot(plan) Error in MEPlot.design(plan) : The design obj must be of a type containing FrF2 or pb. IAPlot(plan) Error in IAPlot.design(plan) : The design obj must be of a type containing FrF2 or pb. Only adding a fake factor to the plan the plot commands are working *** plan - FrF2(8, 4, factor.names=c(T,C,K,Q), default.level=c(-,+), randomize = FALSE) Sincerely, Andrea B. On 25 Jun, 2010, at 10:35 AM, Dennis Murphy wrote: Hi: MEPlot, IAPlot and cubePlot come from the FrF2 package; the DanielPlot function is in both package BsMD and FrF2. Try library(FrF2) and then run your code again; it worked for me... If you check the list of functions in BHH2 under HTML help, you'll find that none of the plot functions you used below are found in that package, but they are all found under FrF2. HTH, Dennis On Thu, Jun 24, 2010 at 4:17 PM, Andrea Bernasconi DG andrea.bernasconi...@gmail.com wrote: Hi R HELP, I consider the 2^3 factorial experiment described at page 177 of the book Statistics for Experimenters: Design, Innovation, and Discovery by George E. P. Box, J. Stuart Hunter, William G. Hunter (BHH2). This example use the following data in file BHH2-Data/tab0502.dat at ftp://ftp.wiley.com/ in /sci_tech_med/statistics_experimenters/BHH2-Data.zip run T C K y 1 1 -1 -1 -1 60 2 2 1 -1 -1 72 3 3 -1 1 -1 54 4 4 1 1 -1 68 5 5 -1 -1 1 52 6 6 1 -1 1 83 7 7 -1 1 1 45 8 8 1 1 1 80 Using these data and the R BHH2 package, I was not able to reproduce the very simple results in the BHH2 book. In particular, the following solution will have no meaning since K is categorical: ( plan - lm(y ~ (T+C+K)^2, data = DATA) ) MEPlot(plan) # Main Effects IAPlot(plan) # Interactions Effects DanielPlot(plan) cubePlot(plan, T, C, K) I decided to rebuilt the data using: plan - FrF2(8, 3, factor.names=c(T,C,K), default.level=c(-,+), randomize = FALSE) ( plan - add.response(plan, y) ) giving: T C K y 1 - - - 60 2 + - - 72 3 - + - 54 4 + + - 68 5 - - + 52 6 + - + 83 7 - + + 45 8 + + + 80 class=design, type= full factorial Unfortunately the following plot commands do not work: MEPlot(plan) IAPlot(plan) DanielPlot(plan) The error is: Error in MEPlot.design(plan) : The design obj must be of a type containing FrF2 or pb. Why? If I add a fake factor to the plan the plot commands work, but the solution will have no meaning: plan - FrF2(8, 4, factor.names=c(T,C,K,Q), default.level=c(-,+), randomize = FALSE) ( plan - add.response(plan, y) ) MEPlot(plan) IAPlot(plan) DanielPlot(plan) Sincerely, Andrea B. [[alternative HTML version deleted]] __ R-help@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. Mobile +41 79 621 74 07 [[alternative HTML version deleted]] __ R-help@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.
[R] BBH2 and FrF2 packages
Hi R HELP, I consider the 2^3 factorial experiment described at page 177 of the book Statistics for Experimenters: Design, Innovation, and Discovery by George E. P. Box, J. Stuart Hunter, William G. Hunter (BHH2). This example use the following data in file BHH2-Data/tab0502.dat at ftp://ftp.wiley.com/ in /sci_tech_med/statistics_experimenters/BHH2-Data.zip run T C K y 1 1 -1 -1 -1 60 2 2 1 -1 -1 72 3 3 -1 1 -1 54 4 4 1 1 -1 68 5 5 -1 -1 1 52 6 6 1 -1 1 83 7 7 -1 1 1 45 8 8 1 1 1 80 Using these data and the R BHH2 package, I was not able to reproduce the very simple results in the BHH2 book. In particular, the following solution will have no meaning since K is categorical: ( plan - lm(y ~ (T+C+K)^2, data = DATA) ) MEPlot(plan) # Main Effects IAPlot(plan) # Interactions Effects DanielPlot(plan) cubePlot(plan, T, C, K) I decided to rebuilt the data using: plan - FrF2(8, 3, factor.names=c(T,C,K), default.level=c(-,+), randomize = FALSE) ( plan - add.response(plan, y) ) giving: T C K y 1 - - - 60 2 + - - 72 3 - + - 54 4 + + - 68 5 - - + 52 6 + - + 83 7 - + + 45 8 + + + 80 class=design, type= full factorial Unfortunately the following plot commands do not work: MEPlot(plan) IAPlot(plan) DanielPlot(plan) The error is: Error in MEPlot.design(plan) : The design obj must be of a type containing FrF2 or pb. Why? If I add a fake factor to the plan the plot commands work, but the solution will have no meaning: plan - FrF2(8, 4, factor.names=c(T,C,K,Q), default.level=c(-,+), randomize = FALSE) ( plan - add.response(plan, y) ) MEPlot(plan) IAPlot(plan) DanielPlot(plan) Sincerely, Andrea B. [[alternative HTML version deleted]] __ R-help@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.
[R] Which is the easiest (most elegant) way to force aov to treat numerical variables as categorical ?
Hi R help, Hi R help, Which is the easiest (most elegant) way to force aov to treat numerical variables as categorical ? Sincerely, Andrea Bernasconi DG PROBLEM EXAMPLE I consider the latin squares example described at page 157 of the book: Statistics for Experimenters: Design, Innovation, and Discovery by George E. P. Box, J. Stuart Hunter, William G. Hunter. This example use the data-file /BHH2-Data/tab0408.dat from ftp://ftp.wiley.com/ in /sci_tech_med/statistics_experimenters/BHH2-Data.zip. The file tab0408.dat contains following DATA: DATA driver cars additive y 1 11A 19 2 21D 23 3 31B 15 4 41C 19 5 12B 24 6 22C 24 7 32D 14 8 42A 18 9 13D 23 10 23A 19 11 33C 15 12 43B 19 13 14C 26 14 24B 30 15 34A 16 16 44D 16 Now summary( aov(MODEL, data=DATA) ) Df Sum Sq Mean Sq F value Pr(F) cars 1 12.8 12.800 0.8889 0.3680 driver 1 115.2 115.200 8. 0.0179 * additive 3 40.0 13.333 0.9259 0.4634 Residuals 10 144.0 14.400 --- Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 This results differ from book result at p 159, since cars and driver are treated as numerical variables by aov. BRUTE FORCE SOLUTION Manually transforming cars and driver into categorical variables, I obtain the correct result: DATA_AB driver cars additive y 1 D1 C1A 19 2 D2 C1D 23 3 D3 C1B 15 4 D4 C1C 19 5 D1 C2B 24 6 D2 C2C 24 7 D3 C2D 14 8 D4 C2A 18 9 D1 C3D 23 10 D2 C3A 19 11 D3 C3C 15 12 D4 C3B 19 13 D1 C4C 26 14 D2 C4B 30 15 D3 C4A 16 16 D4 C4D 16 summary( aov(MODEL, data=DATA_AB) ) Df Sum Sq Mean Sq F value Pr(F) cars 3 24 8.000 1.5 0.307174 driver 3216 72.00013.5 0.004466 ** additive 3 40 13.333 2.5 0.156490 Residuals6 32 5.333 --- Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 QUESTION Which is the easiest (most elegant) way to force driver and cars from DATA to be treated as categorical variables by aov? More generally, which is the easiest way to force aov to treat numerical variables as categorical ? Sincerely, Andrea Bernasconi DG PROBLEM EXAMPLE I consider the latin squares example described at page 157 of the book: Statistics for Experimenters: Design, Innovation, and Discovery by George E. P. Box, J. Stuart Hunter, William G. Hunter. This example use the data-file /BHH2-Data/tab0408.dat from ftp://ftp.wiley.com/ in /sci_tech_med/statistics_experimenters/BHH2-Data.zip. The file tab0408.dat contains following DATA: DATA driver cars additive y 1 11A 19 2 21D 23 3 31B 15 4 41C 19 5 12B 24 6 22C 24 7 32D 14 8 42A 18 9 13D 23 10 23A 19 11 33C 15 12 43B 19 13 14C 26 14 24B 30 15 34A 16 16 44D 16 Now summary( aov(MODEL, data=DATA) ) Df Sum Sq Mean Sq F value Pr(F) cars 1 12.8 12.800 0.8889 0.3680 driver 1 115.2 115.200 8. 0.0179 * additive 3 40.0 13.333 0.9259 0.4634 Residuals 10 144.0 14.400 --- Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 This results differ from book result at p 159, since cars and driver are treated as numerical variables by aov. BRUTE FORCE SOLUTION Manually transforming cars and driver into categorical variables, I obtain the correct result: DATA_AB driver cars additive y 1 D1 C1A 19 2 D2 C1D 23 3 D3 C1B 15 4 D4 C1C 19 5 D1 C2B 24 6 D2 C2C 24 7 D3 C2D 14 8 D4 C2A 18 9 D1 C3D 23 10 D2 C3A 19 11 D3 C3C 15 12 D4 C3B 19 13 D1 C4C 26 14 D2 C4B 30 15 D3 C4A 16 16 D4 C4D 16 summary( aov(MODEL, data=DATA_AB) ) Df Sum Sq Mean Sq F value Pr(F) cars 3 24 8.000 1.5 0.307174 driver 3216 72.00013.5 0.004466 ** additive 3 40 13.333 2.5 0.156490 Residuals6 32 5.333 --- Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 QUESTION Which is the easiest (most elegant) way
Re: [R] Which is the easiest (most elegant) way to force aov to treat numerical variables as categorical ?
I think I found the solution ! cc-factor(cars) dd-factor(driver) MODEL-y~cc+dd+additive summary(aov(MODEL,data=DATA)) On 14 Jun, 2010, at 2:52 PM, Andrea Bernasconi DG wrote: Hi R help, Hi R help, Which is the easiest (most elegant) way to force aov to treat numerical variables as categorical ? Sincerely, Andrea Bernasconi DG PROBLEM EXAMPLE I consider the latin squares example described at page 157 of the book: Statistics for Experimenters: Design, Innovation, and Discovery by George E. P. Box, J. Stuart Hunter, William G. Hunter. This example use the data-file /BHH2-Data/tab0408.dat from ftp://ftp.wiley.com/ in /sci_tech_med/statistics_experimenters/BHH2-Data.zip. The file tab0408.dat contains following DATA: DATA driver cars additive y 1 11A 19 2 21D 23 3 31B 15 4 41C 19 5 12B 24 6 22C 24 7 32D 14 8 42A 18 9 13D 23 10 23A 19 11 33C 15 12 43B 19 13 14C 26 14 24B 30 15 34A 16 16 44D 16 Now summary( aov(MODEL, data=DATA) ) Df Sum Sq Mean Sq F value Pr(F) cars 1 12.8 12.800 0.8889 0.3680 driver 1 115.2 115.200 8. 0.0179 * additive 3 40.0 13.333 0.9259 0.4634 Residuals 10 144.0 14.400 --- Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 This results differ from book result at p 159, since cars and driver are treated as numerical variables by aov. BRUTE FORCE SOLUTION Manually transforming cars and driver into categorical variables, I obtain the correct result: DATA_AB driver cars additive y 1 D1 C1A 19 2 D2 C1D 23 3 D3 C1B 15 4 D4 C1C 19 5 D1 C2B 24 6 D2 C2C 24 7 D3 C2D 14 8 D4 C2A 18 9 D1 C3D 23 10 D2 C3A 19 11 D3 C3C 15 12 D4 C3B 19 13 D1 C4C 26 14 D2 C4B 30 15 D3 C4A 16 16 D4 C4D 16 summary( aov(MODEL, data=DATA_AB) ) Df Sum Sq Mean Sq F value Pr(F) cars 3 24 8.000 1.5 0.307174 driver 3216 72.00013.5 0.004466 ** additive 3 40 13.333 2.5 0.156490 Residuals6 32 5.333 --- Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 QUESTION Which is the easiest (most elegant) way to force driver and cars from DATA to be treated as categorical variables by aov? More generally, which is the easiest way to force aov to treat numerical variables as categorical ? Sincerely, Andrea Bernasconi DG PROBLEM EXAMPLE I consider the latin squares example described at page 157 of the book: Statistics for Experimenters: Design, Innovation, and Discovery by George E. P. Box, J. Stuart Hunter, William G. Hunter. This example use the data-file /BHH2-Data/tab0408.dat from ftp://ftp.wiley.com/ in /sci_tech_med/statistics_experimenters/BHH2-Data.zip. The file tab0408.dat contains following DATA: DATA driver cars additive y 1 11A 19 2 21D 23 3 31B 15 4 41C 19 5 12B 24 6 22C 24 7 32D 14 8 42A 18 9 13D 23 10 23A 19 11 33C 15 12 43B 19 13 14C 26 14 24B 30 15 34A 16 16 44D 16 Now summary( aov(MODEL, data=DATA) ) Df Sum Sq Mean Sq F value Pr(F) cars 1 12.8 12.800 0.8889 0.3680 driver 1 115.2 115.200 8. 0.0179 * additive 3 40.0 13.333 0.9259 0.4634 Residuals 10 144.0 14.400 --- Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 This results differ from book result at p 159, since cars and driver are treated as numerical variables by aov. BRUTE FORCE SOLUTION Manually transforming cars and driver into categorical variables, I obtain the correct result: DATA_AB driver cars additive y 1 D1 C1A 19 2 D2 C1D 23 3 D3 C1B 15 4 D4 C1C 19 5 D1 C2B 24 6 D2 C2C 24 7 D3 C2D 14 8 D4 C2A 18 9 D1 C3D 23 10 D2 C3A 19 11 D3 C3C 15 12 D4 C3B 19 13 D1 C4C 26 14 D2 C4B 30 15 D3 C4A 16 16 D4 C4D 16 summary( aov(MODEL, data=DATA_AB) ) Df Sum Sq Mean Sq F
[R] What is the best way to plots surfaces in 3 dimensions?
Hi R help, What is the best way to plots surfaces in 3 dimensions? I also have the following availability problem with plot3d and scatterplot3d, and wireframe: install.packages(scatterplot3d) Warning: unable to access index for repository http://cran.ch.r-project.org/bin/macosx/leopard/contrib/2.10 Warning message: In getDependencies(pkgs, dependencies, available, lib) : package scatterplot3d is not available install.packages(plot3d) Warning: unable to access index for repository http://cran.ch.r-project.org/bin/macosx/leopard/contrib/2.10 Warning message: In getDependencies(pkgs, dependencies, available, lib) : package plot3d is not available wireframe(volcano, shade = TRUE, + aspect = c(61/87, 0.4), + light.source = c(10,0,10)) Error: could not find function wireframe Sincerely, Andrea [[alternative HTML version deleted]] __ R-help@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.
[R] How to obtain the coefficients from a summary of aov ?
Dear Madame, Dear Sir, I am able to obtain the coefficients from a 'summary' of 'lm', but NOT from a 'summary' of 'aov'. The following example shows my steps. ## Initialize rm(list = ls()) # remove (almost) everything in the working environment utils::data(npk, package=MASS) # get data model - yield ~ block + N*P*K ## Using lm npk.lm - lm(model, npk) ( s.npk.lm - summary(npk.lm) ) ... Estimate Std. Error t value Pr(|t|) (Intercept) 54.8750 0.8021 68.415 2e-16 *** block11.7125 1.3893 1.233 0.24131 block21.6792 0.8021 2.093 0.05822 . block3 -1.8229 0.5672 -3.214 0.00744 ** ... s.npk.lm$coef[block1,Pr(|t|)] # this works [1] 0.2413061 ## Using aov npk.aov - aov(model, npk) ( s.npk.aov - summary(npk.aov) ) ... Df Sum Sq Mean Sq F value Pr(F) block5 343.29 68.659 4.4467 0.015939 * N1 189.28 189.282 12.2587 0.004372 ** P1 8.40 8.402 0.5441 0.474904 ... s.npk.aov$coef[block,Pr(F)] # this does NOT works ... NULL ... How to obtain the coefficients from a 'summary' of 'aov' ? In advance, I thank you very much for your eventual answer. Sincerely, Andrea Bernasconi mobile: +41 79 621 74 07 URL: http://web.me.com/andrea.bernasconi.dg/Andrea_Bernasconi_DG_home_page/HOME.html __ R-help@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.
Re: [R] How to obtain the coefficients from a summary of aov ?
Thank you David, but how to get the value of 0.015939 present in s.npk.aov, and not given by s.npk.aov$coef[block,Pr(F)] ? On the other, the procedure to extract coefficients from a summary of lm or aov should be the same. Andrea On 21 Apr, 2010, at 3:20 PM, David Winsemius wrote: On Apr 21, 2010, at 8:37 AM, Andrea Bernasconi DG wrote: Dear Madame, Dear Sir, I am able to obtain the coefficients from a 'summary' of 'lm', but NOT from a 'summary' of 'aov'. The following example shows my steps. ## Initialize rm(list = ls()) # remove (almost) everything in the working environment @#$%^*() DON'T DO THAT luckily I left off the l when I copied and pasted but otherwise this would have trashed my workspace. utils::data(npk, package=MASS) # get data model - yield ~ block + N*P*K ## Using lm npk.lm - lm(model, npk) ( s.npk.lm - summary(npk.lm) ) ... Estimate Std. Error t value Pr(|t|) (Intercept) 54.8750 0.8021 68.415 2e-16 *** block11.7125 1.3893 1.233 0.24131 block21.6792 0.8021 2.093 0.05822 . block3 -1.8229 0.5672 -3.214 0.00744 ** ... s.npk.lm$coef[block1,Pr(|t|)] # this works [1] 0.2413061 ## Using aov npk.aov - aov(model, npk) str(npk.aov) npk.aov$coefficients (Intercept) block2 block3 block4 block5 block6 N1 51.825 3.425 6.750 -3.900 -3.500 2.325 9.850 P1 K1 N1:P1 N1:K1 P1:K1N1:P1:K1 0.417 -1.917 -3.767 -4.700 0.567 NA Or reading the help pages one might have tried, although I will admit that the differences in parametrization confounded my efforts at describing a linear combination of those results to create the simpler result offered above: ?model.tables model.tables(npk.aov, effects) Tables of effects block block 1 2 3 4 5 6 -0.850 2.575 5.900 -4.750 -4.350 1.475 N N 0 1 -2.8083 2.8083 P P 0 1 0.5917 -0.5917 K K 0 1 1.9917 -1.9917 N:P P N 0 1 0 -0.9417 0.9417 1 0.9417 -0.9417 N:K K N 0 1 0 -1.175 1.175 1 1.175 -1.175 P:K K P 01 0 0.14167 -0.14167 1 -0.14167 0.14167 ( s.npk.aov - summary(npk.aov) ) ... Df Sum Sq Mean Sq F value Pr(F) block5 343.29 68.659 4.4467 0.015939 * N1 189.28 189.282 12.2587 0.004372 ** P1 8.40 8.402 0.5441 0.474904 ... s.npk.aov$coef[block,Pr(F)] # this does NOT works ... NULL ... How to obtain the coefficients from a 'summary' of 'aov' ? In advance, I thank you very much for your eventual answer. Sincerely, Andrea Bernasconi mobile: +41 79 621 74 07 URL: http://web.me.com/andrea.bernasconi.dg/Andrea_Bernasconi_DG_home_page/HOME.html __ R-help@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. David Winsemius, MD West Hartford, CT mobile: +41 79 621 74 07 __ R-help@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.
