#Chapitre 2
#page 22
set.seed(314)
x=rnorm(50)
Y=function(x){x^2}
par(oma=rep(0,4));par(mar=c(3, 3, 1, 1)+0.1,mgp=c(2,1,0))
plot(Y,xlab="X",from=-3, to=3,lwd=2,col="blue")
points(x,Y(x),col="red")

#page 23
data(anscombe)
with(anscombe,cor(x1,y1))
## [1] 0.8164205
with(anscombe,cor(x2,y2))
## [1] 0.8162365
with(anscombe,cor(x3,y3))
## [1] 0.8162867
with(anscombe,cor(x4,y4))
## [1] 0.8165214
old.par <- par(no.readonly = TRUE)

par(mar=c(3,3,2,1))
par(mgp=c(2,1,0))

#page 24
layout(matrix(1:4,nrow=2,ncol=2,byrow=TRUE))
with(anscombe,plot(x1,y1,main=substitute(rho(x1,y1) == x, 
     list(x = format(with(anscombe,cor(x1,y1)), digits = 3, nsmall = 3)))))
with(anscombe,plot(x2,y2,main=substitute(rho(x2,y2) == x, 
     list(x = format(with(anscombe,cor(x2,y2)), digits = 3, nsmall = 3)))))
with(anscombe,plot(x3,y3,main=substitute(rho(x3,y3) == x, 
     list(x = format(with(anscombe,cor(x3,y3)), digits = 3, nsmall = 3)))))
with(anscombe,plot(x4,y4,main=substitute(rho(x4,y4) == x, 
     list(x = format(with(anscombe,cor(x4,y4)), digits = 3, nsmall = 3)))))

par(old.par)
layout(1)

#page 25
library(mvtnorm)
help(package="mvtnorm")
dmvnorm(c(0,0), c(0,0), diag(2), log=FALSE)
## [1] 0.1591549
#page 26
#Si le package n'est pas installe, enlever le commentaire
#puis executer la commande ci-dessous.
#if(!require("lattice")){install.packages("lattice")}
library(lattice)
# N2(c(0,0),I_2)
g <- expand.grid(x=seq(-2,2,0.05),y=seq(-2,2,0.05))
g$z <- dmvnorm(x=cbind(g$x,g$y),mean=c(0,0), sigma=diag(2),log=FALSE)
wireframe(z ~ x*y, data = g,colorkey = TRUE,drape=TRUE)

#page 27
# N2(c(0,0),matrix(c(1,0.75,0.75,1),byrow=T,nrow=2))
var <- matrix(c(1,0.75,0.75,1),byrow=T,nrow=2)
g <- expand.grid(x = seq(-2,2,0.05), y = seq(-2,2,0.05))
g$z <- dmvnorm(x=cbind(g$x,g$y),mean=c(0,0), sigma=var, log=FALSE)
wireframe(z ~ x*y, data = g,colorkey = TRUE,drape=TRUE)

#Si le package n'est pas installe, enlever le commentaire
#puis executer la commande ci-dessous.
#if(!require("rgl")){install.packages("rgl")}
library(rgl)
# N2(c(0,0),matrix(c(1,0,0,1),byrow=T,nrow=2))
g <- expand.grid(x = seq(-4,4,0.05), y = seq(-4,4,0.05))
g$z <- dmvnorm(x=cbind(g$x,g$y),mean=c(0,0),
               sigma=diag(2), log=FALSE)
g2z <- matrix(g$z*5000,byrow=T, nrow=length(seq(-4,4,0.05)))
g2x <- 10*(1:nrow(g2z))
g2y <- 10*(1:ncol(g2z))
zlim <- range(g2y)
zlen <- zlim[2]-zlim[1]+1
colorlut <- terrain.colors(zlen) # table des couleurs
col <- colorlut[ g2y-zlim[1]+1 ] # couleur par point
open3d()
## glX 
##   1
surface3d(g2x, g2y, g2z, color=col, back="lines")
#Pour sauvegarder le graphique enlever le commentaire de la commande ci-dessous
# snapshot3d("rgl_norm01.png")

# N2(c(0,0),matrix(c(1,0.75,0.75,1),byrow=T,nrow=2))
g <- expand.grid(x = seq(-4,4,0.05), y = seq(-4,4,0.05))
g$z <- dmvnorm(x=cbind(g$x,g$y),mean=c(0,0), 
               sigma=matrix(c(1,0.75,0.75,1), byrow=T,nrow=2), log=FALSE)
g2z <- matrix(g$z*5000,byrow=T, nrow=length(seq(-4,4,0.05)))
g2x <- 10*(1:nrow(g2z))
g2y <- 10*(1:ncol(g2z))
zlim <- range(g2y)
zlen <- zlim[2]-zlim[1]+1
colorlut <- terrain.colors(zlen) # table des couleurs
col <- colorlut[ g2y-zlim[1]+1 ] # couleur par point
open3d()
## glX 
##   2
surface3d(g2x, g2y, g2z, color=col, back="lines")
#Pour sauvegarder le graphique enlever le commentaire de la commande ci-dessous
# snapshot3d("rgl_norm02.png")

#page 30
cor(x,Y(x))
## [1] -0.3902511
#page 32
#Si le package n'est pas installe, enlever le commentaire
#puis executer la commande ci-dessous.
#if(!require("ModStatR")){install.packages("ModStatR")}
library(ModStatR)
#La fonction Gauss2F1 a ete integree `a la bibliotheque ModStatR
ModStatR::Gauss2F1
## function (a, b, c, x) 
## {
##     args <- lapply(as.list(match.call())[-1L], eval, parent.frame())
##     names <- if (is.null(names(args))) 
##         character(length(args))
##     else names(args)
##     dovec <- names %in% vectorize.args
##     do.call("mapply", c(FUN = FUN, args[dovec], MoreArgs = list(args[!dovec]), 
##         SIMPLIFY = SIMPLIFY, USE.NAMES = USE.NAMES))
## }
## <bytecode: 0x7fdcbda1ea98>
## <environment: 0x7fdcbda2aff8>
(Gauss2F1(1/2,1/2,(length(x)-2)/2,1-cor(x,Y(x))^2))*cor(x,Y(x))
## [1] -0.3938385
#La fonction Gauss2F1gsl a ete integree `a la bibliotheque ModStatR
ModStatR::Gauss2F1gsl
## function (a, b, c, x) 
## {
##     args <- lapply(as.list(match.call())[-1L], eval, parent.frame())
##     names <- if (is.null(names(args))) 
##         character(length(args))
##     else names(args)
##     dovec <- names %in% vectorize.args
##     do.call("mapply", c(FUN = FUN, args[dovec], MoreArgs = list(args[!dovec]), 
##         SIMPLIFY = SIMPLIFY, USE.NAMES = USE.NAMES))
## }
## <bytecode: 0x7fdcbd042c18>
## <environment: 0x7fdcbd0512e0>
(Gauss2F1gsl(1/2,1/2,(length(x)-2)/2,1-cor(x,Y(x))^2))*cor(x,Y(x))
## [1] -0.3938385
#page 33
library(mvtnorm)
sigma <- matrix(c(4,2.5,2.5,3), ncol=2)
cov2cor(sigma)
##           [,1]      [,2]
## [1,] 1.0000000 0.7216878
## [2,] 0.7216878 1.0000000
cor_res = NULL; cor_sb_res = NULL
cor_sb_gsl_res = NULL; cor_sbapprox_res = NULL
set.seed(2718)
for(iii in 1:1000){
  simul_XY <- rmvnorm(n=20, mean=c(1,2), sigma=sigma)
  X=simul_XY[,1]
  Y=simul_XY[,2]
  cor_res = c(cor_res,cor(X,Y))
  cor_sb_res = c(cor_sb_res,(Gauss2F1(1/2,1/2,
                (length(X)-2)/2,1-cor(X,Y)^2))*cor(X,Y))
  cor_sb_gsl_res = c(cor_sb_gsl_res,(Gauss2F1gsl(1/2,1/2,
                (length(X)-2)/2,1-cor(X,Y)^2))*cor(X,Y))
  cor_sbapprox_res = c(cor_sbapprox_res,cor(X,Y)+
                cor(X,Y)*(1-cor(X,Y)^2)/(length(X)-2))
}

#Si le package n'est pas installe, enlever le commentaire
#puis executer la commande ci-dessous.
#if(!require("yarrr")){install.packages("yarrr")}
library(yarrr)
## Loading required package: jpeg
## Loading required package: BayesFactor
## Loading required package: coda
## Loading required package: Matrix
## ************
## Welcome to BayesFactor 0.9.12-4.2. If you have questions, please contact Richard Morey (richarddmorey@gmail.com).
## 
## Type BFManual() to open the manual.
## ************
## Loading required package: circlize
## ========================================
## circlize version 0.4.8
## CRAN page: https://cran.r-project.org/package=circlize
## Github page: https://github.com/jokergoo/circlize
## Documentation: http://jokergoo.github.io/circlize_book/book/
## 
## If you use it in published research, please cite:
## Gu, Z. circlize implements and enhances circular visualization 
##   in R. Bioinformatics 2014.
## ========================================
## yarrr v0.1.6. Citation info at citation('yarrr'). Package guide at yarrr.guide()
## Email me at Nathaniel.D.Phillips.is@gmail.com
long_res = cbind(data.frame(res=c(cor_res,cor_sb_res,
  cor_sb_gsl_res,cor_sbapprox_res)),type=rep(1:4,rep(5000,4)))
pirateplot(res~type,data=long_res)
abline(h=cov2cor(sigma)[2,1],lwd=2,col="red")

mean_res = c(cov2cor(sigma)[2,1],with(long_res,
  tapply(res,type,mean)))
names(mean_res) <- c("true value","cor_res","cor_sb_res",
  "cor_sb_gsl_res","cor_sbapprox_res")
