#Chapitre 10
require(BioStatR)
## Le chargement a nécessité le package : BioStatR
#page 403
foret<-rep(1:3,c(10,10,10))
hauteur<-c(23.4,24.4,24.6,24.9,25,26.2,26.1,24.8,25.5,25.8,18.9,21.1,21.1,
  22.1,22.5,23.5,22.7,21.3,22.2,21.7,22.5,22.9,23.7,24,24,24.5,24.3,24.2,
  23.4,23.9)
foret<-factor(foret)
arbre<-data.frame(foret,hauteur)
rm(foret)
rm(hauteur)
arbre
##    foret hauteur
## 1      1    23.4
## 2      1    24.4
## 3      1    24.6
## 4      1    24.9
## 5      1    25.0
## 6      1    26.2
## 7      1    26.1
## 8      1    24.8
## 9      1    25.5
## 10     1    25.8
## 11     2    18.9
## 12     2    21.1
## 13     2    21.1
## 14     2    22.1
## 15     2    22.5
## 16     2    23.5
## 17     2    22.7
## 18     2    21.3
## 19     2    22.2
## 20     2    21.7
## 21     3    22.5
## 22     3    22.9
## 23     3    23.7
## 24     3    24.0
## 25     3    24.0
## 26     3    24.5
## 27     3    24.3
## 28     3    24.2
## 29     3    23.4
## 30     3    23.9
moyennes<-tapply(arbre$hauteur,arbre$foret,mean)
moyennes
##     1     2     3 
## 25.07 21.71 23.74
#page 404
variances<-tapply(arbre$hauteur,arbre$foret,var)
variances
##         1         2         3 
## 0.7356667 1.5565556 0.4026667
#page 405
moy.g<-mean(arbre$hauteur)
moy.g
## [1] 23.50667
mean(moyennes)
## [1] 23.50667
#page 406
plot(arbre$foret,arbre$hauteur)
points(1:3,moyennes,pch="@")
abline(h=moy.g)

pdf("ch11fig101.pdf")
plot(arbre$foret,arbre$hauteur)
points(1:3,moyennes,pch="@")
abline(h=moy.g)
dev.off()
## quartz_off_screen 
##                 2
#page 409
options(contrasts=c("contr.sum","contr.poly"))
modele1<-lm(hauteur~foret,data=arbre)
anova(modele1)
## Analysis of Variance Table
## 
## Response: hauteur
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## foret      2 57.265 28.6323  31.874 7.809e-08 ***
## Residuals 27 24.254  0.8983                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
modele1_aov<-aov(hauteur~foret,data=arbre)
summary(modele1_aov)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## foret        2  57.26  28.632   31.87 7.81e-08 ***
## Residuals   27  24.25   0.898                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#page 410
options(contrasts=c("contr.sum","contr.poly"))

#page 411
residus<-residuals(modele1)
shapiro.test(residus)
## 
##  Shapiro-Wilk normality test
## 
## data:  residus
## W = 0.95202, p-value = 0.1914
length(residus)
## [1] 30
#En plus : les r\'esidus des deux mod\`eles sont \'egaux
all(residuals(modele1)==residuals(modele1_aov))
## [1] TRUE
#page 413
bartlett.test(residus~foret,data=arbre)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  residus by foret
## Bartlett's K-squared = 3.8798, df = 2, p-value = 0.1437
coef(modele1)
## (Intercept)      foret1      foret2 
##   23.506667    1.563333   -1.796667
#En plus : les coefficients des deux mod\`eles sont \'egaux
all(coef(modele1)==coef(modele1_aov))
## [1] TRUE
#page 414
-sum(coef(modele1)[2:3])
## [1] 0.2333333
dummy.coef(modele1)
## Full coefficients are 
##                                                 
## (Intercept):      23.50667                      
## foret:                   1          2          3
##                  1.5633333 -1.7966667  0.2333333
#En plus : fonctionne aussi avec le mod\`ele aov et introduction de la 
#fonction model.tables
dummy.coef(modele1_aov)
## Full coefficients are 
##                                                 
## (Intercept):      23.50667                      
## foret:                   1          2          3
##                  1.5633333 -1.7966667  0.2333333
model.tables(modele1_aov)
## Tables of effects
## 
##  foret 
## foret
##       1       2       3 
##  1.5633 -1.7967  0.2333
if(!("granova" %in% rownames(installed.packages()))){
  install.packages("granova")}
