Simuler des données avec {simstudy} sous R

La démarche générale de {simstudy} est la suivante :

1. Bien décrire les caractéristiques de chaque variable (loi de répartition, moyenne, dépendance vis à vis d'une autre variable...) avec defData().

2. Générer un jeu de données en lui imposant une taille avec genData().

Son argument dist permet de définir la loi de répartition.

1. Distributions continues :

• "normal" : Distribution normale (paramètres : moyenne, variance (écart type au carré)).

• "gamma" : Distribution Gamma (paramètres : forme, échelle).

• "beta" : Distribution Beta (paramètres : alpha, beta).

• "uniform" : Distribution uniforme continue (paramètres : bornes inférieure et supérieure).

• "exponential" : Distribution exponentielle (paramètre : taux).

• "lognormal" : Distribution log-normale (paramètres : moyenne logarithmique, écart-type logarithmique).

• "uniformInt" : Distribution uniforme discrète (paramètres : bornes inférieure et supérieure pour des entiers).

2. Distributions discrètes :

• "binary" : Distribution binaire (Bernoulli) avec une probabilité de succès (paramètre : probabilité de 1).

• "poisson" : Distribution de Poisson (paramètre : lambda).

• "negBinomial" : Distribution binomiale négative (paramètres : nombre d'échecs, probabilité de succès).

• "categorical" : Distribution catégorielle (paramètre : une chaîne de probabilités séparées par des points-virgules, par exemple, "0.3;0.5;0.2" pour une variable à 3 catégories).

• "uniformInt" : Génère un nombre entier tiré d'une distribution uniforme discrète (paramètres : minimum, maximum).

3. Distributions spécifiques :

• "nonrandom" : Génère une variable déterministe selon une formule donnée (pas de variance aléatoire).

• "mixture" : Crée une variable selon un mélange de plusieurs distributions. Il faut définir les proportions de chaque composante.

Prenons un exemple. Je veux simuler un jeu de données :

• Id, identifiant unique de chaque ligne.

• Variétés de pommes, "Gala", "Granny" ou "Fuji"

• Couleur, une variable qui dit la couleur dominante de la pomme (toutes rouges, sauf les Granny, vertes)

• Poids, poids de pommes fluctuant autour de 175 g avec SD de 10, mais renseigne toi car le poids doit varier selon la variété).

• Hauteur, ... de trognon qui soit corrélée au poids (valeurs crédibles).

• Diamètre qui corresponde grosso modo à une formule calculée en fonction du poids et de la hauteur.

Créons les variables Variété et Id :

def <- defData(varname = "Variété", dist = "categorical", formula = "1/3;1/3;1/3", id = "Id")

On peut aussi ajouter la variable couleur à def

def <- defData(def, varname = "Couleur", dist = "nonrandom", formula = "ifelse(Variété == 2, 'vert','rouge') ")

Je l'ai ajouté à la variété numéro 2 (Granny) car les catégories de pommes n'ont pas encore de nom.

A ce stade, on pourrait ajouter d'autres variables, mais l'on va générer le jeu de données (150 valeurs) afin de nommer les variétés avant de rajouter des variables.

Générer les données

dt <- genData(150, def) 

Remplaçons les catégories de variétés par leurs noms réels.

dt[, Variété := factor(Variété , levels = c(1,2,3), labels = c("Gala", "Granny", "Fuji"))]

Ajoutons la variable normale poids telle qu'elle change en fonction de la variété.

defW <- defCondition(condition = "Variété == 'Gala'",

formula = 160, variance = 64, dist = "normal") 

defW <- defCondition(defW, condition = "Variété == 'Granny'",

formula = 180, variance = 100, dist = "normal") 

defW <- defCondition(defW, condition = "Variété == 'Fuji'",

formula = 200, variance = 144, dist = "normal") 

dt <- addCondition(defW, dt, newvar = "Poids") 

Définir la hauteur du trognon corrélée au poids

def <- defDataAdd(varname = "Hauteur", formula = "0.2 * Poids", variance = 4, dist = "normal") 

# Ajouter Hauteur aux données existantes

dt <- addColumns(def, dt) 

# Ajouter la définition de Diamètre sans écraser Hauteur 

#Nouveau def que l'on n'ajoute pas à def

def <- defDataAdd(varname = "Diamètre", formula = "0.5 * sqrt(Poids) + 0.1 * Hauteur", variance = 4, dist = "normal")

# Ajouter la variable Diamètre aux données existantes

dt <- addColumns(def, dt) 

# Ajouter une variable de traitement avec trtAssign

dt<- trtAssign(dt, nTrt = 2, grpName = "traitement", balanced = TRUE) 

Si l'on souhaite générer des données corrélées, {simstudy} s'y prête parfaitement.

On peut lui préciser la corrélation entre deux variables (rho) ou lui glisser une matrice de corrélation.

Générer des variables corrélées en précisant un coefficient de corrélation

# Exemple Grain et Fourrage ont une corrélation de 0.3

dt <- genCorData(200, # Nb de lignes

      mu = c(41.5,23), # Moyennes

      sigma = c(3.3,1.94), # Ecart-types

      rho = 0.3, # r de Pearson

      corstr = "cs", # compound symmetry" (toutes les paires ont la même corrélation), sinon "ind" pour indépendante (non corrélées)

cnames = c("Grain","Fourrage")) # Noms de variables

cor.test(dt$Grain, dt$Fourrage)

# Exemple 2 - 3 variables corrélées à 0.4

dt <- genCorData(1000, mu = c(4, 12, 3), sigma = 3, rho = 0.4, corstr = "cs", cnames = c("x0",

    "x1", "x2"))

dt <- genCorData(1000,

mu = c(4, 12, 3),

sigma = c(1, 2, 3),

corMatrix = C)

dt

Générer deux variables corrélées en précisant le coefficient de corrélation

En précisant une matrice de corrélation, on peut aisément générer un jeu de données corrélées. 

Cette fonctionnalité rejoint ce que l'on peut faire plus classiquement avec mvrnorm().

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) 

