---
title: "Analyzing Golf Ball Backspin"
author: "Tyler Wiederich"
date: "Fall 2022"
output: 
  word_document: default
  pdf_document: default
  html_document: default
---

```{r setup, include=FALSE}
options(continue = " ")
knitr::opts_chunk$set(prompt = TRUE, comment = NA, background = "white",
tidy = TRUE, tidy.opts = list(width.cutoff = 60, blank = FALSE))
```

# Introduction

In the previous lab, we analyzed the carry distance of various golf ball types via the difference of two means. In this lab, we will look at the backspin (measured in RPM) for these same types. The data were collected as a part of the experiment by Mark Crossfield. We will analyze this data using ANOVA methodology, and we will gain experience with the ``emmeans`` R package.

* Backspin is important for a golfer because it allows the golfer to control the ball as it lands. This is especially important when taking shots with clubs called irons. The data collected are shots hit using a 7-iron club.

* For this example, we will assume that more backspin is more desirable.

The corresponding data can be found in the 7IronSpin.csv file. Below is a preview of the data set.

```{r dataview}
# make sure you have this file in the directory of your .rmd file
spin <- read.csv("7IronSpin.csv")

head(spin)
```

**Questions**

* What is the factor? What are the factor levels?

* Identify the response variable.

* What is the experimental unit?

* Describe the design of this experiment. What type of design is this? Are there any drawbacks of this design? What do you think could be done to improve the design of the experiment?

Let's take a look at a plot of the spin rates obtained by Mark Crossfield.

```{r visual, echo=FALSE, fig.width=7, fig.align='center'}
par(mfrow = c(1,1))
boxplot(formula = Backspin ~ Ball, data = spin, main = "Box and Dot Plot", ylab = "Backspin (RPM)", xlab = "Ball Type", pars = list(outpch = NA), col = NA, cex.axis = 0.85)

#seed set for same jitter
set.seed(4422)
stripchart(x = spin$Backspin ~ spin$Ball, lwd = 2, col = "red", method = "jitter", vertical = T, pch = 1, add = T)
```

**Question** - Do you think there are differences among the means? Remember, a plot like this is meant to give an initial impression of the data relative to the research hypotheses. Focus on the shifts of the points, the variability of the points, and the location of the medians.

\newpage

# Analysis of Variance

Recall that our goal is to determine whether there is a difference in mean backspin levels among different ball types.

ANOVA methods allow us to test the following set of hypotheses 
$$\begin{align}
H_0 :\,\, &\mu_1 = \mu_2 = \ldots = \mu_t \\
H_a :\,\, & \text{At least one pair of means are not equal}
\end{align}$$
where $t$ is the number of treatments. For our problem, we can express the full set of hypotheses as 
$$\begin{align}
H_0:\,\, &\mu_{QStar}=\mu_{QStarT}=\mu_{Soft}=\mu_{ZStar}=\mu_{ZStarX} \\
H_a:\,\, &\text{At least one pair of means are not equal}
\end{align}$$
We will conduct this hypothesis test using $\alpha = 0.05$.

In the class notes, we discussed the test statistic method and the p-value method for the evaluation of the aforementioned hypothesis test. 

* For the test statistic method, we use $F_{obs}=\frac{MST}{MSE}$ where $MST$ and $MSE$ represent the mean square for the treatments and the mean square for the error, respectively. Recall that we reject $H_0$ when $F_{obs} > F_{\alpha,t-1,N-t}$ where $N$ denotes the total number of observations and $t$ denotes the number of treatments.

* For the p-value method, we calculate the $P(X>F_{obs})$ where $X$ is a random variable from an F distribution with $\nu_1 = t-1$ and $\nu_2 = N-t$ degrees of freedom and $F_{obs}$ is calculated as above. We reject $H_0$ when the p-value is less than $\alpha$.

Let's proceed with these calculations piece by piece first, and then we will verify our calculations using built in R functions. We will need to calculate $MST$ and $MSE$. Recall, the formulas are
$$\begin{align}
MST &= \frac{\sum\limits_{i=1}^{t}n_i(\bar{y}_{i\cdot}-\bar{y}_{\cdot\cdot})^2}{t-1} \\[8pt]
MSE &= \frac{\sum\limits_{i=1}^{t}\sum\limits_{j=1}^{n_i}(y_{ij}-\bar{y}_{i\cdot})^2}{N-t}
\end{align}$$