library(granova)
## Le chargement a nécessité le package : car
## Le chargement a nécessité le package : carData
granova.1w(arbre$hauteur,arbre$foret)

## $grandsum
##     Grandmean        df.bet       df.with        MS.bet       MS.with 
##         23.51          2.00         27.00         28.63          0.90 
##        F.stat        F.prob SS.bet/SS.tot 
##         31.87          0.00          0.70 
## 
## $stats
##   Size Contrast Coef Wt'd Mean  Mean Trim'd Mean Var. St. Dev.
## 2   10         -1.80     21.71 21.71       21.82 1.56     1.25
## 3   10          0.23     23.74 23.74       23.87 0.40     0.63
## 1   10          1.56     25.07 25.07       25.10 0.74     0.86
pdf("chap10fig102.pdf")
print(granova.1w(arbre$hauteur,arbre$foret))
## $grandsum
##     Grandmean        df.bet       df.with        MS.bet       MS.with 
##         23.51          2.00         27.00         28.63          0.90 
##        F.stat        F.prob SS.bet/SS.tot 
##         31.87          0.00          0.70 
## 
## $stats
##   Size Contrast Coef Wt'd Mean  Mean Trim'd Mean Var. St. Dev.
## 2   10         -1.80     21.71 21.71       21.82 1.56     1.25
## 3   10          0.23     23.74 23.74       23.87 0.40     0.63
## 1   10          1.56     25.07 25.07       25.10 0.74     0.86
dev.off()
## quartz_off_screen 
##                 2
#page 416
if(!("granovaGG" %in% rownames(installed.packages()))){
  install.packages("granovaGG")}
library(granovaGG)
## Le chargement a nécessité le package : ggplot2
granovagg.1w(arbre$hauteur,arbre$foret)
## 
## By-group summary statistics for your input data (ordered by group means)
##   group group.mean trimmed.mean contrast variance standard.deviation group.size
## 2     2      21.71        21.82    -1.80     1.56               1.25         10
## 3     3      23.74        23.87     0.23     0.40               0.63         10
## 1     1      25.07        25.10     1.56     0.74               0.86         10
## 
## Below is a linear model summary of your input data
## 
## Call:
## lm(formula = score ~ group, data = owp$data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.8100 -0.4550  0.0750  0.5425  1.7900 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  23.5067     0.1730 135.844  < 2e-16 ***
## group1        1.5633     0.2447   6.388 7.66e-07 ***
## group2       -1.7967     0.2447  -7.342 6.75e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9478 on 27 degrees of freedom
## Multiple R-squared:  0.7025, Adjusted R-squared:  0.6804 
## F-statistic: 31.87 on 2 and 27 DF,  p-value: 7.809e-08

pdf("chap10fig103.pdf")
print(granovagg.1w(arbre$hauteur,arbre$foret))
## 
## By-group summary statistics for your input data (ordered by group means)
##   group group.mean trimmed.mean contrast variance standard.deviation group.size
## 2     2      21.71        21.82    -1.80     1.56               1.25         10
## 3     3      23.74        23.87     0.23     0.40               0.63         10
## 1     1      25.07        25.10     1.56     0.74               0.86         10
## 
## Below is a linear model summary of your input data
## 
## Call:
## lm(formula = score ~ group, data = owp$data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.8100 -0.4550  0.0750  0.5425  1.7900 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  23.5067     0.1730 135.844  < 2e-16 ***
## group1        1.5633     0.2447   6.388 7.66e-07 ***
## group2       -1.7967     0.2447  -7.342 6.75e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9478 on 27 degrees of freedom
## Multiple R-squared:  0.7025, Adjusted R-squared:  0.6804 
## F-statistic: 31.87 on 2 and 27 DF,  p-value: 7.809e-08
dev.off()
## quartz_off_screen 
##                 2
#page 419
modele2<-aov(hauteur~foret,data=arbre)
model.tables(modele2)
## Tables of effects
## 
##  foret 
## foret
##       1       2       3 
##  1.5633 -1.7967  0.2333
TukeyHSD(modele2)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = hauteur ~ foret, data = arbre)
## 
## $foret
##      diff        lwr        upr     p adj
## 2-1 -3.36 -4.4109317 -2.3090683 0.0000000
## 3-1 -1.33 -2.3809317 -0.2790683 0.0110330
## 3-2  2.03  0.9790683  3.0809317 0.0001544
plot(TukeyHSD(modele2))