mean_res
##       true value          cor_res       cor_sb_res   cor_sb_gsl_res 
##        0.7216878        0.7205181        0.7224218        0.7224218 
## cor_sbapprox_res 
##        0.7240996
#page 34
all.equal(cor_sb_res, cor_sb_gsl_res)
## [1] TRUE
#page 36
#Si le package n'est pas installe, enlever le commentaire
#puis executer la commande ci-dessous.
#if(!require("pwr")){install.packages("pwr")}
library(pwr)
pwr.r.test(r=0.30,n=50,sig.level=0.05,alternative="two.sided")
## 
##      approximate correlation power calculation (arctangh transformation) 
## 
##               n = 50
##               r = 0.3
##       sig.level = 0.05
##           power = 0.5715558
##     alternative = two.sided
pwr.r.test(r=0.30,n=50,sig.level=0.05,alternative="greater")
## 
##      approximate correlation power calculation (arctangh transformation) 
## 
##               n = 50
##               r = 0.3
##       sig.level = 0.05
##           power = 0.6911395
##     alternative = greater
#page 37
pwr.r.test(r=0.3,power=0.80,sig.level=0.05, alternative="two.sided")
## 
##      approximate correlation power calculation (arctangh transformation) 
## 
##               n = 84.07364
##               r = 0.3
##       sig.level = 0.05
##           power = 0.8
##     alternative = two.sided
#page 40
#Si le package n'est pas installe, enlever le commentaire
#puis executer la commande ci-dessous.
#if(!require("BioStatR")){install.packages("BioStatR")}
library(BioStatR)
data("Quetelet")
str(Quetelet)
## 'data.frame':    66 obs. of  3 variables:
##  $ sexe  : Factor w/ 2 levels "f","h": 2 1 1 1 1 1 1 2 1 2 ...
##  $ poids : int  60 57 51 55 50 50 48 72 52 64 ...
##  $ taille: int  170 169 172 174 168 161 162 189 160 175 ...
table(Quetelet$sexe)
## 
##  f  h 
## 25 41
#page 41
#Si le package n'est pas installe, enlever le commentaire
#puis executer la commande ci-dessous.
#if(!require("GGally")){install.packages("GGally")}
library(GGally)
## Loading required package: ggplot2
## 
## Attaching package: 'ggplot2'
## The following object is masked from 'package:yarrr':
## 
##     diamonds
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
pairquet <- ggpairs(Quetelet)
print(pairquet)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Quet_H <- subset(Quetelet,subset=sexe=="h") 
nrow(Quet_H)
## [1] 41
#page 42
#Si le package n'est pas installe, enlever le commentaire
#puis executer la commande ci-dessous.
#if(!require("MVN")){install.packages("MVN")}
library(MVN)
## sROC 0.1-2 loaded
result = mvn(data = Quet_H[,-1], mvnTest = "mardia",
               univariateTest = "SW", univariatePlot = "histogram",
               multivariatePlot = "qq", multivariateOutlierMethod =
               "adj", showOutliers = TRUE, showNewData = TRUE)

result$multivariateNormality
##              Test        Statistic           p value Result
## 1 Mardia Skewness 3.58579754080577 0.464953104910953    YES
## 2 Mardia Kurtosis 1.19851769719842  0.23071553760644    YES
## 3             MVN             <NA>              <NA>    YES
#page 43
result$univariateNormality
##           Test  Variable Statistic   p value Normality
## 1 Shapiro-Wilk   poids      0.9728    0.4238    YES   
## 2 Shapiro-Wilk  taille      0.9737    0.4511    YES
result$multivariateOutliers
##    Observation Mahalanobis Distance Outlier
## 31          31               27.078    TRUE
## 12          12               19.530    TRUE
Quet_H[c("12","31"),]
##    sexe poids taille
## 12    h    72    164
## 31    h    77    200
#page 44
result = mvn(data = Quet_H[,-1], mvnTest = "mardia",
               univariateTest = "SW", univariatePlot = "histogram",
               multivariatePlot = "persp", multivariateOutlierMethod =
                 "adj", showOutliers = TRUE)

result = mvn(data = Quet_H[,-1], mvnTest = "mardia", 
             univariateTest = "SW", univariatePlot = "histogram", 
             multivariatePlot = "contour", multivariateOutlierMethod = 
               "adj", showOutliers = TRUE)

#page 45
result = mvn(data = Quet_H[,-1], mvnTest = "mardia",
             univariateTest = "SW", univariatePlot = "box",
             multivariatePlot = "qq", multivariateOutlierMethod =
               "adj", showOutliers = TRUE)
## Warning in uniPlot(data, type = univariatePlot): Box-Plots are based on
## standardized values (centered and scaled).

result = mvn(data = Quet_H[,-1], mvnTest = "mardia",
             univariateTest = "SW", univariatePlot = "scatter",
             multivariatePlot = "qq", multivariateOutlierMethod =
               "adj", showOutliers = TRUE)

cor.test(Quet_H[,"poids"],Quet_H[,"taille"])
## 
##  Pearson's product-moment correlation
## 
## data:  Quet_H[, "poids"] and Quet_H[, "taille"]
## t = 4.6441, df = 39, p-value = 3.823e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.3540921 0.7641071
## sample estimates:
##       cor 
## 0.5967362
#page 46
corobs = cor.test(Quet_H[,"poids"],Quet_H[,"taille"])$estimate
tanh(atanh(corobs)-qnorm(.975)/sqrt(41-3))
##       cor 
## 0.3540921
tanh(atanh(corobs)+qnorm(.975)/sqrt(41-3))
##       cor 
## 0.7641071
tanh(atanh(corobs)-qnorm(.975)/sqrt(41-3))-corobs/(2*(41-1))
##       cor 
## 0.3466329
tanh(atanh(corobs)+qnorm(.975)/sqrt(41-3))-corobs/(2*(41-1))
##       cor 
## 0.7566479
#page 47
#Si le package n'est pas installe, enlever le commentaire
#puis executer la commande ci-dessous.
#if(!require("MBESS")){install.packages("MBESS")}
library(MBESS)
ci.R(corobs,1,41-2)
## $Lower.Conf.Limit.R
## [1] 0.3498668
## 
## $Prob.Less.Lower
## [1] 0.025
## 
## $Upper.Conf.Limit.R
## [1] 0.7594818
## 
## $Prob.Greater.Upper
## [1] 0.025
sqrt(41-2)*corobs/sqrt(1-corobs^2)
##     cor 
## 4.64412
2*(1-pt(sqrt(41-2)*abs(corobs)/sqrt(1-corobs^2),41-2))
##          cor 
## 3.823038e-05
#page 48
#La fonction rho est disponible dans le package ModStatR
ModStatR::rho
## function(x, y, indices){
##   restest <- cor.test(x[indices], y[indices],
##                       method = c("pearson"))
##   return(restest$estimate)
## }
## <bytecode: 0x7fdc8c80ecf0>
## <environment: namespace:ModStatR>
#Si le package n'est pas installe, enlever le commentaire
#puis executer la commande ci-dessous.
#if(!require("boot")){install.packages("boot")}
library(boot)
## 
## Attaching package: 'boot'
## The following object is masked from 'package:lattice':
## 
##     melanoma
boot.rho_H <- boot(Quet_H[,"poids"], y=Quet_H[,"taille"], rho,
                    R=10000)

boot.ci(boot.rho_H)
## Warning in boot.ci(boot.rho_H): bootstrap variances needed for studentized
## intervals
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = boot.rho_H)
## 
## Intervals : 
## Level      Normal              Basic         
## 95%   ( 0.4213,  0.7594 )   ( 0.4216,  0.7593 )  
## 
## Level     Percentile            BCa          
## 95%   ( 0.4342,  0.7719 )   ( 0.4137,  0.7532 )  
## Calculations and Intervals on Original Scale
#page 49
rho_0 <- 0.78
n <- 41

sqrt(n-3)*(atanh(corobs)-atanh(rho_0))
##       cor 
## -2.202592
2*(1-pnorm(abs(sqrt(n-3)*(atanh(corobs)-atanh(rho_0)))))
##        cor 
## 0.02762352
sqrt(n-3)*(atanh(corobs)-atanh(rho_0)-rho_0/(2*(n-1)))
##       cor 
## -2.262695
2*(1-pnorm(abs(sqrt(n-3)*(atanh(corobs)-
                            atanh(rho_0)-rho_0/(2*(n-1))))))
##       cor 
## 0.0236545
#page 50
#Si le package n'est pas installe, enlever le commentaire
#puis executer la commande ci-dessous.
#if(!require("ModStatR")){install.packages("ModStatR")}
library(ModStatR)
#La fonction corrdist est disponible dans le package ModStatR
ModStatR::corrdist
## function (rho, rho_0, n) {
##   y = (n-2)*gamma(n-1)*(1-rho_0^2)^((n-1)/2)*(1-rho^2)^((n-4)/2)
##   y = y/(sqrt(2*pi)*gamma(n-1/2)*(1-rho_0*rho)^(n-3/2))
##   y = y*Gauss2F1(1/2,1/2,(2*n-1)/2,(rho_0*rho+1)/2)
##   return(y)
## }
## <bytecode: 0x7fdcebac4fa8>
## <environment: namespace:ModStatR>
#La fonction corrdistapprox est disponible dans le package ModStatR
ModStatR::corrdistapprox
## function (rho, rho_0, n) {
##   y = (n-2)*gamma(n-1)*(1-rho_0^2)^((n-1)/2)*(1-rho^2)^((n-4)/2)
##   y = y/(sqrt(2*pi)*gamma(n-1/2)*(1-rho_0*rho)^(n-3/2))
##   y = y*(1+1/4*(rho_0*rho+1)/(2*n-1)+9/32*(rho_0*rho+1)^2/
##            (2*n-1)/(2*n+1))
##   return(y)
## }
## <bytecode: 0x7fdceba72758>
## <environment: namespace:ModStatR>
#La fonction corrdistapprox2 est disponible dans le package ModStatR
ModStatR::corrdistapprox2
## function (rho, rho_0, n) {
##   y = (n-2)*gamma(n-1)*(1-rho_0^2)^((n-1)/2)*(1-rho^2)^((n-4)/2)
##   y = y/(sqrt(2*pi)*gamma(n-1/2)*(1-rho_0*rho)^(n-3/2))
##   y = y*(
##     1+1/4*(rho_0*rho+1)/(2*n-1)+ 
##     9/32*(rho_0*rho+1)^2/(2*n+1)/(2*n-1)+ 
##     75/128*(rho_0*rho+1)^3/((2*n+3)*(4*n^2-1))+ 
##     3675/2048*(rho_0*rho+1)^4/((2*n+3)*(2*n+5)*(4*n^2-1))+ 
##     59535/8192*(rho_0*rho+1)^5/((2*n+3)*(2*n+5)*(2*n + 7)*(4*n^2-1))
##   )
##   return(y)
## }
## <bytecode: 0x7fdceba23520>
## <environment: namespace:ModStatR>
print(corrdist(.5,0,20))
## [1] 0.167112
print(corrdistapprox(.5,0,20))
## [1] 0.1671105
print(corrdistapprox2(.5,0,20))
## [1] 0.167112
#page 51
integrate(corrdist, lower=-1, upper=corobs, rho_0=rho_0, n=n)
## 0.01175932 with absolute error < 3.6e-10
2*integrate(corrdist, lower=-1, upper=corobs, rho_0=rho_0,
            n=n)$value
