This worksheet is for extra practice on the last week of lecture material. It is not attached to a tutorial, and thus my solutions will be available immediately (but you are still encouraged to try this worksheet first, and then look at my solutions to see how you did).
Packages
library(tidyverse)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.4 ✔ readr 2.1.6
✔ forcats 1.0.1 ✔ stringr 1.6.0
✔ ggplot2 4.0.1 ✔ tibble 3.3.0
✔ lubridate 1.9.4 ✔ tidyr 1.3.1
✔ purrr 1.2.0
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(broom)
I use broom later, but you may not need to.
Writing a function to do wide \(t\)-test
The way we know how to run a two-sample \(t\)-test is to arrange the data “long”: that is, to have one column with all the data in it, and a second column saying which of the two groups each data value belongs to. However, sometimes we get data in two separate columns, and we would like to make a function that will run the \(t\)-test on data in this form.
As an example, suppose we have data on student test scores for some students who took a course online (asynchronously), and some other students who took the same course in-person. The students who took the course online scored 32, 37, 35, 28; the students who took the course in-person scored 35, 31, 29, 25. (There were actually a lot more students than this, but these will do to test with.)
Enter these data into two vectors called online and classroom respectively.
If you are careful you can copy and paste my values rather than typing them in again.
\(\blacksquare\)
Describe what happens when you pass each of your two vectors into enframe.
Try it and see, then talk about what you got:
online_df <-enframe(online)online_df
classroom_df <-enframe(classroom)classroom_df
When you feed enframe a vector, it creates a dataframe with two columns. The first is called name and just numbers the observations; the second is called value and contains the actual values that were in the original vector. So in these cases, we get two four-row dataframes (because each of the vectors had four observations in it).
I gave these names because we are going to use them again in a moment.
Extra: enframe also works with “named vectors”. You might have seen this in the context of quantile, which by default gives you a five-number summary with the percentiles labelled as the names of the result:
q <-quantile(classroom)q
0% 25% 50% 75% 100%
25 28 30 32 35
enframe(q)
If your vector has names, enframe puts them in the name column; if it does not (as in our case above), it numbers the observations.
\(\blacksquare\)
The function bind_rows glues together two (or more) dataframes (most commonly, with the same column names), one above the other. Its inputs are as many dataframes as you want to give it, and in addition an input .id which is the name of a column that will identify which dataframe each row came from. Use all of this to glue together the two dataframes you made with enframe.
As you read this question, you’ll see why it was sensible to give the dataframes names. I’m going to use the ones I created:
d <-bind_rows(online_df, classroom_df, .id ="group")d
Choose your own name for the column identifying which dataframe each value came from. These are teaching methods, but since we are going to be writing a function with these ideas, it’s better to use a more “generic” name like group.
I named this one as well.
Extra: Note that the groups are numbered, but that the group labels are actually text, so that something like making a boxplot out of this dataframe will in fact work just fine:
ggplot(d, aes(x = group, y = value)) +geom_boxplot()
Contrast this with the two-sample one about three- and four-bedroom houses, where the number of bedrooms was actually a number, and so we had to do some extra work to get the boxplot to come out right.
\(\blacksquare\)
Explain briefly why the dataframe you just created is suitable for running a two-sample \(t\)-test on, and run a two-sided (Welch) two-sample \(t\)-test.
The work we just did to put the two input vectors one above each other means that we now have all the data values in one column, with a separate column indicating which group they belong to (the one I called group). The column called name is not of relevance to us, and we just ignore it.
Hence, the layout is now exactly right for running t.test on:
t.test(value ~ group, data = d)
Welch Two Sample t-test
data: value by group
t = 1.0498, df = 5.9776, p-value = 0.3344
alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
95 percent confidence interval:
-3.998992 9.998992
sample estimates:
mean in group 1 mean in group 2
33 30
\(\blacksquare\)
Write a function that takes two vectors as input. Call them x and y. The function should run a two-sample (two-sided, Welch) \(t\)-test on the two vectors as input and return the output from the \(t\)-test.
Solution
The reason I had you do the previous questions was so that you have written everything that you need to put into your function. (This is a good strategy: get it working in ordinary code first. Then, when you write your function, there’s less to go wrong.)
t_test_wide <-function(x, y) { d1 <-enframe(x) d2 <-enframe(y) d <-bind_rows(d1, d2, .id ="group")t.test(value ~ group, data = d)}
The function you are writing will apply to any two columns of input, not just test scores in two different classes, and future-you will appreciate knowing the roles of things in your function, ready for when future-you needs to change something about how the function works. Hence, it’s a good idea to give them generic names. Make sure that your inputs to t.test reflect the column names that are in what I called d (that is, probably value that came from enframe and whatever you called the input to .id).
R style, and thus best practice here, is to have the last line of the function be the value that is returned from the function (here, the output from t.test, all of it). There is no need to wrap it in return; it works (unlike in Python) to do it without, and is better style in R.
Extra: Also, if you are thinking about Python, you might be wondering about the indentation. Python requires the indentation, but R doesn’t (because of the curly brackets around the function). Having said that, though, it’s good style to do the indentation as I did it (and also benefits future-you when you come back to your code in six months and have forgotten how you did it): indent each line of the function once (two spaces, which is also the default TAB in R Studio). If you have any pipelines in your code, indent them again so that future-you can see what belongs to the pipeline and what does not.
As an example of my last point, you could have written your function like this:
This time I indented the t.test line two more spaces to show that it belongs to the pipeline that starts with bind_rows. The data = . inside the t.test means “use the output dataframe from bind_rows as the input dataframe to t.test”.
I actually don’t like the second way as much: it might be more “efficient” in that you don’t create temporary variables, but it tries to do everything at once and it’s much more likely that the reader (which might be future-you) will get confused by what is going on. I favour having one line do one thing at a time: in t_test_wide, it’s very clear that I am doing this:
turn input vector x into a dataframe
turn input vector y into a dataframe
glue those dataframes together one above the other
run the t-test on that dataframe
with one line of code corresponding to each of those things. If you find your eyes glazing over as you read some code, it is likely that the coder is trying to be clever rather than trying to be clear.
\(\blacksquare\)
Test your function on the same data you used earlier, and verify that you get the same P-value.
Solution
Use classroom and online as the two inputs to the function. They can be either way around, since the test is two-sided:
t_test_wide(online, classroom)
Welch Two Sample t-test
data: value by group
t = 1.0498, df = 5.9776, p-value = 0.3344
alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
95 percent confidence interval:
-3.998992 9.998992
sample estimates:
mean in group 1 mean in group 2
33 30
My second version also works:
t_test_wide_2(online, classroom)
Welch Two Sample t-test
data: value by group
t = 1.0498, df = 5.9776, p-value = 0.3344
alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
95 percent confidence interval:
-3.998992 9.998992
sample estimates:
mean in group 1 mean in group 2
33 30
The same P-value indeed, and everything else the same also.
If the P-value is not the same (or you get an error from your function), you will need to go back and fix up your function until it works. Debugging is a big part of coding. Most of us don’t get things right the first time. Within limits,1 it doesn’t matter how long it takes, as long as you get there by the time your work is due.
\(\blacksquare\)
Modify your function to return just the P-value from the \(t\)-test, as a number. Hint: you will need to find out how to get the P-value out of the t.test output. You might find this page useful.
Solution
This means finding out how to get just the P-value out of the t.test output. It’s there on the output, but how do we get it out?
The page I linked to in the hint actually tells you all you need to know:
it mentions t.test in the very first paragraph, which should motivate you to keep reading, even though we have only used broom for regressions so far.
look for the example of tidy on a regression problem. The paragraph below that tells you how to get things out of the tidy output.
the section on Hypothesis Testing shows you how tidy works for a two-sample \(t\)-test, and demonstrates that the column you’ll want is called p.value.
Putting those together says that what you need to do is this:
save the \(t\)-test output
run tidy on it, and save that
use $ to pull out the P-value.
When we were playing around before, we made a dataframe d like this:
d
That leads to this (displaying results for checking):
t_d <-t.test(value ~ group, data = d)t_d
Welch Two Sample t-test
data: value by group
t = 1.0498, df = 5.9776, p-value = 0.3344
alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
95 percent confidence interval:
-3.998992 9.998992
sample estimates:
mean in group 1 mean in group 2
33 30
t_d_tidy <-tidy(t_d)t_d_tidy
t_d_tidy$p.value
[1] 0.3343957
and that can go into our function, minus the displaying of results, with the last thing being what is returned from it:
Extra: another approach, which turns out to be slightly simpler (but requires more thinking to get there) is to find out all the things that are output from t.test. To do that, run your initial t.test again, but once again save the result:
The third of those is the one we want. Get it out with a dollar sign, as if it were a column of a dataframe:2
t_d$p.value
[1] 0.3343957
This, I think, also works (in a rather more tidyverse fashion):3
t_d %>%pluck("p.value")
[1] 0.3343957
This function names has quietly shown up a couple of times in this course (once in combination with power.t.test and once when I showed you what augment did in regression). That leads to this kind of function:
This is one line shorter than before, because it avoids the use of tidy, but it meant that we had to find out that t.test has an item p.value that we can access. The output from tidy is (by design) more predictable.
Yours might be different again from those, but if it achieves the same kind of thing and it works, I am good with it.
The page I linked to in the hint is a good introduction to the thinking behind broom. I found it by searching for “r package broom”.
\(\blacksquare\)
Test your modified function and demonstrate that it does indeed return only the P-value (and not the rest of the \(t\)-test output).
Solution
I have two functions, so I have two tests (but you will only have one):
t_test_wide(online, classroom)
[1] 0.3343957
t_test_wide_names(online, classroom)
[1] 0.3343957
Success.
\(\blacksquare\)
Modify your function to allow any inputs that t.test accepts. Demonstrate that your modified function works by obtaining a pooled \(t\)-test for the test score data, with a one-sided alternative that the online students had a higher mean mark than the classroom students.
Solution
This is the business about passing input arguments on elsewhere without having to intercept them one by one, using ...:
t_test_wide <-function(x, y, ...) { d1 <-enframe(x) d2 <-enframe(y) d <-bind_rows(d1, d2, .id ="group") t_d <-t.test(value ~ group, data = d, ...) t_d_tidy <-tidy(t_d) t_d_tidy$p.value}
Put the ... at the end of the inputs to the function, and again at the end of the inputs to t.test (since that’s where anything unrecognized needs to be sent).
To test: when I constructed the dataframe with all the scores in it, I put the first input first, so if I input the online students first, my alternative hypothesis is that their mean is greater than for the classroom students, so the alternative will need to be greater:
t_test_wide(online, classroom, alternative ="greater",var.equal =TRUE)
[1] 0.1671256
Our function does not have inputs alternative or var.equal, so they get passed on to t.test which does.
The P-value is less than it was before, which says that if we had an a priori reason (before looking at the data) to suspect that the online students would score higher if anything, the evidence of this kind of difference is stronger than the two-sided test indicates.
If your two input vectors are the other way around, your alternative will need to be less, because the first input to bind_rows is the one that gets numbered 1, and thus is the “first” group as far as alternative is concerned:
Extra: the useful comparison is really with a one-sided Welch test. The Welch test is the default, obtained by not passing on a value for var.equal:
t_test_wide(online, classroom, alternative ="greater")
[1] 0.1671979
and this says that there is extremely little difference from the pooled test, meaning that the spreads were not actually unequal enough to worry about.
Code-wise, it doesn’t matter how many inputs are actually caught by the ...: anything that is not an input to your function will be passed on, no matter how many things that is.
Extra extra: the ... mechanism is how tidyverse functions like select or summarize work that accept any number of inputs. The bind_rows function we used earlier is structured in the same way:
All of its “own” inputs have a dot on the front of their names (only .id here, which you might remember is the name of the column that distinguishes which of the original dataframes a row comes from).
\(\blacksquare\)
Footnotes
There is never actually an infinite amount of time to do something, so you do have to know enough to get reasonably close reasonably quickly. The time limit on assignments is there for that reason.↩︎
It is actually an element of a list, but dataframes at their core are actually lists of columns, so the extraction process is the same.↩︎
pluck pulls anything out of anything else; it doesn’t have to be a dataframe.↩︎
