You are expected to complete this assignment on your own: that is, you may discuss general ideas with others, but the writeup of the work must be entirely your own. If your assignment is unreasonably similar to that of another student, you can expect to be asked to explain yourself.
If you run into problems on this assignment, it is up to you to figure out what to do. The only exception is if it is impossible for you to complete this assignment, for example a data file cannot be read. (There is a difference between you not knowing how to do something, which you have to figure out, and something being impossible, which you are allowed to contact me about.)
You must hand in a rendered document that shows your code, the output that the code produces, and your answers to the questions. This should be a file with .html on the end of its name. There is no credit for handing in your unrendered document (ending in .qmd), because the grader cannot then see whether the code in it runs properly. After you have handed in your file, you should be able to see (in Attempts) what file you handed in, and you should make a habit of checking that you did indeed hand in what you intended to, and that it displays as you expect.
Hint: render your document frequently, and solve any problems as they come up, rather than trying to do so at the end (when you may be close to the due date). If your document will not successfully render, it is because of an error in your code that you will have to find and fix. The error message will tell you where the problem is, but it is up to you to sort out what the problem is.
64 subjects were given some simple tasks to do, and the number of tasks completed in 10 minutes was recorded. This is shown in the column response in our data. Before attempting the tasks, each subject was randomly assigned to either receive or not receive a drug (in xdrug), and also to receive or not receive electroshock therapy (in xshock). The researchers were interested in whether either the drug or the electroshock or the combination of both affected the number of tasks the subject could solve. (The subjects were all fully informed of the risks, and agreed to sign consent forms.) The data are in http://ritsokiguess.site/datafiles/drugshock.csv.
As usual. My dataframe is called task (call yours whatever you like):
my_url <- "http://ritsokiguess.site/datafiles/drugshock.csv"
task <- read_csv(my_url)Rows: 64 Columns: 3
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): xdrug, xshock
dbl (1): response
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.taskNote that the categories of the two categorical variables were not just yes and no, but variations on yes and no that you can tell apart.
Extra: nothing challenging for you, but I had some work to do to get the data into this form for you (mostly done while waiting at Service Ontario to get my health card renewed). The data came from here:
data("Shkdrug", package = "BSDA")
ShkdrugI think this was originally intended for a one-way analysis treating the four treatment combinations as levels of a single factor, but it looks like a two-factor analysis to me. Let’s see what those treatments are:
Shkdrug %>% count(treatment)It looks like Drug and NoDrug, but sometimes NoDrug is NoDg. It also looks like Shock and NoShock, except that it isn’t that either. So we have some tidying up to do. First, split up the treatment column into two columns, one for the drug and one for electroshock, which is done with separate_wider_delim:
Shkdrug %>% separate_wider_delim(treatment,
names = c("xdrug", "xshock"),
delim = "/") -> task
taskI called my new dataframe task, and for some reason called the new columns xdrug and xshock. Possibly I was distracted by waiting for my number to come up. Anyway. The next thing is to sort out those category levels. There are several possibilities for both xdrug and xshock, and we need to handle them all, so case_when is better than ifelse:
task %>%
mutate(xdrug = case_when(xdrug == "Drug" ~ "drug_yes",
xdrug == "NoDg" ~ "drug_no",
xdrug == "NoDrug" ~ "drug_no"
)) %>%
mutate(xshock = case_when(xshock == "NoS" ~ "shock_no",
xshock == "Shk" ~ "shock_yes",
xshock == "S" ~ "shock_yes")) -> task
taskcase_when always looks better if you take five seconds to line things up. If you don’t, it makes your reader wonder what else you are sloppy about.
I gave the new categories names starting with drug_ and shock_ rather than just yes and no, so as to be less confusing for you later. Let’s check that we have something sensible now:
task %>% count(xdrug, xshock)There were supposed to be 16 subjects in each subcategory, so it looks good. This is what I saved for you.
With one quantitative and two categorical variables, a grouped boxplot is the standard one:
ggplot(task, aes(x = xdrug, y = response, fill = xshock)) + geom_boxplot()
The other categorical variable, xshock, goes as fill and is distinguished by colour. Using colour instead of fill will outline the boxes with the appropriate colour, and is equally good.
You might already be suspecting an interaction, since there appears to be an effect of electroshock for the subjects that did not get the drug, but not for the subjects that did get the drug. Seeing the effect of one factor depending on the level of the other is a hallmark of an interaction.
I asked you to put xdrug on the \(x\)-axis just to give the grader something more uniform to assess. In actual fact, you could just as well have xshock on the \(x\)-axis and xdrug as fill.
Again, put xdrug on the \(x\)-axis. For this graph, you need the mean value of response for each combination of xdrug and xshock, which is a group-by and summarize. Then plot the means as points joined by lines for each value of xshock, which is done by setting xshock as both colour and group:
task %>% group_by(xdrug, xshock) %>%
summarize(mean_res = mean(response)) %>%
ggplot(aes(x = xdrug, y = mean_res, colour = xshock, group = xshock)) + geom_point() + geom_line()`summarise()` has grouped output by 'xdrug'. You can override using the
`.groups` argument.
One point for finding the means, two for making the graph from them (the latter is the main thing here). Other approaches:
d say, and then use that dataframe as the first input to ggplotstat_summary in lecture, but it’s in the slides, and if you can make it work to get this graph, I’m good with that.The two lines on the interaction plot do not seem at all parallel, so we should expect to see an interaction. To decide whether it is going to be significant, remember that the interaction plot does not show variability, and to get a sense of how much variability there is, we need to go back to the boxplot. The two boxes on the left look definitely different and the two on the right definitely do not, so we expect to see significance somewhere. (To be more precise, an interaction means that the effect of xshock depends on the level of xdrug, and what we just said means that it definitely does appear to.)
Include an interaction as the first step of any two-way ANOVA, and see if it is significant:
task.1 <- aov(response ~ xdrug * xshock, data = task)
summary(task.1) Df Sum Sq Mean Sq F value Pr(>F)
xdrug 1 1.00 1.00 0.568 0.45398
xshock 1 39.06 39.06 22.189 1.51e-05 ***
xdrug:xshock 1 18.06 18.06 10.260 0.00218 **
Residuals 60 105.62 1.76
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1With a P-value of 0.002, the interaction effect is definitely significant, confirming what we saw from the graphs.
Extra: The strong significance of xshock is not relevant; it is actually significant because subjects who had the electroshock treatment got fewer of the tasks done on average than those who did not, but by studying the interaction with simple effects (next), we get a more nuanced conclusion than that.
xshock for each level of xdrug, stating your conclusions in the context of the data.A simple effects analysis is appropriate because the interaction is significant, and this will help us to understand why it was significant. One point.
Two points for each of the simple effects analyses, one for each level of xdrug. In each case, extract the level of xdrug you care about, run an aov predicting response from xshock, and display the results. You might need to find out the two levels of xdrug: you can either go back to the data and look, or you can count them:
task %>% count(xdrug)Then do the two analyses, for example like this:
task %>%
filter(xdrug == "drug_no") %>%
aov(response ~ xshock, data = .) %>%
summary() Df Sum Sq Mean Sq F value Pr(>F)
xshock 1 55.12 55.12 39.61 6.15e-07 ***
Residuals 30 41.75 1.39
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1For the subjects who did not get the drug, there is a significant effect of xshock on response (P-value 0.0000006), and you can check the plots to see what kind of effect it was: if a subject did not get the drug, they did significantly more tasks if they didn’t get electroshock either.
And then the other one:
task %>%
filter(xdrug == "drug_yes") %>%
aov(response ~ xshock, data = .) %>%
summary() Df Sum Sq Mean Sq F value Pr(>F)
xshock 1 2.00 2.000 0.939 0.34
Residuals 30 63.87 2.129 With a P-value of 0.34, this is not significant, so that if a subject did get the drug, then electroshock (or not) made no difference to the number of tasks they completed.
These conclusions, you note, are entirely what you would have guessed from the graphs, but now you have P-values to support them with. Two points for each filter, aov, summary, and conclusion in words in the context of the data.
I mentioned the effect of xshock in our first aov, the one with the interaction in it. Now you see that inferring an overall effect of xshock would have been a mistake: yes, subjects who received the electroshock treatment completed fewer tasks on average overall, but the reason for that was that the subjects who received electroshock completed a lot fewer tasks if they did not get the drug, and it took the simple effects analysis to enable us to see that. Hence, trying to interpret a main effect when the interaction is significant is a mistake, and doing so will likely mislead you.
Comments: I did these in one pipeline each, because that seemed to be simplest, and I knew I was only doing one thing each time (running aov and displaying the results). There is no reason to do Tukey here (see next part); if there was, it would be better to save the result of the aov each time, look at the summary of it, and then if significant look at the TukeyHSD of it. There are a couple of technical details in the way I did it, both concerning the fact that in a pipeline there is “the thing that came out of the previous step” that is taken as the first input of the next step. In aov, what came out of the filter is going to be the dataframe used in data, which is not (by default) first, so I used the special . notation to indicate “here is where the result of the previous step goes”. In the next line, the output of aov is the input to summary, which is how a pipeline works, so I don’t have to say what it is called (it doesn’t have a name beyond . anyway).
For you, feel free to save the result of the filter and/or the result of the aov and use them as inputs to the next thing. As long as you get the same two ANOVA tables as I did by what looks like sensible code, I am happy.
Extra 1: As far as I can tell, there is no particular reason for this dataset to condition on xdrug and look for an effect of xshock. It would be perfectly good to do it the other way around. In the auto noise example in lecture, we cared more about the effect of filter type and the engine size was like a “block” (in experimental design terms, a factor that makes a difference but that we don’t really care about; putting it in the analysis will reduce the error sum of squares and thus give us better tests for the things we really do care about). In that context, it makes more sense to condition on the block and say “in this block, this happened”. Here, though, you can imagine that the researchers were interested in both factors, and so there was no special reason to do it this way around. I asked you to do it like this to be consistent with your graphs (which I was also more specific with than I might have been), so that you could find out why your significant simple effect was significant by looking straight at your graphs. I invite you to switch everything around: make xshock the \(x\) on your graphs, and condition on xshock for your simple effects, and see what difference it makes to your conclusions.
Extra 2: I did some digging and found out that the model formula input of aov is called formula. It seems that if you specify that by name, it will take data = from the output of the previous line, and so you no longer need the data = .:
task %>%
filter(xdrug == "drug_no") %>%
aov(formula = response ~ xshock) %>%
summary() Df Sum Sq Mean Sq F value Pr(>F)
xshock 1 55.12 55.12 39.61 6.15e-07 ***
Residuals 30 41.75 1.39
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The logic, I think, is that the first two inputs to aov are formula and data in that order; since we specified formula by name, the first remaining input to aov is now data, so that’s where the output from the filter goes.
Because the xshock factor only has two levels, so if you find a significant difference due to xshock, you already know that its two levels are different in terms of response without any further analysis being needed.
I suppose the bit up to the first comma is all you need, but it is better to say something like the rest of it as well, to show that you understand.
Extra: if you do it anyway (as you can find out for yourself easily enough if you are not sure), this is what you get:
task %>%
filter(xdrug == "drug_no") %>%
aov(response ~ xshock, data = .) %>%
TukeyHSD() Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = response ~ xshock, data = .)
$xshock
diff lwr upr p adj
shock_yes-shock_no -2.625 -3.476797 -1.773203 6e-07you get exactly the same P-value as you got before, so you have really gained nothing. You get a confidence interval for the difference in means, sure, but even that tells you not much more than your graph does (for subjects that did not receive the drug, the number of tasks completed was less for the subjects that received electroshock).