## [1] 0.02351864
plot(seq(-1,1,.01), corrdist(seq(-1,1,.01), rho_0, n) ,type="l",
       xlab="Pearson correlation", ylab="Density")
polygon(x=c(seq(-1,corobs, length=100), corobs, 0),
          y=c(sapply(seq(-1,corobs, length=100),
          corrdist, rho_0=rho_0, n=n), 0,0), col="grey")

#page 52
#La fonction ref.cor.test est disponible dans le package ModStatR
ModStatR::ref.cor.test
## function(corobs, rho_0, n){
##   if(corobs<rho_0){
##     return(2*integrate(corrdist, lower=-1, upper=corobs, 
##                        rho_0=rho_0, n=n)$value)
##   }
##   if(corobs>rho_0){
##     return(2*integrate(corrdist, lower=corobs, upper=1,
##                        rho_0=rho_0, n=n)$value)
##   }
##   if(corobs==rho_0){
##     return(1)
##   }
## }
## <bytecode: 0x7fdcd9f15ce8>
## <environment: namespace:ModStatR>
ref.cor.test(corobs=corobs,rho_0=rho_0,n=n)
## [1] 0.02351864
Quet_F <- subset(Quetelet,subset=sexe=="f")
nrow(Quet_F)
## [1] 25
#Si le package n'est pas installe, enlever le commentaire
#puis executer la commande ci-dessous.
#if(!require("jmuOutlier")){install.packages("jmuOutlier")}
library(jmuOutlier)
set.seed(1133)
perm.cor.test(Quet_F[,"poids"],Quet_F[,"taille"])
## [[1]]
## [1] "Permutation correlation test.  Method is pearson"
## 
## [[2]]
## [1] "p-value was estimated based on 20000 simulations."
## 
## $alternative
## [1] "two.sided"
## 
## $p.value
## [1] 0.0033
#page 53
#Si le package n'est pas installe, enlever le commentaire
#puis executer la commande ci-dessous.
#if(!require("boot")){install.packages("boot")}
library(boot)
boot.rho_F <- boot(Quet_F[,"poids"], y=Quet_F[,"taille"], rho,
                     R=10000)
boot.ci(boot.rho_F)
## Warning in boot.ci(boot.rho_F): bootstrap variances needed for studentized
## intervals
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = boot.rho_F)
## 
## Intervals : 
## Level      Normal              Basic         
## 95%   ( 0.3525,  0.7542 )   ( 0.3765,  0.7770 )  
## 
## Level     Percentile            BCa          
## 95%   ( 0.3439,  0.7443 )   ( 0.2760,  0.7160 )  
## Calculations and Intervals on Original Scale
#Complement en ligne, calcul exact du test de permutatio avec un tres petit echantillon
library(mvtnorm)
set.seed(1116)
sigma <- matrix(c(4,3,3,3), ncol=2)
petit_ech <- rmvnorm(n=8, mean=c(1,2), sigma=sigma)
petit_ech
##            [,1]     [,2]
## [1,]  1.8792081 2.477339
## [2,]  2.2013862 3.317487
## [3,]  1.5618803 2.526577
## [4,]  1.6926084 2.537411
## [5,]  4.9256143 3.866493
## [6,]  3.2565043 3.884740
## [7,]  2.1189746 3.725687
## [8,] -0.5696987 1.727317
cor(petit_ech[,1],petit_ech[,2])
## [1] 0.8562032
#Si le package n'est pas installe, enlever le commentaire
#puis executer la commande ci-dessous.
#if(!require("combinat")){install.packages("combinat")}
library(combinat)
## 
## Attaching package: 'combinat'
## The following object is masked from 'package:utils':
## 
##     combn
spcor <- sapply(permn(petit_ech[,1]), y=petit_ech[,2], 
                method="pearson", cor)
mean(abs(spcor)>=abs(cor(petit_ech[,1], petit_ech[,2])))
## [1] 0.0005704365
#page 56
Mes5_red_gv <- Mesures5[Mesures5$espece=="glycine violette",-5]

sumMasse=sum(Mes5_red_gv$masse)
sumTaille=sum(Mes5_red_gv$taille)
sumGraines=sum(Mes5_red_gv$graines)
sumMasseSec=sum(Mes5_red_gv$masse_sec)

#page 57
MatCov=matrix(NA,4,4,dimnames = list(c("masse","taille",
       "graines","masse_sec"),c("masse","taille","graines",
       "masse_sec")))
MatCov[1,2] <- MatCov[2,1] <- 1/(nrow(Mes5_red_gv)-1)*
  (sum(Mes5_red_gv$masse*Mes5_red_gv$taille)-
       1/nrow(Mes5_red_gv)*sumMasse*sumTaille)
MatCov[1,3] <- MatCov[3,1] <- 1/(nrow(Mes5_red_gv)-1)*
  (sum(Mes5_red_gv$masse*Mes5_red_gv$graines)-
       1/nrow(Mes5_red_gv)*sumMasse*sumGraines)
MatCov[1,4] <- MatCov[4,1] <- 1/(nrow(Mes5_red_gv)-1)*
  (sum(Mes5_red_gv$masse*Mes5_red_gv$masse_sec)-
       1/nrow(Mes5_red_gv)*sumMasse*sumMasseSec)
MatCov[2,3] <- MatCov[3,2] <- 1/(nrow(Mes5_red_gv)-1)*
  (sum(Mes5_red_gv$taille*Mes5_red_gv$graines)-
       1/nrow(Mes5_red_gv)*sumTaille*sumGraines)
MatCov[2,4] <- MatCov[4,2] <- 1/(nrow(Mes5_red_gv)-1)*
  (sum(Mes5_red_gv$taille*Mes5_red_gv$masse_sec)-
       1/nrow(Mes5_red_gv)*sumTaille*sumMasseSec)
MatCov[3,4] <- MatCov[4,3] <- 1/(nrow(Mes5_red_gv)-1)*
  (sum(Mes5_red_gv$graines*Mes5_red_gv$masse_sec)-
       1/nrow(Mes5_red_gv)*sumGraines*sumMasseSec)
MatCov[1,1] <- 1/(nrow(Mes5_red_gv)-1)*
  (sum(Mes5_red_gv$masse^2)-1/nrow(Mes5_red_gv)*
       sumMasse^2)
MatCov[2,2] <- 1/(nrow(Mes5_red_gv)-1)*
  (sum(Mes5_red_gv$taille^2)-1/nrow(Mes5_red_gv)*
       sumTaille^2)
MatCov[3,3] <- 1/(nrow(Mes5_red_gv)-1)*
  (sum(Mes5_red_gv$graines^2)-1/nrow(Mes5_red_gv)*
       sumGraines^2)
MatCov[4,4] <- 1/(nrow(Mes5_red_gv)-1)*
  (sum(Mes5_red_gv$masse_sec^2)-1/nrow(Mes5_red_gv)*
       sumMasseSec^2)
MatCov
##               masse    taille   graines masse_sec
## masse     111.60654 50.942597 11.949643  37.56280
## taille     50.94260 24.376987  5.583896  17.16522
## graines    11.94964  5.583896  1.906169   4.15513
## masse_sec  37.56280 17.165221  4.155130  13.06861
#page 58
cov(Mes5_red_gv)
##               masse    taille   graines masse_sec
## masse     111.60654 50.942597 11.949643  37.56280
## taille     50.94260 24.376987  5.583896  17.16522
## graines    11.94964  5.583896  1.906169   4.15513
## masse_sec  37.56280 17.165221  4.155130  13.06861
matsisj = outer(sqrt(diag(MatCov)),sqrt(diag(MatCov)),"*")
matsisj
##               masse    taille   graines masse_sec
## masse     111.60654 52.159671 14.585641  38.19087
## taille     52.15967 24.376987  6.816645  17.84862
## graines    14.58564  6.816645  1.906169   4.99109
## masse_sec  38.19087 17.848623  4.991090  13.06861
MatCov/matsisj
##               masse    taille   graines masse_sec
## masse     1.0000000 0.9766664 0.8192745 0.9835545
## taille    0.9766664 1.0000000 0.8191560 0.9617112
## graines   0.8192745 0.8191560 1.0000000 0.8325095
## masse_sec 0.9835545 0.9617112 0.8325095 1.0000000
cor(Mes5_red_gv)
##               masse    taille   graines masse_sec
## masse     1.0000000 0.9766664 0.8192745 0.9835545
## taille    0.9766664 1.0000000 0.8191560 0.9617112
## graines   0.8192745 0.8191560 1.0000000 0.8325095
## masse_sec 0.9835545 0.9617112 0.8325095 1.0000000
Matcor = cor(Mes5_red_gv)
#Si le package n'est pas installe, enlever le commentaire
#puis executer la commande ci-dessous.
#if(!require("corrplot")){install.packages("corrplot")}
library(corrplot)
## corrplot 0.84 loaded
corrplot(Matcor)

#page 59
corrplot(Matcor, method="circle", title="M\u00e9thode = circle",
         mar = c(0, 0, 4, 0))

corrplot(Matcor, method="square", title="M\u00e9thode = square",
         mar = c(0, 0, 4, 0))

#page 60
corrplot(Matcor, method="ellipse", title="M\u00e9thode = ellipse",
         mar = c(0, 0, 4, 0))

corrplot(Matcor, method="number", title="M\u00e9thode = number",
         mar = c(0, 0, 4, 0))

corrplot(Matcor, method="shade", title="M\u00e9thode = shade", 
         mar = c(0, 0, 4, 0))

corrplot(Matcor, method="color", title="M\u00e9thode = color", 
         mar = c(0, 0, 4, 0))

#page 61
corrplot(Matcor, method="pie", title="M\u00e9thode = pie", 
         mar = c(0, 0, 4, 0))

corrplot(Matcor, method="pie", type="full", title="Type = full", 
         mar = c(0, 0, 4, 0))

corrplot(Matcor, method="pie", type="lower", title="Type = lower", 
         mar = c(0, 0, 4, 0))

corrplot(Matcor, method="pie", type="upper", title="Type = upper", 
         mar = c(0, 0, 4, 0))

