#- Simpsons paradox: https://en.wikipedia.org/wiki/Simpson%27s_paradox
#- una tendencia que aparece en varios grupos de datos desaparece cuando estos grupos se combinan y en su lugar aparece la tendencia contraria para los datos agregados.

#- Ejemplos con R:
#- https://paulvanderlaken.com/2017/09/27/simpsons-paradox-two-hr-examples-with-r-code/
#- https://solomonkurz.netlify.app/post/2019-10-09-individuals-are-not-small-groups-i-simpson-s-paradox/
#- https://solomonkurz.netlify.app/post/2019-10-14-individuals-are-not-small-groups-ii-the-ecological-fallacy/
#- http://svmiller.com/blog/2020/01/illustrate-correlation-fallacies-limitations-in-r/
#- https://rpubs.com/shampjeff/blog_post_2  (con datos de pingüinos)
#- https://www.pipinghotdata.com/posts/2021-10-11-estimating-correlations-adjusted-for-group-membership/


#- 1er ejemplo: con datos de pingüinos  ----------------------------------------
#- https://www.pipinghotdata.com/posts/2021-10-11-estimating-correlations-adjusted-for-group-membership/

library(tidyverse)
#- cojo datos de pingüinos y arreglo un poco 
df <- palmerpenguins::penguins %>% 
  dplyr::select(species, bill_depth_mm, body_mass_g) %>% 
  tidyr::drop_na()

p <- ggplot(data = df, 
       mapping = aes(x = bill_depth_mm, y = body_mass_g, color = species)) +
     geom_point()
p 

p +  facet_wrap(. ~ species, ncol = 3) 


p <- p + geom_smooth(method = "lm", se = FALSE) 

p

p + geom_smooth(method = "lm", se = FALSE, color = "black") 



#- creo df_ALL donde no se diferencia x especies de pingüino
df_ALL <- df %>% 
  dplyr::mutate(species = "ALL")

#- junto df y df_ALL
df_ok <- bind_rows(df, df_ALL)

#- para q se vean usamos geom_jitter()
p <- ggplot(data = df_ok, 
            mapping = aes(x = bill_depth_mm, y = body_mass_g, color = species)) +
     geom_jitter()
p 

p <- ggplot(data = df_ok, 
            mapping = aes(x = bill_depth_mm, y = body_mass_g, color = species)) +
     geom_point() + 
     geom_smooth(method = "lm", se = FALSE) 
p

p +  facet_wrap(vars(species), ncol = 4)     
p +  facet_wrap(vars(species), nrow = 2, ncol = 2)     


#- venga vamos a arreglar el orden de species para q "ALL" vaya al final
str(df_ok)  #- species es character
df_ok <- df_ok %>% mutate(species = forcats::as_factor(species))
str(df_ok)  #- species ahora es un factor
levels(df_ok$species)  #- se nos ha ordenado de chiripa (xq ALL estaba en mayusculas)
#- tendríamos q haber hecho esto
df_ok <- df_ok %>% mutate(species = forcats::fct_relevel(species, "ALL", after = 3))


#- el gráfico final con "ALL" al final -----------------------------------------
p <- ggplot(data = df_ok, 
            mapping = aes(x = bill_depth_mm, y = body_mass_g, color = species)) +
     geom_point() + 
     geom_smooth(method = "lm", se = FALSE) +
     facet_grid(cols = vars(species)) 
p

p +  facet_wrap(vars(species), nrow = 2, ncol = 2)     





#- 2º ejemplo: datos de school admissions to University of California, Berkeley ------------------------------------
#- code for: https://solomonkurz.netlify.app/post/2019-10-09-individuals-are-not-small-groups-i-simpson-s-paradox/
library(tidyverse)

UCBAdmissions <- as.data.frame(UCBAdmissions)
df <- UCBAdmissions %>% 
  as_tibble() %>% 
  pivot_wider(id_cols = c(Dept, Gender),names_from = Admit, values_from = Freq) %>% 
  mutate(total = Admitted + Rejected)

df_table <- df %>% 
  group_by(Gender) %>% 
  summarise(percent_admitted = (100 * sum(Admitted) / sum(total)) %>% round(digits = 1))

df_a <- df %>%  
  mutate(dept = str_c("department ", Dept)) %>% 
  pivot_longer(cols = Admitted:Rejected,
               names_to = "admit",
               values_to = "n") 
  
df_a %>% 
  ggplot(aes(x = Gender, y = n, fill = admit)) +
  geom_col(position = "dodge") +
  scale_fill_viridis_d(option = "B", end = .8) +
  xlab(NULL) +
  theme(panel.grid = element_blank()) +
  facet_wrap(~dept)

df %>%
  mutate(p = Admitted / total) %>% 
  ggplot(aes(x = Dept, y = p)) +
  geom_hline(yintercept = .5, color = "white") +
  geom_point(aes(color = Gender, size = total),
             position = position_dodge(width = 0.3)) +
  scale_color_manual(NULL, values = c("red3", "blue3")) +
  scale_y_continuous("admission probability", limits = 0:1) +
  xlab("department") +
  theme(panel.grid = element_blank())

