Ce qui précède correspond à une liste de questions posées par les étudiants pendant le TP (groupe B1 de 4M015 2017). Comme il y a beaucoup de possibilités pour faire des codes jolis, je trouve qu’il est plus raisonnable d’apprendre les détails petit à petit. Je sélectionne les problèmes qui sont (sans doute) intéressants pour qu’on puisse avancer ensemble.

N’hésitez pas à m’envoyer un email pour poser des questions.

source:

1 Remarques sur les TP

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")
Figure 1: Sélection de variables par BIC

par(mfrow = c(1,3))
plot(reg, scale = "Cp")
plot(reg, scale = "r2")
plot(reg, scale = "adjr2")
Figure 2: Sélection de variables par AIC, R2 et R2 ajusté

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()
Figure 3: Modèle proposé par BIC

2 Programmation avec R Notebook et Rmarkdown

R Notebook est un éditeur de .Rmd dans Rstudio.

2.1 Stucture de .Rmd

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

2.1.1 Exemples de chunks

2.1.1.1 Calculs, importations des données, autres manipulations sans sortie

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

2.1.1.2 Les plots

sans template

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

avec template

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

dans les commentaires, on utilise

`r fig$cap('label')`

pour la référence. J’ai donné un exemple dans le template.

2.1.1.3 Fonction avec resultat(print)

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

2.1.2 Pour aller plus loin

Vous pouvez essayer de regarder cet article écrit par l’auteur de Rmarkdown. Il est possible de cliquer directement sur le rouage pour changer les paramètres des chunks.

3 Remarques sur les plots

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)
Figure 4: Rotations de xlab

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'))