STAC32 Assignment 1

Packages

library(tidyverse)

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.

198 children

Body measurements for 198 children were sampled, and can be found in http://datafiles.ritsokiguess.site/kids198.csv. There are four columns: Height (inches), Weight (pounds), Age (months), Sex (identify as female or male).

(2 points) Read the data from the file into a dataframe, and display (at least some of) the dataframe.

This is a .csv file, as you can tell from the filename, so read_csv will read it in. These are straightforward because you don’t need to say anything beyond “it’s a .csv file”:

my_url <- "http://datafiles.ritsokiguess.site/kids198.csv"
children <- read_csv(my_url)

Rows: 198 Columns: 4
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): Sex
dbl (3): Height, Weight, Age

ℹ 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.

The output from read_csv tells you how many rows and columns there are, so that you can verify that they are as expected (something you should always do).

Give your dataframe a name that says what it contains. children or kids is fine; so is body_measurements with the words in the variable name separated by an underscore so that your reader can easily see what it is. Something like my_data is a bad name, because both you and your reader will soon forget what the data represent.

If you instead use the old-fashioned read.csv with a dot, you have revealed that you are not paying attention to how we do things in class, and you should expect to lose points.

To view (some of) the dataframe, this is simplest:

children

If you have set things up right (from “Assignment 0”) your reader will be able to page through the dataframe to check that things are as expected.

Another good approach that you might have seen is this one:

glimpse(children)

Rows: 198
Columns: 4
$ Height <dbl> 67.8, 63.0, 50.1, 55.7, 63.2, 48.8, 63.8, 61.3, 61.1, 54.7, 68.…
$ Weight <dbl> 166, 93, 54, 69, 115, 52, 108, 89, 118, 80, 139, 129, 55, 113, …
$ Age    <dbl> 210, 144, 119, 130, 157, 102, 198, 155, 199, 134, 186, 162, 122…
$ Sex    <chr> "male", "female", "male", "female", "male", "male", "female", "…

which displays the columns going down and the rows going across (that is, “transposed” from the usual view), but it equally well enables your reader to check that the values you read in make sense.

Extra: If you wish, you can verify that it really is a .csv file by downloading the file and taking a look at it. It will probably open in your spreadsheet application, which suggests where it came from, or if you are able to view the file itself, you will see that the top of it looks like this:

Height,Weight,Age,Sex
67.8,166,210,male
63,93,144,female
50.1,54,119,male
55.7,69,130,female
63.2,115,157,male
48.8,52,102,male
63.8,108,198,female

that is to say, literally values separated by commas.

(2 points) How do you know that your dataframe has the correct number of observations and variables?

There should be one row for each of the 198 observations and one column for each of the 4 variables.

To see whether that’s what we have: read the output from read_csv (which directly tells you how many rows and columns you have and whether they are text, as Sex is, or numerical, as the others are). Or, when you look at your dataframe by typing its name, the output tells you how many rows there are and you can see there are four columns. Or, if you used glimpse, you get the number of rows and columns directly at the top of the output.

Make the connection between rows and observations, and between columns and variables, and then say how many of each of them you have (and that they are correct from the information given in the question).

(3 points) Make a suitable plot of age and height.

These are both quantitative variables, so the right plot is a scatterplot. Think about which is an outcome or response variable: it seems that height depends on age rather than the other way around, so put height on the $y$-axis. Don’t forget that the variable names have Capital Letters:

ggplot(children, aes(x = Age, y = Height)) + geom_point()

Extra: it is better not to add a regression line to this plot. Here’s what happens if you do:

ggplot(children, aes(x = Age, y = Height)) + 
  geom_point() + geom_smooth(method = "lm")

`geom_smooth()` using formula = 'y ~ x'

What you see here is that the points:

on the right are mostly below the line (between ages af 200 to 225 months)
just to the left of those (175–200 months) are mostly above the line
on the left end (100–125 months) are mostly below the line.

That is, the trend is not linear.¹