pdf("chap10fig104.pdf")
plot(TukeyHSD(modele2))
dev.off()
## quartz_off_screen 
##                 2
#En plus : export des graphiques en niveaux de gris et aux formats .png ou .ps
png("chap10fig102.png")
granova.1w(arbre$hauteur,arbre$foret)
## $grandsum
##     Grandmean        df.bet       df.with        MS.bet       MS.with 
##         23.51          2.00         27.00         28.63          0.90 
##        F.stat        F.prob SS.bet/SS.tot 
##         31.87          0.00          0.70 
## 
## $stats
##   Size Contrast Coef Wt'd Mean  Mean Trim'd Mean Var. St. Dev.
## 2   10         -1.80     21.71 21.71       21.82 1.56     1.25
## 3   10          0.23     23.74 23.74       23.87 0.40     0.63
## 1   10          1.56     25.07 25.07       25.10 0.74     0.86
dev.off()
## quartz_off_screen 
##                 2
postscript("chap10fig102.ps")
granova.1w(arbre$hauteur,arbre$foret)
## $grandsum
##     Grandmean        df.bet       df.with        MS.bet       MS.with 
##         23.51          2.00         27.00         28.63          0.90 
##        F.stat        F.prob SS.bet/SS.tot 
##         31.87          0.00          0.70 
## 
## $stats
##   Size Contrast Coef Wt'd Mean  Mean Trim'd Mean Var. St. Dev.
## 2   10         -1.80     21.71 21.71       21.82 1.56     1.25
## 3   10          0.23     23.74 23.74       23.87 0.40     0.63
## 1   10          1.56     25.07 25.07       25.10 0.74     0.86
dev.off()
## quartz_off_screen 
##                 2
pdf("chap10fig102bw.pdf",colormodel="gray")
granova.1w(arbre$hauteur,arbre$foret)
## $grandsum
##     Grandmean        df.bet       df.with        MS.bet       MS.with 
##         23.51          2.00         27.00         28.63          0.90 
##        F.stat        F.prob SS.bet/SS.tot 
##         31.87          0.00          0.70 
## 
## $stats
##   Size Contrast Coef Wt'd Mean  Mean Trim'd Mean Var. St. Dev.
## 2   10         -1.80     21.71 21.71       21.82 1.56     1.25
## 3   10          0.23     23.74 23.74       23.87 0.40     0.63
## 1   10          1.56     25.07 25.07       25.10 0.74     0.86
dev.off()
## quartz_off_screen 
##                 2
postscript("chap10fig102bw.ps",colormodel="gray")
granova.1w(arbre$hauteur,arbre$foret)
## $grandsum
##     Grandmean        df.bet       df.with        MS.bet       MS.with 
##         23.51          2.00         27.00         28.63          0.90 
##        F.stat        F.prob SS.bet/SS.tot 
##         31.87          0.00          0.70 
## 
## $stats
##   Size Contrast Coef Wt'd Mean  Mean Trim'd Mean Var. St. Dev.
## 2   10         -1.80     21.71 21.71       21.82 1.56     1.25
## 3   10          0.23     23.74 23.74       23.87 0.40     0.63
## 1   10          1.56     25.07 25.07       25.10 0.74     0.86
dev.off()
## quartz_off_screen 
##                 2
png("chap10fig103.png")
granovagg.1w(arbre$hauteur,arbre$foret)
## 
## By-group summary statistics for your input data (ordered by group means)
##   group group.mean trimmed.mean contrast variance standard.deviation group.size
## 2     2      21.71        21.82    -1.80     1.56               1.25         10
## 3     3      23.74        23.87     0.23     0.40               0.63         10
## 1     1      25.07        25.10     1.56     0.74               0.86         10
## 
## Below is a linear model summary of your input data
## 
## Call:
## lm(formula = score ~ group, data = owp$data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.8100 -0.4550  0.0750  0.5425  1.7900 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  23.5067     0.1730 135.844  < 2e-16 ***
## group1        1.5633     0.2447   6.388 7.66e-07 ***
## group2       -1.7967     0.2447  -7.342 6.75e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9478 on 27 degrees of freedom
## Multiple R-squared:  0.7025, Adjusted R-squared:  0.6804 
## F-statistic: 31.87 on 2 and 27 DF,  p-value: 7.809e-08
dev.off()
## quartz_off_screen 
##                 2
postscript("chap10fig103.ps")
granovagg.1w(arbre$hauteur,arbre$foret)
## 
## By-group summary statistics for your input data (ordered by group means)
##   group group.mean trimmed.mean contrast variance standard.deviation group.size
## 2     2      21.71        21.82    -1.80     1.56               1.25         10
## 3     3      23.74        23.87     0.23     0.40               0.63         10
## 1     1      25.07        25.10     1.56     0.74               0.86         10
## 
## Below is a linear model summary of your input data
## 
## Call:
## lm(formula = score ~ group, data = owp$data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.8100 -0.4550  0.0750  0.5425  1.7900 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  23.5067     0.1730 135.844  < 2e-16 ***
## group1        1.5633     0.2447   6.388 7.66e-07 ***
## group2       -1.7967     0.2447  -7.342 6.75e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9478 on 27 degrees of freedom
## Multiple R-squared:  0.7025, Adjusted R-squared:  0.6804 
## F-statistic: 31.87 on 2 and 27 DF,  p-value: 7.809e-08
## Warning in grid.Call.graphics(C_segments, x$x0, x$y0, x$x1, x$y1, x$arrow): la
## semi-transparence n'est pas supportée sur ce périphérique : reporté seulement
## une fois par page
dev.off()
## quartz_off_screen 
##                 2
pdf("chap10fig103bw.pdf",colormodel="gray")
granovagg.1w(arbre$hauteur,arbre$foret)
## 
## By-group summary statistics for your input data (ordered by group means)
##   group group.mean trimmed.mean contrast variance standard.deviation group.size
## 2     2      21.71        21.82    -1.80     1.56               1.25         10
## 3     3      23.74        23.87     0.23     0.40               0.63         10
## 1     1      25.07        25.10     1.56     0.74               0.86         10
## 
## Below is a linear model summary of your input data
## 
## Call:
## lm(formula = score ~ group, data = owp$data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.8100 -0.4550  0.0750  0.5425  1.7900 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  23.5067     0.1730 135.844  < 2e-16 ***
## group1        1.5633     0.2447   6.388 7.66e-07 ***
## group2       -1.7967     0.2447  -7.342 6.75e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9478 on 27 degrees of freedom
## Multiple R-squared:  0.7025, Adjusted R-squared:  0.6804 
## F-statistic: 31.87 on 2 and 27 DF,  p-value: 7.809e-08
dev.off()
## quartz_off_screen 
##                 2
postscript("chap10fig103bw.ps",colormodel="gray")
granovagg.1w(arbre$hauteur,arbre$foret)
## 
## By-group summary statistics for your input data (ordered by group means)
##   group group.mean trimmed.mean contrast variance standard.deviation group.size
## 2     2      21.71        21.82    -1.80     1.56               1.25         10
## 3     3      23.74        23.87     0.23     0.40               0.63         10
## 1     1      25.07        25.10     1.56     0.74               0.86         10
## 
## Below is a linear model summary of your input data
## 
## Call:
## lm(formula = score ~ group, data = owp$data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.8100 -0.4550  0.0750  0.5425  1.7900 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  23.5067     0.1730 135.844  < 2e-16 ***
## group1        1.5633     0.2447   6.388 7.66e-07 ***
## group2       -1.7967     0.2447  -7.342 6.75e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9478 on 27 degrees of freedom
## Multiple R-squared:  0.7025, Adjusted R-squared:  0.6804 
## F-statistic: 31.87 on 2 and 27 DF,  p-value: 7.809e-08
## Warning in grid.Call.graphics(C_segments, x$x0, x$y0, x$x1, x$y1, x$arrow): la
## semi-transparence n'est pas supportée sur ce périphérique : reporté seulement
## une fois par page
dev.off()
## quartz_off_screen 
##                 2
#page 426
#Exercice 10.1
#1)
options(contrasts=c(unordered="contr.sum", ordered="contr.poly"))