corrplot(Matcor, method="pie", diag=FALSE, title="Diag = FALSE", 
         mar = c(0, 0, 4, 0))

#page 62
corrplot(Matcor, method="pie", diag=FALSE, addCoef.col="orange")
corrplot(Matcor, method="pie", diag=FALSE, addCoef.col="orange", 
         number.cex=.75)

#Si le package n'est pas installe, enlever le commentaire
#puis executer la commande ci-dessous.
#if(!require("RColorBrewer")){install.packages("RColorBrewer")}
library(RColorBrewer)
corrplot(Matcor, method="pie", diag=FALSE, col=brewer.pal(n=8, 
    name="PuOr"))

#page 63
corrplot(Matcor, method="pie", diag=FALSE, tl.srt=45)

#Si le package n'est pas installe, enlever le commentaire
#puis executer la commande ci-dessous.
#if(!require("ggcorrplot")){install.packages("ggcorrplot")}
library(ggcorrplot)
ggcorrplot(Matcor)

ggcorrplot(Matcor, method="circle", type ="lower")

#page 64
library(GGally)
ggcorr(Mes5_red_gv)

#page 65
Mes5_red_gv = subset(Mesures5[,-5],subset=Mesures5$espece==
                       "glycine violette")
mvn(data = Mes5_red_gv[,c("masse","taille","masse_sec")],
      mvnTest = "mardia")$multivariateNormality
##              Test        Statistic            p value Result
## 1 Mardia Skewness 235.670663794595 5.552943391373e-45     NO
## 2 Mardia Kurtosis 18.8407839472993                  0     NO
## 3             MVN             <NA>               <NA>     NO
#page 66
Mes5_red_lr = subset(Mesures5[,-5],subset=Mesures5$espece==
                       "laurier rose")
nrow(Mes5_red_lr)
## [1] 72
mvn(data = Mes5_red_lr[,c("masse","taille","masse_sec")],
    mvnTest = "mardia", univariateTest = "SW",
    univariatePlot = "histogram", multivariateOutlierMethod =
      "adj", showOutliers = TRUE)$multivariateNormality

##              Test         Statistic           p value Result
## 1 Mardia Skewness  9.68412009567582 0.468629772370085    YES
## 2 Mardia Kurtosis -1.49587064286049 0.134687367872201    YES
## 3             MVN              <NA>              <NA>    YES
#page 67
tousalafois=cor.mtest(Mes5_red_lr[,c("masse","taille",
                                     "masse_sec")])
tousalafois$p
##              [,1]         [,2]         [,3]
## [1,] 0.000000e+00 2.516155e-20 2.928949e-10
## [2,] 2.516155e-20 0.000000e+00 3.913295e-04
## [3,] 2.928949e-10 3.913295e-04 0.000000e+00
#page 68
tousalafois$p<.05/3
##      [,1] [,2] [,3]
## [1,] TRUE TRUE TRUE
## [2,] TRUE TRUE TRUE
## [3,] TRUE TRUE TRUE
mvn(data = Mes5_red_gv[,c("masse","taille")], mvnTest =
      "mardia")$multivariateNormality
##              Test        Statistic              p value Result
## 1 Mardia Skewness 43.1097790148276 9.81908756021404e-09     NO
## 2 Mardia Kurtosis  3.4626491853959 0.000534885297369447     NO
## 3             MVN             <NA>                 <NA>     NO
mvn(data = Mes5_red_gv[,c("masse","masse_sec")], mvnTest =
      "mardia")$multivariateNormality
##              Test        Statistic              p value Result
## 1 Mardia Skewness 207.002931056865 1.17221434869874e-43     NO
## 2 Mardia Kurtosis 25.0465246005361                    0     NO
## 3             MVN             <NA>                 <NA>     NO
#page 69
mvn(data = Mes5_red_gv[,c("taille","masse_sec")], mvnTest =
      "mardia")$multivariateNormality
##              Test        Statistic              p value Result
## 1 Mardia Skewness 57.4923156600727 9.75288619274453e-12     NO
## 2 Mardia Kurtosis 6.38914935544299  1.6681100944993e-10     NO
## 3             MVN             <NA>                 <NA>     NO
library(ModStatR)
# La fonction perm.cor.mtest a ete integree `a la bibliotheque ModStatR
perm.cor.mtest
## function(mat, alternative= "two.sided", 
##                           method= "pearson", num.sim= 20000, ...) 
##   {
##   mat <- as.matrix(mat)
##   n <- ncol(mat)
##   p.mat <- cor.mat <- matrix(NA, n, n)
##   diag(p.mat) <- 0
##   diag(cor.mat) <- 1
##   for (i in 1:(n-1)) {
##     for (j in (i + 1):n) {
##       tmp <- jmuOutlier::perm.cor.test(x = mat[, i], y = mat[, j], 
##                            alternative=alternative, method=method, 
##                            num.sim=num.sim)
##       p.mat[i, j] <- p.mat[j, i] <- tmp$p.value
##       cor.mat[i, j] <- cor.mat[j, i] <- cor(x = mat[, i],
##                                             y = mat[, j], method=method, ...)
##     }
##   }
##   list(p = p.mat, cor = cor.mat)
## }
## <bytecode: 0x7fdcdad2de10>
## <environment: namespace:ModStatR>
permtoutalafois <- perm.cor.mtest(Mes5_red_gv)
permtoutalafois
## $p
##      [,1] [,2] [,3] [,4]
## [1,]    0    0    0    0
## [2,]    0    0    0    0
## [3,]    0    0    0    0
## [4,]    0    0    0    0
## 
## $cor
##           [,1]      [,2]      [,3]      [,4]
## [1,] 1.0000000 0.9766664 0.8192745 0.9835545
## [2,] 0.9766664 1.0000000 0.8191560 0.9617112
## [3,] 0.8192745 0.8191560 1.0000000 0.8325095
## [4,] 0.9835545 0.9617112 0.8325095 1.0000000
#page 70
permtoutalafois$p<.05/6
##      [,1] [,2] [,3] [,4]
## [1,] TRUE TRUE TRUE TRUE
## [2,] TRUE TRUE TRUE TRUE
## [3,] TRUE TRUE TRUE TRUE
## [4,] TRUE TRUE TRUE TRUE
cor(Mes5_red_gv)
##               masse    taille   graines masse_sec
## masse     1.0000000 0.9766664 0.8192745 0.9835545
## taille    0.9766664 1.0000000 0.8191560 0.9617112
## graines   0.8192745 0.8191560 1.0000000 0.8325095
## masse_sec 0.9835545 0.9617112 0.8325095 1.0000000
#page 71
nrow(na.omit(Mes5_red_lr[,c("taille","masse_sec")]))
## [1] 69
ref.cor.mtest(Mes5_red_lr[,c("masse","taille","masse_sec")],0.7)
## $p
##             [,1]         [,2]         [,3]
## [1,] 0.000000000 0.0036869954 0.6190419077
## [2,] 0.003686995 0.0000000000 0.0005226388
## [3,] 0.619041908 0.0005226388 0.0000000000
## 
## $cor
##           [,1]      [,2]      [,3]
## [1,] 1.0000000 0.8407223 0.6705079
## [2,] 0.8407223 1.0000000 0.4150282
## [3,] 0.6705079 0.4150282 1.0000000
## 
## $n
##      [,1] [,2] [,3]
## [1,]   72   72   69
## [2,]   72   72   69
## [3,]   69   69   69
#page 72
rescormest <- cor.mtest(Mes5_red_lr[,c("masse","taille","masse_sec")])
rescormest$lowCI
##           [,1]      [,2]      [,3]
## [1,] 1.0000000 0.7563713 0.5156606
## [2,] 0.7563713 1.0000000 0.1977747
## [3,] 0.5156606 0.1977747 1.0000000
rescormest$uppCI
##           [,1]      [,2]      [,3]
## [1,] 1.0000000 0.8975717 0.7829391
## [2,] 0.8975717 1.0000000 0.5934179
## [3,] 0.7829391 0.5934179 1.0000000
#page 73
rescormestbonf <- cor.mtest(Mes5_red_lr[,c("masse","taille", 
                  "masse_sec")],conf.level = 1-0.05/3)
rescormestbonf$lowCI
##           [,1]      [,2]      [,3]
## [1,] 1.0000000 0.7331158 0.4753710
## [2,] 0.7331158 1.0000000 0.1459423
## [3,] 0.4753710 0.1459423 1.0000000
rescormestbonf$uppCI
##           [,1]      [,2]      [,3]
## [1,] 1.0000000 0.9072638 0.8027662
## [2,] 0.9072638 1.0000000 0.6269343
## [3,] 0.8027662 0.6269343 1.0000000
Matcor_lr=cor(Mes5_red_lr[,c("masse","taille","masse_sec")],
                use="pairwise.complete.obs")
par(ask = FALSE)
for (i in seq(0.1, 0, -0.005)) {
  tmp <- cor.mtest(Mes5_red_lr[,c("masse","taille","masse_sec")],
                     conf.level = 1-i)
  corrplot(Matcor_lr, p.mat = tmp$p, low = tmp$lowCI,
             upp = tmp$uppCI, pch.col = "red", sig.level = i, plotCI =
               "rect", cl.pos = "n", title = substitute(alpha == x,
                                                          list(x = format(i, digits = 3, nsmall = 3))),
             mar = c(0, 0, 4, 0))
  Sys.sleep(0.15)
}

#page 75
#La fonction rho.mult est disponible dans le package ModStatR
ModStatR::rho.mult
## function(mat, indices){
##   tempcor <- as.vector(cor(x = mat[indices, ], method="pearson", 
##                            use="pairwise.complete.obs"))[outer(1:ncol(mat),1:ncol(mat),">")]
##   return(tempcor)
## }
## <bytecode: 0x7fdcdf7b76f0>
## <environment: namespace:ModStatR>
library(boot)
boot.mcor <- boot(Mes5_red_lr[,c("masse","taille","masse_sec")],
                    rho.mult, R=10000)