You might have been wondering what happens if you use geom_smooth without the method = "lm". The answer to that is “try it and see”:

ggplot(children, aes(x = Age, y = Height)) + 
  geom_point() + geom_smooth()

`geom_smooth()` using method = 'loess' and formula = 'y ~ x'

This adds a “smooth trend” to the plot, one that is a smoothed-out version of the data, so that you can see what kind of trend you have: in this case, a linear upward trend until age about 175 months, then it levels off after that.²

I haven’t shown you the plain geom_smooth in lecture yet (at least, not as I write this), but you might have experimented with it, and it is a reasonable thing to ask about, or find out about yourself (citing your source, such as this one. Scroll down to the “six common patterns” and click the Run button for each one.)

(2 points) Does the trend that you see on your graph make practical sense? Explain briefly.

The trend is upward, saying that older children are taller. This makes sense because children grow, at least until a certain point, so we would have expected older children to be taller. The “certain point” appears to be something like 180 months (15 years).

The best answer says that children get taller as they get older, but the growth levels off at a certain age (and say something about what that age is, as far as you can judge it from your graph).

Extra: the fact that height levels off is the reason why a linear trend is not a good thing to add to your graph: a linear trend keeps on going up, and the data points don’t.

(2 points) Make a suitable plot of age, height, and sex.

Sex is categorical, so to your previous scatterplot add Sex as colour:

ggplot(children, aes(x = Age, y = Height, colour = Sex)) + 
  geom_point()

Extra: if you found out about geom_smooth, you can add that here too. What happens this time is that you get a separate smooth trend for each Sex (the thing in colour):

ggplot(children, aes(x = Age, y = Height, colour = Sex)) + 
  geom_point() + geom_smooth()

`geom_smooth()` using method = 'loess' and formula = 'y ~ x'

(2 points) What do you learn from your plot of Question 5 that you did not see from your plot of Question 3? Explain briefly.

This must be something to do with differences between sexes.

The plot of Question 5 shows that boys end up taller than girls (on the right, the blue dots are mostly above the red ones). However, it does not appear to be the case that boys are taller than girls at all ages, because the red and blue dots at ages less than about 175 months (about 15 years) are all mixed up. What appears to happen is that the girls stop growing taller sooner than the boys do: at about 160 months for girls, at about 175 months for boys.

Extra: This is perhaps clearer on the graph with the smooth trends: height grows with age in about the same way for boys and girls up to about 160 months, at which point the girls stop growing taller but the boys keep going. I’m not sure I agree with the smooth trend for the boys, although it is actually true that almost all of the boys of age over 200 months are about 70 inches tall, and a fair number of the boys between 175 and 200 months are shorter than that.

Appendectomies

The appendix is a small organ in the lower right abdomen. It helps balance good and bad bacteria in the gut and supports the immune system. Source. The most common problem affecting the appendix is appendicitis. This occurs when an infection causes swelling and irritation. Appendicitis is a medical emergency, and immediate surgery is needed to remove the appendix (an appendectomy).

A 1992 study investigated the cost of a “typical” appendectomy at each of 59 hospitals in North Carolina. Each hospital was also classified as one of “metropolitan”, “regional”, or “rural”; this classification is in the variable region. The data are in the file http://datafiles.ritsokiguess.site/appendectomy.txt.

(1 point) Use your web browser to examine the data file. What separates one data value from the next?

This is what I see (or some of it, at least):

The cost and the region are separated by a colon (:).

I showed you a screenshot to convince you that you and I are looking at the same thing; you don’t need that. All you need is to say “the data values are separated by a colon”, or something equivalent to that.

(2 points) Read in the data and store the values in a suitably-named dataframe. Display at least some of your dataframe.

This one requires a little bit of thought. This is not a .csv and the data values are not separated by spaces, so we need to use read_delim and tell it what the data values are separated by:

my_url <- "http://datafiles.ritsokiguess.site/appendectomy.txt"
appendectomy <- read_delim(my_url, ":")

