twins.R

library(tidyverse)

# simulation --------------------------------------------------------------
add_noise <- function(x, sd, scale = TRUE){
  out <- x + rnorm(length(x), sd = sd)
  if(scale)out <- scale(out)
  out
}
n <- 100
twins_m <-
  tibble(
    h = rnorm(n),
    type = "m",
    family_id = str_c("m", seq_len(n)),
    x1 = add_noise(h, 1),
    x2 = add_noise(h, 1),
    z1 = add_noise(h, 1),
    z2 = add_noise(h, 1)
  )

twins_d <- tibble(
  h = rnorm(n),
  type = "d",
  family_id = str_c("d", seq_len(n)),
  x1 = add_noise(h, 2),
  x2 = add_noise(h, 2),
  z1 = add_noise(h, 2),
  z2 = add_noise(h, 2)
)

twins <- bind_rows(twins_m, twins_d)
rm(twins_d, twins_m, n)

# reshape1 ----------------------------------------------------------------
# the simulated data needs to be shaped as your data
# just ignore it

twins <- pivot_longer(twins, matches("[xz][12]")) %>% 
  mutate(var = str_extract(name, "^[xz]"),
         individual_id = str_c(str_extract(name, "[12]$"), family_id)) %>% 
  select(-name) %>% 
  pivot_wider(names_from = var, values_from = value)


# differences -------------------------------------------------------------
# here I calculate the differences between twins
twin_diff <- function(x){
  if(length(x) != 2L)stop("The number of twins is not two!")
  abs(x[[1]] - x[[2]])
}
twins_diff <- twins %>%
  group_by(family_id) %>% 
  mutate(x_diff = twin_diff(x), z_diff = twin_diff(z), .groups = "drop") %>% 
  ungroup()


# line-plot ---------------------------------------------------------------
twins <- twins %>% 
  mutate(twin_nr = as.character(rep(1:2, NROW(twins)/2)))

ggplot(twins, aes(twin_nr, x, group = family_id)) +
  geom_line(alpha = .3) +
  facet_wrap(~type) +
  theme_minimal() +
  scale_x_discrete(expand = c(0.05, 0.05))

# box-plot ----------------------------------------------------------------

# I don't really like boxplots
ggplot(twins_diff, aes(type, x_diff)) +
  geom_boxplot() +
  theme_minimal()

# Violine Plots say more about the distribution
ggplot(twins_diff, aes(type, x_diff)) +
  geom_violin() +
  geom_boxplot(width = .1, alpha = 0, fill = NA) +
  theme_minimal()