boot.mcor
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = Mes5_red_lr[, c("masse", "taille", "masse_sec")], 
##     statistic = rho.mult, R = 10000)
## 
## 
## Bootstrap Statistics :
##      original        bias    std. error
## t1* 0.8407223 -0.0009012164  0.02889353
## t2* 0.6705079 -0.0055613930  0.06953622
## t3* 0.4150282 -0.0056937945  0.09105317
#page 76
#La fonction boot.mcor.ic est disponible dans le package ModStatR
ModStatR::boot.mcor.ic
## function(mat, boot.mcor.res, conflevel = 0.95){
##   indboot.ci=function(i,type, conf = 0.95){
##     if(type=="perc") 
##       return(boot.ci(boot.mcor.res,index=i,type=type,conf=conf)$perc)
##     if(type=="bca")
##       return(boot.ci(boot.mcor.res,index=i,type=type,conf=conf)$bca)
##   }
##   bootperc <- matrix(unlist(lapply(1:ncol(mat),indboot.ci,type= 
##                                      "perc",conf=conflevel)),ncol=5,byrow=TRUE) 
##   bootabc <- matrix(unlist(lapply(1:ncol(mat),indboot.ci,type= 
##                                     "bca",conf=conflevel)),ncol=5,byrow=TRUE) 
##   cor.ic.percentile.low=matrix(NA,ncol(mat),ncol(mat)) 
##   diag(cor.ic.percentile.low) <- 1 
##   cor.ic.percentile.low[outer(1:ncol(mat),1:ncol(mat),">")]<-bootperc[,4] 
##   cor.ic.percentile.low[outer(1:ncol(mat),1:ncol(mat),"<")]<-bootperc[,4] 
##   cor.ic.percentile.up=matrix(NA,ncol(mat),ncol(mat)) 
##   diag(cor.ic.percentile.up) <- 1 
##   cor.ic.percentile.up[outer(1:ncol(mat),1:ncol(mat),">")]<-bootperc[,5] 
##   cor.ic.percentile.up[outer(1:ncol(mat),1:ncol(mat),"<")]<-bootperc[,5] 
##   cor.ic.BCa.low=matrix(NA,ncol(mat),ncol(mat)) 
##   diag(cor.ic.BCa.low) <- 1 
##   cor.ic.BCa.low[outer(1:ncol(mat),1:ncol(mat),">")]<-bootabc[,4] 
##   cor.ic.BCa.low[outer(1:ncol(mat),1:ncol(mat),"<")]<-bootabc[,4] 
##   cor.ic.BCa.up=matrix(NA,ncol(mat),ncol(mat)) 
##   diag(cor.ic.BCa.up) <- 1 
##   cor.ic.BCa.up[outer(1:ncol(mat),1:ncol(mat),">")]<- bootabc[,5] 
##   cor.ic.BCa.up[outer(1:ncol(mat),1:ncol(mat),"<")]<-bootabc[,5] 
##   return(list(cor.ic.percentile.low=cor.ic.percentile.low, 
##               cor.ic.percentile.up=cor.ic.percentile.up,cor.ic.BCa.low= 
##                 cor.ic.BCa.low,cor.ic.BCa.up=cor.ic.BCa.up)) 
## }
## <bytecode: 0x7fdc9b09f150>
## <environment: namespace:ModStatR>
boot.mcor.ic.res <- boot.mcor.ic(Mes5_red_lr[,c("masse","taille","masse_sec")],boot.mcor)
boot.mcor.ic.res
## $cor.ic.percentile.low
##           [,1]      [,2]      [,3]
## [1,] 1.0000000 0.7767109 0.5133714
## [2,] 0.7767109 1.0000000 0.2173860
## [3,] 0.5133714 0.2173860 1.0000000
## 
## $cor.ic.percentile.up
##           [,1]      [,2]      [,3]
## [1,] 1.0000000 0.8897671 0.7843293
## [2,] 0.8897671 1.0000000 0.5754015
## [3,] 0.7843293 0.5754015 1.0000000
## 
## $cor.ic.BCa.low
##           [,1]      [,2]      [,3]
## [1,] 1.0000000 0.7743566 0.5150255
## [2,] 0.7743566 1.0000000 0.2277593
## [3,] 0.5150255 0.2277593 1.0000000
## 
## $cor.ic.BCa.up
##           [,1]      [,2]      [,3]
## [1,] 1.0000000 0.8885415 0.7854277
## [2,] 0.8885415 1.0000000 0.5812631
## [3,] 0.7854277 0.5812631 1.0000000
corrplot(Matcor_lr, low = boot.mcor.ic.res$cor.ic.percentile.low,
         upp = boot.mcor.ic.res$cor.ic.percentile.up, pch.col = "red",
         plotCI = "rect", cl.pos = "n", mar = c(0, 0, 4, 0),
         title = "IC Bootstrap percentile 95%")

#page 77
corrplot(Matcor_lr, low = boot.mcor.ic.res$cor.ic.BCa.low, 
         upp = boot.mcor.ic.res$cor.ic.BCa.up, pch.col = "red", 
         plotCI = "rect", cl.pos = "n", mar = c(0, 0, 4, 0), 
         title = "IC Bootstrap BCa 95%")

boot.mcor.ic.res.bonf <- boot.mcor.ic(Mes5_red_lr[,c("masse",
   "taille","masse_sec")], boot.mcor, conflevel = 1-0.05/3)
boot.mcor.ic.res.bonf
## $cor.ic.percentile.low
##           [,1]      [,2]      [,3]
## [1,] 1.0000000 0.7610604 0.4708826
## [2,] 0.7610604 1.0000000 0.1709946
## [3,] 0.4708826 0.1709946 1.0000000
## 
## $cor.ic.percentile.up
##           [,1]      [,2]      [,3]
## [1,] 1.0000000 0.8982317 0.8072041
## [2,] 0.8982317 1.0000000 0.6047321
## [3,] 0.8072041 0.6047321 1.0000000
## 
## $cor.ic.BCa.low
##           [,1]      [,2]      [,3]
## [1,] 1.0000000 0.7590688 0.4751161
## [2,] 0.7590688 1.0000000 0.1817004
## [3,] 0.4751161 0.1817004 1.0000000
## 
## $cor.ic.BCa.up
##           [,1]      [,2]      [,3]
## [1,] 1.0000000 0.8973285 0.8093626
## [2,] 0.8973285 1.0000000 0.6177995
## [3,] 0.8093626 0.6177995 1.0000000
corrplot(Matcor_lr, mar = c(0, 0, 4, 0),
         low = boot.mcor.ic.res.bonf$cor.ic.percentile.low,
         upp = boot.mcor.ic.res.bonf$cor.ic.percentile.up,
         pch.col = "red", plotCI = "rect", cl.pos = "n", title =
         "IC Bootstrap percentile 95%\navec correction de Bonferroni")

#page 78
corrplot(Matcor_lr, mar = c(0, 0, 4, 0),
         low = boot.mcor.ic.res.bonf$cor.ic.BCa.low,
         upp = boot.mcor.ic.res.bonf$cor.ic.BCa.up,
         pch.col = "red", plotCI = "rect", cl.pos = "n", title =
           "IC Bootstrap BCa 95%\navec correction de Bonferroni")

#page 81
sigmaZ=cov(Mes5_red_gv)
sigmaZ
##               masse    taille   graines masse_sec
## masse     111.60654 50.942597 11.949643  37.56280
## taille     50.94260 24.376987  5.583896  17.16522
## graines    11.94964  5.583896  1.906169   4.15513
## masse_sec  37.56280 17.165221  4.155130  13.06861
sigma11=sigmaZ[1,1]
sigma11
## [1] 111.6065
sigma21=sigmaZ[1,2:4]
sigma21
##    taille   graines masse_sec 
##  50.94260  11.94964  37.56280
sigma22=sigmaZ[2:4,2:4]
sigma22
##              taille  graines masse_sec
## taille    24.376987 5.583896  17.16522
## graines    5.583896 1.906169   4.15513
## masse_sec 17.165221 4.155130  13.06861
#page 82
sigma22inverse=solve(sigma22)
sigma22inverse
##               taille    graines  masse_sec
## taille     0.5544019 -0.1196444 -0.6901493
## graines   -0.1196444  1.7350573 -0.3945076
## masse_sec -0.6901493 -0.3945076  1.1084420
corM=sqrt(t(sigma21)%*%sigma22inverse%*%sigma21/sigma11)
corM
##           [,1]
## [1,] 0.9900271
sqrt(summary(lm(masse~taille+graines+masse_sec,
                data=Mes5_red_gv))$r.squared)
## [1] 0.9900271
#page 84
library(MBESS)
Expected.R2(.5, 10, 5)
## [1] 0.7536454
Variance.R2(.5, 10, 5)
## [1] 0.02680003
#page 86
Mes5_red_lr = subset(Mesures5[,-5],subset=Mesures5$espece==
                       "laurier rose")
Rmultobs_lr=sqrt(summary(lm(masse~taille+masse_sec,
                              data=Mes5_red_lr))$r.squared)

n=nrow(na.omit(Mes5_red_lr[,c("masse","taille","masse_sec")]))
n
## [1] 69
#page 87
p=3
Fobs=(n-p)/(p-1)*Rmultobs_lr^2/(1-Rmultobs_lr^2)
Fobs
## [1] 151.3025
mvn(data = Mes5_red_lr[,c("masse","taille","masse_sec")], 
    mvnTest = "mardia")$multivariateNormality
##              Test         Statistic           p value Result
## 1 Mardia Skewness  9.68412009567582 0.468629772370085    YES
## 2 Mardia Kurtosis -1.49587064286049 0.134687367872201    YES
## 3             MVN              <NA>              <NA>    YES
n-p
## [1] 66
p-1
## [1] 2
pf(Fobs,p-1,n-p, lower.tail = FALSE)
## [1] 2.230659e-25
fstat <- summary(lm(masse~taille+masse_sec, data= 
                      Mes5_red_lr))$fstatistic
fstat
##    value    numdf    dendf 
## 151.3025   2.0000  66.0000
#page 88
pf(fstat[1], fstat[2], fstat[3], lower.tail=FALSE)
##        value 
## 2.230659e-25
summary(lm(masse~taille+masse_sec,data=Mes5_red_lr))
## 
## Call:
## lm(formula = masse ~ taille + masse_sec, data = Mes5_red_lr)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.98597 -0.29118 -0.06216  0.32882  1.02173 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.44846    0.37067  -3.908 0.000222 ***
## taille       0.33645    0.02876  11.700  < 2e-16 ***
## masse_sec    1.16109    0.16935   6.856 2.92e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4265 on 66 degrees of freedom
##   (3 observations deleted due to missingness)
## Multiple R-squared:  0.8209, Adjusted R-squared:  0.8155 
## F-statistic: 151.3 on 2 and 66 DF,  p-value: < 2.2e-16
.Machine$double.eps
## [1] 2.220446e-16
format.pval(pf(fstat[1L], fstat[2L], fstat[3L], lower.tail=FALSE))
## [1] "< 2.22e-16"
#page 89
library(MBESS)
ss.power.R2(Population.R2=.3^2, Specified.N=15, alpha.level=.05, 
            p=2)
