Simuler des données avec R
Simulations de données sous R afin de contrôler la pertinence de son protocole AVANT expérience.
L'essentiel de cette page
Il peut être nécessaire de simuler des valeurs sous R pour tester son raisonnement : c'est même absolument nécessaire en science où l'on devrait toujours prépublier son raisonnement à l'avance (éviter ainsi le biais de sélection ou le biais de jugement a posteriori) et tester sa solidité à l'avance de son expérience sur des données simulées avant de mettre la main au porte-monnaie.
Voici ci-dessous quelques fonctions basiques pour simuler ses données.
On notera en partie 2 la fonction "maison" range_to_sd() qui va permettre de retrouver l'écart-type théorique d'une population lorsqu'on a seulement un intervalle théorique. Très pratique pour inventer des valeurs collant à la théorie.
Point important : il existe des méthodes pour contrôler que les échantillons ont une bonne taille et vérifier l'effet détectable avec un échantillon.
1- Simuler des valeurs avec des fonctions classiques
Créer une série de valeurs
# Créer une liste de valeurs de 1 à 150
x <- 1:150
Créer une liste de 10000 valeurs suivant la loi normale, moyenne 20, écart-type 3.
x <- rnorm(10000,20,3)
Afficher l'ensemble des objets présents dans la console "R".Créer une liste de 10000 valeurs suivant la loi de poisson et autour d'un événement de probabilité maximale de 2.
x <- rpois(10000,2)
Créer une liste de 1000 valeurs suivant une loi uniforme, minimum 10, maximum 90.
x <- runif(1000,10,90)
Simuler un échantillon à partir des données d'un vecteur
# Prenons un vecteur x
x <- runif(1000,10,90)
# Voici un échantillon de x de 20 valeurs
echantillon <- sample(x,20) ; echantillon
Mettre les valeurs d'une liste dans un ordre aléatoire
x <- sample(x)
# Tirer au sort avec sample avec un biais de sélection
# Simuler un échantillon de 20 individus issus d'une population à 30% de filles
sample(c("Fille","Garçon"),20,replace=T,prob=c(0.3,0.7))
# Répéter 1 et 2, 3 fois de suite chacun dans 2 ordres différents
rep(1:2,each=3) # 1 1 1 2 2 2
rep(1:2, times = 3) # 1 2 1 2 1 2
Exemple : simulons une variable poids qui varient en fonction d'une variable catégorielle.
Approche 1 - rudimentaire
cat <- rep(c("A","B","T"),each=20) ; cat
poids <- c(rnorm(20,50,10),rnorm(20,55,10),rnorm(20,52,10))
Approche 2 - automatisable grâce à Map() qui applique une formation à plusieurs vecteurs en parallèles
# Définir les paramètres
n <- 20 # Individus par groupe ; moyennes <- c(50, 52, 55) ; ecart_types <- c(10, 12, 8)
cat <- rep(c("A", "B", "T"), each=n)
# Générer les poids en utilisant Map pour appliquer la fonction avec des arguments multiples
poids <- unlist(Map(function(m, sd) rnorm(n, m, sd), moyennes, ecart_types))
boxplot(poids~cat)
2- Simuler des valeurs normales à partir d'un intervalle
Je sais que la production de maïs dans certaines conditions varie entre 33 et 50 quintaux à l'hectare. Bien, simulons alors les rendements de différentes parcelles ? Mais avec quelle moyenne ? Quel écart-type ? La fonction range_to_sd() va nous aider.
# Calculer la moyenne
(33+50)/2 # La moyenne est de 41.5
# La fonction range_to_sd() - code à copier-coller sans le modifier
range_to_sd <- function(min, max, conf.level = 0.99) {
z <- qnorm((1 + conf.level) / 2)
return((max - min) / (2 * z))
}
# Calculer l'écart-type théorique
range_to_sd(33,50) # écart-type de 3.3
# Il ne reste qu'à simuler une population
population <- rnorm(10000,41.5,3.3)
quantile(population,p=c(.005,0.995)) # On voit bien que l'on retrouve les rendements théoriques dans 99% des cas
3- Simuler des données corrélées à une variable
Je souhaite simuler une variable corrélée à une ou plusieurs variables déjà existantes en fixant la moyenne et l'écart-type ?
Je peux utiliser le fonction cor.e() dérivée de la fonction du même nom proposée par Maxime Lopez.
Cliquer sur ce lien pour dérouler le code de la fonction cor.e() à copier-coller
cor.e <- function(dataframe, cor.f, meani=0, sdi=1){
# dataframe le dataframe des variables initiales
# cor.f la liste des corrélations souhaitées entre la variables initiales et la variable crée z
# meani et sdi la moyenne et l'ecart type de la variable créée
mat.ini <- as.matrix(dataframe)
cor1 <- cor(mat.ini)
z <- rnorm(nrow(mat.ini), mean=meani, sd=sdi)
z <- z*sqrt(sdi^2/var(z))
z <- z+(meani-mean(z))
m1 <- cbind(mat.ini, z)
m2 <- scale(m1)
c1 <- chol(var(m2))
m3 <- m2 %*% solve(c1)
cor2 <- cbind(rbind(cor1, z=cor.f), z=c(cor.f,1))
m4 <- m3 %*% chol(cor2)
m5 <- sweep( m4, 2, attr(m2, 'scaled:scale'), '*')
m5 <- sweep( m5, 2, attr(m2, 'scaled:center'), '+')
return(m5[,ncol(m5)])
}
Exemple : j'ai une variable x1 d'une moyenne de 20 et écart-type de 3.
Je souhaite simuler x2 d'une moyenne de 10 et d'écart-type 1 corrélée à x1 à 0.3.
Allons-y :
x1 <- rnorm(1000,20,3)
x2 <- cor.e(x1,0.3,10,1)
cor(x1,x2)
mean(x2)
sd(x2)
Comment simuler des variables corrélées entre elles en fixant les moyennes et écart-type de ces variables.
Facile ! Suivons la méthode ci-dessous :
library(MASS)
# k : nombre de variables souhaité
k <- 2
# Moyenne de chaque variable - ici 40 et 20
mu <- c(40,20)
# Ecart-type de chaque variable - ici 4 et 2
sd <- c(4, 2)
# Corrélation entre les variables - ici 0.4
cor <- 0.4
# Simuler une matrice de covariance
M <- matrix(cor, nrow=k, ncol=k) * outer(sd, sd)
diag(M) <- sd^2
# Générer des données
nobs <- 100000
data <- mvrnorm(nobs, mu = mu, Sigma = M)
# Vérifier les moyennes et les variances
colMeans(data)
apply(data, 2, sd)
cor(data)
x <- data[,1] # Ma première variable, moyenne 40, écart-type 4
y <- data[,2] # Ma deuxième variable, moyenne 20, écart-type 2
4- Répétitions de tests statistiques, extraire la p-value de ses données simulées afin de contrôler que notre expérience sera discriminante dans 95% des cas et d'adapter ses tailles d'échantillons
Remarque : pour extraire la p-value d'un test statistique, on met $p.value à la fin.
Pour m.test() de {KefiR}, il faut ajouter les arguments return, verbose et plot :
pval = m.test(x,y,return=FALSE,verbose=FALSE, plot=FALSE, boot=FALSE)
# Commande pour pouvoir exécuter le test G.
# library(DescTools)
# Taille d'échantillon envisagée
n <- 20
# Liste des p-value à l'issu de 1000 simulations
pval <- c()
for (i in 1:1000) {
# Simuler mes variables.
# ici j'en simule 2 qualitatives
# mais je peux aussi simuler des quantitatives avec rnorm()
# cat 1, catégorielle A, B ou C
conditions <- c( rep("A",n),rep("B",n),rep("C",n))
# cat 2, catégorielle germé ou pas germé avec prop différentes selon condition
germination <- c(
sample(c("Germé","Pas germé"),n,replace=T,prob=c(.3,.7)),
sample(c("Germé","Pas germé"),n,replace=T,prob=c(.4,.6)),
sample(c("Germé","Pas germé"),n,replace=T,prob=c(.3,.7)))
# Test stat : n'importe lequel marche selon les variable
# test G, Khi², anova, ancova, student, corrélation et comparaison de moyennes
# Remarque, si vous utilsez m.test() de{KefiR},il renvoie directement une p-value en faisant le code suivant (à adapter)
# m.test(quant,cal,return=FALSE,verbose=FALSE)
# Récupérer la pvalue de la simulation
pval_temp <- GTest(table(conditions,germination))$p.value
# L'ajouter à la listes des p-values
pval <- c(pval,pval_temp)
}
#hist(pval)
# Quantile 95% - cette valeur permet de vérifier
# que 95% des p-value sont inférieures à 0.05
# cette valeur doit donc être < 0.05
quantile(pval,prob=0.95)