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

Tamano de muestra y variables independientes

La estructura de la población y la edad de la menarquia se simulan de la misma forma que en el ejemplo anterior.

Peso del infante

Al igual que se simulo la edad a la menarquia, se simula el peso promedio (kg; media = 9.67, DE = 1.17) al final de la lactancia (Terry et al., 2009: tabla 1) para cada madre:

\[\nu_j \sim \mathcal{N}(\nu, \varsigma_0)\] Donde j va de 1 al número de madres y \(\nu_j\) es el peso promedio al final de la lactancia para las hijas de la j-ésima madre. A continuación, se simula el peso de las hijas de cada madre al término de la lactancia.

\[infante_{ij} \sim \mathcal{N}(\nu_j, \varsigma)\]

Donde i va de 1 a \(n_j\), e \(infante_{ij}\) es el peso al término de la lactancia de la i-ésima hija de la j-ésima madre.

# Madre 
infante_mom <- rnorm(N_mom,9.67,1.17/sqrt(N_mom)) 

# Hijas
infante <- vector(mode="numeric",length = sum(n_off))
infante[1:n_off[1]] <- rnorm(n_off[1],infante_mom[1],1.17)
cum_n_off <- cumsum(n_off)
for(j in 2:N_mom){
  infante[(cum_n_off[j-1]+1):cum_n_off[j]] <- rnorm(n_off[j],infante_mom[j],1.17) 
}

Variables dependientes

Estatura

En esta simulación, además de tomar en cuenta que existen diferencias entre familias en la relación entre la menarquia y la estatura, también se simulan las diferencias entre familias en la inversión materna durante la lactancia (Figura 5.8). Utilizando los mismos parámetros que en la sección anterior (Onland-Mored et al., 2005), el modelo de la relación entre la estatura (cm), la edad a la menarquia (años) y el peso del infante (kg) es el siguiente:

\[estatura_{ij} = \alpha_j + \beta_{infante_j}*infante_{ij} + \beta_j*menarquia_{ij} + e_{ij}\\ \begin{bmatrix} \alpha_j\\\beta_{infante_j}\\ \beta_j \end{bmatrix} \sim \mathcal{MVN}\left( \begin{bmatrix} \alpha_{0}\\ \beta_{infante_0} \\ \beta_{0} \end{bmatrix}, S\right)\]

alfa.0 <- 161.5  # cm
beta.0 <- 0.31
sigma_e <- 0.1
beta_i.0 <- -0.15
e <- rnorm(sum(n_off),0,sigma_e)  # error aleatorio
rho <- matrix(c(1, 0.3, 0.0, 0.3, 1, -0.3, 0, -0.3, 1), nrow = 3) #relacion entre infante y menarquia
s_alpha <- 0.1 
s_beta <- 0.4 
s_beta_i <- 0.1
Sigmas <- c(s_alpha, s_beta, s_beta_i)
factor_var <- 0.1 # representa la variabilidad ambiental, entre mas grande mayor variacion ambiental
S <- factor_var * diag(Sigmas) %*% rho %*% diag(Sigmas) # matriz de covarianza para los par?metros

# Simulaci?n de los par?metros para cada madre
alfa_beta_beta_i <- mvrnorm(N_mom, mu=c(alfa.0,beta.0,beta_i.0),Sigma=S) 

