#Chapitre 10
require(BioStatR)
## Loading required package: BioStatR
#page 402
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 403
variances<-tapply(arbre$hauteur,arbre$foret,var)
variances
##         1         2         3 
## 0.7356667 1.5565556 0.4026667
#page 404
moy.g<-mean(arbre$hauteur)
moy.g
## [1] 23.50667
mean(moyennes)
## [1] 23.50667
#page 405
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 408
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 409
options(contrasts=c("contr.sum","contr.poly"))

#page 410
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ésidus des deux modèles sont égaux
all(residuals(modele1)==residuals(modele1_aov))
## [1] TRUE
#page 412
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èles sont égaux
all(coef(modele1)==coef(modele1_aov))
## [1] TRUE
#page 413
-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èle 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)
## Loading required package: car
## Loading required 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 415
if(!("granovaGG" %in% rownames(installed.packages()))){
  install.packages("granovaGG")}
library(granovaGG)
## Loading required 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
## 2     2      21.71        21.82    -1.80     1.56               1.25
## 3     3      23.74        23.87     0.23     0.40               0.63
## 1     1      25.07        25.10     1.56     0.74               0.86
##   group.size
## 2         10
## 3         10
## 1         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