# Your turn 25. Let's analyze a dataset, `ToothGrowth` built-into R. It's results from a study on the effect of vitamin C on tooth growth in guinea pigs. The researchers were interested in the effects of two variables:

• `supp`: type of Vitamin C given to the guinea pigs (**OJ** = orange juice; **VC** = a liquid supplement)

• `dose`: the amount of Vitamin C given (**low**, **med**, or **high**).

This gives us a total of 6 treatment groups (*OJ-low, OJ-med, OJ-high, VC-low, VC-med, VC-high*). A total of 60 guinea pigs were randomly assigned, 10 to each of those groups.

Summary statistics and visualizations are displayed below.

**(a)** Evaluate the conditions necessary to conduct an AxB ANOVA. Are you concerned about any of the assumptions being violated?

**(b)** Based on the interaction plot, do you think we'll conclude that a significant interaction effect exists? Explain.

**(c)** Conduct an AxB (Supplement x Dose) ANOVA and draw a conclusion about the interaction effect.

**(d)** If the interaction effect is **not statistically significant**, then write a conclusion regarding the main effects of the type of supplement and dose. If the interaction effect **is statistically significant**, split your data by `dose` and test for differences between the **OJ** and **VC** treatments.

This table shows the means (and standard deviations) of the 10 guinea pigs in each combination of treatments. . | Orange Juice | AbsorbAcid | Total ------:| ------------:| -------------:| ---------------:| Low | 13.23 (4.460) | 07.98 (2.747) | **10.605 (4.4998)** Med | 22.70 (3.911) | 16.77 (2.515) | **19.735 (4.4154)** High | 26.06 (2.655) | 26.14 (4.798) | **26.100 (3.7742)** **Total** | **20.66 (6.606)** | **16.96 (8.266)** | **N=60, 18.813 (7.6493)** ```{r} # Density plot of each group densityplot(~len | dose*supp, data=ToothGrowth, main="Tooth Growth by Dose and Supplement", xlab="Tooth Growth", layout=c(3,2)) # Side-by-side boxplots bwplot(supp ~ len | dose, data=ToothGrowth, ylab="supplement", xlab="Tooth Growth", main="Tooth Growth by Dose and Supplement", layout=(c(1,3))) # Interaction plot with(ToothGrowth, interaction.plot(x.factor=dose, trace.factor=supp, response=len, fun=mean, type="b", legend=T, ylab="mean tooth length", main="Interaction Plot", pch=c(1,19))) # Levene's test leveneTest(len ~ supp * dose, data=ToothGrowth) ```

```{r, eval=FALSE} ########## Here is where the ########## ########## lab report begins ########## ``` **(a) Evaluate the conditions necessary to conduct an AxB ANOVA. Are you concerned about any of the assumptions being violated?** **Type your answer here**

**(b) Based on the interaction plot, do you think we'll conclude that a significant interaction effect exists? Explain.** **Type your answer here**

**(c) Conduct an AxB (Supplement x Dose) ANOVA and draw a conclusion about the interaction effect.** ```{r} # Construct your anova model and store it with a name # Remember, the model should look like y ~ x1 + x2 + x1x2 # so your syntax should be: # name_of_model <- aov(y ~ x1 + x2 + x1x2, data=name_of_data) # The variables are len, supp, and dose. The data is ToothGrowth # Replace the XXXXX with the name of the model you just constructed # This will summarize the results and print an interaction plot anova(XXXXX) plotModel(XXXXX) ```

**Type any conclusion you can make from your analysis. Did you find a significant interaction?**

**(d) If the interaction effect *is statistically significant*, split your data by `dose` and test for differences between the *OJ* and *VC* treatments.** ```{r} # Split the data with dplyr # I'll go ahead and split the data into 3 subsets: # (1) lowdose, (2) meddose, (3) highdose lowdose <- ToothGrowth %>% filter(dose=="low") meddose <- ToothGrowth %>% filter(dose=="med") highdose <- ToothGrowth %>% filter(dose=="high") # Run a t-test for each of the subsets # Replace the XXXXX values with the appropriate variables t.test(XXXXX ~ XXXXX, data=lowdose, var.equal=FALSE, alternative = c("greater")) t.test(XXXXX ~ XXXXX, data=meddose, var.equal=FALSE, alternative = c("greater")) t.test(XXXXX ~ XXXXX, data=highdose, var.equal=FALSE, alternative = c("less")) # I didn't show this in class, but you could run a bunch of post-hoc tests # instead of splitting the data. Here's Tukey's test for all 6 # combinations of treatments. Just replace the name_of_model TukeyHSD(name_of_model, conf.level=.95) # Bonferroni amd Holm's corrections would work, too. # Here I test the lowdose means with the Bonferroni adjustment with(lowdose, pairwise.t.test(len, supp, p.adjust.method="bonferroni")) # Here I test the meddose means with Holm's method with(meddose, pairwise.t.test(len, supp, p.adjust.method="holm")) ```

26. This time, analyze that same `ToothGrowth` dataset using a randomization-based equivalent to the AxB ANOVA. Test the interaction effect; then, split the data by dose and run randomization-based tests to compare pairs of means.

```{r} # Follow along and see if you understand this code. # I'll use the broom package library(broom) # First, we'll re-run the model aovmodel <- aov(len ~ supp * dose, data=ToothGrowth) # Now we'll use the "tidy" command to store those results tidy(aovmodel) # Store the 3rd statistic (F-value for interaction = 4.8465) obs_F_interaction<- tidy(aovmodel)$statistic[3] # Now let's run 5,000 randomizations, each time storing the F-value # We will shuffle the interaction term random_interact <- do(5000) * tidy(aov(len ~ supp + dose + shuffle(supp):shuffle(dose), data=ToothGrowth))$statistic[3] # We'll plot a histogram of the F values from these 10,000 randomizations histogram(~random_interact, data = random_interact, xlim=c(-.1, 7), xlab = "F-statistic for interaction", width=.05, col="steelblue", border="white") # p-value tally(~ (random_interact >= obs_F_interaction), format="prop", data=random_interact) # With our split data, let's run a randomized t-test test.stat <- diffmean(len ~ supp, data=lowdose) rnd.stat <- do(5000) * diffmean(len ~ shuffle(supp), data=lowdose) p <- tally(~ (diffmean <= test.stat), format="prop", data=rnd.stat) p ```