#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