#page 427
#2)
variete<-rep(1:6,c(5,5,5,5,5,5))
vitamine<-c(93.6,95.3,96,93.7,96.2,95.3,96.9,95.8,97.3,97.7,94.5,97,97.8,97,
98.3,98.8,98.2,97.8,97.2,97.9,94.6,97.8,98,95,98.9,93.2,94.4,93.8,95.6,94.8)
variete<-factor(variete)
exo1<-data.frame(variete,vitamine)
modele1<-aov(vitamine~variete,data=exo1)
residus1<-residuals(modele1)
shapiro.test(residus1)
## 
##  Shapiro-Wilk normality test
## 
## data:  residus1
## W = 0.9563, p-value = 0.2485
length(residus1)
## [1] 30
bartlett.test(residus1~variete,data=exo1)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  residus1 by variete
## Bartlett's K-squared = 5.6023, df = 5, p-value = 0.3469
#page 428
#3)
modele1
## Call:
##    aov(formula = vitamine ~ variete, data = exo1)
## 
## Terms:
##                 variete Residuals
## Sum of Squares   45.836    38.512
## Deg. of Freedom       5        24
## 
## Residual standard error: 1.266754
## Estimated effects may be unbalanced
summary(modele1)
##             Df Sum Sq Mean Sq F value  Pr(>F)   
## variete      5  45.84   9.167   5.713 0.00131 **
## Residuals   24  38.51   1.605                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#page 429
#4)
granovagg.1w(vitamine,group=variete)
## 
## By-group summary statistics for your input data (ordered by group means)
##   group group.mean trimmed.mean contrast variance standard.deviation group.size
## 6     6      94.36        94.33    -1.92     0.85               0.92          5
## 1     1      94.96        95.00    -1.32     1.54               1.24          5
## 2     2      96.60        96.67     0.32     1.03               1.01          5
## 5     5      96.86        96.93     0.58     3.73               1.93          5
## 3     3      96.92        97.27     0.64     2.14               1.46          5
## 4     4      97.98        97.97     1.70     0.34               0.58          5
## 
## The following groups are likely to be overplotted
##   group group.mean contrast
## 5     5      96.86     0.58
## 3     3      96.92     0.64
## 
## Below is a linear model summary of your input data
## 
## Call:
## lm(formula = score ~ group, data = owp$data)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -2.420 -0.795  0.150  0.925  2.040 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  96.2800     0.2313 416.298  < 2e-16 ***
## group1       -1.3200     0.5172  -2.552  0.01748 *  
## group2        0.3200     0.5172   0.619  0.54189    
## group3        0.6400     0.5172   1.238  0.22785    
## group4        1.7000     0.5172   3.287  0.00311 ** 
## group5        0.5800     0.5172   1.122  0.27316    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.267 on 24 degrees of freedom
## Multiple R-squared:  0.5434, Adjusted R-squared:  0.4483 
## F-statistic: 5.713 on 5 and 24 DF,  p-value: 0.001311

