1.6 Non-Gaussian phenotypes

To simulate non-Gaussian data, we can specify a link function and a family as arguments to simulate_population(). Underneath the predictors are being simulated as multivariate normal (on the latent scale), and then the resulting phenotype is transformed (onto the expected scale) and then binomial or Poisson sampling is applied (the observed scale).

Here is an example to simulate Poisson distributed data: \[ y_i \sim Poisson(\lambda_i) \] \[ \lambda_i = exp( \beta_0 + \boldsymbol{x}_{i} \boldsymbol{\beta} + \epsilon_i ) \] \[ \boldsymbol{x}_i \sim \mathcal{N}(\boldsymbol{\mu}_x, \Sigma_x) \] \[ \epsilon_i \sim \mathcal{N}(0,\sigma^2_\epsilon) \]

The only change in the code that is needed is the addition of the link and family arguments.

squid_data <- simulate_population(
  parameters = list(
    observation = list(
      names = c("temperature","rainfall"),
      beta = c(0.2,0.1)
    ),
    residual = list(
      mean = 1.75,
      vcov = 0.2
    )
  ),
  n = 2000,
  family = "poisson", 
  link = "log"
)

data <- get_population_data(squid_data)
head(data)

##   y temperature    rainfall  residual squid_pop
## 1 0   0.2490819  0.18108147 0.8249614         1
## 2 3  -0.3059614 -0.01854387 1.0974162         1
## 3 3  -1.9080617 -0.10657783 1.4348215         1
## 4 5  -0.3827834 -0.63354609 2.0413375         1
## 5 5   2.1639215  0.02263793 0.8399427         1
## 6 4  -0.2962397 -0.14766836 1.6809449         1

plot(table(data$y), ylab="Frequency", xlab="z")

glm(y ~ temperature + rainfall, data, family="poisson")

## 
## Call:  glm(formula = y ~ temperature + rainfall, family = "poisson", 
##     data = data)
## 
## Coefficients:
## (Intercept)  temperature     rainfall  
##      1.8545       0.1857       0.1124  
## 
## Degrees of Freedom: 1999 Total (i.e. Null);  1997 Residual
## Null Deviance:       5169 
## Residual Deviance: 4537  AIC: 11550

Available families are ‘gaussian’, ‘poisson’ or ‘binomial’ and link functions ‘identity’, ‘log’, ‘inverse’, ‘sqrt’, ‘logit’, ‘probit’.