---
title: "Experimental Design in R"
author: "Joschka Schwarz"
---

```{r}
#| include: false
source(here::here("R/setup-ggplot2-tie.R"))
options(dplyr.summarise.inform = FALSE)
```

Experimental design is a crucial part of data analysis in any field, whether you work in business, health or tech. If you want to use data to answer a question, you need to design an experiment! In this course you will learn about basic experimental design, including block and factorial designs, and commonly used statistical tests, such as the t-tests and ANOVAs. You will use built-in R data and real world datasets including the CDC NHANES survey, SAT Scores from NY Public Schools, and Lending Club Loan Data. Following the course, you will be able to design and analyze your own experiments!

# 1. Introduction to Experimental Design

An introduction to key parts of experimental design plus some power and sample size calculations.

## Intro to experimental design

Theory. Coming soon ...
<p class="dc-cookie-banner-text">                DataCamp and our partners use cookies and similar technologies to improve your learning experience, offer data science content relevant to your interests, improve the site and to show more relevant advertisements. You can change your mind at any time.            



**1. Intro to experimental design**

Hi, my name is Kaelen Medeiros and I’m a data scientist who works in the technology and health industries. I’m here to teach you about experimental design in R. An experiment starts with a question. The experiment involves collecting data with the question in mind and will include analyzing the data to seek an answer. In this course, we’ll focus on asking good questions - in statistical language, formulating clear hypotheses, design the data collection process, and the analysis of collected data.

**2. Steps of an experiment**

The three high-level steps of an experiment are planning, design, and analysis. For planning, you start with your hypothesis -- your question, or even a series of questions. What are you hoping to answer? What is the population of interest, those to whom it applies? What will your dependent variable be, the outcome, which hopefully can be measured to answer the question? What are your independent or explanatory variables, the variables you think may explain the dependent variable? Design questions may naturally follow from planning questions. Choosing a design might entail knowing you want to study different variables and the possible interaction effects of those variables, so you choose a factorial design. Then, if your dependent variable will be a yes/no answer, you know you’re going to be dealing with some kind of logistic regression when you get to analysis.We’re using open data in this course, so we don't know the original experimental design, but that’s okay. The data we'll use throughout the course has been cleaned and altered by me as if it was collected as part of our experiments.

**3. Experimental timeline**

Here's an example timeline for an experiment I designed at a past job. There's no standard for timelines, however!

**4. Key components of an experiment**

The three key components of an experiment include randomization, replication, and blocking. All of these are done to assess variability across our study population, meaning we're looking to explain the variation in an outcome by the different explanatory variables. In order to keep bias low, we need to randomize, create a replicable experiment, and sometimes block. Say we want to test who can score the most free-throws in 5 minutes. We select as one group the high school basketball team, and for the other an English class; this is an incredibly biased experiment right out of the gate. Any conclusions we draw will be based on our poor selection of groups and lack of randomization.

**5. Randomization**

Randomization is a key tenet of any experiment. Randomization helps ensure that variability in outcome due to outside factors that we're not studying in an experiment are evenly distributed among treatment groups. One example of randomization is double-blind medical trials, where neither the patient nor the researcher knows if the patient is receiving treatment or not. The patient is randomized by a third party into one of the two groups.

**6. Replication**

Replication is the idea that we need to repeat our experiment in order to assess what variability looks like. We certainly can't analyze data if we test a drug's efficacy on only one patient; maybe it worked for them, but how do we know if it will also work for other people? Our experiment must be designed to be easily replicable to get a wide range of outcomes.

**7. Blocking**

Blocking is a technique to help control variability. A classic example is to test the effects of a drug on male and female patients, blocking by sex, to account for treatment variability, in this case, known differences in the drug reactions of male and female bodies.

**8. Let's practice!**

Now, let's do some example exercises to demonstrate all three of these key experimental design concepts.

## A basic experiment

Let's dive into experimental design. Note that all of these concepts will be covered in more detail in the next video, "Hypothesis Testing."

`ToothGrowth` is a built-in R dataset from a study that examined the effect of three different doses of Vitamin C on the length of the odontoplasts, the cells responsible for teeth growth in 60 guinea pigs, where tooth length was the measured outcome variable.

Built-in data can be loaded with the `data()` function. The dataset will be loaded as a data frame with the same name passed as an argument to `data()`. For example, you can load the famous `iris` dataset using `data("iris")`.

If you wanted to conduct a two-sided t-test with the famous `mtcars` dataset, it would look like this, where `x` is the outcome in question, `alternative` is set to `"two.sided"`, and `mu` is the value you're testing to see if the mean of `mpg` is *not* equal to.

```{r}
data(mtcars)

t.test(x = mtcars$mpg, alternative = "two.sided", mu = 40)
```

Suppose you know that the average length of a guinea pigs odontoplasts is 18 micrometers. Conduct a two-sided t-test on the ToothGrowth dataset. Here, a two-sided t-test will check to see if the mean of `len` is *not equal* to 18.

**Steps**

1. Load the `ToothGrowth` dataset with the `data()` function (like shown above.) Run this line and then type `ToothGrowth` into the console to be sure it loaded.
2. Use `t.test()` to test if the `len` variable is not equal to 18 micrometers (a two-sided t-test.) There is an example above that you should model your code after.

```{r}
# Load the ToothGrowth dataset
data("ToothGrowth")

# Perform a two-sided t-test
t.test(x = ToothGrowth$len, alternative = "two.sided", mu = 18)
```

Excellent job! Given the high p-value, we fail to reject the null hypothesis that the mean of `len` is equal to 18. That is, we don't have evidence that it is different from 18 micrometers. P-values and hypothesis testing will be covered in more detail in the next video.

## Randomization

Randomization of subjects in an experiment helps spread any variability that exists naturally between subjects evenly across groups. For ToothGrowth, an example of effective randomization would be randomly assigning male and female guinea pigs into different experimental groups, ideally canceling out any existing differences that naturally exist between male and female guinea pigs.

In the experiment that yielded the `ToothGrowth` dataset, guinea pigs were randomized to receive Vitamin C either through orange juice or ascorbic acid, indicated in the dataset by the `supp` variable. It's natural to wonder if there is a difference in tooth length by supplement type - a question that a t-test can also answer!

Starting with this exercise, you should use `t.test()` and other modeling functions with formula notation: 

`t.test(outcome ~ explanatory_variable, data = dataset)`

This can be read: "test `outcome` by `explanatory_variable` in my `dataset`." The default test for `t.test()` is a two-sided t-test.

You no longer have to explicitly declare `dataset$outcome`, because the `data` argument is specified.

**Steps**

1. Conduct the proper test to determine if there is a difference in tooth length based on supplement type, and save the results as an object `ToothGrowth_ttest`.
2. Load the `broom` package.
3. Tidy the `ToothGrowth_ttest` with `tidy()`. This will print the results to the console.

```{r}
# Perform a t-test
ToothGrowth_ttest <- t.test(len ~ supp, data = ToothGrowth)

# Load broom
library(broom)

# Tidy ToothGrowth_ttest
tidy(ToothGrowth_ttest)
```

Nice job! Given the p-value of around 0.06, there seems to be no evidence to support the hypothesis that there's a difference in mean tooth length by supplement type, or, more simply, that there is no difference in mean tooth length by supplement type. Generally in most experiments, any p-value above 0.05 will offer no evidence to support the given hypothesis.

## Replication

Recall that replication means you need to conduct an experiment with an adequate number of subjects to achieve an acceptable statistical power. Sample size and power will be discussed in more detail later in this chapter.

Let's examine the `ToothGrowth` dataset to make sure they followed the principle of replication. We'll use `dplyr` to do some exploratory data analysis (EDA). The data has already been loaded for you.

When using `dplyr` functions, we can utilize the pipe operator, `%>%`, to chain functions together. An example using `mtcars`:

```{r}
# Package
library(dplyr)

data(mtcars)

mtcars %>%
    count(cyl)
```

`count()` groups `mtcars` by `cyl` and then counts how many there are of each number of cylinders.

**Steps**

1. Load the `dplyr` package.
2. Use `count()` to determine how many guinea pigs were given each supplement and dose.

```{r}
# Load dplyr
library(dplyr)

# Count number of observations for each combination of supp and dose
ToothGrowth %>% 
    count(supp, dose)
```

Great job! The researchers seem to have tested each combination of `supp` and `dose` on 10 subjects each, which is low, but was deemed adequate for this experiment.

## Blocking

Though this is not true, suppose the supplement type is actually a nuisance factor we'd like to control for by blocking, and we're actually only interested in the effect of the dose of Vitamin C on guinea pig tooth growth. 

If we block by supplement type, we create groups that are more similar, in that they'll have the same supplement type, allowing us to examine only the effect of dose on tooth length.

We'll use the `aov()` function to examine this. `aov()` creates a linear regression model by calling `lm()` and examining results with `anova()` all in one function call. To use `aov()`, we'll still need functional notation, as with the randomization exercise, but this time the formula should be `len ~ dose + supp` to indicate we've blocked by supplement type. (We'll cover `aov()` and `anova()` in more detail in the next chapter.)

`ggplot2` is loaded for you.

**Steps**

1. Make a boxplot to visually examine if the tooth length is different by `dose`. `dose` has been converted to a factor variable for you.
2. Use `aov()` to detect the effect of `dose` and `supp` on `len`. Save as a model object called `ToothGrowth_aov`.
3. Examine `ToothGrowth_aov` with `summary()` to determine if dose has a significant effect on tooth length.

```{r}
# Create a boxplot with geom_boxplot()
ggplot(ToothGrowth, aes(x = dose, y = len)) + 
    geom_boxplot()

# Create ToothGrowth_aov
ToothGrowth_aov <- aov(len ~ dose + supp, data = ToothGrowth)

# Examine ToothGrowth_aov with summary()
summary(ToothGrowth_aov)
```

Congrats! You have just designed your first Randomized Complete Block Design (RCBD) experiment. We'll learn more about this type of experiment in Chapter 3. Given the very small observed p-value for `dose`, it appears we have evidence to support the hypothesis that mean `len` is different by `dose` amount.

## Hypothesis testing

Theory. Coming soon ...
<p class="dc-cookie-banner-text">                DataCamp and our partners use cookies and similar technologies to improve your learning experience, offer data science content relevant to your interests, improve the site and to show more relevant advertisements. You can change your mind at any time.            



**1. Hypothesis testing**

A very important concept in experimental design is the formation and testing of a hypothesis, or your central research question. For the Tooth Growth dataset we worked with previously, the hypothesis concerned the effect of different doses and administration methods of Vitamin C on the length of tooth growth in the guinea pig. Let's dig in a little more and look at how to create a research hypothesis.

**2. Breaking down hypothesis testing:**

There are really two hypotheses that are grouped together: the null and alternative hypotheses.The null hypothesis is exactly what it sounds like, and the implications change depending on what you're testing. For example, in the tooth growth experiment, the null hypothesis is: "There is no effect of vitamin C dosage or administration type on guinea pig tooth growth."There's some nuance involved in the alternative hypothesis, and its construction will help lead you to the correct test. If you're testing if the mean is only less than or greater than a value (like you did in the first exercise), it's a one-sided test. If you're testing that it's not equal to some number, that's a two sided test. The one/two sided rule applies both to if you're testing one or two groups' means.Recall when we conducted a two sided test to determine if the mean length of tooth growth was not equal to 18. The p-value was 0.4135, so at the 0.05 significance level, we fail to reject the null hypothesis. We have no strong evidence to suggest the mean is not equal to 18.

**3. Power and sample size**

Directly related to hypothesis testing is the idea of power. Power is the probability that the test correctly rejects the null hypothesis when the alternative hypothesis is true. One "golden rule" in statistics is to aim to have 80% power in your experiments, which you'll need an adequate sample size to achieve.Effect size, in the context of power analysis, is a standardized measure of the difference you're trying to detect, calculated as the difference between group means divided by the pooled standard deviation of the data. It's easier to detect a larger difference in means than a smaller one. Sample size is important in experiments. In general, as sample size increases, power increases; you've collected more information, so you're more capable of examining your hypotheses.You need two of these three pieces to calculate the other: if you have a given power and effect size, you can generally calculate sample size. Let's review an example.

