ANOVA suite 18_1 SSP : predictions

Hulk a bu un jus d’ANOVA petit; il s’est transformer et devenu tres fort !

SIMU DE VECTORS PARAMETRICS 𝛔 et 𝞵

library(gplots)

## 
## Attachement du package : 'gplots'

## L'objet suivant est masqué depuis 'package:stats':
## 
##     lowess

library(knitr)
set.seed(123)
V1=rnorm(19,6.1,1.5)#parameters table line1...
V2=rnorm(20,6.6,1.8)#parameters table line1...
V3=rnorm(15,7.7,2)
V4=rnorm(17,8.7,1.5)
#make a long vector
V1234=c(V1,V2,V3,V4)
Myf=c(rep("0mat",19),rep("1er",20),rep("2er",15),rep("3er",17))
Myf=factor(Myf)#gl is eq : old fx gl Splus
simu181=data.frame("measure"=V1234,"Groupsj"=Myf)
summary(simu181)

##     measure       Groupsj  
##  Min.   : 3.150   0mat:19  
##  1st Qu.: 6.137   1er :20  
##  Median : 7.172   2er :15  
##  Mean   : 7.303   3er :17  
##  3rd Qu.: 8.400            
##  Max.   :12.038

kable(head(simu181,10))

measure	Groupsj
5.259286	0mat
5.754734	0mat
8.438062	0mat
6.205763	0mat
6.293932	0mat
8.672598	0mat
6.791374	0mat
4.202408	0mat
5.069721	0mat
5.431507	0mat

save(simu181,file= "simu181.Rda")
load("simu181.Rda")
getwd()

## [1] "C:/R/SSP"

A1=simu181
summary(A1)

##     measure       Groupsj  
##  Min.   : 3.150   0mat:19  
##  1st Qu.: 6.137   1er :20  
##  Median : 7.172   2er :15  
##  Mean   : 7.303   3er :17  
##  3rd Qu.: 8.400            
##  Max.   :12.038

plotmeans(A1$measure~A1$Groupsj,col=c(1,2,3,4),cex=3)

Verifions juste l’homogeneite de la variance (within Group.j)

library(car)

## Le chargement a nécessité le package : carData

leveneTest(A1$measure~A1$Groupsj)#P(T\t>>0.05% OK

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3  0.4056 0.7495
##       67

Pourquoi peut on faire des predictions avec un ANOVA

