1.1 Importation des packages

On a 2 méthodes pour importer les packages :

library(nom de package)
require(nom de package)

La différence est que require = try + library. C’est-à-dire, require a une sortie de type logical.

print(require(pacakge_n_existe_pas))

Alors, on conclut qu’on utilise toujours library pour importer un package existant. On n’utilise require que dans le cas :

if(!require(package))  install.packages('package')

Si on ne veut utiliser qu’une ou deux fonctions ou variables d’un package, il faut eviter (si possible) l’utilisation de library. En revanche, on peut directement écrire

package::fonction()
package::variable

qui ressemble beaucoup à

data.frame$nom_de_colone

1.2 Indexation par []

Il y a 3 types de index par [] dans R : int, logical, char

• int : les positions
• logical : les conditions
• char : les noms (par ex. les noms de colonne dans un dataframe)

1.3 Manipuler les données avec le package dplyr

En fait, dans la vraie vie, il est très rare qu’on manipule directement les données (par ex. data.frame), car ce n’est pas assez efficace si on a plein de données à traiter. Je vous donne une introduction au package dplyr, qui est assez facile à comprendre. Pour aller plus loin, il vaut mieux lire un peu la documentation officielle.

1.4 Un exemple d’utilisation de apply() pour TP4 et TP5

Ici, je vous donne un exemple de l’utilisation de apply() en respectant les exos des TP4 et TP5

rttrans <- function(n,q,theta){
  z = rt(n,q)
  x = z + theta 
        return(x)
}

rttrans_K <- function(n,q,theta,K){
  mat = vector()
  for (i in 1:K){
    mat = cbind(mat,rttrans(n,q,theta))
  }
    return(mat)
}

estimation_K <- function(mat, method){
  if (method == 'mean'){
    return(apply(mat,2,mean))
  }else if (method == 'median'){
    return(apply(mat,2,median))
  }else if (method == 'mean_trim'){
    return(apply(mat,2,function(mat) mean(mat,trim = 0.2)))
  }else{
    print('Wrong input for method.')
}
}

risques <- function(theta,hat.theta){
  return(mean((hat.theta - theta)^2))
}

######################################################################

n = 20
q = 5
theta = 1
K = 50

estimation_mean_trim <- estimation_K(rttrans_K(n,q,theta,K),'mean_trim')
risques(theta,estimation_mean_trim)

1.5 Sélection de variables par leaps

X = seq(-10,10,by = 0.1)
n = length(X)
sigma = 4

epsilon = rnorm(n,0,sigma)

f = function(x){
  return(0.02*x^3 + 0.5*x^2 + 5*sin(x))
}

Y = f(X) + epsilon

######################################################################

df_data <- data.frame(Y,X,X^2, X^3, X^4, X^5, X^6, sin(X), cos(X), exp(X))
model = lm(Y~., data = df_data)

######################################################################

library(leaps)
reg <- regsubsets(Y~., data = df_data, nbest = 1)
plot(reg, scale = "bic")

par(mfrow = c(1,3))
plot(reg, scale = "Cp")
plot(reg, scale = "r2")
plot(reg, scale = "adjr2")

model_bic <- lm(Y~cbind(X^2,X^3,sin(X)))
Y_predict = predict(model_bic,data.frame(X))
plot(X, f(X), type = 'l', lwd = 2, xlab = 'X', ylab = 'Y',col = 'darkred')
points(X, Y, type = 'b', lwd = 1, col = 'grey', pch = 1, lty = 2)
lines(X, Y_predict, lwd = 2)
legend('topleft', c('fonction à estimer','observations','estimation'), col = c('darkred','grey',1), lwd = c(2,1,2), pch = c(-1,1,-1), lty = c(1,2,1))
grid()

2.1 Stucture de .Rmd

• metadata (configuration)
• texte
• chunk de codes
• texte
• chunk de codes …

```{r nom1, include=FALSE}
seq(1, 10)
```

```{r nom2, echo=FALSE}
plot(rnorm(1000))
```

```{r nom3, echo=FALSE, fig.cap=fig$cap('label','description')}
plot(rnorm(1000))
```

`r fig$cap('label')`

```{r echo=FALSE}
df <- read.table(...)
summary(data.frame(df))
```

3.1 Rotations de xlab

# méthode 1
par(mfrow = c(1,2))
boxplot(cars, las = 3)

# méthode 2
boxplot(cars, xaxt = "n",  xlab = "")
axis(1, labels = FALSE)
text(x =  seq_along(names(cars)), y = par("usr")[3] - 10, srt = 55, adj = 1,
     labels = names(cars), xpd = TRUE)

On remarque que c’est beaucoup plus facile avec le package ggplot2.

3.2 legend() en dehors d’un plot

# on va faire un plot de 2 nuages des points de dim 2
A <- data.frame(x=rnorm(100, 20, 2), y=rnorm(100, 20, 2))
B <- data.frame(x=rnorm(100, 21, 1), y=rnorm(100, 21, 1))

# ajuster la marge
# le default est c(5,4,4,2) + 0.1
# regarder ?par qu'est-ce qu'il représente
par(mar=c(5.1, 4.1, 4.1, 7.1), xpd=TRUE)

plot(y ~ x, A, ylim=range(c(A$y, B$y)), xlim=range(c(A$x, B$x)), pch=1, col='grey')
points(y ~ x, B, pch=2, col='pink')

legend("topright", inset=c(-0.15,0), legend=c("A","B"), pch=c(1,2), title="Group", col=c('grey','pink'))

4.1 Devoir 1

Le devoir 1 est principalement un exo sur la rédaction et l’utilisation de Rmarkdown.

4.1.3 Gaussien ou pas ?

Les hypothèses de type gaussien sont souvent assurées par le TCL. Donc, tracer l’histogramme d’un échantillon de taille très limitée n’est pas forcément pertinant. Dans le problème de l’anorexie, il y a environ 30 données dans chaque groupe. Je vous donne un petit exemple intéressant.

On commence par simuler un échantillon de taille 1000, qui suit une loi normale centrée réduite.

ech <- rnorm(1000)
hist(ech)

Ensuite, on tire un échantillon de taille 30 parmi les 1000 observations qu’on a simulées.

par(mfrow = c(3,3))
replicate(9, hist(sample(ech, size = 30, replace = FALSE)))

Qu’est-ce qu’on observe ? Qu’est-ce qu’on peut conclure ?