3 Multi-response Models
We can simulate multiple response variables, that covary at different hierarchical levels. For example, we might have two phenotypes, body mass and behaviour, that covary at the between individual and within individual (residual) levels. In the case of such a simple random effects model, we have a covariance matrix at each level:
\[ \boldsymbol{y}_{ij} = \boldsymbol{\beta}_0 + \boldsymbol{u}_j + \boldsymbol{\epsilon}_{ij} \] \[ \boldsymbol{u}_i \sim \mathcal{N}(0, \Sigma_u) \] \[ \boldsymbol{\epsilon}_{i} \sim \mathcal{N}(0, \Sigma_{\epsilon}) \] \[ \Sigma_u = \begin{bmatrix} \sigma^2_{u_1} & \sigma_{u_1u_2} \\ \sigma_{u_1u_2} & \sigma^2_{u_2} \end{bmatrix} , \Sigma_{\epsilon} = \begin{bmatrix} \sigma^2_{\epsilon_1} & \sigma_{\epsilon_1\epsilon_2} \\ \sigma_{\epsilon_1\epsilon_2} & \sigma^2_{\epsilon_2} \end{bmatrix} \]
We can indicate that there are multiple phenotypes within the parameter list using the n_response argument. If we have \(q\) response variables, vcov now needs to either be a vector of length \(q\) (the variances, if we assume the covariances are 0) or a \(q*q\) covariance matrix. So below, we simulate 2 response variables, with a covariance between them at both individual and residual levels.
squid_data <- simulate_population(
data_structure=make_structure(structure = "individual(100)",repeat_obs=10),
n_response = 2,
parameters=list(
individual = list(
vcov = matrix(c(1,0.5,0.5,1),nrow=2,ncol=2,byrow=TRUE)
),
residual = list(
vcov = matrix(c(1,0.5,0.5,1),nrow = 2,ncol = 2,byrow=TRUE)
)
)
)We haven’t specified any names (of responses or predictors). By default, simulate_population() will add a number to the default name to indicate which reponse variable it refers to, so here we have y1 and y2 for the response variable, and individual_effect1 and individual_effect2 etc.
data <- get_population_data(squid_data)
head(data)## y1 y2 individual_effect1 individual_effect2 residual1
## 1 1.7154548 -0.3400219 1.917971 0.7785752 -0.2025162
## 2 0.5628202 0.3404365 1.917971 0.7785752 -1.3551508
## 3 2.2598837 0.1039745 1.917971 0.7785752 0.3419127
## 4 1.2400569 1.7172877 1.917971 0.7785752 -0.6779141
## 5 0.8401039 0.2945199 1.917971 0.7785752 -1.0778671
## 6 2.6972366 0.0814206 1.917971 0.7785752 0.7792656
## residual2 individual squid_pop
## 1 -1.1185971 1 1
## 2 -0.4381387 1 1
## 3 -0.6746007 1 1
## 4 0.9387125 1 1
## 5 -0.4840553 1 1
## 6 -0.6971546 1 1We can name the response variables easily, by giving the response_name argument a vector of names
squid_data <- simulate_population(
data_structure=make_structure(structure = "individual(100)",repeat_obs=10),
n_response = 2,
response_name = c("body_mass","behaviour"),
parameters=list(
individual = list(
vcov = matrix(c(1,0.5,0.5,1),nrow=2,ncol=2,byrow=TRUE)
),
residual = list(
vcov = matrix(c(1,0.5,0.5,1),nrow = 2,ncol = 2,byrow=TRUE)
)
)
)
data <- get_population_data(squid_data)
head(data)## body_mass behaviour individual_effect1 individual_effect2 residual1
## 1 -1.7790297 -2.1313181 -0.2306292 0.2132262 -1.5484005
## 2 0.9989910 1.5914308 -0.2306292 0.2132262 1.2296201
## 3 0.5519600 0.4977197 -0.2306292 0.2132262 0.7825892
## 4 0.8959494 0.9204513 -0.2306292 0.2132262 1.1265785
## 5 0.5997486 0.8608984 -0.2306292 0.2132262 0.8303778
## 6 -1.4619544 -0.8689692 -0.2306292 0.2132262 -1.2313253
## residual2 individual squid_pop
## 1 -2.3445443 1 1
## 2 1.3782046 1 1
## 3 0.2844935 1 1
## 4 0.7072251 1 1
## 5 0.6476722 1 1
## 6 -1.0821955 1 1