library(MASS)
library(ggplot2)
library(gridExtra)

Tamano de muestra

Se simula una población de n=51 mujeres. Supongamos que estas mujeres pertenecen a diferentes familias que estan agrupadas por madre. Digamos, por ejemplo, que hay 10 madres y cada madre tiene un numero diferente de hijas. La anterior agrupación tiene como objetivo simular el escenario biológico en el que el fenotipo individual varia menos entre hijas de la misma madre que entre hijas de diferentes madres.

set.seed(1986)
N_mom <- 10 #Numero de madres
n_off <- c(6,4,5,7,6,3,7,5,4,4) #Numero de hijas
mother <- rep(1:N_mom, times = n_off) #Indice de la madre

Variables independientes

Menarquia

Siguiendo al artículo de Terry y colaboradores (2009), se simula la edad a la menarquia promedio (años; media = 12.5, DE = 1.72) para cada madre a partir del siguiente modelo:

\[\mu_j \sim \mathcal{N}(\mu, \sigma_0)\] Donde j va de 1 al número de madres y \(\mu_j\) es la edad promedio a la menarquia de las hijas de la j-ésima madre. A continuación, se simula la edad a la menarquia de las hijas de cada madre.

\[menarquia_{ij} \sim \mathcal{N}(\mu_j, \sigma)\]

Donde i va de 1 a \(n_j\), y \(menarquia_{ij}\) es la edad a la menarquia de la i-ésima hija de la j-ésima madre.

#madre
menar_mom <- rnorm(N_mom, 12.5, 1.72/sqrt(N_mom))

#hijas
menarquia <- vector(mode="numeric",length = sum(n_off))
menarquia[1:n_off[1]] <- rnorm(n_off[1],menar_mom[1], 1.72)
cum_n_off <- cumsum(n_off)
for(j in 2:N_mom){
  menarquia[(cum_n_off[j-1]+1):cum_n_off[j]] <- rnorm(n_off[j],menar_mom[j],1.72) 
}

Variables dependientes

Estatura

Se simula la asociación que nos interesa, esto es, la relación entre la estatura (cm) y la edad a la menarquia sin hacer explícita la variación en la inversión hecha por diferentes madres (Figura 5.8). Siguiendo el artículo de Onland-Mored y colaboradores (2005) se escogieron los parámetros que se utilizan en el siguiente modelo:

\[ estatura_{ij} = \alpha + \beta*menarquia_{ij} + e_{ij} \]

Donde \(estatura_{ij}\) es la estatura de la i-ésima hija de la j-ésima madre y \(e_{ij}\) es el error aleatorio.

# Parametros
alfa <- 161.5  # cm
beta <- 0.31
sigma_e <- 0.2
e <- rnorm(sum(n_off),0,sigma_e)  # error aleatorio

# La estatura en funcion de la edad a la menarquia
estatura <- alfa + beta*menarquia + e

# Base de datos de las variables creadas
d <- data.frame(estatura = estatura,
                menarquia = menarquia,
                mother = as.factor(mother))

Visualización