Lichti/STAT 503 – Statistical Methods for Biology Modified: 2019‐13‐03 1 STAT 503 – Statistical Methods for Biology Tutorial 5: Probability distributions and simulation Tutorial 5 introduces R functions for probability distributions and random simulation (what is sometimes called Monte Carlo simulation). These functions include the probability mass or density functions for different distributions, as well as the corresponding cumulative distribution functions, quantile functions, and functions that simulate random draws from a random variable. Working through this tutorial, you will learn: 1) The nomenclature of probability distribution functions in R 2) How to obtain values from the cumulative distribution function and mass function for a probability distribution given values for its parameters. 3) How to find the quantiles of a probability distribution. 4) How to simulate data from a variety of distributions. 5) About the sample() function, which allows you to randomly sample from a dataset 6) How to use the replicate() function to repeat a simulation experiment 7) How to use sample() and replicate() together for bootstrap resampling 8) (OPTIONAL) How to set up a simulation to help plan data collection 9) (OPTIONAL) How to simulate the mechanism underlying the normal distribution. You will also get more practice using mutate() and generating plots, and will use a simulation experiment to explore the origins of the normal distribution. Please attempt all of the exercises and ask for help if you encounter any problems or confusion. 1. Getting started In this tutorial, we will use the dplyr and ggplot2 packages. Remember that when you attach the dplyr package, it may mask a few items from the base and stats packages. This is okay. There is no data to load with this tutorial. 2. Probability distributions in R R contains probability and simulation functions for many different distributions, and all of distribution functions follow a standard naming convention: there is a one-letter prefix that tells you what the function does, followed by a root name that identifies the name of the distribution. For instance, the root name for the normal distribution in R is norm, and the binomial distribution is binom. R has four different functions for each distribution. These are introduced below using the binomial distribution, but the syntax is similar for other distributions as well. Lichti/STAT 503 – Statistical Methods for Biology Modified: 2019‐13‐03 2 2.1 An introduction using the binomial distribution For each distribution, R has four different functions, which are listed for the binomial distribution below: pbinom() – the cumulative distribution function. The prefix “p” is for “probability.” It is used for the cdf instead of the mass function because the cdf can be used to find probabilities for any distribution, whether it is discrete or continuous. dbinom() – the “density” function. For discrete distributions, the “d” function returns the value of the probability mass function, which is conceptually analogous to the density function for a continuous distribution. qbinom() – the “quantile” function. This function returns the value of the variable that corresponds to a given quantile of the distribution. rbinom() – the “random” simulation function. The random simulator draws random numbers from the distribution. We will use it quite a bit later in this tutorial. To see how each of these functions works, we’ll set up an example using the binomial distribution with sample size 6 and probability of success 0.3. First, we’ll calculate the probability mass function and cumulative distribution function manually. Then we’ll see how to get the same results using the R functions listed above. 2.1.1 Probability functions – Start by loading dplyr. It is not required for the distributions, but we will use it to organize our work. Once it is loaded, create a data frame with 1 column, named x, containing the integers from 0 to 6 (see code below). Next, add a column for the probability mass function. We’ll call this column px_manual. You can add it to bin_example using the mutate() function: library(dplyr) n <- 6 p <- 0.3 bin_example <- data.frame(x = 0:n) bin_example <- mutate( bin_example, px_manual = choose(n, x) * p^x * (1 - p)^(n - x) ) (if you need a refresher on mutate(), see Tutorial 2, Section 3.5). In the code above, I defined objects for and , and set them equal to 6 and 0.3, respectively. I could have put 6 and 0.3 directly into the code wherever I need them, but defining them as objects this way means that if I want to change the values, I only need to edit the code in one place. This minimizes my opportunity to make mistakes. For example, I cannot accidentally change in the equation for the mass function, but forget to change it in the original definition for the variable . Lichti/STAT 503 – Statistical Methods for Biology Modified: 2019‐13‐03 3 As a check to make sure that we've done everything correctly, we can get the sum of . It should be 1 (and it is, more or less – if you actually test the equality, R will return FALSE; it is off by about 3 10 due to the limitations of computer precision): > bin_example x px_manual 1 0 0.117649 2 1 0.302526 3 2 0.324135 4 3 0.185220 5 4 0.059535 6 5 0.010206 7 6 0.000729 > sum(bin_example$px_manual) [1] 1 Next, we will use the cumsum() function to calculate the cumulative distribution function (cdf). Recall that the cdf of a random variable is defined as and describes the probability that the random variable will take a value less than or equal to . The name cumsum stands for cumulative sum; the function adds the values of a vector up in order of appearance, returning a running summation: bin_example <- mutate( bin_example, Fx_manual = cumsum(px_manual) ) Running this and looking at the result, we have: > bin_example x px_manual Fx_manual 1 0 0.117649 0.117649 2 1 0.302526 0.420175 3 2 0.324135 0.744310 4 3 0.185220 0.929530 5 4 0.059535 0.989065 6 5 0.010206 0.999271 7 6 0.000729 1.000000 Notice that in the first line, px_manual and Fx_manual are equal, but in the second row, Fx_manual is the sum of the first two values in px_manual, the third value of Fx_manual is the sum of the first three px_manual values, etc. Now we have values for both and that we’ve calculated by hand. Let’s compare these to the equivalent functions from R. We’ll add the values from both R functions in one mutate() call: bin_example <- mutate( bin_example, px_R = dbinom(x = x, size = n, prob = p), Fx_R = pbinom(q = x, size = n, prob = p) ) Lichti/STAT 503 – Statistical Methods for Biology Modified: 2019‐13‐03 4 The arguments in these functions have different names (x for a generic variable value in dbinom(), versus q for a quantile in pbinom()), but are really the same thing. They both represent the value of the variable at which you want to evaluate the function. This is always true: the first argument for any probability function in R that starts with the prefix d is x, and the first argument for any function starting with the prefix p is q. After the first argument come one or more arguments that specify the parameters for the particular distribution. These will be different for different distributions. For the binomial, the parameters are the sample size ( ) and the probability of success ( ). Looking at the results for the binomial distribution, we see that the values perfectly match our manual calculations: > bin_example x px_manual Fx_manual px_R Fx_R 1 0 0.117649 0.117649 0.117649 0.117649 2 1 0.302526 0.420175 0.302526 0.420175 3 2 0.324135 0.744310 0.324135 0.744310 4 3 0.185220 0.929530 0.185220 0.929530 5 4 0.059535 0.989065 0.059535 0.989065 6 5 0.010206 0.999271 0.010206 0.999271 7 6 0.000729 1.000000 0.000729 1.000000 So, if we get the same answer, why do we need dbinom() or pbinom() First, they’re convenient. If you are working with large datasets, they may also run more quickly and do a better job of managing your computer’s memory than manual calculations do. Finally, although we can easily perform the manual calculations for the binomial, we cannot do this with other distributions. 2.1.2 Bad arguments – You may be wondering what happens if you put in an impossible value for , , or . Try different inputs and find out. For example, > dbinom(x = -2, size = 6, prob = 0.3) [1] 0 > dbinom(x = 2, size = -1, prob = 0.3) [1] NaN Warning message: In dbinom(x = 2, size = -1, prob = 0.3) : NaNs produced > dbinom(x = 2, size = 6.3, prob = 0.3) [1] NaN Warning message: In dbinom(x = 2, size = 6.3, prob = 0.3) : NaNs produced > dbinom(x = 2, size = 6, prob = 2) [1] NaN Warning message: In dbinom(x = 2, size = 6, prob = 2) : NaNs produced In the first row, I gave R a value of that falls outside the domain of the distribution, and it correctly returned a probability density of 0. In the others, I gave it nonsense parameter values, and it returned NaN, which stands for “Not a Number,” and is treated by R like a missing value. Lichti/STAT 503 – Statistical Methods for Biology Modified: 2019‐13‐03 5 2.1.3 The lower.tail and log arguments – Remember that p functions return the cumulative distribution function for a variable, . All cdf functions in R have a lower.tail argument that has a default value of TRUE. Changing it to lower.tail = FALSE will return 1 . For example, > pbinom(q = 3, size = 6, prob = 0.3) [1] 0.92953 > pbinom(q = 3, size = 6, prob = 0.3, lower.tail = FALSE) [1] 0.07047 In addition, both the probability (p) and density (d) functions have a log argument, which defaults to FALSE. If changed to log = TRUE, the natural log of the probability will be returned. This is especially useful if you need to program a maximum likelihood estimation. > pbinom(q = 3, size = 6, prob = 0.3, log = TRUE) [1] -0.0730762 > exp(-0.0730762) [1] 0.92953 2.1.4 Quantile functions – Recall that to find any quantile of a random variable, you can set and then solve for , where 0 1. This implies that the quantile function is the inverse of the cumulative distribution function. When we call pbinom(), our first argument is a value for , and we get a the smallest value of for which back. For qbinom() and all quantile functions in R, the first argument is the probability that we want the quantile for, and we get back. For example, # p returned from pbinom() above > qbinom(p = 0.92953, size = 6, prob = 0.3) [1] 3 The relationship between the quantiles and cumulative distribution function is illustrated in Figure 1. Note that, because the binomial is a discrete random variable, we may get the same result for several different values of , > qbinom(p = seq(0.86,0.95,0.01), size = 6, prob = 0.3) [1] 3 3 3 3 3 3 3 4 4 4 and the level of resolution will increase if is bigger, > qbinom(p = seq(0.86,0.95,0.01), size = 600, prob = 0.3) [1] 192 193 193 194 194 195 196 197 198 199 Lichti/STAT 503 – Statistical Methods for Biology Modified: 2019‐13‐03 6 Figure 1: The curves show cumulative distribution functions for two populations: ～ 100, 0.2 (black) and ～ 100, 0.7 (red). The horizontal dotted line corresponds to a cumulative probability of 0.9. Ninety-percent (90%) of the values for the first population fall at or below 25, which is the value of where the black curve crosses the horizontal line. Similarly, 90% of the second population falls at or below 76. These values can be obtained by running qbinom(0.9, size = 100, prob = 0.2) for the first population, or qbinom(0.9, size = 100, prob = 0.7) for the second. pbinom(25, size = 100, prob = 0.20) reverses the process, returning the value of the y-axis where the black line crosses the black, vertical, dashed line (the actual value returned is 0.9125). Code for this plot may be found in the Appendix. 2.1.5 Random number generation and simulation – The function rbinom() returns a vector of random draws from ～ , . We can use simulations like this in many ways. For instance, after we have fit a model, we can use simulations to help us visualize its implications or make predictions. We can also use simulations to help us determine whether a particular model is appropriate for a given dataset, and we can use simulations to help plan studies. Finally, we can use simulation to help understand statistical concepts. The first argument for all random generation functions is named n and specifies the number of draws that you want to generate. For example, > rbinom(n = 5, size = 6, prob = 0.3) [1] 3 0 2 3 1 Lichti/STAT 503 – Statistical Methods for Biology Modified: 2019‐13‐03 7 gives me a vector of 5 numbers drawn from the binomial distribution with sample size 6 and probability of success 0.3. Since these are random numbers, you will probably not get the same vector that I got. 2.2 The normal distribution The probability, density, quantile, and random generation functions for other distributions work very similarly to pbinom(), dbinom(), qbinom(), and rbinom(). The only difference, other than the distribution name, is that each distribution will have its own parameter values. For the normal distribution, these are the mean (mean) and standard deviation (sd). Thus, > dnorm(3, mean = 1.2, sd = 1) [1] 0.07895016 evaluates the probability density function, , for the normal distribution with a mean of 1.2 and a standard deviation of 1 at 3. The corresponding cumulative distribution function, > pnorm(3, mean = 1.2, sd = 1) [1] 0.9640697 returns the probability that a value drawn at random from ～ 1.2, 1 would be less than or equal to 3 (in other words, it returns the value of ). The values returned by the three functions are summarized in Figure 2. Figure 2: Illustration of the normal density, probability, and quantile functions, dnorm(), pnorm(), and qnorm(), for ~ 1.2, 1 . Code for the figure is in the Appendix. pnorm(3, 1.2, 1) is the area of the shaded region dnorm(3, 1.2, 1) is the height of the curve at x = 3qnorm(0.9641, 1.2, 1) is the upper bound of the shaded region 0.0 0.1 0.2 0.3 0.4 -4 -2 0 2 4 6 x f x Lichti/STAT 503 – Statistical Methods for Biology Modified: 2019‐13‐03 8 The quantile function works as it did with the binomial distribution. Applying it, we can see that the empirical rule (95% of the outcomes within 2 standard deviations and 68% within 1 standard deviation of the mean) is only a rough approximation: > qnorm(c(0.025, 0.16, 0.5, 0.84, 0.975), mean=0, sd=1) [1] -1.9599640 -0.9944579 0.0000000 0.9944579 1.9599640 We can simulate values from the normal distribution using rnorm(). For example, > x <- rnorm(100000, mean = 1.2, sd = 1) > qplot(x, col=I(‘black’), fill=I(‘cornflowerblue’)) creates the histogram in Figure 3. All of the normal distribution functions in R use the default parameter values mean = 0 and sd = 1. Therefore, rnorm(1000) draws 1000 random values from the standard normal distribution. Figure 3: Histogram of 100,000 random draws from the normal distribution with a mean of 1.2 and a standard deviation of 1. 0 3000 6000 9000 12000 -2 0 2 4 6 x Lichti/STAT 503 – Statistical Methods for Biology Modified: 2019‐13‐03 9 2.3 Other distributions R provides functions for a large number of other distributions. Those that we may encounter in this class include the Poisson distribution, which is used to model the number of events that occur in some interval or unit (e.g., the number of tornados in a year, or the number of trees in a hectare of land), the distribution (sometimes spelled out as chi-square), which describes the distribution of sample sums-of-squares (the numerator of the sample variance) for draws from a standard normal distribution, the Student’s t distribution, and the F-distribution. Brief descriptions for several of these other distributions appear in Table 1. Although they are here primarily for future reference, I encourage you to spend some time playing with the simulation functions to see how each one behaves under different parameterizations (e.g., generate 10,000 random draws, and then draw a histogram, as I did for the normal distribution above). Table 1: Some important statistical distributions and their R representations. Distribution R name Parameters Application Uniform unif min – the minimum possible value max – the ma
ximum possible value The uniform distribution is helpful for randomizing sample selection and treatments, and for simulation Poisson pois lambda – the expected number of events in a given interval A discrete distribution that models count data Student’s t t df – degrees of freedom ( 1, inherited from the sample variance) The distribution of the sample mean given a normal population and unknown variance Chi-square chisq df – degrees of freedom Chi-square models the distribution of a sum-of-squares It appears in several statistical tests. F f df1 – numerator degrees of freedom df2 – denominator degrees of freedom F is the ratio of two chi-square random variables. Along with Student’s t, it is the main distribution used in ANOVA and regression. Lichti/STAT 503 – Statistical Methods for Biology Modified: 2019‐13‐03 10 3. Random sampling from a fixed set of values: sample() In addition to functions that allow random simulation from a known distribution, base R includes the sample() function, which provides a simple mechanism for randomly sampling from a particular set of values, either with or without replacement. 3.1 Random sampling without replacement Suppose that we have a vector of integers from 1 to 10: x <- 1:10 > x [1] 1 2 3 4 5 6 7 8 9 10 To simulate a simple random sample from x with sample size = 3, we can call, > sample(x, size = 3) [1] 9 4 6 > sample(x, size = 3) [1] 7 4 9 (this is a random simulation, so your results will vary). Each time we call sample(x, size = 3), a new set of 3 values will be randomly sampled from the vector x. This will work on non-numeric values as well. For example: > sample(letters, size = 3) [1] “q” “f” “r” 3.2 Permutation If the size argument is set equal to the length of x, or if it is omitted, then the result will be a vector that contains all of the values in x, but in a new, randomly selected order: > sample(x, size = length(x)) [1] 6 3 7 2 10 8 9 1 5 4 > sample(x) [1] 5 4 6 10 9 1 8 7 2 3 This type of sampling, in which we select a sample with the same size as the original set of values, and do so without replacement so that the result is a random reordering of the original set, is called permutation. 3.3 Sampling with replacement By default, sample() draws values without replacement, so attempting to set size to a value greater than the length of x will generate an error. However, setting the optional argument replace = TRUE will instruct sample() to draw values with replacement. When sampling with replacement, we may or may not get duplicate values when size length(x), > sample(x, size = 3, replace = TRUE) Lichti/STAT 503 – Statistical Methods for Biology Modified: 2019‐13‐03 11 [1] 8 1 2 > sample(x, size = 3, replace = TRUE) [1] 1 2 2 but duplicates are guaranteed if size > length(x): > sample(x, size = 10, replace = TRUE) [1] 5 5 4 1 3 8 1 4 4 4 > sample(x, size = 13, replace = TRUE) [1] 7 4 10 10 3 10 9 10 2 2 1 10 8 3.4 Probability-weighted samples All of the examples thus far have used simple random samples either with or without replacement. Although we will not use it in this class, the sample() function also allows you to set the probability or frequency associated with of each element of the vector that is being sampled. For instance, > sample(1:4, size = 20, replace = TRUE, prob = (1:4)/10) [1] 7 4 10 10 3 10 9 10 2 2 1 10 8 Here, I’ve instructed sample() to give me a vector of 20 values, where each element is an integer from 1 to 4, but the value 1 only has a 10% chance of occurring, whereas 2 has a 20% chance, 3 has a 30% chance, and 4 has a 40% chance. This can be seen by evaluating the argument given to the prob argument: > (1:4)/10 [1] 0.1 0.2 0.3 0.4 In the example that we’ve just looked at, I specified a proper probability mass function (the values sum to 1), but this is not strictly necessary. For instance, suppose we have data on ABO blood types in 5000 individuals, where 130 are type O, 2654 are type A, 1842 are type AB, and 374 are type B, and we would like to simulate a draw of 20 individuals from this population. One option would be to create a vector of 5000 values and take a simple random sample: abo <- c( rep("O", 130), rep("A", 2654), rep("AB", 1842), rep("B", 374) ) > sample(abo, size = 20) # note that replace = FALSE by default [1] “A” “AB” “AB” “B” “A” “A” “AB” “O” “O” “AB” “A” “AB” [13] “A” “A” “AB” “AB” “A” “AB” “A” “A” Alternatively, we can get the same random simulation process by calling, sample( c(“O”,”A”,”AB”,”B”), # values size = 20, # desired sample size replace = TRUE, # required b/c size > length(x) prob = c(130, 2654, 1842, 374) # frequencies in population ) Lichti/STAT 503 – Statistical Methods for Biology Modified: 2019‐13‐03 12 3.5 Example applications for sample() The sample function can be useful in a variety of situations. For example, as we have just seen, it can be used to simulate sampling from a known population. This allows us to get a preview of what our data might look like. As demonstrated in lecture, sample() can be used with the replicate() function (see Section 3.6) to see how well a particular estimation method might be affected by sample size or other aspects of sample design. 3.5.1 Randomizing treatment assignments – In experiments, sample() can be used to randomly assign individuals to treatments. For example, suppose that we are conducting a drug study in which we are looking at four different dosage levels, one of which is a control. Therefore, we have four treatments. For the sake of argument, suppose that we also have 100 subjects that have been selected to participate in the study, and that they have each been assigned a unique ID number from 1 to 100. The following code will randomly assign each individual to a treatment: treatment_levels <- c("control", "group 1", "group 2", "group 3") data <- data.frame( id = sample(1:100), treatment = factor(treatment_levels)) ) write.csv(data, file="treatment_assignments.csv", row.names=FALSE) The first statement creates a character vector that identifies each of the four treatment levels. Then, in the second statement, a data frame is created that has subject ID in the first column and a factor indicating treatment group in the second column. However, sample() has been called on the ID numbers to randomly permute them (see Section 3.2). The vector treatment_levels is converted into a factor, and is also recycled so that there are a total of 100 rows in the data frame (note that this would not work if 100/4 was not a whole number; it would produce an error). Finally, write.csv() outputs the result to a CSV file (you would need to use the standard procedures to make sure the file is written to the directory where you want it). If desired, arrange() could be called to sort the data frame by treatment and ID before saving it. 3.5.2 Bootstrap resampling – One common use of the sample() function is in bootstrap resampling. The theory underlying this methodology is that, assuming that you have a representative sample of independent units from a population, you can obtain an approximate confidence interval for almost any population parameter by repeatedly resampling from the observed data to build an empirical representation (i.e., a histogram) of the sampling distribution for the parameter's estimator. To generate a single bootstrap sample, we begin with some data. Here. I'll simulate some "data" that consist of 25 observations from a normally distributed variable with mean 34.9 and standard deviation of 2.71 (feel free to supply your own units): x <- rnorm(25, mean = 34.9, sd = 2.71) Lichti/STAT 503 – Statistical Methods for Biology Modified: 2019‐13‐03 13 > x [1] 30.68825 31.84423 33.18473 34.61263 33.14514 38.19584 [7] 40.99194 34.96162 32.17812 31.61524 40.35282 33.64944 [13] 30.15063 35.29006 36.36103 36.80707 41.95605 34.12107 [19] 30.75378 36.62187 37.20257 30.07368 33.98218 33.51193 [25] 31.98624 Once we have the data, we can sample from it with replacement to generate a new, simulated data vector that represents a plausible alternative dataset – that is, one that we could reasonably imagine having drawn from the population instead of the data that we really saw, given our actual data. xboot <- sample(x, replace = TRUE) > xboot [1] 30.68825 34.96162 37.20257 34.61263 35.29006 34.61263 [7] 35.29006 34.12107 31.84423 36.36103 33.18473 40.35282 [13] 35.29006 34.96162 34.61263 36.80707 32.17812 30.07368 [19] 34.61263 34.12107 34.61263 36.36103 33.98218 33.51193 [25] 36.62187 (compare these values to the original values for x, above). Taking the mean and standard deviation of the bootstrap sample gives us one new “observation” of the sample mean and sample standard deviation. The statistics calculated from the bootstrap sample are similar to the original values, but not identical to them: > mean(x) [1] 34.56953 > sd(x) [1] 3.334223 > mean(xboot) [1] 34.65073 > sd(xboot) [1] 2.140218 By repeating this process many times, we can build a frequency distribution for plausible values of the estimators, and . The easiest way to accomplish this is the replicate() function, as demonstrated in the next section. 4. Repetition of a random simulation: replicate() In most situations, it is not very useful to simulate a set a values only one time. For example, by itself, the single bootstrap sample generated in Section 3.5.2 does not provide us with very much useful information. To generate confidence intervals, we would need to repeat the process of drawing a new bootstrap sample many times, obtain estimates of the mean and standard deviation for each one, and then find the empirical quantiles from those means and standard deviations. The replicate() function makes it very simple to code this kind of exercise. However, before seeing how to use replicate(), it is helpful to see how this might be done if we didn’t have replicate() available. Lichti/STAT 503 – Statistical Methods for Biology Modified: 2019‐13‐03 14 4.1 The for loop A standard approach to this sort of task that works in almost every computer language (including R) is the for loop. As an example, look at the following code and its output: for(i in 1:3){ print(i + 1) } [1] 2 [1] 3 [1] 4 Notice that in the output, we have three separate 1-element vectors, equal to 1 + 1, 2 + 1, and 3 + 1. The for loop specifies a variable called an iterator, which I’ve named i, that takes on each value in a specified vector (here, the integers from 1 to 3), one at a time. Each time the iterator is assigned a new value, the block of code inside the curly braces gets executed. In this code, the first time through the loop, i is set to 1, so R prints i + 1 = 1 + 1 = 2. The second time through, it uses i = 2, and so on until it reaches the end of the last possible value of i. It is worth noting that the words “for” and “in” are reserved words in R. You cannot assign objects to these names. To generate a set of bootstrap samples using the for loop, we first create a list object to store our results, and then run the bootstrap sample code from Section 3.5.2 inside the loop, with a few minor modifications: xboot <- list() # an empty list named xboot for(i in 1:3){ xboot[[i]] <- sample(x, replace = TRUE) } > length(xboot) [1] 3 In the code above, we repeat the process of sampling x with replacement 3 times, and append each of the resulting vectors of 25 simulated observations to the list xboot. Notice the double square brackets used to index xboot; this notation is necessary for the indexing to work properly with lists. The result is a list of 3 vectors, each with 25 elements: > xboot [[1]] [1] 38.19584 36.36103 34.96162 34.96162 34.12107 36.36103 [7] 31.61524 30.07368 33.51193 40.35282 31.61524 30.07368 [13] 31.61524 31.84423 30.07368 33.14514 34.12107 33.51193 [19] 31.98624 35.29006 36.36103 34.61263 37.20257 40.35282 [25] 33.64944 [[2]] [1] 33.98218 38.19584 34.61263 35.29006 36.62187 36.62187 [7] 37.20257 31.98624 40.35282 36.80707 33.14514 38.19584 [13] 38.19584 34.61263 31.84423 34.12107 40.35282 37.20257 [19] 30.07368 34.96162 36.62187 36.80707 34.61263 35.29006 [25] 36.36103 Lichti/STAT 503 – Statistical Methods for Biology Modified: 2019‐13‐03 15 [[3]] [1] 31.84423 34.61263 30.15063 33.51193 30.07368 33.64944 [7] 30.68825 34.61263 40.99194 37.20257 31.98624 35.29006 [13] 34.12107 36.36103 40.35282 41.95605 33.18473 33.98218 [19] 35.29006 31.98624 30.75378 35.29006 31.84423 34.96162 [25] 32.17812 To make xboot more compact, we can rearrange it into a matrix using the simplify2array() function. This will give us a 25 3 matrix in which each column represents one of our 3 replicated samples: xboot <- simplify2array(xboot) > xboot [,1] [,2] [,3] [1,] 38.19584 33.98218 31.84423 [2,] 36.36103 38.19584 34.61263 [3,] 34.96162 34.61263 30.15063 [4,] 34.96162 35.29006 33.51193 [5,] 34.12107 36.62187 30.07368 [6,] 36.36103 36.62187 33.64944 [7,] 31.61524 37.20257 30.68825 [8,] 30.07368 31.98624 34.61263 [9,] 33.51193 40.35282 40.99194 [10,] 40.35282 36.80707 37.20257 [11,] 31.61524 33.14514 31.98624 [12,] 30.07368 38.19584 35.29006 [13,] 31.61524 38.19584 34.12107 [14,] 31.84423 34.61263 36.36103 [15,] 30.07368 31.84423 40.35282 [16,] 33.14514 34.12107 41.95605 [17,] 34.12107 40.35282 33.18473 [18,] 33.51193 37.20257 33.98218 [19,] 31.98624 30.07368 35.29006 [20,] 35.29006 34.96162 31.98624 [21,] 36.36103 36.62187 30.75378 [22,] 34.61263 36.80707 35.29006 [23,] 37.20257 34.61263 31.84423 [24,] 40.35282 35.29006 34.96162 [25,] 33.64944 36.36103 32.17812 And finally, we can use the apply() function to calculate the mean and standard deviation for each column (for the means, colMeans() is actually faster, but we’ll use apply() here). apply() executes a named function (FUN) on each row (MARGIN = 1) or column (MARGIN = 2) of a matrix: xboot_means <- apply(xboot, MARGIN = 2, FUN = mean) xboot_sds <- apply(xboot, MARGIN = 2, FUN = sd) > xboot_means [1] 34.23884 35.76285 34.27505 > xboot_sds [1] 2.894935 2.487915 3.205830 4.2 replicate() If all we want to do is repeat the same random process many times, then the replicate() function simplifies our code considerably. It takes two arguments: n, which specifies the number of times that code should be replicated, and a block of Lichti/STAT 503 – Statistical Methods for Biology Modified: 2019‐13‐03 16 code to replicate. If the code includes more than one line, it must be enclosed in curly braces ({}). Using replicate, only one line of code is required to generate 3 bootstrap samples and reformat them to a matrix: xboot2 <- replicate(3, sample(x, replace = TRUE)) > xboot2 [,1] [,2] [,3] [1,] 33.64944 33.51193 36.36103 [2,] 37.20257 35.29006 37.20257 [3,] 33.51193 30.15063 31.98624 [4,] 36.80707 36.80707 32.17812 [5,] 30.15063 36.80707 30.75378 [6,] 30.15063 33.18473 33.64944 [7,] 33.14514 37.20257 36.62187 [8,] 30.75378 30.15063 30.68825 [9,] 36.36103 33.51193 41.95605 [10,] 31.61524 33.14514 36.62187 [11,] 36.80707 31.61524 36.36103 [12,] 31.84423 30.68825 34.12107 [13,] 31.61524 30.68825 37.20257 [14,] 30.75378 36.80707 41.95605 [15,] 33.98218 31.84423 33.14514 [16,] 32.17812 30.07368 36.36103 [17,] 31.98624 33.18473 34.12107 [18,] 33.51193 31.61524 35.29006 [19,] 30.07368 36.80707 33.51193 [20,] 33.64944 34.12107 41.95605 [21,] 33.18473 38.19584 38.19584 [22,] 40.99194 36.62187 33.98218 [23,] 34.61263 40.99194 36.80707 [24,] 40.35282 38.19584 34.96162 [25,] 36.80707 33.14514 31.61524 xboot2_means <- apply(xboot2, MARGIN = 2, FUN = mean) xboot2_sds <- apply(xboot2, MARGIN = 2, FUN = sd) > xboot2_means [1] 33.82794 34.17429 35.50429 > xboot2_sds [1] 3.027593 3.027248 3.215012 In this example, replicate() allows us to write cleaner, more readable code, but with more complex processes, it can also produce large performance gains. In R, for loops are inefficient and should only ever be used as a last resort. Whenever possible, take advantage of R’s native vectorization (it’s behavior of executing a function elementwise on a vector), or use functions like replicate() that are designed to process data efficiently. 4.3 Generating a bootstrap confidence interval Bootstrap resampling is typically used to generate approximate confidence intervals when other methods (e.g., Student’s t distribution or maximum likelihood) won’t work. However, to complete this section, we will calculate the bootstrap interval fo