Ateliers R Trucs et Astuces

# Définir les paramètres
n = 10 # Taille d'échantillon par groupe
b0 = 100 # moyenne du groupe 1 (référence)
b1 = 10 # différence de moyenne du groupe 2 vàv de la référence (contraste)
sigma = 10 # écart type

# Simuler les donnée
dfsim = data.frame(Id = 1:(n*2),
                       groupe = c(rep("1", n),
                                  rep("2", n)))
dfsim$xi = ifelse(dfsim$groupe == "1", 0, 1)
dfsim = dfsim %>% 
  mutate(y = rnorm(n(), b0 + b1 * xi, sigma))

# Stocker le modèle et explorer les coefficients
mod <- lm(y ~ groupe, dfsim)
summary(mod)

# Tester l'hypothèse (b0 + b1) = 5 x b0
hypotheses(mod, "(b1 + b2) = 5 * b1")

iter = 100
n = c(5:200)
b0 = 100
b1 = 10
sigma = 10
counter = 0 # compteur à incrémenter de 1 à chaque itération

out = matrix(NA, nrow = iter*length(n), ncol = 7)
colnames(out) = c("iter", # numéro de simulation
                  "n", # nombre de mesures par groupe
                  "b1_sim", # différence simulée
                  "b1_est", # différence estimée
                  "p", # p-value
                  "diff_est", # erreur de mesure (est - sim)
                  "perc_diff" # % erreur (est - sim) / sim * 100
)

# Barre de progression
pb <- progress_bar$new(
  format = "  downloading [:bar] :percent eta: :eta",
  total = nrow(out), clear = FALSE, width= 60)

set.seed(42)

for (i in 1:iter) {
  for (reps in 1:length(n)) {
    pb$tick()
    Sys.sleep(1 / nrow(out))
    
    counter = counter + 1
    
    # Génération du jeu de données simulées
    dfsim = data.frame(Id = 1:(n[reps]*2),
                       groupe = c(rep("1", n[reps]),
                                  rep("2", n[reps])))
    dfsim$xi = ifelse(dfsim$groupe == "1", 0, 1)
    dfsim = dfsim %>% 
      mutate(y = rnorm(n(), b0 + b1 * xi, sigma))
    
    # modèle statistique
    lm.sim = lm(y ~ groupe, dfsim)
    b1_est = as.numeric(coef(lm.sim)[2]) # valeur de b1 estimée
    p = p_value(lm.sim)$p[2] # p-value pour b1
    diff_est = b1_est - b1
    perc_diff = diff_est / b1 * 100
    
    # stockage des données à la ième ligne de la matrice 'out'
    out[counter, 1] = i # numéro de simulation
    out[counter, 2] = n[reps] # nombre de mesures par groupe
    out[counter, 3] = b1 # différence simulée
    out[counter, 4] = b1_est # différence estimée
    out[counter, 5] = p # p-value
    out[counter, 6] = diff_est # erreur de mesure (est - sim)
    out[counter, 7] = perc_diff # % erreur (est - sim) / sim * 100
  }
  
}

out %>% 
  data.frame() %>% 
  ggplot(aes(x = n, y = b1_est, fill = n)) +
  geom_point(alpha = .9, size = 2, shape = 21, color = "black") +
  geom_hline(yintercept = b1, linetype = "dashed", size = 1.5) +
  scale_fill_viridis(option = "H", direction = -1) +
  xlab("Taille d'echantillon par groupe") + 
  ylab("b1")

out = as.data.frame(out)
out = out %>% 
  mutate(S = ifelse(b1_est < 0, 1, 0),
           M = ifelse(abs(b1_est - b1_sim) / b1_sim > .2,
                      1, 0))

p1 = out %>% 
  ggplot(aes(x = n, y = b1_est, fill = factor(M))) +
  geom_point(alpha = .9, size = 2, shape = 21) +
  geom_hline(yintercept = b1, linetype = "dashed", size = 1.5) +
  scale_fill_manual(values = c("white", "red")) +
  xlab("Taille d'echantillon par groupe") + 
  ylab("b1") +
  theme(legend.position = "none") +
  labs(title = "Erreur de type M",
       subtitle = "Estimés > 20 % d'erreur")

