Re: [R] Levene Test with R
Hi, library(car) ?levene.test At 16:36 5/7/2007, along zeng wrote: Hi All, is there Levene' test in R ? If not ,Could you give me some advice about Levene test with R? Thanks a lot! I am waiting for yours. __ R-help@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 and provide commented, minimal, self-contained, reproducible code. __ R-help@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 and provide commented, minimal, self-contained, reproducible code.
[R] help on the use of ldBand
Hi R Users, I am trying to use the ldBand package. Together with the package, I have downloaded the ld98 program (version for windows) as indicated in the help page on ldBand. I did it, but obtained an error message Error in (head + 1):length(w) : Argument NA/NaN when I copied the help examples, so it seems that a conection between R and ld98 is not well performed in my computer. Did I put ld98.exe in the wrong place? If so, where should I put it? Thanks a lot in advance, Berta Ibañez Beroiz __ R-help@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 and provide commented, minimal, self-contained, reproducible code.
[R] Fwd: Re: aov y lme
Dear Prof. Ripley and Christoph, thank you very much for your comments. You have helped me a lot. Thanks, Tomas Goicoa Dear Prof. Ripley Thank you for your email. Yes, this is of course the correct syntax to save us the extra calculation. And I forgot the lower.tail = FALSE for pf() in my example to obtain the p-value. Thank you for the corrections and Best regards, Christoph Buser Prof Brian Ripley writes: On Mon, 22 Jan 2007, Christoph Buser wrote: Dear Tomas You can produce the results in Montgomery Montgomery with lme. Please remark that you should indicate the nesting with the levels in your nested factor. Therefore I recreated your data, but used 1,...,12 for the levels of batch instead of 1,...,4. purity-c(1,-2,-2,1,-1,-3,0,4, 0,-4, 1,0, 1,0,-1,0,-2,4,0,3, -3,2,-2,2,2,-2,1,3,4,0,-1,2,0,2,2,1) suppli-factor(c(rep(1,12),rep(2,12),rep(3,12))) batch-factor(c(rep(1:4,3), rep(5:8,3), rep(9:12,3))) material-data.frame(purity,suppli,batch) As you remarked you can use aov summary(material.aov-aov(purity~suppli+suppli:batch,data=material)) Df Sum Sq Mean Sq F value Pr(F) suppli2 15.056 7.528 2.8526 0.07736 . suppli:batch 9 69.917 7.769 2.9439 0.01667 * Residuals24 63.333 2.639 --- Signif. codes: 0 $,1rx(B***$,1ry(B 0.001 $,1rx(B**$,1ry(B 0.01 $,1rx(B*$,1ry(B 0.05 $,1rx(B.$,1ry(B 0.1 $,1rx(B $,1ry(B 1 Remark that the test of suppli is not the same as in Montgomery. Here it is wrongly tested against the Residuals and I don't think so: aov() is doing the correct thing for the model specified. The aov() model you want is probably aov(purity ~ suppli + Error(suppli:batch), data=material) and this gives summary(.Last.value) Error: suppli:batch Df Sum Sq Mean Sq F value Pr(F) suppli 2 15.056 7.528 0.969 0.4158 Residuals 9 69.917 7.769 Error: Within Df Sum Sq Mean Sq F value Pr(F) Residuals 24 63.333 2.639 you should perform the calculate the test with: (Fsuppi - summary(material.aov)[[1]][1,Mean Sq]/ summary(material.aov)[[1]][2,Mean Sq]) pf(Fsuppi, df1 = 2, df2 = 9) You want the other tail. -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@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 and provide commented, minimal, self-contained, reproducible code.
[R] aov y lme
Dear R user, I am trying to reproduce the results in Montgomery D.C (2001, chap 13, example 13-1). Briefly, there are three suppliers, four batches nested within suppliers and three determinations of purity (response variable) on each batch. It is a two stage nested design, where suppliers are fixed and batches are random. y_ijk=mu+tau_i+beta_j(nested in tau_i)+epsilon_ijk Here are the data, purity-c(1,-2,-2,1, -1,-3, 0,4, 0,-4, 1, 0, 1,0,-1,0, -2,4,0,3, -3,2,-2,2, 2,-2,1,3, 4,0,-1,2, 0,2,2,1) suppli-factor(c(rep(1,12),rep(2,12),rep(3,12))) batch-factor(rep(c(1,2,3,4),9)) material-data.frame(purity,suppli,batch) If I use the function aov, I get material.aov-aov(purity~suppli+suppli:batch,data=material) summary(material.aov) Df Sum Sq Mean Sq F value Pr(F) suppli2 15.056 7.528 2.8526 0.07736 . suppli:batch 9 69.917 7.769 2.9439 0.01667 * Residuals24 63.333 2.639 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 and I can estimate the variance component for the batches as (7.769- 2.639)/3=1.71 which is the way it is done in Montgomery, D. I want to use the function lme because I would like to make a diagnosis of the model, and I think it is more appropriate. Looking at Pinheiro and Bates, I have tried the following, library(nlme) material.lme-lme(purity~suppli,random=~1|suppli/batch,data=material) VarCorr(material.lme) Variance StdDev suppli =pdLogChol(1) (Intercept) 1.563785 1.250514 batch = pdLogChol(1) (Intercept) 1.709877 1.307622 Residual2.638889 1.624466 material.lme Linear mixed-effects model fit by REML Data: material Log-restricted-likelihood: -71.42198 Fixed: purity ~ suppli (Intercept) suppli2 suppli3 -0.417 0.750 1.583 Random effects: Formula: ~1 | suppli (Intercept) StdDev:1.250514 Formula: ~1 | batch %in% suppli (Intercept) Residual StdDev:1.307622 1.624466 Number of Observations: 36 Number of Groups: suppli batch %in% suppli 312 From VarCorr I obtain the variance component 1.71, but I am not sure if this is the way to fit the model for the nested design. Here, I also have a variance component for suppli and this is a fixed factor. Can anyone give me a clue? [[alternative HTML version deleted]] __ R-help@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 and provide commented, minimal, self-contained, reproducible code.
[R] (no subject)
Dear R user, I am trying to reproduce the results in Montgomery D.C (2001, chap 13, example 13-1). Briefly, there are three suppliers, four batches nested within suppliers and three determinations of purity (response variable) on each batch. It is a two stage nested design, where suppliers are fixed and batches are random. y_ijk=mu+tau_i+beta_j(nested in tau_i)+epsilon_ijk Here are the data, purity-c(1,-2,-2,1, -1,-3, 0,4, 0,-4, 1, 0, 1,0,-1,0, -2,4,0,3, -3,2,-2,2, 2,-2,1,3, 4,0,-1,2, 0,2,2,1) suppli-factor(c(rep(1,12),rep(2,12),rep(3,12))) batch-factor(rep(c(1,2,3,4),9)) material-data.frame(purity,suppli,batch) If I use the function aov, I get material.aov-aov(purity~suppli+suppli:batch,data=material) summary(material.aov) Df Sum Sq Mean Sq F value Pr(F) suppli2 15.056 7.528 2.8526 0.07736 . suppli:batch 9 69.917 7.769 2.9439 0.01667 * Residuals24 63.333 2.639 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 and I can estimate the variance component for the batches as (7.769- 2.639)/3=1.71 which is the way it is done in Montgomery, D. I want to use the function lme because I would like to make a diagnosis of the model, and I think it is more appropriate. Looking at Pinheiro and Bates, I have tried the following, library(nlme) material.lme-lme(purity~suppli,random=~1|suppli/batch,data=material) VarCorr(material.lme) Variance StdDev suppli =pdLogChol(1) (Intercept) 1.563785 1.250514 batch = pdLogChol(1) (Intercept) 1.709877 1.307622 Residual2.638889 1.624466 material.lme Linear mixed-effects model fit by REML Data: material Log-restricted-likelihood: -71.42198 Fixed: purity ~ suppli (Intercept) suppli2 suppli3 -0.417 0.750 1.583 Random effects: Formula: ~1 | suppli (Intercept) StdDev:1.250514 Formula: ~1 | batch %in% suppli (Intercept) Residual StdDev:1.307622 1.624466 Number of Observations: 36 Number of Groups: suppli batch %in% suppli 312 From VarCorr I obtain the variance component 1.71, but I am not sure if this is the way to fit the model for the nested design. Here, I also have a variance component for suppli and this is a fixed factor. Can anyone give me a clue? [[alternative HTML version deleted]] __ R-help@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 and provide commented, minimal, self-contained, reproducible code.
[R] split-plot multiple comparisons
Dear R user, I am new with split-plot designs and I have problems with multiple comparisons. This data correspond to an split-plot experiment with two replications (bloque).(Hoshmand, 2006 pp 138). Briefly, the whole-plot factor is Nitrogen concentration (nitrogeno) and the subplot factor is the variety of corn (hibrido). The aim is to determine if major differences in yield response to nitrogen fertilization exists among 5 hibrids of corn. yield-c(130,150,170,165, 125,150,160,165, 110,140,155,150, 115,140,160,140, 115,170,160,170, 135,170,190,185, 150,160,180,200, 135,155,165,175, 130,150,175,170, 145,180,195,200) hibrido-factor(rep(c(rep(P3747,4),rep(P3732,4),rep(Mol17,4),rep(A632,4),rep(LH74,4)),2)) bloque-factor(c(rep(I,20),rep(II,20))) nitrogeno-factor(rep(c(0,70,140,210),10)) maiz-data.frame(yield,hibrido,bloque,nitrogeno) names(maiz)-c(yield,hibrido,bloque,nitrogeno) options(contrasts=c(contr.helmert,contr.poly)) I fit the model maiz.aov-aov(yield~hibrido+nitrogeno+nitrogeno*hibrido+Error(bloque/nitrogeno)) and I try the orthogonal contrasts summary(maiz.aov,split=list(hibrido=list(c1=1,c2=2,c3=3,c4=4))) Error: bloque Df Sum Sq Mean Sq F value Pr(F) Residuals 1 4100.6 4100.6 Error: bloque:nitrogeno Df Sum Sq Mean Sq F value Pr(F) nitrogeno 3 12051.9 4017.3 42.756 0.005825 ** Residuals 3 281.994.0 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Error: Within Df Sum Sq Mean Sq F valuePr(F) hibrido 4 2466.25 616.56 20.5521 3.820e-06 *** hibrido: c11 1501.56 1501.56 50.0521 2.629e-06 *** hibrido: c21 438.02 438.02 14.6007 0.001504 ** hibrido: c31 301.04 301.04 10.0347 0.005968 ** hibrido: c41 225.63 225.63 7.5208 0.014458 * hibrido:nitrogeno 12 863.75 71.98 2.3993 0.051992 . hibrido:nitrogeno: c1 3 454.69 151.56 5.0521 0.011893 * hibrido:nitrogeno: c2 3 72.40 24.13 0.8044 0.509595 hibrido:nitrogeno: c3 3 276.04 92.01 3.0671 0.058027 . hibrido:nitrogeno: c4 3 60.62 20.21 0.6736 0.580663 Residuals 16 480.00 30.00 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 I understand that the line hibrido: c1 means that the first and second levels of hibrido are statistically significant, but I am not sure about the interpretation of the lines hibrido:nitrogeno: c1. Does it means that the first and second levels of hibrido are statistically significant according to the levels of the variable nitrogeno? If this is so, how can I guess the levels of nitrogeno in which they are different? I also would like to answer questions such as the following for a particular level of nitrogeno, is the mean of the first level of hibrido different than the second (third, fourth, fifth) ? Given a level of hibrido, is its mean different deppending of the level of nitrogeno? Is the mean of yield within the first level of hibrido an nitrogeno different than the mean of yield within the second level of hibrido and the third of nitrogeno? Thank you very much, Tomás Goicoa [[alternative HTML version deleted]] __ R-help@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 and provide commented, minimal, self-contained, reproducible code.
Re: [R] effects in ANCOVA
Dear Chuck, thank you very much indeed. I was looking for that and I could not find it. Cheers Tomas __ R-help@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 and provide commented, minimal, self-contained, reproducible code.
[R] effects in ANCOVA
Dear R users, I am trying to fit the following ANCOVA model in R2.4.0 Y_ij=mu+alpha_i+beta*(X_ij-X..)+epsilon_ij Particularly I am interested in obtaining estimates for mu, and the effects alpha_i I have this data (from the book Applied Linear Statistical Models by Neter et al (1996), page 1020) y-c(38,43,24,39,38,32,36,38,31,45,27,21,33,34,28) x-c(21,34,23,26,26,29,22,29,30,28,18,16,19,25,29) xmean-x-mean(x) grupo-factor(rep(c(1,2,3),5)) datos-data.frame(y,xmean,grupo) and I have done the following modelo-lm(y~xmean+grupo,data=datos) summary(modelo) Call: lm(formula = y ~ xmean + grupo, data = datos) Residuals: Min 1Q Median 3Q Max -2.4348 -1.2739 -0.3363 1.6710 2.4869 Coefficients: Estimate Std. Error t value Pr(|t|) (Intercept) 39.8174 0.8576 46.432 5.66e-14 *** xmean 0.8986 0.1026 8.759 2.73e-06 *** grupo2 -5.0754 1.2290 -4.130 0.00167 ** grupo3 -12.9768 1.2056 -10.764 3.53e-07 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 1.873 on 11 degrees of freedom Multiple R-Squared: 0.9403, Adjusted R-squared: 0.9241 F-statistic: 57.78 on 3 and 11 DF, p-value: 5.082e-07 In the book, the estimates are From this book, mu=33.8 alpha_1=6.017 alpha_2=0.942 beta=0.899 Is it possible to obtain this estimates and their standard errors from the model I fitted? Thanks in advance Tomas Goicoa [[alternative HTML version deleted]] __ R-help@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 and provide commented, minimal, self-contained, reproducible code.
[R] Problems in batch mode
Dear R-users, I am trying to run some simulations in batch mode. In an older version of the program, I used rterm --vsize=100M --nsize=5000K --restore --save input file output file, however, in the new version R 2.2.0 , the parameters vsize and nsize are ignored. I can use the command memory.limit to increase memory, but I am not sure if this corresponds to vsize and nsize. Does anyone know how can I set these quantities now? __ R-help@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
