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)

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)

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

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)