pdf("chap10fig105.pdf")
granovagg.1w(vitamine,group=variete)
## 
## By-group summary statistics for your input data (ordered by group means)
##   group group.mean trimmed.mean contrast variance standard.deviation group.size
## 6     6      94.36        94.33    -1.92     0.85               0.92          5
## 1     1      94.96        95.00    -1.32     1.54               1.24          5
## 2     2      96.60        96.67     0.32     1.03               1.01          5
## 5     5      96.86        96.93     0.58     3.73               1.93          5
## 3     3      96.92        97.27     0.64     2.14               1.46          5
## 4     4      97.98        97.97     1.70     0.34               0.58          5
## 
## The following groups are likely to be overplotted
##   group group.mean contrast
## 5     5      96.86     0.58
## 3     3      96.92     0.64
## 
## Below is a linear model summary of your input data
## 
## Call:
## lm(formula = score ~ group, data = owp$data)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -2.420 -0.795  0.150  0.925  2.040 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  96.2800     0.2313 416.298  < 2e-16 ***
## group1       -1.3200     0.5172  -2.552  0.01748 *  
## group2        0.3200     0.5172   0.619  0.54189    
## group3        0.6400     0.5172   1.238  0.22785    
## group4        1.7000     0.5172   3.287  0.00311 ** 
## group5        0.5800     0.5172   1.122  0.27316    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.267 on 24 degrees of freedom
## Multiple R-squared:  0.5434, Adjusted R-squared:  0.4483 
## F-statistic: 5.713 on 5 and 24 DF,  p-value: 0.001311
dev.off()
## quartz_off_screen 
##                 2
#page 431
#6)
Tukey1 <- TukeyHSD(modele1, conf.level = 0.95)
Tukey1
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = vitamine ~ variete, data = exo1)
## 
## $variete
##      diff        lwr         upr     p adj
## 2-1  1.64 -0.8371495  4.11714947 0.3468489
## 3-1  1.96 -0.5171495  4.43714947 0.1803668
## 4-1  3.02  0.5428505  5.49714947 0.0107386
## 5-1  1.90 -0.5771495  4.37714947 0.2058535
## 6-1 -0.60 -3.0771495  1.87714947 0.9733815
## 3-2  0.32 -2.1571495  2.79714947 0.9985151
## 4-2  1.38 -1.0971495  3.85714947 0.5310572
## 5-2  0.26 -2.2171495  2.73714947 0.9994551
## 6-2 -2.24 -4.7171495  0.23714947 0.0926651
## 4-3  1.06 -1.4171495  3.53714947 0.7697394
## 5-3 -0.06 -2.5371495  2.41714947 0.9999996
## 6-3 -2.56 -5.0371495 -0.08285053 0.0399329
## 5-4 -1.12 -3.5971495  1.35714947 0.7278111
## 6-4 -3.62 -6.0971495 -1.14285053 0.0017510
## 6-5 -2.50 -4.9771495 -0.02285053 0.0470143
#page 432
#4)
if(!("multcomp" %in% rownames(installed.packages()))){
  install.packages("multcomp")}