**4. Power and sample size calculations**

We're going to calculate power and sample size using the pwr package in this course, but there are built-in power functions you can use if you're more familiar with the concepts.Let's look at calculating power for an ANOVA, or Analysis of Variance test. This is good prep for when we execute these later in the course. The pwr.anova.test() function takes 5 arguments, of which one must be entered as NULL so it can be calculated. k is the number of groups in the comparison, n is the number of observations per group, f is the effect size, then you have to enter a significance level and a power.So to calculate power for a test with three groups, 20 people per group, with an effect size of 0.2 and a significance level of 0.05, you would enter this code. Calculating it returns a power of 0.25--not great! We probably can't detect that small of an effect size with so few people in each group.

**5. Let's practice!**

Let's go to the exercises and review one- and two-sided tests, plus examine the pwr package to calculate power and sample size.

## One sided vs. Two-sided tests

Recall in the first exercise that we tested to see if the mean of the guinea pigs' teeth in `ToothGrowth` was not equal to 18 micrometers. That was an example of a two-sided t-test: we wanted to see if the mean of `len` is some other number on either side of 18.

We can also conduct a one-sided t-test, explicitly checking to see if the mean is less than or greater than 18. Whether to use a one- or two-sided test usually follows from your research question. Does an intervention cause longer tooth growth? One-sided, greater than. Does a drug cause the test group to lose more weight? One-sided, less than. Is there a difference in mean test scores between two groups of students? Two-sided test.

The `ToothGrowth` data has been loaded for you.

**Steps**

1. Test to see if the mean of the length variable of `ToothGrowth` is *less than* 18.


```{r}
# Less than
t.test(x = ToothGrowth$len,
       alternative = "less",
       mu = 18)
```

2. Test to see if the mean of the length variable of `ToothGrowth` is *greater than* 18.


```{r}
# Greater than
t.test(x = ToothGrowth$len,
       alternative = "greater",
       mu = 18)
```

Excellent! It turns out the mean of `len` is actually very close to 18, so neither of these tests tells us much about the mean of tooth length.

## pwr package Help Docs exploration

The `pwr` package has been loaded for you. Use the console to look at the documentation for the `pwr.t.test()` function. The list of arguments are specialized for a t-test and include the ability to specify the alternative hypothesis. 

If you'd like, take some time to explore the different `pwr` package functions and read about their inputs.

What does a call to any `pwr.*()` function yield?

> *Question*
> ---
> ???<br>
> <br>
> ⬜ A vector containing the thing to be calculated.<br>
> ⬜ A data frame of the power, sample size, and other inputs.<br>
> ✅ An object of class "power.htest".<br>
> ⬜ An object of class "integer".<br>

A call to any `pwr.*()` function returns an object of class `power.htest`, which can then be manipulated in the same way as many different R objects.

## Power & Sample Size Calculations

One key part of designing an experiment is knowing the required sample size you'll need to be able to test your hypothesis. 

The `pwr` package provides a handy function, `pwr.t.test()`, which will calculate that for you. However, you do need to know your desired significance level (often 0.05), if the test is one- or two-sided, if the data is from one sample, two samples, or paired, the effect size, and the power. Some of this information will be given to you or can be reasoned from the design. 

A power or sample size calculation is usually different each time you conduct one, and the details of the calculation strongly depend on what kind of experiment you're designing and what your end goals are.

**Steps**

1. Load the `pwr` package.
2. Calculate power using an effect size of 0.35, a sample size of 100 in each group, and a significance level of 0.10.

```{r}
# Load the pwr package
library(pwr) 

# Calculate power
pwr.t.test(n = 100,
           d = 0.35,
           sig.level = 0.1,
           type = "two.sample",
           alternative = "two.sided",
           power = NULL)
```

3. Calculate the sample size needed with an effect size of 0.25, a significance level of 0.05, and a power of 0.8.

```{r}
# Calculate sample size
pwr.t.test(n = NULL,
           d = 0.25, 
           sig.level = 0.05, 
           type = "one.sample", 
           alternative = "greater", 
           power = 0.8)
```

Because this sample size calculation was for a one-sided test, we only need 101 subjects, not 101 in each group. As you design experiments in the future, the `pwr` package includes functions for calculating power and sample size for a variety of different tests, including ANOVA (more on that in the next chapter!)

# 2. Basic Experiments

Explore the Lending Club dataset plus build and validate basic experiments, including an A/B test.

## ANOVA, single and multiple factor experiments

Theory. Coming soon ...
<p class="dc-cookie-banner-text">                DataCamp and our partners use cookies and similar technologies to improve your learning experience, offer data science content relevant to your interests, improve the site and to show more relevant advertisements. You can change your mind at any time.            



**1. ANOVA, single and multiple factor experiments**

Now, let's explore a few key concepts in experimental design: the ANOVA test for 3+ groups, single and multiple factor experiments and the importance of completely randomized design. I'll also introduce the open dataset we'll be using for this chapter, which is data from Lending Club, a loan company.

**2. ANOVA**

So far, we've done some basic comparative experiments where we examined the difference in mean between two groups (or two time periods of the same group) using a t-test. What do we do if we have more than two groups to compare?The answer is an ANOVA test, which stands for Analysis of Variance. It allows us to compare the means of three or more groups, though there's a bit of a catch--we'll only know if at least one of the means is different from the others, but not which mean specifically.The test is still informative and can be implemented in R in a few different ways. In the first, you build a model object with the lm() or glm() functions and then call anova() on that model object. The second implementation, the aov() function, calls lm() internally and both builds the model object for you and conducts the ANOVA test.

**3. Single factor experiments**

A single factor experiment is, like the model example from the last slide, an experiment with one possible explanatory variable. In this example, model 1 is a linear regression model with some outcome variable y and explanatory factor variable x (a single factor).Ideally, a single factor experiment also has a completely randomized design, which means there's some structure in your experiment: if applicable, subjects are randomly assigned to the treatment or control group. A classic, textbook example of a single factor completely randomized design is testing cotton fabric strength. You can test the tensile strength of different cotton fabrics in a random order, and all that differs is the percent cotton in the fabric.

**4. Multiple factor experiments**

A multiple factor experiment expands on the single factor experiment and includes multiple possible explanatory factor variables that may be influencing the outcome variable. This might be an experiment that takes into account not just how much Vitamin C a guinea pig was given, but also the delivery method (a callback to our tooth growth example from chapter 1.)

**5. Intro to Lending Club data**

Across this chapter, we'll use an open dataset from the loan company Lending Club, as downloaded from Kaggle. It's a fairly large dataset with about 890,000 observations and 75 variables, so often in the exercises you'll work with a subset of this dataset. The outcome we'll be most interested in is the amount of loan funded. We'll test different explanatory variables that may influence the amount plus analyze an A/B test using this data.

**6. Let's practice!**

Let's do some exploratory data analysis on the Lending Club data and begin to explore some single and multiple factor experiments.

## Exploratory Data Analysis (EDA) Lending Club

A sample of 1500 observations from the Lending Club dataset has been loaded for you and is called `lendingclub`. Let's do some EDA on the data, in hopes that we'll learn what the dataset contains. We'll use functions from `dplyr` and `ggplot2` to explore the data.

**Steps**

1. Use `glimpse()` from `dplyr` to look at the variable names, types, and some observations.


```{r}
# Package
library(readr)

# Data
lendingclub <- read_csv("data/lendclub.csv")

# Examine the variables with glimpse()
glimpse(lendingclub)
```

2. Pipe the dataset into `summarize()` and find the median loan amount, mean interest rate, and mean annual income. No need to name these outputs anything.


```{r}
# Find median loan_amnt, mean int_rate, and mean annual_inc with summarize()
lendingclub %>% summarize(median(loan_amnt), 
                          mean(int_rate), 
                          mean(annual_inc))
```

3. Use `ggplot2` to build a bar chart of counts of responses to the `purpose` variable. The axes have been flipped for you using `coord_flip()` so the labels are easier to read.


```{r}
# Examine the variables with glimpse()
glimpse(lendingclub)

# Find median loan_amnt, mean int_rate, and mean annual_inc with summarize()
lendingclub %>% summarize(median(loan_amnt), mean(int_rate), mean(annual_inc))

# Use ggplot2 to build a bar chart of purpose
ggplot(data=lendingclub, aes(x = purpose)) + 
  geom_bar() +
  coord_flip()
```

4. Run the code to use `recode()`, creating the `purpose_recode` variable, which takes `purpose` and pares it down to a more manageable number of levels; we're creating 4 (`debt_related`, `big_purchase`, `home_related`, `life_change`.)


```{r}
# Examine the variables with glimpse()
glimpse(lendingclub)

# Find median loan_amnt, mean int_rate, and mean annual_inc with summarize()
lendingclub %>% summarize(median(loan_amnt), mean(int_rate), mean(annual_inc))

# Use ggplot2 to build a bar chart of purpose
ggplot(data=lendingclub, aes(x = purpose)) + 
  geom_bar() +
  coord_flip()

# Use recode() to create the new purpose_recode variable
lendingclub$purpose_recode <- lendingclub$purpose %>% recode( 
        "credit_card" = "debt_related", 
        "debt_consolidation" = "debt_related",
        "medical" = "debt_related",
        "car" = "big_purchase", 
        "major_purchase" = "big_purchase", 
        "vacation" = "big_purchase",
        "moving" = "life_change", 
        "small_business" = "life_change", 
        "wedding" = "life_change",
        "house" = "home_related", 
        "home_improvement" = "home_related")
```

You can see that the original `purpose` variable had quite a few levels, which were very detailed. By using `recode()` here, you created `purpose_recode`, which has a much more manageable 4 general levels that describe the purpose for people's loans.

## How does loan purpose affect amount funded?

In the last exercise, we pared the `purpose` variable down to a more reasonable 4 categories and called it `purpose_recode`. As a data scientist at Lending Club, we might want to design an experiment where we examine how the loan purpose influences the amount funded, which is the money actually issued to the applicant. 

Remember that for an ANOVA test, the null hypothesis will be that all of the mean funded amounts are equal across the levels of `purpose_recode`. The alternative hypothesis is that at least one level of `purpose_recode` has a different mean. We will not be sure which, however, without some post hoc analysis, so it will be helpful to know how ANOVA results get stored as an object in R.

**Steps**

1. Use `lm()` to look at how the `purpose_recode` variable affects `funded_amnt`. Save the model as an object called `purpose_recode_model`.
2. Use`summary()` to examine `purpose_recode_model`. These are the results of the linear regression.
3. Call `anova()` on `purpose_recode_model`. Save as an object called `purpose_recode_anova`. Print it to the console by typing it.  
4. Finally, examine the class of `purpose_recode_anova`.

```{r}
# Build a linear regression model, purpose_recode_model
purpose_recode_model <- lm(funded_amnt ~ purpose_recode, data = lendingclub)

# Examine results of purpose_recode_model
summary(purpose_recode_model)

# Get anova results and save as purpose_recode_anova
purpose_recode_anova <- anova(purpose_recode_model)

# Print purpose_recode_anova
purpose_recode_anova

# Examine class of purpose_recode_anova
class(purpose_recode_anova)
```

Based on the very small p-value, `purpose_recode_anova`'s results indicate that there is evidence to support the hypothesis that the mean loan amounts are different for at least one combination of `purpose_recode`'s levels. You also saw that the ANOVA results are saved as a data frame, which is nice in case you need to access results later. Loans aren't issued in a vacuum, however, and it's likely that more than just purpose influences the amount funded.

## Which loan purpose mean is different?