Rows: 59 Columns: 2
── Column specification ────────────────────────────────────────────────────────
Delimiter: ":"
chr (1): region
dbl (1): cost

ℹ 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.

appendectomy

The name appendectomy is a long one to type (but see the Extra below), so I have no objection if you use something shorter that still gets the message across. I would be fine with app as a minimum.

There are, as you should verify for yourself, 59 rows (hospitals) and 2 columns (variables). If you mess it up (eg., try to use read_csv), you will probably get only one column with all the data as text, and the fact that there is only one column is your warning sign that there is a better way to do it.

This way also works, though you will need to think about why:

appendectomy <- read_delim(my_url)

Rows: 59 Columns: 2
── Column specification ────────────────────────────────────────────────────────
Delimiter: ":"
chr (1): region
dbl (1): cost

ℹ 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.

What happens here is that read_delim guesses what is separating the data values by looking at the first few lines of the file. You can tell what it guessed by looking at the Delimiter: in the message output by read_delim; in this case, it guessed “colon” and it was right.

If you do it this way, you will need to say why it still worked even though you didn’t specify the delimiter. Something like “read_delim figured out that the data values were separated by a colon, as shown in the Delimiter: line of the read_delim output”.

Extra: you don’t have to shy away from long descriptive names for your variables. Start typing the name of your variable (the default is a minimum of three characters), pause, and R Studio will bring up a completion box. Use the up and down arrows to select the thing you want.

(2 points) Make a suitable graph of the two variables in your dataframe.

One quantitative variable cost and one categorical variable region, so a boxplot is called for:

ggplot(appendectomy, aes(x = region, y = cost)) + geom_boxplot()

In this course, the boxes on a boxplot go up and down, so the x and y need to be this way around. (Later in the course, you can exercise a bit of artistic judgment, but even then, you will need to have a reason for your boxes to go across rather than up and down.)

Extra: one valid reason would be that the labels on the $x$-axis are so long that they overwrite each other. Here, if we had (say) six regions with names as long as metropolitan, that might be in danger of happening. In that case, what you can do is to interchange the role of x and y:

ggplot(appendectomy, aes(y = region, x = cost)) + geom_boxplot()

There is typically more room horizontally on the $y$ axis for labels than horizontally on the $x$ axis (imagine if you had 20 groups instead of 3). On the other hand, the groups are usually explanatory (the region a hospital is in might influence the cost of an appendectomy, rather than the other way around) so, other things equal, the categorical variable is logically best put on the $x$-axis, unless you have a good reason for moving it.

(2 points) What does your plot tell you about how the cost of an appendectomy compares across the three regions on average?

The median cost (of a “typical” appendectomy) is almost the same across the three regions, or the median cost is slightly higher in rural hospitals (with metropolitan and regional hospitals being about the same on average).

Make sure it is clear that:

you know the lines across the boxes on a boxplot are the median and not the mean. “Average” is the umbrella term for any of the measures of centre such as mean, median, or even the mode.
if you want to say that the cost is higher in rural hospitals on average, make sure you qualify it with a word like “slightly” to make it clear that any difference is small. A hospital administrator or someone in medical insurance (the target audience for an analysis of these data) is unlikely to be much interested in very small differences in median like these.

As I asked the question, I am not interested in spread (the height of the boxes) or in any outliers (such as the rural hospital with very expensive appendectomies). So you don’t need to talk about those.

Allergy meds

In people who get seasonal allergies, various medications are popular, but often these medications come with side-effects such as insomnia or drowsiness. In a study, allergy sufferers were randomly allocated to one of three allergy medications: Seldane-D, Pseudoephedrine or a placebo, and were asked to note the first “adverse event” (side-effect) that they noted (insomnia, drowsiness, headache). The data are in http://datafiles.ritsokiguess.site/allergy_meds.txt, with the column event containing the side-effect that was experienced by that person.

(2 points) Examine the data and appropriately read the data into a dataframe.