p2 = out %>% 
  ggplot(aes(x = n, y = b1_est, fill = factor(S))) +
  geom_point(alpha = .9, size = 2, shape = 21) +
  geom_hline(yintercept = b1, linetype = "dashed", size = 1.5) +
  scale_fill_manual(values = c("white", "red")) +
  xlab("Taille d'echantillon par groupe") + 
  ylab("b1") +
  theme(legend.position = "none") +
  labs(title = "Erreur de type S")

p1 + p2

set.seed(42)
penguins %>% 
  drop_na() %>% 
  group_by(sex, species) %>% 
  slice_sample(n = 5) %>% 
  ggplot(aes(y = flipper_length_mm, x = sex, group = species)) +
  geom_point() +
  facet_wrap(~species) +
  ylab("Flipper length (mm)")

set.seed(42)
df <- penguins %>% 
  drop_na() %>% 
  group_by(sex, species) %>% 
  slice_sample(n = 5) %>% 
  summarise(mean_fl = mean(flipper_length_mm, na.rm = T),
            sd_fl = sd(flipper_length_mm, na.rm = T)) %>% 
  mutate(N = 5) %>%
  pivot_wider(names_from = sex, values_from = c(mean_fl, sd_fl)) %>% 
  mutate(Delta_mean = mean_fl_male - mean_fl_female) %>% 
  select(species, N, 
         mean_fl_female, sd_fl_female,
         mean_fl_male, sd_fl_male,
         Delta_mean) %>% 
  rename(Species = species,
         # Females = mean_fl_female,
         # Males = mean_fl_male,
         Delta = Delta_mean)
df %>% 
  tt(digits = 3) %>% 
  setNames(c("", "", "mean", 
             "sd", "mean", "sd", "")) %>% 
  group_tt(j = list(
    "Species" = 1,
    "N" = 2,
    "Females" = 3:4,
    "Males" = 5:6,
    "Delta" = 7))

# Generer le dataframe
df = crossing(n = 1:n, Sex, Sp) %>% 
  mutate(Sex_bool = ifelse(Sex == "Males", 1, 0),
         Sp_bool_1 = ifelse(Sp == "Ch", 1, 0),
         Sp_bool_2 = ifelse(Sp == "Gen", 1, 0))

# Générer les données
set.seed(42)
df = df %>% 
  mutate(flipper_length = rnorm(n(), 
                      b0 + 
                        b1 * Sex_bool +
                        b2 * Sp_bool_1 + 
                        b3 * Sp_bool_2 +
                        b4 * Sex_bool * Sp_bool_1 +
                        b5 * Sex_bool * Sp_bool_2,
                      sigma))

# Graphique
df %>% 
  ggplot(aes(x = Sex, y = flipper_length, color = Sp, group = Sp)) +
  geom_point(size = 4, alpha = .6) +
  geom_smooth(method = "lm") +
  facet_wrap(~Sp) +
  scale_color_manual(values =
                       c("darkorange","darkorchid","cyan4")) +
  theme(legend.position = "none")

df_parms = data.frame(Parameter = 
                        c("(Intercept)", "SexMales",
                          "SpCh", "SpGen", "SexMales:SpCh",
                          "SexMales:SpGen"),
                      Coefficient = c(b0, b1, b2, b3, b4, b5))

parms <- model_parameters(lm.penguins) %>% data.frame()
parms %>% 
  ggplot(aes(x = Coefficient, y = Parameter)) +
  geom_pointrange(aes(xmin = CI_low, xmax = CI_high)) +
  geom_point(data = df_parms,
             aes(x = Coefficient, 
                 y = Parameter, 
                 fill = "red"), 
             shape = 23, 
             size = 2.5, 
             alpha = .6) +
  theme(legend.position = "none")

p <- plot_comparisons(lm.penguins, variables = "Sex", condition = "Sp")

df_dim = data.frame(marginaleffects_term_index = rep("Sex", 3),
                    Sp = c("Ad", "Ch", "Gen"),
                    Estimate = c(mu_A_M - mu_A_F,
                                 mu_C_M - mu_C_F,
                                 mu_G_M - mu_G_F))