```{r functions, include=FALSE}
sumSq <- function(x){
  sum((x-mean(x))^2)
}
```

```{r ANOVA}
#We will calculate mean square of treatments first. 
#We need our treatment means, the grand/overall mean, and the sample sizes
save.mean <- aggregate(Backspin~Ball, data = spin, FUN = mean)[,2]
save.n <- aggregate(Backspin~Ball, data = spin, FUN = length)[,2]
tot.mean <- mean(spin$Backspin)
N <- length(spin$Backspin)
t <- length(unique(spin$Ball))

#now we compute mst
mst <- sum(save.n*(save.mean - tot.mean)^2)/(t-1)

#now we compute mse.
#We evaluate the inner summation for each treatment.
#sumSq is a custom function. Please see the markdown document for details
inner <- aggregate(Backspin~Ball, data = spin, FUN = sumSq)[,2]

#now we add those together and divide by N-t
mse <- sum(inner)/(N-t)

#to calculate our observed F statistic we take MST/MSE
F.obs <- mst/mse

#now we need the critical value for our F-statistic
nu1 <- t-1
nu2 <- N-t
alpha <- 0.05
F.crit <- qf(alpha, df1 = nu1, df2 = nu2, lower.tail = F)
F.crit
#now let's calculate the p-value
p.val <- pf(F.obs, df1 = nu1, df2 = nu2, lower.tail = F)
```

Now, let's take a look at the ANOVA table we just computed (the code is hidden).

```{r anovaTable, echo=FALSE}
Df = c(t-1,N-t)
Sum.Sq = c(sum(save.n*(save.mean - tot.mean)^2),sum(inner))
Mean.Sq = round(c(mst, mse))
F.value = c(round(F.obs,2),"")
p.value = c(formatC(p.val, format = "e", digits = 2),"")
data.frame(Df,Sum.Sq, Mean.Sq, F.value, p.value, row.names = c("Treatment", "Error"))
```

Let's take a look at what R produces and check to see if it matches our output. 

```{r ANOVASimple}
spin$Ball <- factor(spin$Ball)
mod.fit <- aov(Backspin ~ Ball, data = spin)
summary(mod.fit)
```

As we can see, the ``aov()`` function in R outputs the same values that we obtained using our own calculations. Now, let's use the results to make a conclusion for the original hypothesis test. 

* For the test statistic, we see that `r round(F.obs,2)` $>$ `r round(F.crit,2)`, therefore, we reject $H_0$ and conclude there is at least one pair of ball type means that are different.

* For the p-value method `r formatC(p.val, format = "e", digits = 2)`$< 0.05$, therefore, we reject $H_0$ and conclude there is at least one pair of ball type means that are different.

\newpage

**Questions**

* Think about what we have concluded thus far. Do you think our analysis is finished?

* What important questions have been left unanswered? How may we go about answering these questions?

\newpage

# Multiple Comparisons

After we reject $H_0$ in an ANOVA setting, we must now locate which population means are different. With most research questions, it is not useful to know only that there exists some difference among the means. More often than not we wish to know *which pairs of means are different*. This is the question that multiple comparisons help us to answer.

The definition of multiple comparisons is quite literal. In statistics, multiple comparisons refers to conducting a set of inferences simultaneously. In our example, we are interested in comparing all pairs of means. We now wish to extend the analysis a step further to reach conclusions about which means are higher and which means are lower.

The easiest way to do this in R is to use the ``pairwise.t.test()`` function. We see the results below.

```{r pairwiseNaive}
pairwise.t.test(x = spin$Backspin, g = spin$Ball, p.adjust.method = "none", alternative = "two.sided")
```

Unsurprisingly, ``pairwise.t.test()`` computes p-values for hypothesis tests involving two independent samples. The tests are for all possible combinations of our factor levels. We can use the entries in the output to see which comparisons yield significant differences. However, there is a problem with the way we have performed this test. Notice, we have not adjusted our p-values for the fact that we are performing multiple tests. Therefore, our experimentwise error rate will be inflated.

Let's take a look at the p-values when using the Bonferroni adjustment.