library(multcomp)
## Le chargement a nécessité le package : mvtnorm
## Le chargement a nécessité le package : survival
## Le chargement a nécessité le package : TH.data
## Le chargement a nécessité le package : MASS
## 
## Attachement du package : 'TH.data'
## L'objet suivant est masqué depuis 'package:MASS':
## 
##     geyser
wht = glht(modele1, linfct = mcp(variete = "Tukey"))
cld(wht)
##    1    2    3    4    5    6 
## "ab" "ac" "bc"  "c" "bc"  "a"
plot(Tukey1)

pdf("chap10fig106.pdf")
plot(Tukey1)
dev.off()
## quartz_off_screen 
##                 2
#page 433
CI <- confint(wht)
fortify(CI)
##      lhs rhs estimate        lwr         upr
## 1  2 - 1   0     1.64 -0.8377859  4.11778586
## 2  3 - 1   0     1.96 -0.5177859  4.43778586
## 3  4 - 1   0     3.02  0.5422141  5.49778586
## 4  5 - 1   0     1.90 -0.5777859  4.37778586
## 5  6 - 1   0    -0.60 -3.0777859  1.87778586
## 6  3 - 2   0     0.32 -2.1577859  2.79778586
## 7  4 - 2   0     1.38 -1.0977859  3.85778586
## 8  5 - 2   0     0.26 -2.2177859  2.73778586
## 9  6 - 2   0    -2.24 -4.7177859  0.23778586
## 10 4 - 3   0     1.06 -1.4177859  3.53778586
## 11 5 - 3   0    -0.06 -2.5377859  2.41778586
## 12 6 - 3   0    -2.56 -5.0377859 -0.08221414
## 13 5 - 4   0    -1.12 -3.5977859  1.35778586
## 14 6 - 4   0    -3.62 -6.0977859 -1.14221414
## 15 6 - 5   0    -2.50 -4.9777859 -0.02221414
ggplot(CI,aes(lhs,estimate,ymin=lwr,ymax=upr))+geom_pointrange()+
  geom_hline(yintercept = 0)