## $Specified.Sample.Size
## [1] 15
## 
## $Actual.Power
## [1] 0.1472873
## 
## $Noncentral.F.Parm
## [1] 1.483516
## 
## $Effect.Size
## [1] 0.0989011
ss.power.R2(Population.R2=.3^2, alpha.level=.05, 
            desired.power=.80, p=2)
## $Necessary.Sample.Size
## [1] 101
## 
## $Actual.Power
## [1] 0.8022579
## 
## $Noncentral.F.Parm
## [1] 9.989011
## 
## $Effect.Size
## [1] 0.0989011
#page 90
Rmultobs_lr
## [1] 0.906061
tanh(atanh(Rmultobs_lr)-qnorm(.975)/sqrt(n-3))
## [1] 0.8521136
tanh(atanh(Rmultobs_lr)+qnorm(.975)/sqrt(n-3))
## [1] 0.9409563
ci.R(Rmultobs_lr,p-1,n-p)
## $Lower.Conf.Limit.R
## [1] 0.847737
## 
## $Prob.Less.Lower
## [1] 0.025
## 
## $Upper.Conf.Limit.R
## [1] 0.9391867
## 
## $Prob.Greater.Upper
## [1] 0.025
#page 91
ci.R2(Rmultobs_lr^2,p-1,n-p)
## $Lower.Conf.Limit.R2
## [1] 0.718658
## 
## $Prob.Less.Lower
## [1] 0.025
## 
## $Upper.Conf.Limit.R2
## [1] 0.8820717
## 
## $Prob.Greater.Upper
## [1] 0.025
#page 92
R_0=0.78
n=69
sqrt(n-3)*(atanh(Rmultobs_lr)-atanh(R_0))
## [1] 3.734649
2*(1-pnorm(abs(sqrt(n-3)*(atanh(Rmultobs_lr)-atanh(R_0)))))
## [1] 0.0001879769
#page 94
m.cov = rbind(c(3,1,1,0), c(1,3,0,1), c(1,0,2,0), c(0,1,0,2))

m.cor.1 = cov2cor(m.cov)
m.cor.1
##           [,1]      [,2]      [,3]      [,4]
## [1,] 1.0000000 0.3333333 0.4082483 0.0000000
## [2,] 0.3333333 1.0000000 0.0000000 0.4082483
## [3,] 0.4082483 0.0000000 1.0000000 0.0000000
## [4,] 0.0000000 0.4082483 0.0000000 1.0000000
#Si le package n'est pas installe, enlever le commentaire
#puis executer la commande ci-dessous.
#if(!require("corpcor")){install.packages("corpcor")}
library("corpcor")
m.pcor.1 = cor2pcor(m.cor.1)
m.pcor.1
##            [,1]       [,2]        [,3]        [,4]
## [1,]  1.0000000  0.4000000  0.43852901 -0.17541160
## [2,]  0.4000000  1.0000000 -0.17541160  0.43852901
## [3,]  0.4385290 -0.1754116  1.00000000  0.07692308
## [4,] -0.1754116  0.4385290  0.07692308  1.00000000
#page 95
m.pcor.2 = cor2pcor(m.cov)
m.pcor.2
##            [,1]       [,2]        [,3]        [,4]
## [1,]  1.0000000  0.4000000  0.43852901 -0.17541160
## [2,]  0.4000000  1.0000000 -0.17541160  0.43852901
## [3,]  0.4385290 -0.1754116  1.00000000  0.07692308
## [4,] -0.1754116  0.4385290  0.07692308  1.00000000
m.cor.2 = pcor2cor(m.pcor.1)
m.cor.2
##               [,1]         [,2]          [,3]         [,4]
## [1,]  1.000000e+00 3.333333e-01  4.082483e-01 6.678397e-16
## [2,]  3.333333e-01 1.000000e+00  3.756598e-16 4.082483e-01
## [3,]  4.082483e-01 2.504399e-16  1.000000e+00 5.828671e-16
## [4,] -2.086999e-16 4.082483e-01 -2.690156e-16 1.000000e+00
z.m.cor.1 <- zapsmall(m.cor.1)
z.m.cor.2 <- zapsmall(m.cor.2)
z.m.cor.1 == z.m.cor.2
##      [,1] [,2] [,3] [,4]
## [1,] TRUE TRUE TRUE TRUE
## [2,] TRUE TRUE TRUE TRUE
## [3,] TRUE TRUE TRUE TRUE
## [4,] TRUE TRUE TRUE TRUE
#Complement en ligne pour illustrer la definition 3
m.precision = solve(m.cov)
m.inv.cor = solve(cov2cor(m.cov))
#Meme resultat `a partir de l'inversion de la matrice de variance-covariance 
#ou `a partir de la matrice de correlation lineaire
-cov2cor(m.precision)
##            [,1]       [,2]        [,3]        [,4]
## [1,] -1.0000000  0.4000000  0.43852901 -0.17541160
## [2,]  0.4000000 -1.0000000 -0.17541160  0.43852901
## [3,]  0.4385290 -0.1754116 -1.00000000  0.07692308
## [4,] -0.1754116  0.4385290  0.07692308 -1.00000000
-cov2cor(m.inv.cor)
##            [,1]       [,2]        [,3]        [,4]
## [1,] -1.0000000  0.4000000  0.43852901 -0.17541160
## [2,]  0.4000000 -1.0000000 -0.17541160  0.43852901
## [3,]  0.4385290 -0.1754116 -1.00000000  0.07692308
## [4,] -0.1754116  0.4385290  0.07692308 -1.00000000
#Ce resultat est le meme que celui de la fonction cor2pcor
zapsmall(m.pcor.1) == zapsmall(-cov2cor(m.precision))
##       [,1]  [,2]  [,3]  [,4]
## [1,] FALSE  TRUE  TRUE  TRUE
## [2,]  TRUE FALSE  TRUE  TRUE
## [3,]  TRUE  TRUE FALSE  TRUE
## [4,]  TRUE  TRUE  TRUE FALSE
#page 96
Mes5_red_lr = subset(Mesures5[,-5],subset=Mesures5$espece== 
                       "laurier rose")
with(Mes5_red_lr, cor.test(masse,masse_sec))
## 
##  Pearson's product-moment correlation
## 
## data:  masse and masse_sec
## t = 7.3977, df = 67, p-value = 2.929e-10
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.5156606 0.7829391
## sample estimates:
##       cor 
## 0.6705079
#page 97
Mes5_red_lr_noNA <- na.omit(Mes5_red_lr[,-3])
sigmaTriplet = cov(Mes5_red_lr_noNA)

sigma11=sigmaTriplet[c(1,3),c(1,3)]
sigma11
##               masse masse_sec
## masse     0.9860401 0.2235081
## masse_sec 0.2235081 0.1126897
sigma21=sigmaTriplet[2,c(1,3),drop=FALSE]
sigma21
##           masse masse_sec
## taille 1.634642  0.275422
sigma22=sigmaTriplet[2,2,drop=FALSE]
sigma22
##          taille
## taille 3.908031
#page 98
par_covar=sigma11-t(sigma21)%*%solve(sigma22)%*%sigma21
par_covar
##               masse  masse_sec
## masse     0.3023059 0.10830523
## masse_sec 0.1083052 0.09327907
cov2cor(par_covar)
##               masse masse_sec
## masse     1.0000000 0.6449619
## masse_sec 0.6449619 1.0000000
res1 = residuals(lm(masse~taille, data=Mes5_red_lr_noNA)) 
res2 = residuals(lm(masse_sec~taille, data=Mes5_red_lr_noNA)) 
cor(res1,res2)
## [1] 0.6449619
corTriplet <- cov2cor(sigmaTriplet)
corTriplet
##               masse    taille masse_sec
## masse     1.0000000 0.8327150 0.6705079
## taille    0.8327150 1.0000000 0.4150282
## masse_sec 0.6705079 0.4150282 1.0000000
#page 99
(corTriplet[3,1]-corTriplet[2,1]*corTriplet[3,2])/
  sqrt(1-corTriplet[2,1]^2)/sqrt(1-corTriplet[3,2]^2)
## [1] 0.6449619
cor2pcor(sigmaTriplet)
##           [,1]       [,2]       [,3]
## [1,] 1.0000000  0.8213989  0.6449619
## [2,] 0.8213989  1.0000000 -0.3488714
## [3,] 0.6449619 -0.3488714  1.0000000
#Si le package n'est pas installe, enlever le commentaire
#puis executer la commande ci-dessous.
#if(!require("ppcor")){install.packages("ppcor")}
library(ppcor)
## Loading required package: MASS
pcor(Mes5_red_lr_noNA)$estimate
##               masse     taille  masse_sec
## masse     1.0000000  0.8213989  0.6449619
## taille    0.8213989  1.0000000 -0.3488714
## masse_sec 0.6449619 -0.3488714  1.0000000
#page 102
pcor(Mes5_red_lr_noNA)
## $estimate
##               masse     taille  masse_sec
## masse     1.0000000  0.8213989  0.6449619
## taille    0.8213989  1.0000000 -0.3488714
## masse_sec 0.6449619 -0.3488714  1.0000000
## 
## $p.value
##                  masse       taille    masse_sec
## masse     0.000000e+00 9.516717e-18 2.915068e-09
## taille    9.516717e-18 0.000000e+00 3.548483e-03
## masse_sec 2.915068e-09 3.548483e-03 0.000000e+00
## 
## $statistic
##               masse    taille masse_sec
## masse      0.000000 11.699880  6.856312
## taille    11.699880  0.000000 -3.024256
## masse_sec  6.856312 -3.024256  0.000000
## 
## $n
## [1] 69
## 
## $gp
## [1] 1
## 
## $method
## [1] "pearson"
pcor(Mes5_red_lr_noNA)$p.value<5/3
##           masse taille masse_sec
## masse      TRUE   TRUE      TRUE
## taille     TRUE   TRUE      TRUE
## masse_sec  TRUE   TRUE      TRUE
#page 103
r_pcor<-pcor(Mes5_red_lr_noNA)
r_pcr_e<-r_pcor$estimate; n<-r_pcor$n; gp<-r_pcor$gp; N<-n-gp-3
tanh(atanh(r_pcr_e)-qnorm(.975)/sqrt(n-gp-3))
##               masse     taille  masse_sec
## masse     1.0000000  0.7249494  0.4804122
## taille    0.7249494  1.0000000 -0.5421967
## masse_sec 0.4804122 -0.5421967  1.0000000
tanh(atanh(r_pcr_e)+qnorm(.975)/sqrt(n-gp-3))
##               masse     taille  masse_sec
## masse     1.0000000  0.8862574  0.7656492
## taille    0.8862574  1.0000000 -0.1204666
## masse_sec 0.7656492 -0.1204666  1.0000000
#page 104
icl<-tanh(atanh(r_pcr_e)-qnorm(.975)/sqrt(N))-r_pcr_e/(2*(N+2))
icu<-tanh(atanh(r_pcr_e)+qnorm(.975)/sqrt(N))-r_pcr_e/(2*(N+2))

