library(MASS)
library(ggplot2)
library(gridExtra)
La estructura de la población y la edad de la menarquia se simulan de la misma forma que en el ejemplo anterior.
A diferencia de la primera simulación, en ésta si se toma en cuenta la variación materna, es decir, este modelo toma en cuenta que existen diferencias entre familias en la relación entre la menarquia y la estatura (i.e. una \(\alpha\) y \(\beta\) diferente por familia). Utilizando los mismos parámetros que anteriormente (Onland-Mored et al., 2005), el modelo del cual se simula la estatura es el siguiente:
\[estatura_{ij} = \alpha_j + \beta_j*menarquia_{ij} + e_{ij}\\ \begin{bmatrix} \alpha_j\\ \beta_j \end{bmatrix} \sim \mathcal{MVN}\left( \begin{bmatrix} \alpha_{0}\\ \beta_{0} \end{bmatrix}, S\right)\]
# Parametros
alfa.0 <- 161.5 # cm
beta.0 <- 0.31
sigma_e <- 0.1
e <- rnorm(sum(n_off),0,sigma_e) # error aleatorio
rho <- matrix(c(1, -0.3, -0.3, 1), nrow = 2)
s_alpha <- 0.2
s_beta <- 0.4
Sigmas <- c(s_alpha, s_beta)
factor_var <- 0.1 # representa la variabilidad ambiental
S <- factor_var * diag(Sigmas) %*% rho %*% diag(Sigmas)
# Simulacion de los parametros para cada madre
alfa_beta <- mvrnorm(N_mom, mu=c(alfa.0,beta.0),Sigma=S)
# Simulacion de la estatura (cm) para las hijas de cada madre
estatura <- vector(mode="numeric",length = sum(n_off))
estatura[1:n_off[1]] <- alfa_beta[1,1] +
alfa_beta[1,2]*menarquia[1:n_off[1]] +
e[1:n_off[1]] # error aleatorio
cum_n_off <- cumsum(n_off)
for(j in 2:N_mom){
estatura[(cum_n_off[j-1]+1):cum_n_off[j]] <- alfa_beta[j,1] +
alfa_beta[j,2]*menarquia[(cum_n_off[j-1]+1):cum_n_off[j]] +
e[(cum_n_off[j-1]+1):cum_n_off[j]] # error aleatorio
}
# Base de datos de las variables creadas
d <- data.frame(estatura = estatura,
menarquia = menarquia,
mother = as.factor(mother))