pdf("chap10fig107.pdf")
print(ggplot(CI,aes(lhs,estimate,ymin=lwr,ymax=upr))+geom_pointrange()+
  geom_hline(yintercept = 0))
dev.off()
## quartz_off_screen 
##                 2
ggplot(aes(lhs,estimate),data=fortify(summary(wht))) + 
  geom_linerange(aes(ymin=lwr,ymax=upr),data=CI) + 
  geom_text(aes(y=estimate+1,label=round(p,3)))+geom_hline(yintercept = 0) +
  geom_point(aes(size=p),data=summary(wht)) +scale_size(trans="reverse")

pdf("chap10fig108.pdf")
ggplot(aes(lhs,estimate),data=fortify(summary(wht))) + 
  geom_linerange(aes(ymin=lwr,ymax=upr),data=CI) + 
  geom_text(aes(y=estimate+1,label=round(p,3)))+geom_hline(yintercept = 0) +
  geom_point(aes(size=p),data=summary(wht)) +scale_size(trans="reverse")
dev.off()
## quartz_off_screen 
##                 2
#page 434
if(!("multcompView" %in% rownames(installed.packages()))){
  install.packages("multcompView")}
library(multcompView)
if(!("plyr" %in% rownames(installed.packages()))){install.packages("plyr")}
library(plyr)
generate_label_df <- function(HSD,flev){
Tukey.levels <- HSD[[flev]][,4]
Tukey.labels <- multcompLetters(Tukey.levels)['Letters']
plot.labels <- names(Tukey.labels[['Letters']])
boxplot.df <- ddply(exo1, flev, function (x) max(fivenum(x$vitamine)) + 0.2)
plot.levels <- data.frame(plot.labels, labels = Tukey.labels[['Letters']],
  stringsAsFactors = FALSE)
labels.df <- merge(plot.levels, boxplot.df, by.x = 'plot.labels', by.y = flev,
  sort = FALSE)
return(labels.df)
}

#page 435
p_base <- ggplot(exo1,aes(x=variete,y=vitamine)) + geom_boxplot() + 
  geom_text(data = generate_label_df(Tukey1, 'variete'), aes(x = plot.labels, 
  y = V1, label = labels))
p_base

pdf("chap10fig109.pdf")
print(p_base)
dev.off()
## quartz_off_screen 
##                 2
#page 436
#Exercice 10.1
#2)
traitement<-rep(1:5,c(7,7,7,7,7))
taux<-c(4.5,2.5,6,4.5,3,5.5,3.5,7.5,3,2.5,4,2,4,5.5,8,6.5,6,3.5,5,
  7,5,2,7.5,4,2.5,5,3.5,6.5,6.5,5.5,6,4.5,4,7,5.5)