Before we examine other factors besides `purpose_recode` that might influence the amount of loan funded, let's examine which means of `purpose_recode` are different. This is the post-hoc test referred to in the last exercise. 

The result of that ANOVA test was statistically significant with a very low p-value. This means we can reject the null hypothesis and accept the alternative hypothesis that at least one mean was different. But which one?

We should use Tukey's HSD test, which stands for Honest Significant Difference. To conduct Tukey's HSD test in R, you can use `TukeyHSD()`:

```{r}
#| eval: false
TukeyHSD(aov_model, "outcome_variable_name", conf.level = 0.9)
```

This would conduct Tukey's HSD test on some `aov_model`, looking at a specific `"outcome_variable_name"`, with a `conf.level` of 90%.

**Steps**

1. Build a model using `aov()` that examines `funded_amnt` by `purpose_recode`. Save it as `purpose_aov`.
2. Use `TukeyHSD()` to conduct the Tukey's HSD test on `purpose_aov` with a confidence level of 0.95. Save as an object called `tukey_output`.
3. Tidy `tukey_output` with `tidy()` from the `broom` package (which has been loaded for you.)

```{r}
# Use aov() to build purpose_aov
purpose_aov <- aov(funded_amnt ~ purpose_recode, data = lendingclub)

# Conduct Tukey's HSD test to create tukey_output
tukey_output <- TukeyHSD(purpose_aov, "purpose_recode", conf.level = 0.95)

# Tidy tukey_output to make sense of the results
tidy(tukey_output)
```

Looking at the p-values for each comparison of the levels of `purpose_recode`, we can see that only a few of the mean differences are statistically significant, for example the differences in the means for the `debt_related` and `big_purchase` loan amounts. In this case, these tiny p-values are most likely to be due to large sample size, and further tests would be required to determine what's actually significant in the case of loans (known as the practical significance.)

## Multiple Factor Experiments

We tested whether the purpose of a loan affects loan amount funded and found that it does. However, we also know that it's unlikely that loans are funded based only on their intended purpose. It's more likely that the company is looking at a holistic picture of an applicant before they decide to issue a loan. 