By “examine” I mean something like “look at the data file in a web browser”. If you do that you’ll see that you have aligned columns, with a variable number of spaces between them. This means that read_table is the thing:

my_url <- "http://datafiles.ritsokiguess.site/allergy_meds.txt"
allergy <- read_table(my_url)


── Column specification ────────────────────────────────────────────────────────
cols(
  event = col_character(),
  medication = col_character()
)

allergy

There were evidently 406 people in the survey, and you see columns of adverse events and of medications. You might have to scroll down a bit to see a different event or medication, but they are there if you go down far enough.

Extra 1: as a general principle, when the name of the data file ends in something other than .csv (or .xlsx or another format that you know how to deal with), it pays to take a look at the data first to get a sense of what you are up against.

Extra 2: read_table reads in data values separated by any amount of whitespace (spaces or tabs), while read_delim reads in data values separated by any character, but only one of that character. If you have something like this:

x:y
1:11
2::10
3:13
4::13

where the data values are separated by colons, but sometimes one colon and sometimes two, you are, at this point in the course, stuck. Here’s what happens if you try to use read_delim:³

colonny_text <- "
x:y
1:11
2::10
3:13
4::13
"
test_data <- read_delim(I(colonny_text), ":")

Warning: One or more parsing issues, call `problems()` on your data frame for details,
e.g.:
  dat <- vroom(...)
  problems(dat)

Rows: 4 Columns: 2
── Column specification ────────────────────────────────────────────────────────
Delimiter: ":"
chr (1): y
dbl (1): x

ℹ 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.

test_data

read_delim infers that there are only two columns, but some of the values in y acquired an extra colon on the front, and the whole column is now text rather than the numbers it should be. In the warning above the output, “parsing issues” means that read_delim had trouble figuring out how many columns there should be:

problems(test_data)

There were two header columns, so read_delim assumed (reasonably) that there would be two columns all the way down. When it found the rows 2::10 and 4::13 with two colons in them, it didn’t know whether you meant there to be two or three columns, and it made the decision that you meant there to be two columns, and the extra colons were actually part of the data.

As an extra to the extra: what you end up doing in this case is to read the data initially as one column (using, say, read_table because there is no whitespace), and then process it afterwards:

test_data <- read_table(I(colonny_text))
test_data

and then (this won’t make much sense to you yet):⁴

test_data %>% 
  separate_wider_regex(`x:y`, c(x = "[^:]*", ":+", y = "[^:]*"))

and then there is a final round of processing to turn those values, which are still text, back into numbers.

Moral: having your data be in a predictable format makes your life much easier!

(3 points) Create a suitable graph of the two variables in this dataset. Hint: look ahead to the next question to see how you might want to arrange things to make that question easier to answer.

The two variables are both categorical, so a good plot is a grouped bar chart. This turns out to be the best way to arrange the graph:

ggplot(allergy, aes(x = medication, fill = event)) + geom_bar(position = "dodge")

The reason this is best is that we want to look within each medication to see which event is most common, so we want to have the three bars for each medication side by side, which is done by having medication as x and the other categorical variable as fill. You might like to think of it this way: the medication is explanatory and the event is response (outcome), in that the medication is causing the event rather than the other way around. When you have a response, that goes best as the fill.

If you have x and fill the other way around:

ggplot(allergy, aes(fill = medication, x = event)) + geom_bar(position = "dodge")

now you have an answer to a question like “out of the people that got a headache, which was the most common medication for them to be on?” This is a (somewhat) reasonable question to want the answer to, but it is not what the next question is asking.

So, the general principle here is that if you have something that is an outcome (response), put that in fill and put the other (explanatory) variable in x.

This is a new idea compared to the grouped boxplot in the lecture notes; there, the categorical variable Sport had a lot of levels (10 different sports), so it was easier to put that on the $x$-axis and have Sex as fill. We ended up in a similar place, though: there, we could say “for each Sport (x), how many athletes were there of each Sex (fill). Another way to think of that is Sport being explanatory and Sex being response, or as”if I know an athlete’s sport, what does that tell me about whether they are male or female?“. In the case of the athletes, knowing the sport sometimes told you a lot about the athlete’s Sex: for example, if the athlete was a netball or water polo player, you knew for certain that the athlete was female or male (respectively).