traitement<-factor(traitement)
exo2<-data.frame(traitement,taux)
modele2<-aov(taux~traitement,data=exo2)
residus2<-residuals(modele2)
shapiro.test(residus2)
## 
##  Shapiro-Wilk normality test
## 
## data:  residus2
## W = 0.97436, p-value = 0.5734
length(residus2)
## [1] 35
bartlett.test(residus2~traitement,data=exo2)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  residus2 by traitement
## Bartlett's K-squared = 3.1361, df = 4, p-value = 0.5353
#page 437
#3)
modele1<-lm(taux~traitement,data=exo2)
anova(modele1)
## Analysis of Variance Table
## 
## Response: taux
##            Df Sum Sq Mean Sq F value Pr(>F)
## traitement  4 19.043  4.7607  1.8687 0.1419
## Residuals  30 76.429  2.5476
#4)
power.anova.test(5,7,19.043,76.42857)
## 
##      Balanced one-way analysis of variance power calculation 
## 
##          groups = 5
##               n = 7
##     between.var = 19.043
##      within.var = 76.42857
##       sig.level = 0.05
##           power = 0.4684833
## 
## NOTE: n is number in each group
#page 438
power.anova.test(groups=5,between.var=19.043,within.var=76.42857,power=.80)
## 
##      Balanced one-way analysis of variance power calculation 
## 
##          groups = 5
##               n = 12.96035
##     between.var = 19.043
##      within.var = 76.42857
##       sig.level = 0.05
##           power = 0.8
## 
## NOTE: n is number in each group
granovagg.1w(taux,group=traitement)
## 
## By-group summary statistics for your input data (ordered by group means)
##   group group.mean trimmed.mean contrast variance standard.deviation group.size
## 2     2       4.07          3.8    -0.76     3.62               1.90          7
## 1     1       4.21          4.2    -0.61     1.65               1.29          7
## 4     4       4.43          4.3    -0.40     4.12               2.03          7
## 5     5       5.57          5.6     0.74     1.12               1.06          7
## 3     3       5.86          5.9     1.03     2.23               1.49          7
## 
## Below is a linear model summary of your input data
## 
## Call:
## lm(formula = score ~ group, data = owp$data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.4286 -1.0714 -0.0714  1.0357  3.4286 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   4.8286     0.2698  17.897   <2e-16 ***
## group1       -0.6143     0.5396  -1.138   0.2639    
## group2       -0.7571     0.5396  -1.403   0.1708    
## group3        1.0286     0.5396   1.906   0.0662 .  
## group4       -0.4000     0.5396  -0.741   0.4643    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.596 on 30 degrees of freedom
## Multiple R-squared:  0.1995, Adjusted R-squared:  0.09272 
## F-statistic: 1.869 on 4 and 30 DF,  p-value: 0.1419

pdf("chap10fig1010.pdf",colormodel="gray")
granovagg.1w(taux,group=traitement)
## 
## By-group summary statistics for your input data (ordered by group means)
##   group group.mean trimmed.mean contrast variance standard.deviation group.size
## 2     2       4.07          3.8    -0.76     3.62               1.90          7
## 1     1       4.21          4.2    -0.61     1.65               1.29          7
## 4     4       4.43          4.3    -0.40     4.12               2.03          7
## 5     5       5.57          5.6     0.74     1.12               1.06          7
## 3     3       5.86          5.9     1.03     2.23               1.49          7
## 
## Below is a linear model summary of your input data
## 
## Call:
## lm(formula = score ~ group, data = owp$data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.4286 -1.0714 -0.0714  1.0357  3.4286 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   4.8286     0.2698  17.897   <2e-16 ***
## group1       -0.6143     0.5396  -1.138   0.2639    
## group2       -0.7571     0.5396  -1.403   0.1708    
## group3        1.0286     0.5396   1.906   0.0662 .  
## group4       -0.4000     0.5396  -0.741   0.4643    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.596 on 30 degrees of freedom
## Multiple R-squared:  0.1995, Adjusted R-squared:  0.09272 
## F-statistic: 1.869 on 4 and 30 DF,  p-value: 0.1419
dev.off()
## quartz_off_screen 
##                 2