p +
  geom_point(data = df_dim,
             aes(x = Sp, 
                 y = Estimate, 
                 fill = "red"), 
             shape = 23, 
             size = 3, 
             alpha = .9) +
  theme(legend.position = "none")

out %>% 
  mutate(sig = case_when(p < 0.05 ~ 1, p >= 0.05 ~ 0)) %>%
  group_by(n) %>% 
  summarise(sig.prop = sum(sig)/n() * 100) %>% 
  ggplot(aes(x = n, y = sig.prop, fill = n)) +
  geom_line() +
  geom_point(alpha = .9, size = 3, shape = 21, color = "black") +
  geom_hline(aes(yintercept=80),linetype="dashed") +
  scale_fill_viridis(option = "H", direction = -1) +
  xlab("Taille d'echantillon par groupe") + 
  ylab("Puissance statistique")

n = 20
Sp = c("Ad", "Ch", "Gen")
b0 = plogis(-.4) # moyenne de survie manchots d'Adélie  = 40 %
b1 = plogis(1) # différence de survie Chinstrap - Adélie = 64 %
b2 = plogis(1.5) # différence de survie Gentoo - Adélie


# Generer le dataframe
df <- crossing(n = 1:n, Sp) %>% 
  mutate(Sp_bool_1 = ifelse(Sp == "Ch", 1, 0),
         Sp_bool_2 = ifelse(Sp == "Gen", 1, 0))

# Générer les données
set.seed(42)
df <- df %>% 
  mutate(logit_p = b0 + 
           b1 * Sp_bool_1 +
           b2 * Sp_bool_2) %>% # Echelle logit
  mutate(p = plogis(logit_p)) %>% # Echelle des probabilités
  mutate(Survie = rbinom(n(), prob = p, size = 1))

glm.binom <- glm(Survie ~ Sp, df, family = "binomial")
model_parameters(glm.binom) # Coefficients

avg_comparisons(glm.binom) # Différence de survie (%)

plot_predictions(glm.binom, by = "Sp")

n = 20
Sp = c("Ad", "Ch", "Gen")
b0 = log(10) # moyenne de fécondité manchots d'Adélie  = 10
b1 = log(2) # différence de fécondité Chinstrap - Adélie = +2
b2 = plogis(4) # différence de fécondité Gentoo - Adélie + 4


# Generer le dataframe
df <- crossing(n = 1:n, Sp) %>% 
  mutate(Sp_bool_1 = ifelse(Sp == "Ch", 1, 0),
         Sp_bool_2 = ifelse(Sp == "Gen", 1, 0))

# Générer les données
set.seed(42)
df <- df %>% 
  mutate(log_lambda = b0 + 
           b1 * Sp_bool_1 +
           b2 * Sp_bool_2) %>% # Echelle log
  mutate(lambda = exp(log_lambda)) %>% # Echelle des comptages
  mutate(Fecondite = rpois(n(), lambda = lambda))

glm.poiss <- glm(Fecondite ~ Sp, df, family = "poisson")
model_parameters(glm.poiss) # Coefficients

avg_comparisons(glm.poiss) # Différence de survie (%)

plot_predictions(glm.poiss, by = "Sp") # Figure

n = 10
sigma = 10
b0 = 100
b1 = 50
Eq = 5 # Région d'équivalence = +/- 5

groupe = c("1", "2")

df = data.frame(Rep = 1:(n[1]*2),
                groupe = c(rep("1", n[1]),
                           rep("2", n[1])))
df$xi = ifelse(df$groupe == "1", 0, 1)

set.seed(42)
df = df %>% 
  mutate(y = rnorm(n(), b0 + b1 * xi, sigma))

lm.eq = lm(y ~ groupe, df)
# pred.b1 = avg_slopes(lm.eq, variables = "groupe")
hyp = hypotheses(lm.eq, 
                 equivalence = 
                   c(- Eq, Eq))
hyp

hyp %>% 
    data.frame() %>% 
    ggplot(aes(y = term, x = estimate)) +
    geom_point() +
    geom_pointrange(aes(xmin = conf.low, xmax = conf.high)) +
    geom_rect(aes(xmin = -Eq, xmax = Eq, ymin = 0, ymax = 2), , fill = "grey", alpha = .6)