I might have shown you in lecture what happens to that plot if you interchange x and fill.

Sometimes it makes more sense to put the categorical variable with more levels on x, and sometimes you have an obvious response variable that belongs in fill, and these can be at odds with each other. You will then need to use your judgment to find a good compromise, bearing in mind your aims with the graph.

(2 points) For each medication, what is the most common side-effect?

Find the tallest bar in each set of three: for the people taking the placebo, the most common side-effect is a headache, but for people taking the other two (real) medications, the most common side-effect is insomnia.

If you did your graph the other way around, with event as x, you have a bit more work to do: you need to look for the tallest red bar (headache) and the tallest green and blue bars (insomnia). You can see that this is harder work than it needs to be; there is value in making choices about your graph that will make your consequent work easier.

Extra: with two categorical variables, you may be familiar with a chi-squared test. That goes like this:

ans <- with(allergy, chisq.test(event, medication))
ans


    Pearson's Chi-squared test

data:  event and medication
X-squared = 37.832, df = 4, p-value = 1.213e-07

(or the equivalent with allergy$event and allergy$medication and without the with.)

This tells you of a significant association: that is, if you know which medication a person was on, that tells you something about which side effect they were likely to have.

It does not, however, tell you where the association is coming from. A way to investigate that is to look at the residuals:

ans$residuals

            medication
event           placebo pseudoephedrine  seldane-d
  drowsiness  2.9079567      -1.6572146 -0.4298012
  headache    2.4234563      -0.6982926 -0.9458606
  insomnia   -3.9189250       1.6060289  1.1191404

This tells you how the observed frequency compares to the expected one (under the null hypothesis that there is no association). When the residual is clearly positive, the observed frequency is bigger than you would expect, and when it’s negative, the observed frequency is smaller than you would expect.

In this case, it seems that insomnia points one way, and the other two events point the other way: a person on placebo is less likely to suffer insomnia and more likely to get one of the other two events, compared to expected, while for a person on either of the two real medications, the opposite is true: insomnia is more likely, and the other two events are less likely, compared to expected.

Another approach is to note that the three largest (in size) residuals are all for placebo, so this medication is different from the two real ones (which, to judge by the residuals, are fairly similar to each other in terms of likely adverse event).

You can decide for yourself if these messages are about the same as the one from your graph. I think the comparison of medications, at least, does make sense: the two real medications have about the same pattern of adverse events (insomnia tallest, drowsiness shortest, and by about the same amount relatively speaking), while the placebo has a much different pattern.

Technical aside: if you have done a chi-squared test by hand, you’ll remember that the procedure goes something like this:

count up the number of observations for each combination of event and medication
for each cell, work out expected frequency (assuming there is no association)
for each cell, work out observed minus expected, square it, divide by expected
add up all of the numbers from the previous step to get the test statistic
look up the test statistic in a chi-squared table to get the P-value.

The residuals here are the square roots of the numbers in the third step above: that is, they are observed minus expected, divided by the square root of expected. This is saying how far apart observed and expected are (like a regression residual), but scaled to account for larger numbers tending to be further apart.

To verify that:

sum(ans$residuals^2)

[1] 37.83248

which is the same (to within rounding) as the chi-squared test statistic above.

Footnotes

Excel will happily plot a straight line and give you an R-squared for any relationship at all, but you know better than that: if it’s not linear, a good statistician won’t add a straight line to the plot.↩︎
You can see that this does a better job of going through the data points than the regression line does.↩︎
You can, if you are careful, read data in from text instead of a file.↩︎
Separate the column called x:y wider, making more columns, specifically a column called x which is characters that are not colons, then one or more colons, then some more characters that are not colons.↩︎