sqrt(n-gp-2)*r_pcr_e/sqrt(1-r_pcr_e^2)
##               masse    taille masse_sec
## masse           Inf 11.699880  6.856312
## taille    11.699880       Inf -3.024256
## masse_sec  6.856312 -3.024256       Inf
r_pcor$statistic
##               masse    taille masse_sec
## masse      0.000000 11.699880  6.856312
## taille    11.699880  0.000000 -3.024256
## masse_sec  6.856312 -3.024256  0.000000
2*(1-pt(sqrt(n-gp-2)*abs(r_pcr_e)/sqrt(1-r_pcr_e^2),n-gp-2))
##                  masse      taille    masse_sec
## masse     0.000000e+00 0.000000000 2.915068e-09
## taille    0.000000e+00 0.000000000 3.548483e-03
## masse_sec 2.915068e-09 0.003548483 0.000000e+00
r_pcor$p.value
##                  masse       taille    masse_sec
## masse     0.000000e+00 9.516717e-18 2.915068e-09
## taille    9.516717e-18 0.000000e+00 3.548483e-03
## masse_sec 2.915068e-09 3.548483e-03 0.000000e+00
library(MBESS)
cis=(sapply(abs(r_pcr_e),ci.R,df.1=1,df.2=N+1, simplify = TRUE))
colnames(cis) <- outer(rownames(r_pcr_e),colnames(r_pcr_e),paste)
cis
##                    masse masse taille masse masse_sec masse masse taille
## Lower.Conf.Limit.R 0.9999976   0.7224134    0.4775546       0.7224134   
## Prob.Less.Lower    0.025       0.025        0.025           0.025       
## Upper.Conf.Limit.R 0.9999976   0.8846597    0.7630296       0.8846597   
## Prob.Greater.Upper 0.025       0.025        0.025           0.025       
##                    taille taille masse_sec taille masse masse_sec
## Lower.Conf.Limit.R 0.9999976     0.1192247        0.4775546      
## Prob.Less.Lower    0.025         0.025            0.025          
## Upper.Conf.Limit.R 0.9999976     0.5391761        0.7630296      
## Prob.Greater.Upper 0.025         0.025            0.025          
##                    taille masse_sec masse_sec masse_sec
## Lower.Conf.Limit.R 0.1192247        0.9999976          
## Prob.Less.Lower    0.025            0.025              
## Upper.Conf.Limit.R 0.5391761        0.9999976          
## Prob.Greater.Upper 0.025            0.025
library(corrplot)
layout(t(1:2))
corrplot(r_pcor$estimate, title = "Corr\u00e9lation partielle", mar = 
           c(0, 0, 4, 0), cl.pos = "r")
corrplot(r_pcor$estimate, low = icl, upp = icu, pch.col = "red", 
         plotCI = "rect", cl.pos = "n", title = 
           "Corr\u00e9lation partielle et IC", mar = c(0, 0, 4, 0))

layout(1)

#page 105
Matcor_lr_noNA=cor(Mes5_red_lr_noNA)
layout(t(1:2))
corrplot(Matcor_lr_noNA, title = "Corr\u00e9lation simple", mar = 
           c(0, 0, 4, 0), cl.pos = "r")
rescormest <- cor.mtest(Mes5_red_lr_noNA)
corrplot(Matcor_lr_noNA, low = rescormest$lowCI, 
         upp= rescormest$uppCI, pch.col = "red", plotCI = "rect", cl.pos= 
           "n", title = "Corr\u00e9lation simple et IC", mar = c(0, 0, 4, 0))

layout(1)