dt <- genCorData(1000,

mu = c(65, 170, 20),

sigma = c(7, 16, 3),

corMatrix = cor_matrix,

      cnames = c("Poids","Taille","IMC à 20 ans"))

dt

# Charger simstudy

library(simstudy)

# Paramètres de la simulation (intégrés directement dans les formules)

nb_periodes <- 100  # Nombre de périodes

# Définir les données de base avec des noms en français

def <- defData(varname = "Temps", dist = "nonrandom", formula = "1:100")

# "tendance" va monter en fonction du temps

def <- defData(def, varname = "tendance", dist = "nonrandom", 

                      formula = "0.5 * Temps")

# "saisonnalite" va monter et descendre avec une périodicité de 12

def <- defData(def, varname = "saisonnalite", dist = "nonrandom", 

                      formula = "10 * sin(2 * 3.14159 * Temps / 12)")  # Remplacement de pi par 3.14159

# "bruit" correspond à des valeurs variants autour de 0

def <- defData(def, varname = "bruit", dist = "normal", 

                      formula = "0", variance = 5^2)

# Valeur est la combinaison de tendance+saisonnalité+bruit

def <- defData(def, varname = "Valeur", dist = "nonrandom", 

                      formula = "tendance + saisonnalite + bruit")

# Générer la série temporelle

donnees <- genData(nb_periodes, def)

# Graphique simple avec plot

plot(donnees$Temps, donnees$Valeur, type = "l", col = "blue", lwd = 2,

     main = "Série Temporelle Simulée avec Tendance, Saisonnalité et Bruit",

     xlab = "Temps", ylab = "Valeur")

Convertir le temps en données calendaires en intégrant la date (mois, année...)

# Ajouter une colonne de dates, en supposant que le 1er jour correspond au 1er janvier 2023

date_debut <- as.Date("2023-01-01")

donnees$Temps <- date_debut + (donnees$Temps - 1) * 30 # (L'expérience par de 2023 point couvre 30 j)

# Afficher les premières lignes pour vérifier

head(donnees)

# Graphique simple avec plot

plot(donnees$Temps, donnees$Valeur, type = "l", col = "red", lwd = 2,

     main = "Série Temporelle Simulée avec Tendance, Saisonnalité et Bruit",

     xlab = "Temps", ylab = "Valeur")