En fait dans oav vous avez la fonction lm (Rgerssion lineaire qui est une droite affine d’équation \(E[Y]=a*X+b\) souvenez vos cours HARMOS).Pous nous statisticien nous appelons l’ordonnée a l’origine 𝛃₀ et 𝛃₁ (–Intercept and slope–)
La condition la plus important pour ce model lm c’est que les résidus e_i (ou l’erreur aléatoire résiduel non prévus du model [connu des econometricien d’exogeneité] doit etre normaux de variance constante et sans correlation.

Qu’est ce qu’un résidu ou erreur _i?

C’est la variation/erreur e_i aléatoire que le model ne peut pas prédire ():

Un exemple du net:

Residus droite de regression

Conditions des résidus N(0,1,i.i.d,no cor)

REGRESSION SUR UN ANOVA

Nous serions tenté de tirer une droite dans le plot Means de l’Anova non? Essayons:

modaov=lm(simu181$measure~simu181$Groupsj)
summary(modaov)

## 
## Call:
## lm(formula = simu181$measure ~ simu181$Groupsj)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.2109 -0.9873 -0.0949  0.9829  4.1935 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          6.3609     0.3571  17.811  < 2e-16 ***
## simu181$Groupsj1er   0.1385     0.4987   0.278  0.78211    
## simu181$Groupsj2er   1.4834     0.5377   2.759  0.00747 ** 
## simu181$Groupsj3er   2.4631     0.5197   4.739 1.16e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.557 on 67 degrees of freedom
## Multiple R-squared:  0.3116, Adjusted R-squared:  0.2808 
## F-statistic: 10.11 on 3 and 67 DF,  p-value: 1.402e-05

#pour nos moyennes de groupes:

kable(tapply(simu181$measure,simu181$Groupsj,mean),caption="mean.j")

mean.j
	x
0mat	6.360942
1er	6.499428
2er	7.844384
3er	8.824073

6.36 est l’intercept donc l’ordonnée a l’origione car 0mat est n’a pas de beta 1 et est inclus dans lintercept( on appel cela ~un contrast)

##sommons les coefficient beta de la regression le grooupe 0 mat etant notre groupe reference [contrast trt]
6.36+0.1385*1

## [1] 6.4985

6.36+1.4*1

## [1] 7.76

6.36+2.4631*1

## [1] 8.8231

Par Magie nous retouvons les moyennes des groupes anova

LA REGRESSION LM

plot(simu181$measure~as.numeric(simu181$Groupsj),col=rainbow(1),pch=20,main="4 groups points and regression line lm")
abline(lm(simu181$measure~as.numeric(simu181$Groupsj)),lwd=2,col4)

#in lm Y is the numerical dependant variable X is the numerical independant variable

Tout comme pour valider nos model en statistiques et c’est valable danstous les model suivant une loi N ouT les résidus doivent etre i.i.d et N(0,1)

AOV181=aov(lm(simu181$measure~simu181$Groupsj))#Anova object
hist(AOV181$residuals,,main="RESIDUS AOV")#pas mal du tout visu

hist(modaov$residuals,col=3,main="RESIDUS LM")

shapiro.test(AOV181$residuals)#P(T>t)>5% Ho accepted donc NORMAL

## 
##  Shapiro-Wilk normality test
## 
## data:  AOV181$residuals
## W = 0.98967, p-value = 0.8316

#ce sont les même !!!? WHY

Comment reconstruire depuis un objets ANOVA cette les points: Predictions + residus = Y_ij

Zerodif=(AOV181$fitted.values+AOV181$residuals)-V1234
kable(head(Zerodif),caption="Y~ij~-Anova")

Y_ij-Anova
x
0
0
0
0
0
0

summary(AOV181) #ANOVA unbalanced

##                 Df Sum Sq Mean Sq F value  Pr(>F)    
## simu181$Groupsj  3  73.51  24.502   10.11 1.4e-05 ***
## Residuals       67 162.36   2.423                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

AOV181$fitted.values

##        1        2        3        4        5        6        7        8 
## 6.360942 6.360942 6.360942 6.360942 6.360942 6.360942 6.360942 6.360942 
##        9       10       11       12       13       14       15       16 
## 6.360942 6.360942 6.360942 6.360942 6.360942 6.360942 6.360942 6.360942 
##       17       18       19       20       21       22       23       24 
## 6.360942 6.360942 6.360942 6.499428 6.499428 6.499428 6.499428 6.499428 
##       25       26       27       28       29       30       31       32 
## 6.499428 6.499428 6.499428 6.499428 6.499428 6.499428 6.499428 6.499428 
##       33       34       35       36       37       38       39       40 
## 6.499428 6.499428 6.499428 6.499428 6.499428 6.499428 6.499428 7.844384 
##       41       42       43       44       45       46       47       48 
## 7.844384 7.844384 7.844384 7.844384 7.844384 7.844384 7.844384 7.844384 
##       49       50       51       52       53       54       55       56 
## 7.844384 7.844384 7.844384 7.844384 7.844384 7.844384 8.824073 8.824073 
##       57       58       59       60       61       62       63       64 
## 8.824073 8.824073 8.824073 8.824073 8.824073 8.824073 8.824073 8.824073 
##       65       66       67       68       69       70       71 
## 8.824073 8.824073 8.824073 8.824073 8.824073 8.824073 8.824073

plot(AOV181$fitted.values,AOV181$residuals)
abline(h=0,col=3)

#on remarque trs peu d'heteroskedasticité
qqnorm(AOV181$residuals)
qqline(AOV181$residuals, col = "red") #NORMAKITY CHECK

Pourquoi ANOVA est couplé avec lm

Tout simplement que une regression sur ANOVA c’est une regression sur les moyennes des group_j [E¦¦X] (ou presque) on les appelent adjusted MEANS: LA suite ds un autre Rpubs

Et si les residus ne sont par N(0,1) Le model est faux!

Hulk est la: DU jus de Kruskall_wallis : Non-parametric alternative to one-way ANOVA test

kruskal.test(A1$measure~A1$Groupsj)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  A1$measure by A1$Groupsj
## Kruskal-Wallis chi-squared = 23.036, df = 3, p-value = 3.97e-05

#Note is a Chisq distributed
pchisq(23.036,3,lower.tail=FALSE)#by hand

## [1] 3.969126e-05

qchisq(0.95,3)

## [1] 7.814728

1-pchisq(7.81,3)#that is a 5% 1-p=lowertailFalse

## [1] 0.05010606