# Simulaci?n 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_beta_i[1,1] +
  alfa_beta_beta_i[1,3] * infante[1:n_off[1]] + 
  alfa_beta_beta_i[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_beta_i[j,1] +
    alfa_beta_beta_i[j,3] * infante[(cum_n_off[j-1]+1):cum_n_off[j]] +
    alfa_beta_beta_i[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 con las variables creadas
d <- data.frame(estatura = estatura,
                menarquia = menarquia,
                infante = infante,
                mother = as.factor(mother))

Vale la pena recalcar que, en esta simulación, el tamaño del infante influye de manera directa en la edad a la menarquia y de manera indirecta en la estatura, lo que puede observarse en la tabla 5.1. El escenario que se intenta simular implica que un tamaño pequeño al término de la lactancia (trayectoria punteada en Figura 5.5) llevará a una edad temprana a la menarquia (i.e. correlación negativa entre lactancia y menarquia). Además, que una menarquia temprana llevará a una baja estatura (i.e. correlación negativa entre la menarquia y la estatura).

estatura menarquia lactancia
estatura 1.0 0.3 0.0
menarquia 0.3 1.0 -0.3
lactancia 0.0 -0.3 1.0

Visualización

library(rethinking)

mod1 <- map2stan(alist(
  estatura ~ dnorm(mu, sigma),
  mu <- a[mother] + bi[mother]*infante + bm[mother]*menarquia,
  a[mother] ~ dnorm(a0, 10),
  bi[mother] ~ dnorm(bi0, 1),
  bm[mother] ~ dnorm(bm0, 1),
  a0 ~ dnorm(160, 10),
  bi0 ~ dnorm(0, 1),
  bm0 ~ dnorm(0, 1),
  sigma ~ dgamma(1, 1)
), data = d)
## In file included from C:/Users/User/Documents/R/win-library/3.5/BH/include/boost/config.hpp:39:0,
##                  from C:/Users/User/Documents/R/win-library/3.5/BH/include/boost/math/tools/config.hpp:13,
##                  from C:/Users/User/Documents/R/win-library/3.5/StanHeaders/include/stan/math/rev/core/var.hpp:7,
##                  from C:/Users/User/Documents/R/win-library/3.5/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
##                  from C:/Users/User/Documents/R/win-library/3.5/StanHeaders/include/stan/math/rev/core.hpp:12,
##                  from C:/Users/User/Documents/R/win-library/3.5/StanHeaders/include/stan/math/rev/mat.hpp:4,
##                  from C:/Users/User/Documents/R/win-library/3.5/StanHeaders/include/stan/math.hpp:4,
##                  from C:/Users/User/Documents/R/win-library/3.5/StanHeaders/include/src/stan/model/model_header.hpp:4,
##                  from file2f9831e0731.cpp:8:
## C:/Users/User/Documents/R/win-library/3.5/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
##  #  define BOOST_NO_CXX11_RVALUE_REFERENCES
##  ^
## <command-line>:0:0: note: this is the location of the previous definition
## In file included from C:/Users/User/Documents/R/win-library/3.5/StanHeaders/include/stan/math/rev/core.hpp:44:0,
##                  from C:/Users/User/Documents/R/win-library/3.5/StanHeaders/include/stan/math/rev/mat.hpp:4,
##                  from C:/Users/User/Documents/R/win-library/3.5/StanHeaders/include/stan/math.hpp:4,
##                  from C:/Users/User/Documents/R/win-library/3.5/StanHeaders/include/src/stan/model/model_header.hpp:4,
##                  from file2f9831e0731.cpp:8:
## C:/Users/User/Documents/R/win-library/3.5/StanHeaders/include/stan/math/rev/core/set_zero_all_adjoints.hpp:14:17: warning: 'void stan::math::set_zero_all_adjoints()' defined but not used [-Wunused-function]
##      static void set_zero_all_adjoints() {
##                  ^
## 
## SAMPLING FOR MODEL 'estatura ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## 
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 2000 [  0%]  (Warmup)
## Iteration:  200 / 2000 [ 10%]  (Warmup)
## Iteration:  400 / 2000 [ 20%]  (Warmup)
## Iteration:  600 / 2000 [ 30%]  (Warmup)
## Iteration:  800 / 2000 [ 40%]  (Warmup)
## Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Iteration: 2000 / 2000 [100%]  (Sampling)
## 
##  Elapsed Time: 13.371 seconds (Warm-up)
##                17.339 seconds (Sampling)
##                30.71 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'estatura ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## 
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
## 
## 
## WARNING: No variance estimation is
##          performed for num_warmup < 20
## 
## Iteration: 1 / 1 [100%]  (Sampling)
## 
##  Elapsed Time: 0 seconds (Warm-up)
##                0 seconds (Sampling)
##                0 seconds (Total)
## 
## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]
precis(mod1,depth=2)
##          Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a[1]   158.31   2.25     154.79     161.73   568    1
## a[2]   160.61   0.71     159.43     161.71   737    1
## a[3]   161.35   0.76     160.20     162.57   825    1
## a[4]   161.35   0.42     160.73     162.01  1000    1
## a[5]   162.01   1.06     160.35     163.65   793    1
## a[6]   162.48   0.76     161.31     163.64   767    1
## a[7]   154.46   5.55     145.88     163.17   699    1
## a[8]   161.05   0.47     160.41     161.84  1000    1
## a[9]   161.62   0.77     160.24     162.65   617    1
## a[10]  162.39   0.66     161.36     163.45   810    1
## bi[1]   -0.09   0.08      -0.23       0.04   559    1
## bi[2]   -0.05   0.06      -0.14       0.04   864    1
## bi[3]   -0.13   0.08      -0.27      -0.01   881    1
## bi[4]   -0.11   0.04      -0.18      -0.05  1000    1
## bi[5]   -0.14   0.08      -0.26      -0.02   856    1
## bi[6]   -0.22   0.05      -0.30      -0.14   884    1
## bi[7]   -0.05   0.05      -0.12       0.04  1000    1
## bi[8]   -0.06   0.06      -0.16       0.02   836    1
## bi[9]   -0.09   0.06      -0.18       0.02   649    1
## bi[10]  -0.19   0.05      -0.27      -0.11   945    1
## bm[1]    0.31   0.11       0.15       0.48   589    1
## bm[2]    0.24   0.03       0.20       0.29  1000    1
## bm[3]    0.23   0.02       0.19       0.27  1000    1
## bm[4]    0.16   0.05       0.09       0.23   884    1
## bm[5]    0.36   0.03       0.31       0.41   894    1
## bm[6]    0.08   0.03       0.04       0.13  1000    1
## bm[7]    0.88   0.38       0.28       1.49   707    1
## bm[8]    0.25   0.03       0.20       0.31  1000    1
## bm[9]    0.44   0.02       0.41       0.48   875    1
## bm[10]   0.28   0.03       0.22       0.33  1000    1
## a0     160.53   2.97     155.81     165.23  1000    1
## bi0     -0.10   0.32      -0.59       0.42  1000    1
## bm0      0.31   0.31      -0.14       0.82  1000    1
## sigma    0.10   0.02       0.08       0.13   322    1
lm(estatura ~ infante + menarquia, data = d)
## 
## Call:
## lm(formula = estatura ~ infante + menarquia, data = d)
## 
## Coefficients:
## (Intercept)      infante    menarquia  
##    165.6208      -0.4599       0.1891