#page 111
library(mvtnorm)
set.seed(1116)
sigma <- matrix(c(4,3,3,3), ncol=2)
petit_ech <- rmvnorm(n=8, mean=c(1,2), sigma=sigma)
cor_spear = cor(petit_ech,method="spearman")
cor_spear
##           [,1]      [,2]
## [1,] 1.0000000 0.8809524
## [2,] 0.8809524 1.0000000
#page 112
cor.test(petit_ech[,1],petit_ech[,2],method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  petit_ech[, 1] and petit_ech[, 2]
## S = 10, p-value = 0.007242
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.8809524
#Si le package n'est pas installe, enlever le commentaire
#puis executer la commande ci-dessous.
#if(!require(pspearman)){install.packages("pspearman")}
library(pspearman)
spearman.test(petit_ech[,1],petit_ech[,2])
## 
##  Spearman's rank correlation rho
## 
## data:  petit_ech[, 1] and petit_ech[, 2]
## S = 10, p-value = 0.007242
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.8809524
2*pspearman(10, 8)
## [1] 0.007242063
#page 113
#Si le package n'est pas installe, enlever le commentaire
#puis executer la commande ci-dessous.
#if(!require(SuppDists)){install.packages("SuppDists")}
library(SuppDists)
2*pSpearman(cor_spear[2,1], 8, lower.tail = FALSE)
## [1] 0.004563492
#Si le package n'est pas installe, enlever le commentaire
#puis executer la commande ci-dessous.
#if(!require(coin)){install.packages("coin")}
library(coin)
## Loading required package: survival
## 
## Attaching package: 'survival'
## The following object is masked from 'package:boot':
## 
##     aml
spearman_test(petit_ech[,1]~petit_ech[,2],distribution= 
                approximate(nresample=99999))
## 
##  Approximative Spearman Correlation Test
## 
## data:  petit_ech[, 1] by petit_ech[, 2]
## Z = 2.3308, p-value = 0.00725
## alternative hypothesis: true rho is not equal to 0
library(combinat)
spcor <- sapply(permn(petit_ech[,1]), y=petit_ech[,2], method= 
                  "spearman", cor)
mean(abs(spcor)>=abs(cor(petit_ech[,1],petit_ech[,2], method= 
                           "spearman")))
## [1] 0.007242063
#page 114
data(Quetelet, package="BioStatR")
with(Quetelet,plot(taille,poids,col=c("red","blue")[sexe],pch= 
                     c(17,15)[sexe], main="Masse en fonction de la taille"))
legend(x="bottomright",legend = c("Femme","Homme"),pch= 
         c(17,15), col=c("red","blue"))

with(Quetelet,spearman_test(taille~ poids, distribution=
                              approximate(nresample=9999)))
## 
##  Approximative Spearman Correlation Test
## 
## data:  taille by poids
## Z = 6.8008, p-value < 1e-04
## alternative hypothesis: true rho is not equal to 0
#page 116
pcor(Mes5_red_lr_noNA, method = "spearman")
## $estimate
##               masse     taille  masse_sec
## masse     1.0000000  0.8030674  0.5460280
## taille    0.8030674  1.0000000 -0.2354607
## masse_sec 0.5460280 -0.2354607  1.0000000
## 
## $p.value
##                  masse       taille    masse_sec
## masse     0.000000e+00 1.750641e-16 1.463849e-06
## taille    1.750641e-16 0.000000e+00 5.324450e-02
## masse_sec 1.463849e-06 5.324450e-02 0.000000e+00
## 
## $statistic
##               masse    taille masse_sec
## masse      0.000000 10.948615  5.294972
## taille    10.948615  0.000000 -1.968231
## masse_sec  5.294972 -1.968231  0.000000
## 
## $n
## [1] 69
## 
## $gp
## [1] 1
## 
## $method
## [1] "spearman"
#page 117
cor(Mes5_red_lr_noNA, method = "spearman")
##               masse    taille masse_sec
## masse     1.0000000 0.8284072 0.6163290
## taille    0.8284072 1.0000000 0.4067113
## masse_sec 0.6163290 0.4067113 1.0000000
pcor(Mes5_red_lr_noNA, method = "spearman")$p.value<5/3
##           masse taille masse_sec
## masse      TRUE   TRUE      TRUE
## taille     TRUE   TRUE      TRUE
## masse_sec  TRUE   TRUE      TRUE
library(corrplot)
layout(t(1:2))
corrplot(r_pcor$estimate, mar = c(0, 0, 4, 0), cl.pos = "r", 
         title = "Corr\u00e9lation partielle\nde Spearman")
corrplot(cor(Mes5_red_lr_noNA, method="spearman"), mar=c(0,0, 
      4, 0), cl.pos = "r", title = "Corr\u00e9lation simple\nde Spearman")

layout(1)

#page 121
library(mvtnorm)
set.seed(1116)
sigma <- matrix(c(4,3,3,3), ncol=2)
petit_ech <- rmvnorm(n=8, mean=c(1,2), sigma=sigma)
cor_kend = cor(petit_ech,method="kendall")
cor_kend
##           [,1]      [,2]
## [1,] 1.0000000 0.7142857
## [2,] 0.7142857 1.0000000
cor.test(petit_ech[,1],petit_ech[,2],method="kendall")
## 
##  Kendall's rank correlation tau
## 
## data:  petit_ech[, 1] and petit_ech[, 2]
## T = 24, p-value = 0.01414
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
##       tau 
## 0.7142857
if(!require(Kendall)){install.packages(Kendall)}
## Loading required package: Kendall
library(Kendall)
Kendall(petit_ech[,1],petit_ech[,2])
## tau = 0.714, 2-sided pvalue =0.018741
#Si le package n'est pas installe, enlever le commentaire
#puis executer la commande ci-dessous.
#if(!require(SuppDists)){install.packages("SuppDists")}
library(SuppDists)
2*pKendall(cor_kend[2,1], 8, lower.tail = FALSE)
## [1] 0.005505952
#page 122
library(combinat)
spcor <- sapply(permn(petit_ech[,1]), y=petit_ech[,2], method= 
                  "kendall", cor)
mean(abs(spcor)>=abs(cor(petit_ech[,1],petit_ech[,2], method= 
                             "kendall")))
## [1] 0.0141369
with(Quetelet,Kendall(taille, poids))
## tau = 0.681, 2-sided pvalue =< 2.22e-16
#page 125
pcor(Mes5_red_lr_noNA, method = "kendall")
## $estimate
##               masse     taille  masse_sec
## masse     1.0000000  0.5977469  0.3946231
## taille    0.5977469  1.0000000 -0.0103710
## masse_sec 0.3946231 -0.0103710  1.0000000
## 
## $p.value
##                  masse       taille    masse_sec
## masse     0.000000e+00 5.684130e-13 1.950263e-06
## taille    5.684130e-13 0.000000e+00 9.004781e-01
## masse_sec 1.950263e-06 9.004781e-01 0.000000e+00
## 
## $statistic
##              masse     taille  masse_sec
## masse     0.000000  7.2078521  4.7585110
## taille    7.207852  0.0000000 -0.1250574
## masse_sec 4.758511 -0.1250574  0.0000000
## 
## $n
## [1] 69
## 
## $gp
## [1] 1
## 
## $method
## [1] "kendall"
#page 126
cor(Mes5_red_lr_noNA, method = "kendall")
##               masse    taille masse_sec
## masse     1.0000000 0.6461236 0.4845355
## taille    0.6461236 1.0000000 0.3061455
## masse_sec 0.4845355 0.3061455 1.0000000
pcor(Mes5_red_lr_noNA, method = "kendall")$p.value<5/3
##           masse taille masse_sec
## masse      TRUE   TRUE      TRUE
## taille     TRUE   TRUE      TRUE
## masse_sec  TRUE   TRUE      TRUE
#page 127
library(corrplot)
layout(t(1:2))
corrplot(r_pcor$estimate, mar = c(0, 0, 4, 0), cl.pos = "r", 
         title = "Corr\u00e9lation partielle\nde Kendall")
corrplot(cor(Mes5_red_lr_noNA, method="kendall"), mar=c(0,0, 
        4, 0), cl.pos = "r", title = "Corr\u00e9lation simple\nde Kendall")

#page 128, exercice 2
read.csv("https://tinyurl.com/y2c68uvw")
##     Maths Sport
## 1   13.12 14.44
## 2   12.43 14.32
## 3   16.26 15.85
## 4   12.61  7.68
## 5    3.54  6.12
## 6   11.50 13.08
## 7    6.43  5.71
## 8   15.79 12.91
## 9    8.02 14.32
## 10  13.86 10.76
## 11  15.57 14.58
## 12  19.43 17.05
## 13  14.76 13.35
## 14  10.03  7.27
## 15  13.74 15.60
## 16  12.01 11.16
## 17   9.06  9.59
## 18  10.26 14.21
## 19   6.86  9.01
## 20  11.41  9.10
## 21  11.96 13.93
## 22  10.19  6.61
## 23  10.19  9.75
## 24   7.83 11.11
## 25   9.63  9.24
## 26  10.67 12.61
## 27  13.10 12.67
## 28  11.99 11.50
## 29  15.23 14.41
## 30  19.42 19.10
## 31  13.29 12.00
## 32   8.85  8.42
## 33  13.34 12.15
## 34  11.61 11.76
## 35   8.72  8.17
## 36   9.03 12.82
## 37   7.27  8.69
## 38  10.82 12.64
## 39  15.25 16.28
## 40  10.78 12.83
## 41  11.60  9.68
## 42  12.43  9.93
## 43  10.60 11.30
## 44  11.98  7.43
## 45   9.78 12.03
## 46   9.88  8.70
## 47  15.15 11.93
## 48   7.54 10.71
## 49  11.30  7.42
## 50  12.20 12.90
## 51  12.88 12.47
## 52  14.99 15.75
## 53  12.25 13.83
## 54  11.79 11.33
## 55  11.66 12.33
## 56  12.17 13.50
## 57   8.61  9.77
## 58   4.70  8.02
## 59   5.99  3.62
## 60  12.26 10.35
## 61  12.01 12.87
## 62  12.89 10.76
## 63  10.34  9.34
## 64   5.48  6.90
## 65  11.99 14.00
## 66   8.93 10.80
## 67  15.30 13.46
## 68   8.04 10.32
## 69  16.82 16.21
## 70  12.63  9.43
## 71  12.87 11.87
## 72  12.64 11.50
## 73   6.32  6.74
## 74   9.03  5.18
## 75  17.12 15.43
## 76   7.84  8.64
## 77  11.18 11.90
## 78  13.56 16.16
## 79  13.00 14.25
## 80  10.90  9.90
## 81  16.33 14.32
## 82   6.83  9.09
## 83  11.81 10.88
## 84  15.98 14.01
## 85  13.38  8.91
## 86   9.65 15.28
## 87  10.54 13.41
## 88  11.38 12.39
## 89  12.37  9.19
## 90  11.38 14.09
## 91  12.29 14.93
## 92  13.26  8.13
## 93  14.65 16.29
## 94  15.52 13.06
## 95  19.95 15.08
## 96  14.76 16.23
## 97   8.95 12.39
## 98   8.64 10.85
## 99  12.04 13.99
## 100  8.60 10.57
## 101 10.81 11.37
## 102 15.34 11.84
## 103 12.18 12.26
## 104 11.40 13.31
## 105 12.21 10.92
## 106 10.31 10.95
## 107 11.06  5.82
## 108  8.03  5.72
## 109 15.61 17.34
## 110 11.01 11.00
## 111  7.92  8.70
## 112 14.04 11.20
## 113  6.64 10.32
## 114 10.28  8.43
## 115 10.52 14.37
## 116  8.91 10.55
## 117 14.98 15.75
## 118 10.21 10.94
## 119 11.78 10.18
#page 137
read.csv("https://tinyurl.com/y2asrzgk")
##     Maths Sport   Age
## 1   13.12 14.44 13.40
## 2   12.43 14.32 14.44
## 3   16.26 15.85 14.17
## 4   12.61  7.68 12.76
## 5    3.54  6.12 11.29
## 6   11.50 13.08 13.38
## 7    6.43  5.71 11.53
## 8   15.79 12.91 13.77
## 9    8.02 14.32 13.54
## 10  13.86 10.76 12.92
## 11  15.57 14.58 14.83
## 12  19.43 17.05 15.24
## 13  14.76 13.35 13.66
## 14  10.03  7.27 11.63
## 15  13.74 15.60 14.37
## 16  12.01 11.16 12.48
## 17   9.06  9.59 11.89
## 18  10.26 14.21 12.58
## 19   6.86  9.01 12.67
## 20  11.41  9.10 11.72
## 21  11.96 13.93 13.47
## 22  10.19  6.61 11.91
## 23  10.19  9.75 12.65
## 24   7.83 11.11 12.78
## 25   9.63  9.24 12.45
## 26  10.67 12.61 12.88
## 27  13.10 12.67 13.44
## 28  11.99 11.50 12.70
## 29  15.23 14.41 14.35
## 30  19.42 19.10 15.38
## 31  13.29 12.00 12.96
## 32   8.85  8.42 11.93
## 33  13.34 12.15 13.84
## 34  11.61 11.76 12.99
## 35   8.72  8.17 10.77
## 36   9.03 12.82 12.86
## 37   7.27  8.69 11.60
## 38  10.82 12.64 13.25
## 39  15.25 16.28 14.15
## 40  10.78 12.83 12.15
## 41  11.60  9.68 13.25
## 42  12.43  9.93 12.73
## 43  10.60 11.30 12.68
## 44  11.98  7.43 12.23
## 45   9.78 12.03 12.92
## 46   9.88  8.70 11.43
## 47  15.15 11.93 13.26
## 48   7.54 10.71 12.55
## 49  11.30  7.42 12.25
## 50  12.20 12.90 13.35
## 51  12.88 12.47 13.20
## 52  14.99 15.75 15.08
## 53  12.25 13.83 13.60
## 54  11.79 11.33 13.58
## 55  11.66 12.33 12.40
## 56  12.17 13.50 13.18
## 57   8.61  9.77 11.71
## 58   4.70  8.02 11.84
## 59   5.99  3.62 10.89
## 60  12.26 10.35 13.29
## 61  12.01 12.87 13.19
## 62  12.89 10.76 13.18
## 63  10.34  9.34 12.45
## 64   5.48  6.90 11.59
## 65  11.99 14.00 13.83
## 66   8.93 10.80 12.24
## 67  15.30 13.46 13.47
## 68   8.04 10.32 12.27
## 69  16.82 16.21 15.31
## 70  12.63  9.43 12.65
## 71  12.87 11.87 12.83
## 72  12.64 11.50 13.13
## 73   6.32  6.74 11.33
## 74   9.03  5.18 11.64
## 75  17.12 15.43 14.51
## 76   7.84  8.64 12.31
## 77  11.18 11.90 12.24
## 78  13.56 16.16 14.57
## 79  13.00 14.25 14.17
## 80  10.90  9.90 12.18
## 81  16.33 14.32 13.88
## 82   6.83  9.09 12.21
## 83  11.81 10.88 13.11
## 84  15.98 14.01 14.16
## 85  13.38  8.91 13.37
## 86   9.65 15.28 13.63
## 87  10.54 13.41 13.91
## 88  11.38 12.39 13.32
## 89  12.37  9.19 13.12
## 90  11.38 14.09 14.10
## 91  12.29 14.93 13.75
## 92  13.26  8.13 12.40
## 93  14.65 16.29 14.28
## 94  15.52 13.06 15.08
## 95  19.95 15.08 15.15
## 96  14.76 16.23 15.15
## 97   8.95 12.39 13.26
## 98   8.64 10.85 12.54
## 99  12.04 13.99 13.63
## 100  8.60 10.57 12.07
## 101 10.81 11.37 12.98
## 102 15.34 11.84 13.89
## 103 12.18 12.26 13.51
## 104 11.40 13.31 13.00
## 105 12.21 10.92 12.91
## 106 10.31 10.95 12.96
## 107 11.06  5.82 12.18
## 108  8.03  5.72 12.02
## 109 15.61 17.34 15.44
## 110 11.01 11.00 12.50
## 111  7.92  8.70 12.31
## 112 14.04 11.20 13.79
## 113  6.64 10.32 13.09
## 114 10.28  8.43 12.63
## 115 10.52 14.37 13.84
## 116  8.91 10.55 12.40
## 117 14.98 15.75 14.14
## 118 10.21 10.94 12.29
## 119 11.78 10.18 12.68