We can examine more than one explanatory factor in a multiple factor experiment. Like our experiments on `ToothGrowth` from Chapter 1, an experimenter can try and control two (or more!) different factors and see how they affect the outcome. We're using open data, so we can't quite control the factors here (they're submitted as someone fills out their loan application), but let's look at how a few other factors affect loan amount funded.

**Steps**

1. Use `aov()` to build a linear model and ANOVA in one step, examining how `purpose_recode` and employment length (`emp_length`) affect the funded amount. Save as an object `purpose_emp_aov` and print the result out.
2. The printed `purpose_emp_aov` does not show p-values, which we might be interested in. Display those by calling `summary()` on the `aov` object.

```{r}
# Use aov() to build purpose_emp_aov
purpose_emp_aov <- aov(funded_amnt ~ purpose_recode + emp_length, data = lendingclub)

# Print purpose_emp_aov to the console
purpose_emp_aov

# Call summary() to see the p-values
summary(purpose_emp_aov)
```

Excellent! You could also perform Tukey's HSD test on this model, but given that `emp_length` has 12 levels, it'll be quite the output. If it was important to the experiment to know the answer, you'd definitely need to look at it.

## Model validation

Theory. Coming soon ...
<p class="dc-cookie-banner-text">                DataCamp and our partners use cookies and similar technologies to improve your learning experience, offer data science content relevant to your interests, improve the site and to show more relevant advertisements. You can change your mind at any time.            



**1. Model validation**

In the last lesson, we built a few linear models with the lm() and aov() functions plus and an anova table with the anova() function. We also did some pre-modeling EDA, though we skipped something pretty crucial that we'll discuss now.

**2. Pre-modeling EDA**

Before modeling you should do some EDA of your data, as in the last lesson. Let's say Lending Club asked you, one of their data scientists, to examine the funded amount of the loan based on verification_status. verification_status is a variable that indicates if the applicant's reported income was somehow verified by Lending Club themselves, verified by another source, or not verified. We looked at the mean and variance in the last lesson, with dplyr code that looked like this:We didn't group by the purpose variable because it hadn't been recoded yet. If we run the second block of code it gives us the median and variance of funded amount.

**3. Pre-modeling EDA continued**

There's more, however! While a boxplot is the kind of graph that non-data scientists don't always respond to (the biggest misconception I hear is people thinking the median is the mean), it's often nice to build one for yourself to see the interquartile range of the variable. We can accomplish this with this ggplot2 code.

**4. Boxplot**

The boxplot here shows no obvious outliers, though there are a few extreme observations on the "not verified" category, represented by the dots in the upper left. We see that "source verified" and "verified" have very similar distributions. This is good news, and you can continue on with modeling.

**5. Post-modeling model validation**

Let's skip ahead a bit and say that you built the ANOVA model for funded amount by verification status and found that the mean funded amount for different verification statuses are significantly different. Furthermore, you did Tukey's HSD test and found that only Verified compared to Source verified is not significantly different from one another. Now what?Now comes post-modeling model validation. This can include looking at different plots, such as a residual or residual versus fitted values plot or a qqplot, testing ANOVA assumptions such as the homogeneity, or the sameness, of variances, or even trying non-parametric alternatives to ANOVA, such as the Kruskal-Wallis test. Non-parametric just means that the test does not assume that the data came from a particular statistical distribution, the way that ANOVA tests assume data is normally distributed.

**6. Post-model validation plots**

The residuals versus fitted plot will show if the model is a good fit if there is a similar scatter pattern for each level of the group variable. We saw that here for this plot. If we saw a different pattern for each level, we could begin to think that there's heteroscedasticity in the residuals, and that the model may not be a great fit. The Normal Q-Q plot should, ideally, show the points scattered around the regression line. One assumption of ANOVA and linear models is that the residuals are normally distributed. If that proves not to be true, your model might not be a good fit, and you may need to try adding explanatory variables or try different modeling techniques.The other two graphs are less commonly discussed, but also have interpretations relevant to your model. A good fit would show in your Scale-Location plot as the residuals increasing with the fitted values, we see that here. The Residuals versus Leverage plot shows which levels are best fitted to the model. Here, the smaller levels seem better fit.

**7. Let's practice!**

Let's dive back in and do the all of the validation steps for the model we already built, now that we know how to use the ANOVA-related functions.

## Pre-modeling EDA

Let's do some EDA with our experiment in mind. Lending Club has now asked you, their data scientist, to examine what effect their Lending Club-assigned loan `grade` variable has on the interest rate, `int_rate`. They're interested to see if the grade they assign the applicant during the process of applying for the loan affects the interest rate ultimately assigned to the applicant during the repayment process. 

The `lendingclub` data has been loaded for you, as has `dplyr` and `ggplot2`.

**Steps**

1. Use `summary()` to look at the `int_rate` variable, and examine its range and interquartile range.

```{r}
# Examine the summary of int_rate
summary(lendingclub$int_rate)
```

2. Using `dplyr`, examine the mean, variance, and median of `int_rate` by `grade`.

```{r}
# Examine the summary of int_rate
summary(lendingclub$int_rate)

# Examine int_rate by grade
lendingclub %>% 
  group_by(grade) %>% 
  summarize(mean = mean(int_rate), var = var(int_rate), median = median(int_rate))
```

3. Use `ggplot2` to create a boxplot of `int_rate` by `grade`.

```{r}
# Examine the summary of int_rate
summary(lendingclub$int_rate)

# Examine int_rate by grade
lendingclub %>% 
  group_by(grade) %>% 
  summarize(mean = mean(int_rate), var = var(int_rate), median = median(int_rate))

# Make a boxplot of int_rate by grade
ggplot(lendingclub, aes(x = grade, y = int_rate)) + 
  geom_boxplot()
```

4. Save a linear model examining this experiment in an object called `grade_aov`. Print the results by calling `summary()`.

```{r}
# Examine the summary of int_rate
summary(lendingclub$int_rate)

# Examine int_rate by grade
lendingclub %>% 
  group_by(grade) %>% 
  summarize(mean = mean(int_rate), var = var(int_rate), median = median(int_rate))

# Make a boxplot of int_rate by grade
ggplot(lendingclub, aes(x = grade, y = int_rate)) + 
  geom_boxplot()

# Use aov() to create grade_aov plus call summary() to print results
grade_aov <- aov(int_rate ~ grade, data = lendingclub)
summary(grade_aov)
```

Excellent job! You can see from the numeric summary and the boxplot that grade seems to heavily influence interest rate. Therefore, the linear model results indicating that `int_rate` is significantly different by `grade` are unsurprising.

## Post-modeling validation plots + variance

In the last exercise, we found that `int_rate` does differ by `grade`. Now we should validate this model, which for linear regression means examining the Residuals vs. Fitted and Normal Q-Q plots.

If you call `plot()` on a model object in R, it will automatically plot both of those plots plus two more. You'll interpret these plots to evaluate model fit. We discussed how to do this in the video.

Another assumption of ANOVA and linear modeling is homogeneity of variance. Homogeneity means "same", and here that would mean that the variance of `int_rate` is the same for each level of `grade`. We can test for homogeneity of variances using `bartlett.test()`, which takes a formula and a dataset as inputs.

**Steps**

1. Run the first line of code with `par()` so the plots will output in a 2 by 2 grid.
2. Call `plot()` on `grade_aov` (which has been created for you) to produce the model diagnostic plots.
3. Test for homogeneity of variances using `bartlett.test()`.

```{r}
# For a 2x2 grid of plots:
par(mfrow = c(2, 2))

# Plot grade_aov
plot(grade_aov)

# Bartlett's test for homogeneity of variance
bartlett.test(int_rate ~ grade, data = lendingclub)
```

Excellent! The residuals on this model are okay, though the residuals on G have a much smaller range than any other level of `grade` (the dots are far less spread out.) The Q-Q plot, however, shows that the residuals are fairly normal. However, given the highly significant p-value from Bartlett's test, the assumption of homogeneity of variances is violated, which is one of the assumptions of an ANOVA model. Therefore, ANOVA might not be the best choice for this experiment. It happens!

## Kruskal-Wallis rank sum test

Given that we found in the last exercise that the homogeneity of variance assumption of linear modeling was violated, we may want to try an alternative. 

One non-parametric alternative to ANOVA is the Kruskal-Wallis rank sum test. For those with some statistics knowledge, it is an extension of the Mann-Whitney U test for when there are more than two groups, like with our `grade` variable. For us, the null hypothesis for this test would be that all of the `int_rate`s have the same ranking by `grade`.

The Kruskal-Wallis rank sum test can be conducted using the `kruskal.test()` function, available in base R. Luckily for you, the use of this function is very similar to using `lm()` or `aov()`: you input a formula and a dataset, and a result is returned.

**Steps**

1. Use `kruskal.test()` to examine whether `int_rate` varies by `grade` when a non-parametric model is employed.

```{r}
# Conduct the Kruskal-Wallis rank sum test
kruskal.test(int_rate ~ grade,
             data = lendingclub)
```

Good job! The low p-value indicates that based on this test, we can be confident in our result, which we found across this experiment, that `int_rate` varies by `grade`.

## A/B testing

Theory. Coming soon ...
<p class="dc-cookie-banner-text">                DataCamp and our partners use cookies and similar technologies to improve your learning experience, offer data science content relevant to your interests, improve the site and to show more relevant advertisements. You can change your mind at any time.            



**1. A/B testing**

A/B tests are an incredibly common type of experiment that a data scientist or data analyst might design. Let's talk about what an A/B test is.

**2. A/B testing**

A/B tests are often discussed in a marketing context and are very common in marketing-related companies, although they're definitely used in other industries as well. One common example of their use is to test customer engagement with different features of a company's website.

**3. Power and sample size in A/B tests**

One important concept that will come back here is power and sample size calculations, as they're pretty crucial in A/B testing. Usually, you'll be calculating sample size, given some power and significance level. Then, you let your A/B test run until you get the sample size you need. You also need to decide an effect size, just as before! You may have guessed by now that an A/B test is simply an application of your basic experimental design knowledge, though they can and do get complicated, so stick with me so you can get good at them.An A/B test changes one thing and one thing only and measures the difference in outcome between these two alternatives.

**4. Lending Club A/B test**

In this section, we'll design and conduct an A/B test with the Lending Club data we've been using. Say Lending Club was interested to see how the color of the website header affected the loan amount, which is how much an applicant asks to borrow. They have a general hypothesis that softer, gentle colors may influence applicants to ask for lower amount of money, perhaps a more reasonable amount that they can more feasibly pay back. They already use a light blue website header, but they've decided to test a second softer, gentler color on the website.Applicants designated Group A were funneled to an application with the light blue existing header while applicants designated group B were funneled to an application with a new mint green header. We'll examine the amount of money the applicant asked for when applying for the loan based on which colored header they were shown -- not a traditional A/B test, which often focuses on a metric such as click-through rate, which in this case would be number of people who applied, but still an interesting test.

**5. Let's practice!**

Let's take everything we've learned in the last two chapters and analyze the Lending Club A/B test.

## Which post-A/B test test?

We'll be testing the mean `loan_amnt`, which is the requested amount of money the loan applicants ask for, based on which color header (green or blue) that they saw on the Lending Club website. 

Which statistical test should we use after we've collected enough data?

> *Question*
> ---
> ???<br>
> <br>
> ⬜ Kruskal-Wallis Rank Sum test<br>
> ✅ T-test<br>
> ⬜ Chi-Square Test<br>
> ⬜ Linear Regression<br>

Excellent! We'll be looking at the mean difference in the amount of loan the applicants asked for, so a t-test is appropriate.

## Sample size for A/B test

We know now that we need to analyze our A/B test results with a t-test after we've collected data. We have two pretty important questions we need to answer before we do that: what's the effect size and what's the sample size required for this test?

In this case, effect size was given to us. Lending Club is looking to detect the relatively small effect size of 0.2. We'll again use the `pwr` package and calculate sample size using an appropriate function to find out how many we'll need to recruit into each group, A and B.

**Steps**

1. Use the correct function from the `pwr` package to calculate the required sample size for each group with `d = 0.2`, a `power` of 0.8, and a 0.05 significance level. Check the `pwr` help docs with `?pwr` if you need help remembering which function to use and what arguments it takes.

```{r}
# Load the pwr package
library(pwr)

# Use the correct function from pwr to find the sample size
pwr.t.test(n = NULL, 
           d = 0.2,
           sig.level = 0.05,
           power = 0.8,
           alternative = "two.sided")
```

Nice! You can see we need about 400 people per group to reach our desired power in this A/B test.

## Basic A/B test

Now that we know the sample size required, and we allowed the experiment to run long enough to get at least 400 people in each group, we can analyze our A/B test. 

Remember that when applicants were using the Lending Club website, they were randomly assigned to two groups, A or B, where A was shown a mint green website header and B was shown a light blue website header. Lending Club was interested to see if website header color choice influenced `loan_amnt`, the amount an applicant asked to borrow.

<img src="http://assets.datacamp.com/production/repositories/1793/datasets/0cc81afa09f43147bf70c484c5f437db215160f6/lendclub2.PNG" alt>

A new dataset, `lendingclub_ab` is available in your workspace. The A/B test was run until there were 500 applicants in each group. Each applicant has been labeled as group A or B. Conduct the proper test to see if the mean of `loan_amnt` is different between the two groups.

**Steps**

1. Create a boxplot of `loan_amnt` by `Group` using `ggplot2`.
2. Conduct the two-sided t-test to test the A/B test results.

```{r}
# Data
lendingclub_ab <- read_rds("data/lendclub_ab.rds")

# Plot the A/B test results
ggplot(lendingclub_ab, aes(x = Group, y = loan_amnt)) + 
  geom_boxplot()

# Conduct a two-sided t-test
t.test(loan_amnt ~ Group, data = lendingclub_ab)
```

Excellent! By looking at both the boxplot and the results of the t-test, it seems that there is no compelling evidence to support the hypothesis that there is a difference the two A/B test groups' mean `loan_amnt`, a result which you would use to help make data-driven decisions at Lending Club.

## A/B tests vs. multivariable experiments

The point of an A/B test is that only one thing is changed and the effect of that change is measured. We saw this with our examples in the video and the last few exercises. On the other hand, a multivariate experiment, such as the `ToothGrowth` experiment from chapter 1, is where a few things are changed (and is similar to a multiple factor experiment, which we covered earlier in this chapter.)

A Lending Club multivariate test can combine all of the explanatory variables we've looked at in this chapter. Let's examine how `Group`, `grade`, and `verification_status` affect `loan_amnt` in the `lendingclub_ab` dataset.

**Steps**

1. Use `lm()` to examine the effect of all three explanatory variables on `loan_amnt`. Save as a model object called `lendingclub_multi`.
2. Examine `lendingclub_multi` with `tidy()` and draw your conclusions.

```{r}
# Build lendingclub_multi
lendingclub_multi <- lm(loan_amnt ~ Group + grade + verification_status, data = lendingclub_ab)

# Examine lendingclub_multi results
tidy(lendingclub_multi)
```

From the results, verification status and having an F `grade` are the factors in this model that have a significant effect on loan amount. Let's move on to the next chapter and conduct more multivariable experiments like this.

# 3. Randomized Complete (& Balanced Incomplete) Block Designs

Use the NHANES data to build a RCBD and BIBD experiment, including model validation and design tips to make sure the BIBD is valid.

## Intro to NHANES and sampling

Theory. Coming soon ...
<p class="dc-cookie-banner-text">                DataCamp and our partners use cookies and similar technologies to improve your learning experience, offer data science content relevant to your interests, improve the site and to show more relevant advertisements. You can change your mind at any time.            



**1. Intro to NHANES and sampling**

In this chapter, we're going to learn how to design and analyze Randomized Complete and Randomized Incomplete Block Design experiments. Before we can do that, however, we'll explore the NHANES dataset and talk about sampling. Let's jump in.

**2. Intro to NHANES dataset**

NHANES, or the National Health and Nutrition Examination Survey, is conducted every two years in the United States to collect information that assists in determination of major diseases and risk factors that are prevalent in the U.S. NHANES combines interviews and physical exams in order to collect demographic, medical, dental, dietary, and general health-related information about participants. A sample is constructed so that the information is collected on individuals who resemble the United States' population. NHANES has existed in some form since the 1960's and since 1999 has been conducted every two years. Accompanying data, codebooks, and analytic guidelines are released by the CDC. The information collected benefits not only government agencies such as the CDC and the Food and Drug Administration, but also researchers interested in risk factors and health conditions in the U.S.

**3. Intro to sampling**

The participants in NHANES aren't chosen randomly; they're sampled according to a scheme designed to represent the U.S. population. People older than 60, African Americans, and Hispanic individuals are oversampled in the collection. Caucasians made up 60% of the U.S. population at the last census, and if they did not over-sample these other groups, we wouldn't end up with enough African American and Hispanic individuals to produce reliable statistics on their health.In this course we're only going to discuss probability sampling methods, or those where some degree of probability is used to select the sampled units. Non-probability sampling includes things like voluntary responses and convenience sampling, which I encourage you to look into yourself. Let's discuss five basic types of sampling.

**4. Sampling**

A simple random sample is the most straightforward: every unit in a population has an equal probability of being sampled. In R, you can implement it with the sample() function. For more information on how to use this function, be sure to read its documentation.Stratified sampling involves splitting your population by some strata variable, such as race, gender, or type, and then taking a simple random sample inside of each stratified group. This is distinct from cluster sampling, where you divide the population into groups called clusters, perhaps every high school in a state, randomly select some number of those clusters, and sample everyone inside the selected clusters. In R, stratified sampling can be carried out simply using group_by() and sample_n() from dplyr. A good way to conduct cluster sampling is the cluster() function in the sampling package.Systematic sampling involves choosing a sample in a systematic way, such as every 5th or 10th or 12.7th (just kidding) unit of the population. It's best implemented in R with a custom function.Multi-stage sampling simply combines one or more of the aforementioned approaches in a logical way.

**5. Let's practice!**

Let's jump in and do some EDA and data cleaning on the NHANES dataset so we can get familiar with using it. We'll also do some resampling of the data to reinforce those concepts.

## NHANES dataset construction

As downloaded from the <a href="https://wwwn.cdc.gov/nchs/nhanes/continuousnhanes/default.aspx?BeginYear=2015">NHANES website</a>, the NHANES datasets are available only as separate .XPT files, a native format to SAS. Luckily for us, the `haven` package exists.

Let's combine the NHANES Demographics, Medical Conditions, and Body Measures datasets, available in their raw .XPT format and accessible through the variables `DEMO_file`, `MCQ_file`, and `BMX_file`. Join all 3 datasets using the `SEQN` variable. A good way to do this is using `Reduce()`, which allows you to combine elements in a helpful way.

The joining code, which is provided for you does the following:


* Creates a list of all 3 datasets (`nhanes_demo`, `nhanes_medical`, `nhanes_bodymeasures`).
* Uses a custom function inside of `Reduce()` to inner join all 3 datasets with the `"SEQN"` variable.
* Saves this as the `nhanes_combined` dataset.
**Steps**

1. Load the `haven` package.
2. Import the three data files with separate calls to `read_xpt()`, where the inputs to these 3 calls to `read_xpt()` are `DEMO_file`, `MCQ_file`, and `BMX_file` and saved as the datasets as `nhanes_demo`, `nhanes_medical`, and `nhanes_bodymeasures`, respectively.
3. Create `nhanes_combined` by merging the 3 datasets you just imported, using the provided code.

```{r}
# Load haven
library(haven)

# Import the three datasets using read_xpt()
nhanes_demo <- read_xpt("data/DEMO_I.XPT")
nhanes_medical <- read_xpt("data/MCQ_I.XPT")
nhanes_bodymeasures <- read_xpt("data/BMX_I.XPT")

# Merge the 3 datasets you just created to create nhanes_combined
nhanes_combined <- list(nhanes_demo, nhanes_medical, nhanes_bodymeasures) %>%
  Reduce(function(df1, df2) inner_join(df1, df2, by = "SEQN"), .)
```

Awesome! Now that we have the NHANES data assembled, let's get to work on EDA & cleaning.

## NHANES EDA

Let's examine our newly constructed dataset with a mind toward EDA. As in the last chapter, it's a good idea to look at both numerical summary measures and visualizations. These help with understanding data and are a good way to find data cleaning steps you may have missed. The `nhanes_combined` dataset has been pre-loaded for you.

Say we have access to NHANES patients and want to conduct a study on the effect of being told by a physician to reduce calories/fat in their diet on weight. This is our treatment; we're pretending that instead of this being a question asked of the patient, we randomly had physicians counsel some patients on their nutrition. However, we suspect that there may be a difference in weight based on the gender of the patient - a blocking factor!

**Steps**

1. Fill in and run the `dplyr` code to find mean weight (`BMXWT`) in kg by our treatment (`MCQ365D`). Is there anything interesting about the `NA` treated patients?
2. Fill in the `ggplot2` code to look at a boxplot of the IQR of patients' weights by the treatment variable.

```{r}
# Fill in the dplyr code
nhanes_combined %>% 
  group_by(MCQ365D) %>% 
  summarize(mean = mean(BMXWT, na.rm = TRUE))

# Fill in the ggplot2 code
nhanes_combined %>% 
  ggplot(aes(as.factor(MCQ365D), BMXWT)) +
  geom_boxplot() +
  labs(x = "Treatment",
       y = "Weight")
```

Great! Now that we have an idea of some of the problems with the data, namely that children weren't given the treatment - that's why we see an `NA` age category. We also have some patients have weights missing, thus the warning that the boxplot throws. It's time for some data cleaning.

## NHANES Data Cleaning

During data cleaning, we discovered that no one under the age of 16 was given the treatment. Recall that we're pretending that the variable that indicates if a doctor has ever advised them to reduce fat or calories in their diet is purposeful nutrition counseling, our treatment. Let's only keep patients who are greater than 16 years old in the dataset.

You also may have noticed that the default settings in `ggplot2` delete any observations with a missing dependent variable, in this case, body weight. One option for dealing with the missing weights, imputation, can be implemented using the `simputation` package. Imputation is a technique for dealing with missing values where you replace them either with a summary statistic, like mean or median, or use a model to predict a value to use.

We'll use `impute_median()`, which takes a dataset and the variable to impute or formula to impute by as arguments. For example, `impute_median(ToothGrowth, len ~ dose)` would fill in any missing values in the variable `len` with the median value for `len` by `dose`. So, if a guinea pig who received a dose of 2.0 had a missing value for the `len` variable, it would be filled in with the median `len` for those guinea pigs with a `dose` of 2.0.

**Steps**

1. Create `nhanes_filter` by using `filter()` to keep anyone older than 16 in the dataset, not including those who are 16. Age is stored in the `RIDAGEYR` variable.
2. Load `simputation`. Use `impute_median()` to fill in the missing observations of `bmxwt` in `nhanes_filter`, grouping by `RIAGENDR`.
3. Recode the `nhanes_final$MCQ365D` variable by setting any observations with a value of 9 to 2 instead. Verify the recoding worked with `count()`.

```{r}
# Filter to keep only those 16+
nhanes_filter <- nhanes_combined %>% filter(RIDAGEYR > 16)

# Load simputation & impute BMXWT by RIAGENDR
library(simputation)
nhanes_final <- impute_median(nhanes_filter, BMXWT ~ RIAGENDR)

# Recode MCQ365D with recode() & examine with count()
nhanes_final$MCQ365D <- recode(nhanes_final$MCQ365D, 
                               `1` = 1,
                               `2` = 2,
                               `9` = 2)
nhanes_final %>% count(MCQ365D)
```

Excellent! Imputation is a powerful tool for dealing with missing data, but should be used with caution, as you can introduce bias into your data if you're not careful how you impute. Now that we have the dataset cleaned, we're ready to learn about RCBDs so we can analyze our experiment.

## Resampling NHANES data

The NHANES data is collected on sampled units (people) specifically selected to represent the U.S. population. However, let's resample the `nhanes_final` dataset in different ways so we get a feel for the different sampling methods.

We can conduct a simple random sample using `sample_n()` from `dplyr`. It takes as input a dataset and an integer of number of rows to sample.

Stratified sampling can be done by combining `group_by()` and `sample_n()`. The function will sample `n` from each of the groups specified in the `group_by()`. 

The `sampling` package's `cluster()` creates cluster samples. It takes in a dataset name, the variable in the set to be used as the cluster variable, passed as a vector with the name as a string (e.g. `c("variable")`), a number of clusters to select, and a method.

**Steps**

1. Use `sample_n()` to select 2500 observations from `nhanes_final` and save as `nhanes_srs`.
2. Create `nhanes_stratified` by using `group_by()` and `sample_n()`. Stratify by `RIAGENDR` and select 2000 of each gender. Confirm that it worked by using `count()` to examine `nhanes_stratified`'s gender variable.
3. Load the `sampling` package. Use `cluster()` to divide `nhanes_final` by `"indhhin2"` into 6 clusters using the `"srswor"` method. Assign to `nhanes_cluster`.

```{r}
# Use sample_n() to create nhanes_srs
nhanes_srs <- nhanes_final %>% sample_n(2500)

# Create nhanes_stratified with group_by() and sample_n()
nhanes_stratified <- nhanes_final %>% group_by(RIAGENDR) %>% sample_n(2000)
nhanes_stratified %>% 
  count(RIAGENDR)

# Load sampling package and create nhanes_cluster with cluster()
library(sampling)
nhanes_cluster <- cluster(nhanes_final, "INDHHIN2", 6, method = "srswor")
```

Excellent! These are some basic sampling methods you can use on your data to create the different kinds of samples that may be necessary in an experiment.

## Randomized Complete Block Designs (RCBD)

Theory. Coming soon ...
<p class="dc-cookie-banner-text">                DataCamp and our partners use cookies and similar technologies to improve your learning experience, offer data science content relevant to your interests, improve the site and to show more relevant advertisements. You can change your mind at any time.            



**1. Randomized Complete Block Designs (RCBD)**

The next design we'll discuss is a Randomized Complete Block Design, often abbreviated to RCBD.

**2. RCBDs**

In the last chapter, we conducted completely randomized experiments. One extension of that type experiment is the RCBD, as described here. We utilize blocking in our experiments when we have a "nuisance" factor, something that may affect the outcome of the experiment, but isn't actually of interest to study in the experiment. Blocking allows for the experimental groups to be more similar. We expect groups like those divided by metal hardness, type of crop, or gender, to have minimal differences inside the group, while across the different groups, we expect there to be much larger differences. Sometimes, the blocking factor can be thought of not as a nuisance factor, but as something you can control for as an experimenter. There's a saying in experimental design - "block what you can, randomize what you can't."Essentially, you create blocks with your nuisance variable or variable you can control, and then run your completely randomized experiment inside each block for your desired number of repetitions.

**3. RCBD workflow**

In most texts or Experimental Design courses, the RCBD is demonstrated with an agricultural example. Say, for instance, you're testing a new fertilizer on your fields. Your outcome will be plant growth.One factor that isn't of terrible interest in this experiment is the type of crop. So you block the fields: one block is corn and one soybeans. Then, inside of each block, you conduct a randomized control trial, allocating the fields randomly to either receive the new fertilizer or to use the standard method of growing the crop. Then, you measure the outcome at the specified time points and compare between the groups. Voila, a RCBD!

**4. agricolae**

While we're on the topic of agriculture, I'd like to introduce the agricolae R package. This package enables you to "draw" some of the different experimental designs possible, including an RCBD. Let's walk through an example of using the package to draw a RCBD. We can create a treatment variable of four treatments, trt, using the R built-in lowercase letters. Then, we create an object rep, and assign it our four planned repetitions. We create our object design dot rcdb with the design dot rcdb function, which takes in our treatment, our rep, a seed, which tells the random number generator where to begin and allows for reproducibility - I'm fond of seed = 42 - and set the serie input, which has to do with how the number blocks are tagged in larger experiments. If we examine the sketch object, we can see how our four repetitions with our four treatments can be randomly repeated in each run.

**5. Let's practice!**

Now that we know the basics of RCDBs and how to draw them with agricolae, let's jump in and try some examples.

## Which is NOT a good blocking factor?

As discussed in the video, the purpose of blocking an experiment is to make the experimental groups more like one another. Groups are blocked by a variable that is known to introduce variability that will affect the outcome of the experiment but is not of interest to study in the experiment itself. 

A rule of thumb in experimental design is often "block what you can, randomize what you cannot", which means you should aim to block the effects you can control for (e.g. sex) and randomize on those you cannot (e.g. smoking status). Variability inside a block is expected to be fairly small, but variability between blocks will be larger.

Which of the following would **NOT** make a good blocking factor for an experiment?

> *Question*
> ---
> ???<br>
> <br>
> ⬜ When studying the effect of radiation treatment on cancer recurrence, treatment hospital type is used a blocking factor.<br>
> ⬜ When studying the effect of four different tips' hardness readings on a metal hardness tester, the known hardness of the metal is used as a block variable.<br>
> ⬜ When testing the effect of a type of fertilizer on plant growth, crop type is used as a blocking factor.<br>
> ✅ When testing the effect of a drug on blood pressure, current pregnancy is used as a blocking factor.<br>

Excellent! Current pregnancy is not really something the experimenter can control, and they're also not likely to have enough subjects pregnant to justify blocking by that factor (unless they sampled for pregnant individuals in the first place!)

## Drawing RCBDs with Agricolae

The `agricolae` package is very helpful when you want to "draw" out the design of an experiment for yourself using R. It can draw many different kinds of experiments, including a randomized complete block design. Here's an example of one:

```{r}
#| eval: false
     [,1] [,2] [,3] [,4]
[1,] "D"  "C"  "A"  "B" 
[2,] "B"  "A"  "D"  "C" 
[3,] "D"  "A"  "B"  "C" 
[4,] "A"  "B"  "D"  "C"
```

In this RCBD, we have 4 blocks (each row of the output). Inside of each block, each treatment "A", "B", "C", and "D" is used, because this is a *complete* design. So if these 4 blocks/rows of the output were four fields of a farmer's, they should give the first field the "D" treatment in the first season, then "C", then "A", then "B".

Let's draw an RCBD design with 5 treatments and 4 blocks, which go in the `r` argument. The `agricolae` package has been loaded for you.

**Steps**

1. Create the object called `designs` using the given code and print it to see all possible designs that `agricolae` can draw.

```{r}
# Package
library(agricolae)

# Create designs using ls()
designs <- ls("package:agricolae", pattern = "design")
print(designs)
```

2. Use `str()` to examine the criteria for `design.rcbd`.

```{r}
# Create designs using ls()
designs <- ls("package:agricolae", pattern = "design")
print(designs)

# Use str() to view design.rcbd's criteria
str(design.rcbd)
```

3. Build the `treats` and `rep` objects. `treats` should be a vector containing the letters A through E, created using `LETTERS[1:5]`. `blocks` should be equal to 4.

```{r}
# Create designs using ls()
designs <- ls("package:agricolae", pattern = "design")
print(designs)

# Use str() to view design.rcbd's criteria
str(design.rcbd)

# Build treats and rep
treats <- LETTERS[1:5]
blocks <- 4
```

4. Create the `my_design_rcbd` object. The `seed` has been set for you, for reproducibility. View the `sketch` part of the object.

```{r}
# Create designs using ls()
designs <- ls("package:agricolae", pattern = "design")
print(designs)

# Use str() to view design.rcbd's criteria
str(design.rcbd)

# Build treats and rep
treats <- LETTERS[1:5]
blocks <- 4

# Build my_design_rcbd and view the sketch
my_design_rcbd <- design.rcbd(treats, r = blocks, seed = 42)
my_design_rcbd$sketch
```

Nice! Now that you have a better idea of what a RCBD looks like, let's try a few examples, including one with the NHANES data we cleaned.

## NHANES RCBD

Recall that our blocked experiment involved a treatment wherein the doctor asks the patient to reduce their fat or calories in their diet, and we're testing the effect this has on weight (`bmxwt`). We plan to block by gender, which in the NHANES dataset is stored as `riagendr`. Recall that blocking is done to create experimental groups that are as similar as possible. Blocking this experiment by gender means that if we observe an effect of the treatment on `bmxwt,` it's more likely that the effect was actually due to the treatment versus the individual's gender.

In your R code, you denote a blocked experiment by using a formula that looks like: `outcome ~ treatment + blocking_factor` in the appropriate modeling function.

`nhanes_final` is available.

**Steps**

1. Use `aov()` to create `nhanes_rcbd`. Recall that the treatment is stored in `mcq365d` and you're testing the outcome `bmxwt`, with the blocking factor `riagendr`.
2. Examine the results of `nhanes_rcbd` with `summary()`.
3. Use `dplyr` functions to examine the mean weights by `mcq365d` and `riagendr`.

```{r}
# Use aov() to create nhanes_rcbd
nhanes_rcbd <- aov(BMXWT ~ MCQ365D + RIAGENDR, data = nhanes_final)

# Check results of nhanes_rcbd with summary()
summary(nhanes_rcbd)

# Print mean weights by mcq365d and riagendr
nhanes_final %>% 
  group_by(MCQ365D, RIAGENDR) %>% 
  summarize(mean_wt = mean(BMXWT, na.rm = TRUE))
```

Nice! It's pretty clear that there truly is a mean difference in weight by gender, so blocking was a good call for this experiment. We also observed a statistically significant effect of the treatment on `bmxwt`, which we hope is actually a result of the treatment. Now that we have the RCBD down, let's tackle Balanced Incomplete Block Designs (BIBD).

## RCBD Model Validation

As we did in the last chapter (and when building any model!) it's a good idea to validate the results. We'll examine the Residuals vs. Fitted and Normal Q-Q plots, though now we'll only see a Constant Leverage plot in place of the other two. A good model has a Q-Q plot showing an approximately normal distribution and no clear patterns across blocks or treatments in the others. 

We can also look at Interaction plots. We hope to see parallel lines, no matter which of the block or the treatment is on the x-axis. If they are, they satisfy a key assumption of the RCBD model called Additivity.

The `nhanes_rcbd` model object from the last exercise has been loaded for you. Examine the results with `summary(nhanes_rcbd)` in the console if you need a refresher.

**Steps**

1. Plot the `nhanes_rcbd` model object, being sure to set up a 2x2 grid of plots beforehand.

```{r}
# Set up the 2x2 plotting grid and plot nhanes_rcbd
par(mfrow = c(2, 2))
plot(nhanes_rcbd)
```

2. Run the code to view the interaction plot between the treatment and gender and observe if the lines are parallel.

```{r}
# Run the code to view the interaction plots
with(nhanes_final, interaction.plot(MCQ365D, RIAGENDR, BMXWT))
```

3. Run the code to view the interaction plot between gender and the treatment (it'll be a little different!) and observe if the lines are parallel.

```{r}
# Run the code to view the interaction plots
with(nhanes_final, interaction.plot(RIAGENDR, MCQ365D, BMXWT))
```

Excellent! The initial diganostic plots show that this model is pretty good but not great - especially at the larger end of the data, the Q-Q plot shows the data might not be normal. The interaction plots show nearly parallel lines, so we can move forward with this model.

## Balanced Incomplete Block Designs (BIBD)

Theory. Coming soon ...
<p class="dc-cookie-banner-text">                DataCamp and our partners use cookies and similar technologies to improve your learning experience, offer data science content relevant to your interests, improve the site and to show more relevant advertisements. You can change your mind at any time.            



**1. Balanced Incomplete Block Designs (BIBD)**

We covered randomized complete block designs, which was a completely randomized experiment executed inside each block in order to control for a "nuisance" factor that you're not interested in studying. What about when you can't fit each treatment in a block? That's an incomplete block design.

**2. Balanced Incomplete Block Designs**

Sometimes it will not be feasible, advisable, or even necessary to test every treatment in each block. In that case, you can possibly build a Balanced Incomplete Block Design. Incomplete block designs do exist and lack the "balance" in each pair of treatments occuring an equal number of times across the experiment. We won't cover these here, but you should know that they're possible.

**3. Is there a BIBD?**

A BIBD isn't always possible, however. You could wing it and plug numbers in and out of the design_bib() function in agricolae, which we'll demonstrate in the exercises, or you can do some quick math to find out if a BIBD is possible.Let t be our number of treatments, k be the number of factor levels of our blocking variable, and r the desired replications of our experiment. To find out if a BIBD is even possible, you can calculate lambda using this formula: replications times number of blocking factors minus one, divided by number of treatments minus 1. If lambda is a whole number, there is a BIBD possible, and if not, there simply isn't one for the given t, k, and r.

**4. BIBD examples (1)**

Say we only have the money for two treatments per block (k) out of four total treatments (t), that is, four blocks and two full replications of the experiment. Even though this may look fully balanced, when you calculate lambda, you quickly see that it equals 2/3, and not a whole number.

**5. Invalid BIBD**

A BIBD isn't possible here; you would need to add more replications of the experiment or a larger block size. In the image version, we can check to see if each pair of treatment occurs in a balanced way. Right away, we can see that AB occurs only once in the first column. The pair AC never occurs in a block. Either one of these methods will tell you this is not a valid BIBD, although calculating the lambda is certainly quicker.

**6. BIBD examples (2)**

In the case of a second example where there are four treatments (t), and you intend to do three treatments per block (k), with three full replications (r), we get a lambda of 2. This BIBD is a go!

**7. Valid BIBD**

Again, we can look at the image to verify that each pair of treatments occurs together twice across the entire experiment. Let's find the number of occurrences of AB, for example. Columns are our blocking factor, so we're checking for pairs column-wise. We see AB in the first and second columns only. AC, for example, occurs in the second and third columns only. If we continued, this would hold true for all pairs of treatments in this study.

**8. Let's practice!**

Don't be too intimidated by BIBDs. Let's jump into some examples of calculating lambda, finding BIBDs using agricolae, and even analyzing a few, including an example with the NHANES data.

## Is a BIBD even possible?

We saw in the video that it's possible a BIBD doesn't exist at all. It's useful to calculate the `lambda` beforehand, and if the result isn't an integer, there isn't a possible BIBD.

Don't bust out the pen and paper yet: the `lambda()` function has been defined for you. It takes as input `t` = number of treatments, `k` = number of treatments per block, and `r` = number of repetitions. It then calculates `lambda` according to this formula:
`r(k - 1) / (t - 1)`.

Try the different combinations of `t`, `k`, and `r` using `lambda()` in your console, and choose the answer which does NOT have a BIBD.

> *Question*
> ---
> ???<br>
> <br>
> ⬜ `t` = 2, `k` = 2, `r` = 2<br>
> ⬜ `t` = 2, `k` = 3, `r` = 4<br>
> ✅ `t` = 3, `k` = 4, `r` = 11<br>
> ⬜ `t` = 12, `k` = 2, `r` = 22<br>

Good job! When designing experiments, you could build your own version of the custom function lambda, or use agricolae to figure it out visually, as we'll see in the next exercise.

## Drawing BIBDs with agricolae

We can also use `agricolae` to draw BIBDs. `design.bib()` takes, at minimum, the treatments (`treats`), an integer `k` corresponding to the number of levels of the blocks, and a `seed` as inputs.

The main thing you should notice about a BIBD is that not every treatment will be used in each block (column) of the output.

From the video and the last exercise, however, you know that sometimes a BIBD isn't valid and that you have to do a little math to be sure your BIBD design is possible. `design.bib()` will return an error message letting you know if a design is not valid.

Let's draw a few BIBDs with `agricolae` so we can see the different warning messages and errors the package provides.

**Steps**

1. Create `my_design_bibd_1` using A, B, and C for the treatments, 4 blocks, and a `seed` of 42.

```{r}
#| error: true
#create my_design_bibd_1
my_design_bibd_1 <- design.bib(LETTERS[1:3], k = 4, seed = 42)
```

2. Create `my_design_bibd_2` using `LETTERS[1:8]` for treatments, 3 blocks, and a `seed` of 42.

```{r}
#| error: true
#create my_design_bibd_2
my_design_bibd_2 <- design.bib(LETTERS[1:8], k = 3, seed = 42)
```

3. Create `my_design_bibd_3` using A, B, C, and D as treatments, 4 blocks, and the same `seed`. Examine the `sketch` of the object.

```{r}
#create my_design_bibd_3
my_design_bibd_3 <- design.bib(LETTERS[1:4], k = 4, seed = 42)
my_design_bibd_3$sketch
```

Nice! You saw two different function errors which help lead you in the right direction with your design, plus one that works. When the design does work, the `sketch` parameter shows the design. The blocks are now the columns, however, unlike with RCBDs.

## BIBD - cat's kidney function

To be sure we truly understand what a BIBD looks like, let's build a dataset containing a BIBD from scratch.

Say we want to test the difference between four different wet foods in cats' diets on their kidney function. Cat food, however, is expensive, so we'll only test 3 foods per block to save some money. Our cats will be blocked by color of cat, as we aren't interested in that as part of our experiment. The outcome will be measured blood creatinine level, an indicator of kidney function and dysfunction in cats and humans alike.

**Steps**

1. The custom function `lambda()` has been loaded for you. Calculate `lambda` with `t` = 4, `k` = 3, and `r` = 3 to make sure a BIBD is possible.
2. Run the code to assemble the dataset. You can see the order in which the food treatments are used in each block.
3. Create `cat_model` with `aov()` according to the description of the experiment above and examine the results with `summary()`. Does type of wet food make a difference on creatinine levels?

```{r}
lambda <- function(t, k, r) {
            return((r*(k-1)) / (t-1))
          }

# Calculate lambda
lambda(4, 3, 3)

# Build the data.frame
creatinine <- c(1.98, 1.97, 2.35, 2.09, 1.87, 1.95, 2.08, 2.01, 1.84, 2.06, 1.97, 2.22)
food <- as.factor(c("A", "C", "D", "A", "B", "C", "B", "C", "D", "A", "B", "D"))
color <- as.factor(rep(c("Black", "White", "Orange", "Spotted"), each = 3))
cat_experiment <- as.data.frame(cbind(creatinine, food, color))

# Create cat_model and examine with summary()
cat_model <- aov(creatinine ~ food + color, data = cat_experiment)
summary(cat_model)
```

Purrfect! It seems there are no differences by type of wet food in kidney function. Hopefully now you can see how a BIBD comes to life!

## NHANES BIBD

Let's jump back into the NHANES data and pretend we have access to NHANES patients ages 18-45. Blocking the adults by race, stored in NHANES as `ridreth1`, we prescribe to our groups either no particular upper body weightlifting regimen, a weightlifting regimen, or a weightlifting regimen plus a prescribed daily vitamin supplement. This information is stored in a variable called `weightlift_treat`.

Those funding the study decide they want it to be a BIBD where only 2 treatments appear in each block. The outcome, arm circumference, is stored as `bmxarmc`. The `nhanes_final` data is loaded for you.

**Steps**

1. Calculate lambda where `t` = 3, `k` = 2, and `r` = 2. Does a BIBD exist here?
2. Create `weightlift_model` and examine the results. The experiment evaluates the outcome `bmxarmc`, where the treatment is stored in `weightlift_treat` and subjects are blocked by `ridreth1`.

```{r}
#| include: false
weightlift_treat <- c(1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1)
```

```{r}
# Data
nhanes_final <- read_rds("data/nhanes_final_18_45.rds")

# Calculate lambda
lambda(3, 2, 2)

# Create weightlift_model & examine results
weightlift_model <- aov(BMXARMC ~ weightlift_treat + RIDRETH1, data = nhanes_final)
summary(weightlift_model)
```

Nice! As it turns out, the weight lifting regimen doesn't seem to have a significant effect on arm circumference when the patient population is blocked by race.

# 4. Latin Squares, Graeco-Latin Squares, & Factorial experiments

Evaluate the NYC SAT scores data and deal with its missing values, then evaluate Latin Square, Graeco-Latin Square, and Factorial experiments.

## Latin squares

Theory. Coming soon ...
<p class="dc-cookie-banner-text">                DataCamp and our partners use cookies and similar technologies to improve your learning experience, offer data science content relevant to your interests, improve the site and to show more relevant advertisements. You can change your mind at any time.            



**1. Latin squares**

Latin Square designs will require us to think back to what we learned about Randomized Controlled Block Designs, as they'll expand on that design. Remember that the goal in RCBDs was to block on a factor to reduce the variance of the outcome due to that factor. Let's look at Latin square designs to see how they expand on that idea.

**2. Latin squares**

Latin square designs are, luckily, not that ominous: they expand on RCBDs by blocking on a second factor. Now, by designing an experiment in this way, we'll be able to control both of the blocking factors' variability and study our actual outcome of interest. Of note -- the treatment and both blocking factors need to have the same number of levels, which is a limitation of the Latin square design. One key assumption of a Latin square experiment is that there is no interaction between any combination of the treatment and blocking factors. If this turns out not to be true, the error term of the model will be inflated or overestimated.One nice thing about Latin square is you can analyze them just as you would a RCBD, by building a linear or other appropriate type of model, then examining it with anova() or using summary to examine the results if you used aov().

**3. Latin square diagram**

You may notice that a Latin Square looks a lot like a RCBD, and you're absolutely right about that. However, in the case of a Latin Square experiment, both the rows and columns of this grid are blocking factors. We know that this is a valid Latin square experiment because all treatments, A through D, appear four times each across the grid, and each treatment appears only once per block. There is only one A, B, C, and D per both row and column. Take a second to verify that for yourself. Visually verifying experimental design starts to feel a lot like playing the puzzle game, Sudoku.

**4. Why is it a Latin square?**

It may be easier to verify quickly that this is a Latin square design if we use shapes to represent each treatment. Personally, I can see a lot faster that there's only one circle per row and column than I can with letters.

**5. Intro to NYC scores**

For this chapter, we'll be using the nyc_scores dataset, which is an open dataset from the city of New York giving the SAT scores for the Reading, Writing, and Math sections from all accredited high schools in the 2014 to 2015 school year. It also includes variables giving the school's borough and stats on enrollment by race and ethnicity, among other variables. We'll use it in this chapter to conduct some educational experiments on a Tutoring Program's effect on the scores. Given that we did not collect this data ourselves, the tutoring program data has been fabricated, and will not be available in the original data should you download it yourself.

**6. Let's practice!**

Let's go to the exercises. We'll begin to get familiar with the nyc_scores dataset and build and analyze a Latin square experiment.

## NYC SAT Scores EDA

Math is a subject the U.S. is consistently behind the rest of the world on, so our experiments will focus on Math score. While the original dataset is an open dataset downloaded from <a href="https://www.kaggle.com/nycopendata/high-schools/data">Kaggle</a>, throughout this chapter I will add a few variables that will allow you to pretend you are an education researcher conducting experiments ideally aimed at raising students' scores, hopefully increasing the likelihood they will be admitted to colleges.

Before diving into analyzing the experiments, we should do some EDA to make sure we fully understand the `nyc_scores` data. In this lesson, we'll do experiments where we block by `Borough` and `Teacher_Education_Level`, so let's examine math scores by those variables. The `nyc_scores` dataset has been loaded for you.

**Steps**

1. Find the mean, variance, and median of `Average_Score_SAT_Math` by `Borough` using `dplyr` methods for EDA as we have used them throughout the course.

```{r}
# Data
# nyc_scores <- read_csv("data/nyc_scores.csv")
nyc_scores <- read_rds("data/nyc_scores.rds")

# Mean, var, and median of Math score by Borough
nyc_scores %>%
    group_by(Borough) %>% 
    summarize(mean = mean(Average_Score_SAT_Math, na.rm = TRUE),
              var = var(Average_Score_SAT_Math, na.rm = TRUE),
              median = median(Average_Score_SAT_Math, na.rm = TRUE))
```

2. Find the mean, variance, and median of `Average_Score_SAT_Math` by `Teacher_Education_Level` using `dplyr` EDA methods.

```{r}
# Mean, var, and median of Math score by Teacher Education Level
nyc_scores %>%
    group_by(Teacher_Education_Level) %>% 
    summarize(mean = mean(Average_Score_SAT_Math, na.rm = TRUE),
              var = var(Average_Score_SAT_Math, na.rm = TRUE),
              median = median(Average_Score_SAT_Math, na.rm = TRUE))
```

3. Find the mean, variance, and median of `Average_Score_SAT_Math` by both `Borough` and `Teacher_Education_Level` using `dplyr` EDA methods.

```{r}
# Mean, var, and median of Math score by both
nyc_scores %>%
    group_by(Borough, Teacher_Education_Level) %>% 
    summarize(mean = mean(Average_Score_SAT_Math, na.rm = TRUE),
              var = var(Average_Score_SAT_Math, na.rm = TRUE),
              median = median(Average_Score_SAT_Math, na.rm = TRUE))
```

Great job! Now that we've examined the data, we can move on to cleaning it, the next important step before analysis.

## Dealing with Missing Test Scores

If we want to use SAT scores as our outcome, we should examine missingness. Examine the pattern of missingness across all the variables in `nyc_scores` using `miss_var_summary()` from the `naniar` package. `naniar` integrates with Tidyverse code styling, including the pipe operator (`%>%`). 

There are 60 missing scores in each subject. Though there are many R packages which help with more advanced forms of imputation, such as `MICE`, `Amelia`, and `mi`, we will continue to use `simputation` and `impute_median()`. 

Create a new dataset, `nyc_scores_2` by imputing Math score by Borough, but note that `impute_median()` returns the imputed variable as type "impute". You'll convert the variable to the numeric in a separate step.

`simputation` and `dplyr` are loaded.

**Steps**

1. Load the `naniar` package.
2. Examine the missingness of variables in `nyc_scores` by piping it to `miss_var_summary()`.

```{r}
# Load naniar
library(naniar)

# Examine missingness with miss_var_summary()
nyc_scores %>% miss_var_summary()
```

3. Create `nyc_scores_2` by imputing the Average Math SAT score by Borough (we're only using Math in our experiments.)

```{r}
# Package
library(mice)

# Examine missingness with md.pattern()
md.pattern(nyc_scores)

# Impute the Math score by Borough
nyc_scores_2 <- simputation::impute_median(dat = nyc_scores, formula = Average_Score_SAT_Math ~ Borough)
```

4. Convert `nyc_scores_2$Average_Score_SAT_Math` to numeric.

```{r}
# Examine missingness with md.pattern()
md.pattern(nyc_scores)

# Impute the Math score by Borough
nyc_scores_2 <- simputation::impute_median(nyc_scores, Average_Score_SAT_Math ~ Borough)

# Convert Math score to numeric
nyc_scores_2$Average_Score_SAT_Math <- as.numeric(nyc_scores_2$Average_Score_SAT_Math)
```

5. Use `dplyr` to examine the median and mean of math score before and after imputation.

```{r}
# Examine missingness with md.pattern()
md.pattern(nyc_scores)

# Impute the Math score by Borough
nyc_scores_2 <- simputation::impute_median(nyc_scores, Average_Score_SAT_Math ~ Borough)

# Convert Math score to numeric
nyc_scores_2$Average_Score_SAT_Math <- as.numeric(nyc_scores_2$Average_Score_SAT_Math)

# Examine scores by Borough in both datasets, before and after imputation
nyc_scores %>% 
  group_by(Borough) %>% 
  summarize(median = median(Average_Score_SAT_Math, na.rm = TRUE), 
              mean = mean(Average_Score_SAT_Math, na.rm = TRUE))
nyc_scores_2 %>% 
  group_by(Borough) %>% 
  summarize(median = median(Average_Score_SAT_Math), 
              mean = mean(Average_Score_SAT_Math))
```

Nice job! Did the median scores change before and after imputation? (Hint: they shouldn't have changed by much, but rounding may have offset them by an integer or two.)

## Drawing Latin Squares with agricolae

We return, once again, to the `agricolae` package to examine what a Latin Square design can look like. Here's an example:

```{r}
#| eval: false
     [,1] [,2] [,3] [,4]
[1,] "B"  "D"  "A"  "C" 
[2,] "A"  "C"  "D"  "B" 
[3,] "D"  "B"  "C"  "A" 
[4,] "C"  "A"  "B"  "D"
```

Since a Latin Square experiment has two blocking factors, you can see that in this design, each treatment appears once in both each row (blocking factor 1) and each column (blocking factor 2).

Look at the help page for `design.lsd()` by typing `?design.lsd` in the console for any help you need designing your Latin Square experiment.

**Steps**

1. Load the `agricolae` package.
2. Create and view the sketch of a Latin Square design, `my_design_lsd`, using treatments A, B, C, D, & E, and a `seed` of 42.

```{r}
# Load agricolae
library(agricolae)

# Design a LS with 5 treatments A:E then look at the sketch
my_design_lsd <- design.lsd(trt = LETTERS[1:5], seed = 42)
my_design_lsd$sketch
```

Perhaps you're thinking to yourself 'This looks a lot like a RCBD'...bingo! It does, but as we know from the video, there are now two blocking factors in a LS design.

## Latin Square with NYC SAT Scores

To execute a Latin Square design on this data, suppose we want to know the effect of our tutoring program, which includes one-on-one tutoring, two small groups, and an in and after-school SAT prep class. A new dataset `nyc_scores_ls` is available that represents this experiment. Feel free to explore the dataset in the console.

We'll block by `Borough` and `Teacher_Education_Level` to reduce their known variance on the score outcome. `Borough` is a good blocking factor because schools in America are funded partly based on taxes paid in each city, so it will likely make a difference in the quality of education.

**Steps**

1. Use `lm()` to test the changes in `Average_Score_SAT_Math` using `nyc_scores_ls`.
2. Tidy `nyc_scores_ls_lm` with the appropriate `broom` function.
3. Examine `nyc_scores_ls_lm` with `anova()`.

```{r}
# Data
nyc_scores_ls <- read_rds("data/nyc_scores_ls.rds")

# Build nyc_scores_ls_lm
nyc_scores_ls_lm <- lm(Average_Score_SAT_Math ~ Tutoring_Program + Borough + Teacher_Education_Level,
                       data = nyc_scores_ls)

# Tidy the results with broom
tidy(nyc_scores_ls_lm)

# Examine the results with anova
anova(nyc_scores_ls_lm)
```

> *Question*
> ---
> At the 0.05 significance level, do we have evidence to believe the tutoring program has an effect on math SAT scores, when blocked by `Borough` and `Teacher_Education_Level`?<br>
> <br>
> ✅ Nope! Given the p-value, we have no reason to reject the null hypothesis.<br>
> ⬜ Yes! Given the p-value, we have reason to believe that the tutoring program had an effect on the Math score.<br>

Excellent! It seems that when we block for `Borough` of the school and `Teacher_Education_Level`, our `Tutoring_Program` isn't having a statistically significant effect on the Math SAT score.

## Graeco-Latin squares

Theory. Coming soon ...
<p class="dc-cookie-banner-text">                DataCamp and our partners use cookies and similar technologies to improve your learning experience, offer data science content relevant to your interests, improve the site and to show more relevant advertisements. You can change your mind at any time.            



**1. Graeco-Latin squares**

Graeco-Latin squares build on Latin squares by adding yet another blocking factor. Let's introduce them quickly so you can get back to the exercises.

**2. Graeco-Latin squares**

This slide might look familiar if you watched the last video, and for good reason; a Graeco-Latin square design simply adds one more blocking factor, but shares the same properties of a Latin square experiment otherwise. We again have to have the same number of treatments as levels in the blocking factors and are operating under the assumption that none of those treatment and blocking factors interact.

**3. Graeco-Latin squares**

Traditionally, Graeco-Latin squares are drawn very similar to Latin squares, with the addition of Greek letters, which indicate the third blocking factors are added next to the Latin letters, which indicate the treatment, as seen in this diagram. For this course, we won't be accessing the Greek alphabet to draw Graeco-Latin squares in R, we'll use the Latin alphabet and numbers. As before, the treatment occurs only once in each row and column. The rows are one blocking factor, the columns another, and the Greek letters the third, so here we notice that A, B, C, and D appear not only once per column and row, but they are also each paired with alpha, beta, gamma, and delta only once.

**4. GLS - explanation**

Again, it may be easier to see why this is a Graeco-Latin square when the letters, representing the treatments, are converted to shapes. You can see that the circles appear only once per row and column, and they're only paired up with each Greek letter once for each of alpha, beta, gamma, and delta, which represent the third blocking factor.

**5. Let's practice!**

Because Graeco-Latin Squares are so similar to Latin Squares, I'm keeping this video short and sweet. Jump in and try the exercises, where we expand the Latin Square experiment from the last lesson.

## NYC SAT Scores Data Viz

In the last lesson, when discussing Latin Squares, we did numerical EDA in the form of looking at means, variances, and medians of the math SAT scores. Another crucial part of the EDA is data visualization, as it often helps in spotting outliers plus gives you a visual representation of the distribution of your variables. 

`ggplot2` has been loaded for you and the `nyc_scores` dataset is available. Create and examine the requested boxplot. How do the medians differ by Borough? How many outliers are present, and where are they mostly present?

**Steps**

1. Create a boxplot of Math SAT scores by `Borough`.
2. Run the code to include a title: `"Average SAT Math Scores by Borough, NYC"`.
3. Change the x- and y-axis labels to read `"Borough (NYC)"` and `"Average SAT Math Scores (2014-15)"`, respectively, using the correct arguments to `labs()`.

```{r}
# Create a boxplot of Math scores by Borough, with a title and x/y axis labels
ggplot(nyc_scores, aes(Borough, Average_Score_SAT_Math)) +
  geom_boxplot() + 
  labs(title = "Average SAT Math Scores by Borough, NYC",
       x = "Borough (NYC)",
       y = "Average SAT Math Score (2014-15)")
```

Beautiful! It's interesting to see the different distribution of scores by Borough and to see that every borough has scores that are outliers, though some more than others.

## Drawing Graeco-Latin Squares with agricolae

As we've seen, `agricolae` provides us the ability to draw all of the experimental designs we've used so far, and they can also draw Graeco-Latin squares. One difference in the input to `design.graeco()` that we haven't seen before is that we'll need to input 2 vectors, `trt1` and `trt2`, which must be of equal length. You can think of `trt1` as your actual treatment and `trt2` as one of your blocking factors. `agricolae` has been loaded for you.

**Steps**

1. Create vectors `trt1` with `LETTERS` A through E and `trt2` with numbers 1 through 5.
2. Make `my_graeco_design` with `design.graeco()`, using and `seed = 42`.
3. Examine the `parameters` and `sketch` of `my_graeco_design`.

```{r}
# Create trt1 and trt2
trt1 <- LETTERS[1:5]
trt2 <- 1:5

# Create my_graeco_design
my_graeco_design <- design.graeco(trt1, trt2, seed = 42)

# Examine the parameters and sketch
my_graeco_design$parameters
my_graeco_design$sketch
```

Superb! You can see that this time the sketch object includes your treatment (the capital letter) and a blocking factor (the number.)

## Graeco-Latin Square with NYC SAT Scores

Recall that our Latin Square exercise in this chapter tested the effect of our tutoring program, blocked by `Borough` and `Teacher_Education_Level`. 

For our Graeco-Latin Square, say we also want to block out the known effect of `Homework_Type`, which indicates what kind of homework the student was given: individual only, small or large group homework, or some combination. We can add this as another blocking factor to create a Graeco-Latin Square experiment.

**Steps**

1. Use `lm()` to test the changes in `Average_Score_SAT_Math` using the `nyc_scores_gls` data.
2. Tidy `nyc_scores_gls_lm` with the appropriate `broom` function.
3. Examine `nyc_scores_gls_lm` with `anova()`.

```{r}
# Data
nyc_scores_gls <- read_rds("data/nyc_scores_gls.rds")

# Build nyc_scores_gls_lm
nyc_scores_gls_lm <- lm(Average_Score_SAT_Math ~ Tutoring_Program + Borough + Teacher_Education_Level + Homework_Type,
                        data = nyc_scores_gls)

# Tidy the results with broom
tidy(nyc_scores_gls_lm)

# Examine the results with anova
anova(nyc_scores_gls_lm)
```

> *Question*
> ---
> At the 0.05 significance level, do we have evidence to believe the tutoring program has an effect on math SAT scores, when blocked by `Borough`, `Teacher_Education_Level`, and `Homework_Type`?<br>
> <br>
> ✅ Nope! Given the p-value, we have no reason to reject the null hypothesis.<br>
> ⬜ Yes! Given the p-value, we have reason to believe that the tutoring program had an effect on the Math score.<br>

Excellent! It seems that when we block for `Borough` of the school and `Teacher_Education_Level`, our `Tutoring_Program` isn't having a statistically significant effect on the Math SAT score.")\n\nsuccess_msg("Bravo! It seems that here, when blocked out by all the other factors, our Tutoring program has no effect on the Math score.

## Factorial experiments

Theory. Coming soon ...
<p class="dc-cookie-banner-text">                DataCamp and our partners use cookies and similar technologies to improve your learning experience, offer data science content relevant to your interests, improve the site and to show more relevant advertisements. You can change your mind at any time.            



**1. Factorial experiments**

Recall that in the Latin squares and Graeco-Latin squares videos I mentioned that one key assumption for both of those designs is no interaction between any combination of the blocking factors and the treatment. That's true, and was also foreshadowing, because there is an experiment type that considers interaction of factor variables called a factorial experiment.

**2. Factorial designs**

A factorial design is one in which two or more factor variables are combined and crossed in an experiment. Each combination of all the levels of all the factors will be considered separately as a factor affecting the outcome, meaning that we'll measure the outcome at every possible combination of all the factors. Let's dig into an example to see how this works.

**3. Factorial example**

Our hypothetical plant growth experiment will require testing all the combinations of levels of high or low water and high or low light. In the diagram, you can see that this includes low water/low light, low water/high light, high water/low light, and high water/high light. Every combination's effect on the outcome will be tested, and you can use Tukey's HSD test to specifically test differences between plant growth for each possible combination.

**4. 2^k factorial experiments**

For this course, we'll only be analyzing 2k factorial experiments, a type where you have some number, k, of two-level factor variables. Usually, a 2k experiment will have factors with a "high" and a "low" level, as we saw in the plant experiment we just discussed. Often, for reasons of cost, the high and low level of factor variables are chosen to be tested because the more levels a factor has, the more combinations begin to appear and the more money the experiment costs. The simplest 2k factor case is 2 factor variables with two levels each. All possible combinations will be tested, but any number, k, of variables can be selected. Often, running a 2k experiment with many 2-level factor variables is a strategy employed to find the factors that are important in your experiment quickly. For example, you run six factor variables and examine results to find what is most affecting the outcome. Keep in mind though, that's 64 different combinations to model and explore!

**5. Let's practice!**

Let's dive into the exercises to see how to design, analyze, and evaluate factorial experiments.

## NYC SAT Scores Factorial EDA

Let's do some more EDA before we dive into the analysis of our factorial experiment. 

Let's test the effect of `Percent_Black_HL`, `Percent_Tested_HL`, and `Tutoring_Program` on the outcome, `Average_Score_SAT_Math`. The `HL` stands for high-low, where a 1 indicates respectively that less than 50% of Black students or that less than 50% of all students in an entire school were tested, and a 2 indicates that greater than 50% of either were tested.

Build a boxplot of each factor vs. the outcome to have an idea of which have a difference in median by factor level (ultimately, mean difference is what's tested.) The `nyc_scores` dataset has been loaded for you.

**Steps**

1. Load `ggplot2`. Create a boxplot of the outcome versus `Tutoring_Program`.

```{r}
# Data
nyc_scores <- read_rds("data/nyc_scores2.rds")

# Load ggplot2
library(ggplot2)

# Build the boxplot for the tutoring program vs. Math SAT score
ggplot(nyc_scores,
       aes(Tutoring_Program, Average_Score_SAT_Math)) + 
       geom_boxplot()
```

2. Using `ggplot2`, create a boxplot of the outcome versus `Percent_Black_HL`.

```{r}
# Build the boxplot for percent black vs. Math SAT score
ggplot(nyc_scores,
       aes(Percent_Black_HL, Average_Score_SAT_Math)) + 
    geom_boxplot()
```

3. Using `ggplot2`, create a boxplot of the outcome versus `Percent_Tested_HL`.

```{r}
# Build the boxplot for percent tested vs. Math SAT score
ggplot(nyc_scores,
       aes(Percent_Tested_HL, Average_Score_SAT_Math)) + 
    geom_boxplot()
```

Excellent! Now, let's move on to the analysis of these factors on the score.

## Factorial Experiment with NYC SAT Scores

Now we want to examine the effect of tutoring programs on the NYC schools' SAT Math score. As noted in the last exercise: the variable `Tutoring_Program` is simply `yes` or `no`, depending on if a school got a tutoring program implemented. For `Percent_Black_HL` and `Percent_Tested_HL`, `HL` stands for high/low. A 1 indicates less than 50% Black students or overall students tested, and a 2 indicates greater than 50% of both.

Remember that because we intend to test all of the possible combinations of factor levels, we need to write the formula like: `outcome ~ factor1 * factor2 * factor3`.

**Steps**

1. Use `aov()` to create a model to test how `Percent_Tested_HL`, `Percent_Black_HL`, and `Tutoring_Program` affect the outcome `Average_Score_SAT_Math`. 
2. Save the outcome as a model object, `nyc_scores_factorial`, and examine this with `tidy()`.

```{r}
# Create nyc_scores_factorial and examine the results
nyc_scores_factorial <- aov(Average_Score_SAT_Math ~ Percent_Tested_HL * Percent_Black_HL * Tutoring_Program, data = nyc_scores)

tidy(nyc_scores_factorial)
```

Whoo! We can see from the results that we can reject the null hypothesis that there is no difference in score based on tutoring program availability. We can also see from the low p-values that there are some interaction effects between the Percent Black and Percent Tested and the tutoring program. Next we need to check the model.

## Evaluating the NYC SAT Scores Factorial Model

We've built our model, so we know what's next: model checking! We need to examine both if our outcome and our model residuals are normally distributed. We'll check the normality assumption using `shapiro.test()`. A low p-value means we can reject the null hypothesis that the sample came from a normally distributed population.

Let's carry out the requisite model checks for our 2^k factorial model, `nyc_scores_factorial`, which has been loaded for you.

**Steps**

1. Test the outcome `Average_Score_SAT_Math` from `nyc_scores` for normality using `shapiro.test()`.
2. Set up a 2 by 2 grid for plots and plot the `nyc_scores_factorial` model object to create the residual plots.

```{r}
# Use shapiro.test() to test the outcome
shapiro.test(nyc_scores$Average_Score_SAT_Math)

# Plot nyc_scores_factorial to examine residuals
par(mfrow = c(2,2))
plot(nyc_scores_factorial)
```

Brilliant! The model appears to be fairly well fit, though our evidence indicates the score may not be from a normally distributed population. Looking at the Q-Q plot, we can see that towards the higher end, the points are not on the line, so we may not be dealing with normality here. If we had more time, we might consider a transformation on the outcome to move towards normality.

## What's next in experimental design

Theory. Coming soon ...
<p class="dc-cookie-banner-text">                DataCamp and our partners use cookies and similar technologies to improve your learning experience, offer data science content relevant to your interests, improve the site and to show more relevant advertisements. You can change your mind at any time.            



**1. What's next in experimental design**

We've reached the end of this course, but there's still much more to consider with respect to designing experiments.

**2. What's next?**

The most immediate example is that there are a huge number of other factorial experiment types out there. We only looked at 2^k factorial experiments, when in fact you can design factorial experiments that include factor variables with more than 2 levels. Of course, all of the combining involved in a lot of factors and factor levels means you might also start to consider a fractional factorial design, where you only consider some fraction of the number of possible combinations.You could continue to expand your knowledge by designing experiments with random factors or could utilize nested or split plot designs, or even a lattice design. There are many other types of designs and plenty of classical texts and other resources if you're interested in deepening your knowledge.What's most important is that design is a valued and integrated part of the experimental process. From A/B tests to studies on people's health, education, or loan fitness, as we considered here, and even beyond, the steps of an experiment should be thought out, the outcome and how to measure it considered, and the possible factors affecting the outcome thought through. In reality, there's always going to be something you'll miss, some unmeasured confounder that introduces noise into your model. Hopefully, what you've picked up throughout this course are some ideas to help reduce that noise and begin to create clear experiments built on the principles of reproducibility, randomization, and replication.

**3. Go forth and design experiments!**

Now it's your turn. Go forth in the world and design experiments!

