index.Rmd

---
title: "Creating simulated data sets in R"
author: "Brad Duthie"
date: "04 February 2026"
output:
  pdf_document: default
  html_document:
    code_folding: show
  word_document: default
bibliography: refs.bib
---

Contents
================================================================================

********************************************************************************

**The ability to simulate data is a useful tool for better understanding statistical analyses and planning experimental designs. These notes illustrate how to simulate data using a variety of different functions in the R programming language, then discuss how data simulation can be used in research. These notes borrow heavily from a Stirling Coding Club session on [randomisation](https://stirlingcodingclub.github.io/randomisation/randomisation_notes.html), and to a lesser extent from a session on [linear models](https://stirlingcodingclub.github.io/linear_modelling/). After working through these notes, the reader should be able to simulate their own data sets and use them to explore data visualisations and statistical analysis. These notes are also available as a  [PDF](https://stirlingcodingclub.github.io/simulating_data/index.pdf).**

********************************************************************************

- [Introduction: Simulating data](#intro)
- [Univariate random numbers](#univariate) 
    - [Random uniform: `runif`](#runif)
    - [Random normal: `rnorm`](#rnorm)
    - [Random poisson: `rpois`](#rpois)
    - [Random binomial: `rbinom`](#rbinom)
- [Random sampling using `sample`](#sample)
    - [Sampling random numbers from a list](#sample_numbers)
    - [Sampling random characters from a list](#sample_characters)
- [Simulating data with known correlations](#mvtnorm)
- [Simulating a full data set](#fullsim)
- [Conclusions](#conclusions)
- [Literature Cited](#refs)

********************************************************************************

<a name="intro">Introduction: Simulating data</a>
================================================================================

The ability generate simulated data is very useful in a lot of research contexts. Simulated data can be used to [better understand statistical methods](https://stirlingcodingclub.github.io/linear_modelling/#makeup), or in some cases to actually [run statistical analyses](https://stirlingcodingclub.github.io/randomisation/randomisation_notes.html#rand_m) (e.g., simulating a null distribution against which to compare a sample). Here I want to demonstrate how to simulate data in R. This can be accomplished with base R functions including [`rnorm`](#rnorm), [`runif`](#runif), [`rbinom`](#rbinom), [`rpois`](#rpois), or `rgamma`; all of these functions sample univariate data (i.e., one variable) from a specified distribution. The function [`sample`](#sample) can be used to sample elements from an R object with or without replacement. Using the [MASS](https://cran.r-project.org/web/packages/MASS/MASS.pdf) library, the [`mvtnorm`](#mvtnorm) function will sample multiple variables with a known correlation structure (i.e., we can tell R how variables should be correlated with one another) and normally distributed errors.

Below, I will first demonstrate how to use some common functions in R for simulating data. Then, I will illustrate how these simulated data might be used to better understand common statistical analyses and data visualisation.


<a name="univariate">Univariate random numbers</a>
================================================================================

Below, I introduce some base R functions that simulate (pseudo)random numbers from a given distribution. Note that most of what follows in this section is a recreation of a similar section in the notes for [randomisation](https://stirlingcodingclub.github.io/randomisation/randomisation_notes.html) analysis in R.

<a name="runif">**Sampling from a uniform distribution**</a>

The `runif` function returns some number (`n`) of random numbers from a [uniform distribution](https://en.wikipedia.org/wiki/Uniform_distribution_(continuous)) with a range from $a$ (`min`) to $b$ (`max`) such that $X\sim\mathcal U(a,b)$ (verbally, $X$ is sampled from a uniform distribution with the parameters $a$ and $b$), where $-\infty < a < b < \infty$ (verbally, $a$ is greater than negative infinity but less than $b$, and $b$ is finite). The default is to draw from a standard uniform distribution (i.e., $a = 0$ and $b = 1$) as done below.

```{r}
rand_unifs_10 <- runif(n = 10, min = 0, max = 1);
```

The above code stores a vector of ten numbers `rand_unifs_10`, shown below. Note that the numbers will be different each time we re-run the `runif` function above.

```{r, echo = FALSE}
print(rand_unifs_10);
```

We can visualise the standard uniform distribution that is generated by plotting a histogram of a very large number of values created using `runif`.

```{r}
rand_unifs_10000 <- runif(n = 10000, min = 0, max = 1);
hist(rand_unifs_10000, xlab = "Random value (X)", col = "grey",
     main = "", cex.lab = 1.5, cex.axis = 1.5);
```

The random uniform distribution is special in some ways. The algorithm for generating random uniform numbers is the starting point for generating random numbers from other distributions using methods such as [rejection sampling](https://en.wikipedia.org/wiki/Rejection_sampling), [inverse transform sampling](https://en.wikipedia.org/wiki/Inverse_transform_sampling), or the [Box Muller](https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform) method [@Box1958]. 

<a name="rnorm">**Sampling from a normal distribution**</a>

The `rnorm` function returns some number (`n`) of randomly generated values given a set mean ($\mu$; `mean`) and standard deviation ($\sigma$; `sd`), such that $X\sim\mathcal N(\mu,\sigma^2)$. The default is to draw from a standard [normal (a.k.a., "Gaussian") distribution](https://en.wikipedia.org/wiki/Normal_distribution) (i.e., $\mu = 0$ and $\sigma = 1$).

```{r}
rand_norms_10 <- rnorm(n = 10, mean = 0, sd = 1);
```

The above code stores a vector of 10 numbers, shown below.

```{r, echo = FALSE}
print(rand_norms_10);
```

We can verify that a standard normal distribution is generated by plotting a histogram of a very large number of values created using `rnorm`.

```{r}
rand_norms_10000 <- rnorm(n = 10000, mean = 0, sd = 1);
hist(rand_norms_10000, xlab = "Random value (X)", col = "grey",
     main = "", cex.lab = 1.5, cex.axis = 1.5);
```

Generating a histogram using data from a simulated distribution like this is often a useful way to visualise distributions, or to see how samples from the same distribution might vary. For example, if we wanted to compare the above distribution with a normal distribution that had a standard deviation of 2 instead of 1, then we could simply sample 10000 new values in `rnorm` with `sd = 2` instead of `sd = 1` and create a new histogram with `hist`. If we wanted to see what the distribution of sampled data might look like given a low sample size (e.g., 10), then we could repeat the process of sampling from `rnorm(n = 10, mean = 0, sd = 1)` multiple times and looking at the shape of the resulting histogram.


<a name="rpois">**Sampling from a poisson distribution**</a>

Many processes in biology can be described by a [Poisson distribution](https://en.wikipedia.org/wiki/Poisson_distribution). A Poisson process describes events happening with some given probability over an area of time or space such that $X\sim Poisson(\lambda)$, where the rate parameter $\lambda$ is both the mean and variance of the Poisson distribution (note that by definition, $\lambda > 0$, and although $\lambda$ can be any positive real number, data are always integers, as with count data). Sampling from a Poisson distribution can be done in R with `rpois`, which takes only two arguments specifying the number of values to be returned (`n`) and the rate parameter (`lambda`).

```{r}
rand_poissons <- rpois(n = 10, lambda = 1.5);
print(rand_poissons);
```

There are no default values for `rpois`. We can plot a histogram of a large number of values to see the distribution when $\lambda = 4.5$ below.


```{r}
rand_poissons_10000 <- rpois(n = 10000, lambda = 4.5);
hist(rand_poissons_10000, xlab = "Random value (X)", col = "grey",
     main = "", cex.lab = 1.5, cex.axis = 1.5);
```

<a name="rbinom">**Sampling from a binomial distribution**</a>

Sampling from a binomial distribution in R with `rbinom` is a bit more complex than using `runif`, `rnorm`, or `rpois`. Like those previous functions, the `rbinom` function returns some number (`n`) of random numbers, but the arguments and output can be slightly confusing at first. Recall that a [binomial distribution](https://en.wikipedia.org/wiki/Binomial_distribution) describes the number of 'successes' for some number of independent trials ($\Pr(success) = p$). The `rbinom` function returns the number of successes after `size` trials, in which the probability of success in each trial is `prob`. For a concrete example, suppose we want to simulate the flipping of a fair coin 1000 times, and we want to know how many times that coin comes up heads ('success'). We can do this with the following code.

```{r}
coin_flips <- rbinom(n = 1, size = 1000, prob = 0.5);
print(coin_flips);
```

The above result shows that the coin came up heads `r coin_flips` times. Note, however, the (required) argument `n` above. This allows the user to set the number of sequences to run. In other words, if we set `n = 2`, then this could simulate the flipping of a fair coin 1000 times once to see how many times heads comes up, then repeating the whole process a second time to see how many times heads comes up again (or, if it is more intuitive, the flipping of two separate fair coins 1000 times).

```{r}
coin_flips_2 <- rbinom(n = 2, size = 1000, prob = 0.5);
print(coin_flips_2);
```

In the above, a fair coin was flipped 1000 times and returned `r coin_flips_2[1]` heads, and then another fair coin was flipped 1000 times and returned `r coin_flips_2[2]` heads. As with the `rnorm` and `runif` functions, we can check to see what the distribution of the binomial function looks like if we repeat this process. Suppose, in other words, that we want to see the distribution of the number of times heads comes up after 1000 flips. We can, for example, simulate the process of flipping 1000 times in a row with 10000 different coins using the code below.

```{r}
coin_flips_10000 <- rbinom(n = 10000, size = 1000, prob = 0.5);
```

I have not printed the above `coin_flips_10000` for obvious reasons, but we can use a histogram to look at the results.

```{r}
hist(coin_flips_10000, xlab = "Random value (X)", col = "grey",
     main = "", cex.lab = 1.5, cex.axis = 1.5);
```

As would be expected, most of the time 'heads' occurs around 500 times out of 1000, but usually the actual number will be a bit lower or higher due to chance. Note that if we want to simulate the results of individual flips in a single trial, we can do so as follows.

```{r}
flips_10 <- rbinom(n = 10, size = 1, prob = 0.5);
```
```{r, echo = FALSE}
print(flips_10);
```

In the above, there are `n = 10` trials, but each trial consists of only a single coin flip (`size = 1`). But we can equally well interpret the results as a series of `n` coin flips that come up either heads (`1`) or tails (`0`). This latter interpretation can be especially useful to write code that randomly decides whether some event will happen (`1`) or not (`0`) with some probability `prob`.


<a name="sample">Random sampling using `sample`</a>
================================================================================

Sometimes it is useful to sample a set of values from a vector or list. The R function `sample` is very flexible for sampling a subset of numbers or elements from some structure (`x`) in R according to some set probabilities (`prob`). Elements can be sampled from `x` some number of times (`size`) with or without replacement (`replace`), though an error will be returned if the `size` of the sample is larger than `x` but `replace = FALSE` (default). 

<a name="sample_numbers">**Sampling random numbers from a list**</a>

To start out simple, suppose we want to ask R to pick a random number from one to ten with equal probability.

```{r}
rand_number_1 <- sample(x = 1:10, size = 1);
print(rand_number_1);
```

The above code will set `rand_number_1` to a randomly selected value, in this case `r rand_number_1`. Because we have not specified a probability vector `prob`, the function assumes that every element in `1:10` is sampled with equal probability. We can increase the `size` of the sample to `10` below.

```{r}
rand_number_10 <- sample(x = 1:10, size = 10);
print(rand_number_10);
```

Note that all numbers from 1 to 10 have been sampled, but in a random order. This is becaues the default is to sample with replacement, meaning that once a number has been sampled for the first element in `rand_number_10`, it is no longer available to be sampled again. To change this and allow for sampling with replacement, we can change the default.

```{r}
rand_number_10_r <- sample(x = 1:10, size = 10, replace = TRUE);
print(rand_number_10_r);
```

Note that the numbers {`r as.numeric(names(which(table(rand_number_10_r) > 1)))`} are now repeated in the set of randomly sampled values above. We can also specify the probability of sampling each element, with the condition that these probabilities need to sum to 1. Below shows an example in which the numbers 1-5 are sampled with a probability of 0.05, while the numbers 6-10 are sampled with a probability of 0.15, thereby biasing sampling toward larger numbers.

```{r}
prob_vec      <- c( rep(x = 0.05, times = 5), rep(x = 0.15, times = 5) );
rand_num_bias <- sample(x = 1:10, size = 10, replace = TRUE, prob = prob_vec);
print(rand_num_bias);
```

Note that `rand_num_bias` above contains more numbers from 6-10 than from 1-5.

<a name="sample_characters">**Sampling random characters from a list**</a>

Sampling characters from a list of elements is no different than sampling numbers, but I am illustrating it separately because I find that I often sample characters for conceptually different reasons. For example, if I want to create a simulated data set that includes three different species, I might create a vector of species identities from which to sample.

```{r}
species <- c("species_A", "species_B", "species_C");
```

This gives three possible categories, which I can now use `sample` to draw from. Assume that I want to simulate the sampling of these three species, perhaps with `species_A` being twice as common as `species_B` and `species_C`. I might use the following code to sample 24 times.

```{r}
sp_sample <- sample(x = species, size = 24, replace = TRUE, 
                    prob = c(0.5, 0.25, 0.25) 
                    );
```

Below are the values that get returned.

```{r, echo = FALSE}
print(sp_sample);
```


<a name="mvtnorm">Simulating data with known correlations</a>
================================================================================

We can generate variables $X_{1}$ and $X_{2}$ that have known correlations $\rho$ with with one another. The code below does this for two standard normal random variables with a sample size of 10000, such that the correlation between them is 0.3.

```{r}
N   <- 10000;
rho <- 0.3;
x1  <- rnorm(n = N, mean = 0, sd = 1);
x2  <- (rho * x1) + sqrt(1 - rho*rho) * rnorm(n = N, mean = 0, sd = 1);
```

Mathematically, these variables are generated by first simulating the sample $x_{1}$ (`x1` above) from a standard normal distribution. Then, $x_{2}$ (`x2` above) is calculated as below,

$x_{2} = \rho x_{1} + \sqrt{1 - \rho^{2}}x_{rand},$

Where $x_{rand}$ is a sample from a normal distribution with the same variance as $x_{1}$. A simple call to the R function `cor` will confirm that the correlation does indeed equal `rho` (with some sampling error).

```{r}
cor(x1, x2);
```


This is useful if we are only interested in two variables, but there is a much more efficient way to generate any number of variables with different variances and correlations to one another. To do this, we need to use the [MASS](https://cran.r-project.org/web/packages/MASS/MASS.pdf) library, which can be installed and loaded as below.

```{r, eval = FALSE}
install.packages("MASS");
library("MASS");
```

```{r, echo = FALSE}
library("MASS");
```

In the [MASS](https://cran.r-project.org/web/packages/MASS/MASS.pdf) library, the function `mvrnorm` can be used to generate any number of variables for a pre-specified covariance structure.

```{r, echo = FALSE}
bumpus <- read.csv("Bumpus_data.csv");
c_bump <- cov(bumpus[,3:5]);
```

Suppose we want to simulate a data set of three measurements from a species of organisms. Measurement 1 ($M_1$) has a mean of $\mu_{M_{1}} =$ `r round(mean(bumpus[,3]), digits = 2)` and variance of $Var(M_1) =$ `r round(c_bump[1, 1], digits = 2)`, measurement 2 ($M_2$) has a mean of $\mu_{M_{1}} =$ `r round(mean(bumpus[,4]), digits = 2)` and variance of $Var(M_2) =$ `r round(c_bump[2, 2], digits = 2)`, and measurement 3 ($M_2$) has a mean of $\mu_{M_{1}} =$ `r round(mean(bumpus[,5]), digits = 2)` and variance of $Var(M_3) =$ `r round(c_bump[3, 3], digits = 2)`. Below is a table summarising.

```{r, echo = FALSE}
library(knitr);
mns <- round(apply(X = bumpus[,3:5], MARGIN = 2, FUN = mean), digits = 2);
vrs <- round(apply(X = bumpus[,3:5], MARGIN = 2, FUN = var), digits = 2);
nms <- c("M1", "M2", "M3");
dfrm <- data.frame(measurement = nms, mean = mns, variance = vrs);
rownames(dfrm) <- NULL;
kable(dfrm);
```


Further, we want the covariance between $M_1$ and $M_2$ to equal $Cov(M_{1}, M_{2}) =$ `r round(c_bump[1, 2], digits = 2)`, the covariance between $M_1$ and $M_3$ to equal $Cov(M_{1}, M_{3}) =$ `r round(c_bump[1, 3], digits = 2)`, and the covariance between $M_2$ and $M_3$ to equal $Cov(M_{2}, M_{3}) =$ `r round(c_bump[2, 3], digits = 2)`. We can put all of this information into a [covariance matrix](https://en.wikipedia.org/wiki/Covariance_matrix) $\textbf{V}$ with three rows and three columns. The diagonal of the matrix holds the variances of each variable, with the off-diagonals holding the covariances (note also that the variance of a variable $M$ is just the variable's covariance with itself; e.g., $Var(M_{1}) = Cov(M_{1}, M_{1})$).

$$
V = \begin{pmatrix}
  Var(M_{1}), & Cov(M_{1}, M_{2}), & Cov(M_{1}, M_{3}) \\
  Cov(M_{2}, M_{1}), & Var(M_{2}), & Cov(M_{2}, M_{3}) \\
  Cov(M_{3}, M_{1}), & Cov(M_{3}, M_{2}), & Var(M_{3}) \\
\end{pmatrix}.
$$

In R, we can create this matrix as follows.

```{r}
matrix_data <- c(12.68, 13.95, 3.07, 13.95, 30.39, 4.70, 3.07, 4.70, 2.18);
cv_mat      <- matrix(data = matrix_data, nrow = 3, ncol = 3, byrow = TRUE);
rownames(cv_mat) <- c("M1", "M2", "M3");
colnames(cv_mat) <- c("M1", "M2", "M3");
```

Here is what `cv_mat` looks like (note that it is symmetrical along the diagonal).

```{r, echo = FALSE}
print(cv_mat);
```

Now we can add the means to a vector in R.

```{r}
mns <- c(159.54, 245.26, 25.52);
```

We are now ready to use the `mvrnorm` function in R to simulate some number `n` of sampled organisms with these three measurements. We use the `mvrnorm` arguments `mu` and `Sigma` to specify the vector of means and covariance matrix, respectively.

```{r}
sim_data <- mvrnorm(n = 40, mu = mns, Sigma = cv_mat);
```

Here are the example data below.

```{r, echo = FALSE}
kable(sim_data);
```

We can check to confirm that the mean values of each column are correct using `apply`.

```{r}
apply(X = sim_data, MARGIN = 2, FUN = mean);
```

And we can check to confirm that the covariance structure of the data is correct using `cov`.

```{r}
cov(sim_data);
```

Note that the values are not exact, but should become closer to the specified values as increase the sample size `n`. We can visualise the data too; for example, we might look at the close correlation between $M_{1}$ and $M_{2}$ using a scatterplot, just as we would for data sampled from the field.

```{r}
par(mar = c(5, 5, 1, 1));
plot(x = sim_data[,1], y = sim_data[,2], pch = 20, cex = 1.25, cex.lab = 1.25,
     cex.axis = 1.25, xlab = expression(paste("Value of ", M[1])),
     ylab = expression(paste("Value of ", M[2])));
```

We could even run an ordination on these simulated data. For example, we could extract the principle components with `prcomp`, then plot the first two PCs to visualise these data. We might, for example, want to compare different methods of ordination using a data set with different, pre-specified properties [e.g., @Minchin1987]. We might also want to use simulated data sets to investigate how different statistical tools perform. I show this in the next section, where I put a full data set together and run [linear models](https://stirlingcodingclub.github.io/linear_modelling/) on it.


<a name="fullsim">Simulating a full data set</a>
================================================================================

Putting everything together, here I will create a data set of three different species from which three different measurements are taken. We can just call these measurements 'length', 'width', and 'mass'. For simplicity, let us assume that these measurements always covary in the same way that we saw with $\textbf{V}$ (i.e., `cv_mat`) above. But let's also assume that we have three species with slightly different mean values. Below is the code that will build a new data set of $N = 20$ samples with four columns: species, length, width, and mass. 

```{r}
N           <- 20;
matrix_data <- c(12.68, 13.95, 3.07, 13.95, 30.39, 4.70, 3.07, 4.70, 2.18);
cv_mat      <- matrix(data = matrix_data, nrow = 3, ncol = 3, byrow = TRUE);
mns_1       <- c(159.54, 245.26, 25.52);
sim_data_1  <- mvrnorm(n = N, mu = mns, Sigma = cv_mat);
colnames(sim_data_1) <- c("Length", "Width", "Mass");
# Below, I bind a column for indicating 'species_1' identity
species     <- rep(x = "species_1", times = 20); # Repeats 20 times
sp_1        <- data.frame(species, sim_data_1);
```

Let us add one more data column. Suppose that we can also sample the number of offspring each organism has, and that the mean number of offspring that an organism has equals one tenth of the organism's mass. To do this, we can use `rpois`, and take advantage of the fact that the argument `lambda` can be a vector rather than a single value. So to get the number of offspring for each organism based on its body mass, we can just insert the mass vector `sp_1$Mass` times 0.1 for `lambda`.

```{r}
offspring   <- rpois(n = N, lambda = sp_1$Mass * 0.1);
sp_1        <- cbind(sp_1, offspring);
```

I have also bound the offspring number to the data set `sp_1`. Here is what it looks like below.

```{r, echo = FALSE}
kable(sp_1);
```

To add two more species, let us repeat the process two more times, but change the expected mass just slightly each time. The code below does this, and puts everything together in a single data set.

```{r}
# First making species 2
mns_2       <- c(159.54, 245.26, 25.52 + 3); # Add a bit
sim_data_2  <- mvrnorm(n = N, mu = mns, Sigma = cv_mat);
colnames(sim_data_2) <- c("Length", "Width", "Mass");
species     <- rep(x = "species_2", times = 20); # Repeats 20 times
offspring   <- rpois(n = N, lambda = sim_data_2[,3] * 0.1);
sp_2        <- data.frame(species, sim_data_2, offspring);
# Now make species 3
mns_3       <- c(159.54, 245.26, 25.52 + 4.5); # Add a bit more
sim_data_3  <- mvrnorm(n = N, mu = mns, Sigma = cv_mat);
colnames(sim_data_3) <- c("Length", "Width", "Mass");
species     <- rep(x = "species_3", times = 20); # Repeats 20 times
offspring   <- rpois(n = N, lambda = sim_data_3[,3] * 0.1);
sp_3        <- data.frame(species, sim_data_3, offspring);
# Bring it all together in one data set
dat <- rbind(sp_1, sp_2, sp_3);
```

Our full data set now looks like the below.

```{r, echo = FALSE}
kable(dat);
```

To summarise, we now have a simulated data set of measurements from three different species, all of which have known variances and covariances of length, width, and mass. Each species has a slightly different mean mass, and for all species, each unit of mass increases the expected number of offspring by 0.1. Because we know these properties of the data for certain, we can start asking questions that might be useful to know about our data analysis. For example, given this covariance structure and these small differences in mass, is a sample size of 20 really enough to even get a significant difference among species masses using an ANOVA? What if we tried to test for differences among masses using some sort of [randomisation approach](https://stirlingcodingclub.github.io/randomisation/randomisation_notes.html#rand_h) Instead? Would this give us more or less power? Let us run an ANOVA to see if the difference between group means (which we know exists) is recovered.

```{r}
aov_result <- aov(Mass ~ species, data = dat);
summary(aov_result);
```

It appears not! What about the relationship between body mass and offspring production that we know exists? Below is a scatterplot of the data for the three different species.

```{r, echo = FALSE}
par(mar = c(5, 5, 1, 1));
s1 <- dat[dat[,1] == "species_1",];
s2 <- dat[dat[,1] == "species_2",];
s3 <- dat[dat[,1] == "species_3",];
plot(x = s1$Mass, y = s1$offspring, pch = 15, cex = 1.25, cex.lab = 1.25,
     cex.axis = 1.25, xlab = "Species mass", ylab = "Species offspring",
     xlim = c(min(dat$Mass) - 1, max(dat$Mass) + 1),
     ylim = c(0, max(dat$offspring)), col = "black");
points(x = s2$Mass, y = s2$offspring, pch = 16, cex = 1.25, col = "red");
points(x = s3$Mass, y = s3$offspring, pch = 17, cex = 1.25, col = "blue");
legend("topleft", pch = 15:17, col = c("black", "red", "blue"), 
       legend = c("Species 1", "Species 2", "Species 3"), cex = 1.25);
```

This looks like there might be a positive relationship, but it is very difficult to determine just from the scatterplot. We can use a generalised linear model to test it with species as a random effect, as we might do if these were data sampled from the field (do not worry about the details of the model here; the key point is that we can use the simulated data with known properties to assess the performance of a statistical test).

```{r}
library(lme4);
mod <- glmer(offspring ~ Mass + (1 | species), data = dat, family = "poisson");
summary(mod);
```

There does not appear to be any effect here either! To get one, it appears that we will need to simulate a larger data set (or a bigger effect size -- or just get lucky when re-simulating a new data set). 

Note that I have run a linear model that might be reasonable given the structure of our data. But the advantage of working with simulated data and knowing for certain what the relationship is between the underlying variables is that we can explore different statistical techniques. For example, we know that our response variable `offspring` is count data, so we are supposed to specify a Poisson error structure using the `family = "poisson"` argument above, right? But what would happen if we just used a normal error structure anyway? Would this really be so bad? Now is the opportunity to test because we *know* what the correct answer is supposed to be! Trying statistical methods that are normally ill-advised can actually be useful here, as it can help us see for ourselves when a technique is bad -- or perhaps when it really is not [e.g., @Ives2015].

<a name="conclusions">Conclusions</a>
================================================================================

Simulating data can be a powerful tool for learning and investigating different statistical analyses. The main benefits of using simulated data are flexibility and certainty. Simulation gives us the flexibility to explore any number of hypotheticals, including different sample sizes, effect sizes, relationships between variables, and error distributions. It also works from a point of certainty; we know what the real relationship is between variables, and what the actual effect sizes are because we define them when generating random samples. So if we want to better understand what would happen if we were unable to sample an important variable in our system, or if we were to use a biased estimator, or if we we were to violate key model assumptions, simulated data is a very useful tool.


<a name="refs">Literature cited</a>
================================================================================