Simuler des données avec R

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

Pour en savoir plus sur les différentes lois de répartitions et la façon de simuler des valeurs selon ces lois, suivez ce lien

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

Pour structurer les répétitions de son tableau, il est particulièrement recommandé d'utiliser expand.grid() :

# J'ai 3 variétés, 6 répétitions et 2 types de sols...

data <- expand.grid(Variete=c("A","B","C"),  Repetition = 1:6, Sol = c("Sol1","Sol2"))

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

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.

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) 

# Avec :

# x1, la variable que l'on souhaite corréler à x2

# 0.3, le coefficient de corrélation recherché entre x1 et x2

# 10, la moyenne de x2

# 1, l'écart-type de x2

cor(x1,x2)

mean(x2)

sd(x2)

Comment simuler des variables corrélées entre elles en fixant les moyennes et écart types de ces variables.
Facile ! Suivons la méthode ci-dessous : il suffit seulement de lui donner en plus une matrice de corrélation.

Cette méthode rejoint de près ce que l'on peut faire avec genCorData() du package dédié à la simulation {simstudy}.

library(MASS)

# Nombre de variables

k <- 3

# Moyenne de chaque variable

mu <- c(65, 170, 20)

# Écart-type de chaque variable

sd <- c(7, 16, 3)

# Corrélations entre les variables

cor_matrix <- matrix(c(  1,  0.4, 0.6,  # Corrélations X1-X1, X1-X2, X1-X3  

0.4, 1,  0.3,  # Corrélations X2-X1, X2-X2, X2-X3 

0.6, 0.3, 1),    # Corrélations X3-X1, X3-X2, X3-X3),

nrow = 3, ncol = 3)# Matrice de covariance

rownames(cor_matrix) <- c("X1-Poids","X2-Taille","X3-IMC")

colnames(cor_matrix) <- rownames(cor_matrix) 

M <- cor_matrix * outer(sd, sd)# Générer des données

nobs <- 1000

data <- mvrnorm(nobs, mu = mu, Sigma = M)

cor(data)

Cette méthode est très efficace et se prête même à des régression autorégressives (Si les colonnes représentent des mesures temporelles ou spatiales, les variables proches dans le temps ou l'espace sont plus corrélées que celles éloignées),  à approfondir...

Attention, si vos corrélations ne sont pas classiques, avec des corrélations de données non normales, vous pouvez faire appel à {Copula}.

1- Simuler deux variables uniformes suivant la relation de corrélation attendue.

Établissons ici une Copule Gaussienne (Dépendance symétrique et essentiellement centrale, pas les queues). Il en existe d'autres telles que la Copule empirique empCopula(), la Copule elliptique tCopula()...

# Charger la bibliothèque nécessaire

library(copula)

# Créer une copule gaussienne avec corrélation 0.4

gauss_cop <- normalCopula(0.4, dim = 2)

# Simuler des données uniformes corrélées

data_copula <- rCopula(100, gauss_cop)

2- Transformer les données uniformes en normales, mais avec variabilité naturelle

Simulons deux variables non normales ayant les caractéristiques attendues (sauf corrélation).

# Générer des tailles de 20 hommes et femmes

tailles <- rnorm(40, mean = c(175,165), sd = c(10,8))

# Générer la longueur des poils

fqce <- rpois(40,c(70,80))

hist(tailles)

hist(fqce)

# Transformer les valeurs générées par la copule en fonction de ces tailles et poids tirés aléatoirement

data_transformed <- data.frame(

  Tailles = quantile(tailles, data_copula[, 1]),  # Appliquer la distribution de taille

  Fréquence = quantile(fqce, data_copula[, 2])    # Appliquer la distribution de poids

)

# Afficher les premières lignes pour vérifier

head(data_transformed)

# Vérifier la corrélation

cor.test(data_transformed$Taille, data_transformed$Fréquence)

# Afficher un graphique

plot(data_transformed$Taille, data_transformed$Fréquence, main = "Taille vs Fréquence cardiaque avec variabilité naturelle")

hist(data_transformed$Fréquence)

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)