```{r bonf}
pairwise.t.test(x = spin$Backspin, g = spin$Ball, p.adjust.method = "bonferroni", alternative = "two.sided")
```

Look at what happened to the p-values. The Bonferroni adjustment ensures that our experimentwise error rate will be no larger than our stated $\alpha$ level. 

While it does not make a difference in this example, it is easy to see that we could very well run into a case where our results are significant (reject means are equal) *without* a multiple comparison adjustment, but they are *not significant* (don't reject means are equal) once we control for type I error appropriately.

Another popular adjustment that controls experimentwise error appropriately is *Tukey's Honestly Significant Difference (HSD)* adjustment. There is a built-in function in R which accepts the object from an ANOVA model and performs the adjustment. 

```{r tukey}
TukeyHSD(x = mod.fit, conf.level = 0.95)
```

The nice thing about this function is it also provides the confidence intervals. You can see the adjusted p values in the last column. 

**Question** - Rank the p-values obtained for the Soft Feel vs. Q Star brands using the three pairwise tests above from most conservative to least conservative (conservative means less likely to reject $H_0$). What do you find?

\newpage

# The ``emmeans`` Package

The [``emmeans``](https://cran.r-project.org/web/packages/emmeans/vignettes/) (Estimated Marginal Means) package is the most comprehensive R package so far for performing multiple comparisons. It has support for many different types of statistical models (beyond ANOVA), dynamic plotting options, and complicated experimental designs (important for STAT 802). Chris believes it will be the standard in the future. 

Let's take a look at some of the functionality of the package. The first thing we have to do is create an ``emmeans`` object using the ANOVA model we fit.

```{r em1}
library(emmeans)
#the specs specifies the predictors or factors over which we want means computed
calc.est <- emmeans(mod.fit, specs = ~Ball)

#by default emmeans displays an optimal amount of digits based on a criteria which you can read about on the vignette;
#we override this with the command below
emm_options(opt.digits = F)

#We can obtain confidence intervals for pairwise comparisons 
confint(contrast(calc.est, method = "pairwise"), level = 0.95)
#if we want to specify an adjustment other than Tukey we add the 'adjust =' argument with the name of the adjustment included

#we can also obtain p-values
test(contrast(calc.est, method = "pairwise"))
```
**Question** - For which pairs of means are there significant differences? How do you know?

Notice that the default adjustment is Tukey's HSD. In fact, six different multiple comparison methods are supported and they can be found in the documentation. Now, let's take a look at some of the graphical options.

We can use the package to obtain visuals of the results. The plotting mechanism utilizes the R package ``ggplot2``, so you must install this package if you wish to create the plots shown. The first plot is the "pairwise P-value plot" shown below

```{r plot1, fig.width=6}
pwpp(calc.est, adjust = "bonferroni")
```

Here, we see the means for each factor level sorted on the y axis. The p-value for each pairwise comparison is denoted by the horizontal position of the vertical line segment along the x axis. The colors and positions of the dots on the vertical line segments indicate the groups being compared. You can view the other plotting options in the documentation.

\newpage

We can also create a plot similar to what we see in the notes, but this is not part of the package.

```{r plot2, echo=FALSE}
plot(factor(c("Qstar","SoftFeel","ZStarXV","ZStar","QStarTour"), levels = c("Qstar","SoftFeel","ZStarXV","ZStar","QStarTour")), y = rep(0,5), xlab = "", ylab = "", yaxt = "n", ylim = c(0.6,10), xaxt = "n", frame.plot = F,)
segments(0.7,2,2.3,2,col = "red", lwd = 5)
segments(2.6,3,3.4,3, col = "blue", lwd = 5)
segments(3.7,4,5.3,4, col = "green", lwd = 5)
text(x = c(1:5), y = rep(1,5), labels = c("QstarTour","ZStar","ZStarXV","SoftFeel","QStar"))
```

**Exploration**

* Perform ANOVA on the carry distance data from the week 10 lab (you may use the ``aov()`` function and you can skip the computations of the mst and mse).

* Regardless of your findings for the ANOVA, use the ``emmeans`` package to find confidence intervals for the pairwise comparisons. Use the Bonferroni adjustment.

* Produce the pairwise p-value plot from the ``emmeans`` package. Explain and interpret the visualization.
