STAC32 Assignment 3

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.

Maximum voluntary ventilation

Maximum Voluntary Ventilation (MVV) is the largest volume of air a person can consciously inhale and exhale over a specific time period, typically 12 to 15 seconds. This measurement is then extrapolated to a volume per minute, expressed in liters per minute (L/min). MVV requires a maximal, sustained effort, reflecting the combined strength and endurance of the respiratory muscles, including the diaphragm and the intercostal muscles. It measures the maximum rate at which the entire breathing apparatus can operate, not just lung size. Source.

It is believed that for healthy young adults, the mean MVV is 110 litres per minute. A random sample of healthy 4th-year university students is taken, with the results in http://datafiles.ritsokiguess.site/mvv.csv. There is one column rate containing the MVV values.

1. (1 point) Read in and display (some of) the data.

With only one column, the “columns” are not actually separated by anything, but the filename has a .csv on the end, so you may as well pretend it’s a CSV file:

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

Rows: 20 Columns: 1
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (1): rate

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

mvv

I called my dataframe mvv because the column in it is called something different.

If you read it in a different way, you should explain why your way will work.

Extra: These are the EXR_C07_S02_13 dataset from the DanielBiostatistics10th package and come from a well-known biostatistics textbook. The data are used as is.

Extra extra: MVV is also the name of a soccer team from Maastricht in the Netherlands. This is the current squad, though you will have to use your Dutch-language skills to work out who plays in which position. Maastricht is in the bottom right-hand corner of the Netherlands, so it is not very surprising that they have a lot of Belgian players. You’ll also note that the team’s colours resemble those of Arsenal.

2. (2 points) Draw a suitable graph of the data.

The obvious choice, with one quantitative variable, is a histogram:

ggplot(mvv, aes(x = rate)) + geom_histogram(bins = 5)

With 20 observations, this is a sensible number of bins. You might need to experiment with your number of bins to get a graph that seems to show the shape. With 6 or 7 bins, the shape looks (with 6 bins) irregular and (with 7 bins) almost flat. Such are the dangers of using a histogram with a smallish sample size. My take is that the irregular shape with 6 bins is mostly or entirely chance.

(added later) It seems that r.datatools gives (for some reason) a different histogram with the same number of bins. Here’s the 5-bin version it gives:

I have let the grader know about this, and so you should be marked fairly on this question even if your histogram is different from mine.

Alternatively, you might look ahead to what you will be using this graph for, and realize that what really matters here is normality. If that’s your thinking, the right graph is a normal quantile plot:

ggplot(mvv, aes(sample = rate)) + stat_qq() + stat_qq_line()

No interpretation needed yet, but if you want to help yourself down the road: my histogram looks more or less normal, and my normal quantile plot is slightly short-tailed but close to normal.

3. (4 points) Run a suitable \(t\)-test to determine whether there is any evidence that the mean MVV (of the population of which these students are a sample) is different from 110. What do you conclude? In your answer, state your null and alternative hypotheses, your P-value, your decision, and a conclusion in the context of the data. Use \(\alpha = 0.05\).

Letting \(\mu\) be the population mean, the null hypothesis is \(H_0: \mu = 110\). We are interested in any difference from 110, so the alternative hypothesis is two-sided: \(H_a: \mu \ne 110\). You can write your hypotheses in words if you prefer: “population mean is equal to 110”, “population mean is not equal to 110”. If you use a symbol like \(\mu\), you must say what it represents.

Now we can run the test. With a two-sided alternative, there is no need to specify an alternative in your code:

t.test(mvv$rate, mu = 110)

    One Sample t-test

data:  mvv$rate
t = 0.12709, df = 19, p-value = 0.9002
alternative hypothesis: true mean is not equal to 110
95 percent confidence interval:
  85.24933 137.95067
sample estimates:
mean of x 
    111.6 

The P-value is 0.90, much greater than 0.05, so we do not reject the null hypothesis (decision). Hence, there is no evidence that the population mean MVV is different from 110.

Extra: in this course, we distinguish between doing a test (to answer the question of whether a particular value, here 110, is believable), which is done as above, getting a P-value, and finding a confidence interval (below), which is used to answer the question of “what is the population mean?” when we do not have a particular value in mind. The reason for doing it this way is that when we do a test, we want to assess the strength of evidence against the null hypothesis as accurately as we can (with a P-value), so that our readers can decide for themselves whether they want to reject the null hypothesis or not. In this case, it is unlikely to make any difference, but if you found a P-value of, say, 0.04, you might have an \(\alpha\) of 0.05 and would reject, but your reader might have \(\alpha = 0.01\), and as long as you actually say what the P-value is, they will know what they would do.

4. (2 points) Obtain and state a 90% confidence interval for the population mean MVV.

Run t.test again, but without the mu and with a conf.level, because the confidence level is not the default 95%:

t.test(mvv$rate, conf.level = 0.90)

    One Sample t-test

data:  mvv$rate
t = 8.8643, df = 19, p-value = 3.535e-08
alternative hypothesis: true mean is not equal to 0
90 percent confidence interval:
  89.83063 133.36937
sample estimates:
mean of x 
    111.6 

The 90% confidence interval is from 89.8 to 133.4 (litres per minute). To decide how many decimal places to quote, look back at the data (whole numbers) and add 1 or at the outside 2 decimals to that. More decimals is claiming accuracy that you do not have.

Some observations about our interval:

• it contains the null-hypothesis mean 110; in fact, 110 is close to the middle of the interval. The confidence interval contains values that we think are “plausible” in the light of the data; since we think 110 is plausible, we have no business rejecting it as a mean.
• it’s rather a long interval, covering what looks like a large range of possible values of mean MVV. It seems that MVV values are rather variable from one person to the next, so it will take a larger sample size to estimate the mean with any sort of accuracy. A consequence of this is that we might not have much power to reject a population mean of 110 (if it is actually wrong). We talk about power of tests later in the course.

Extra: if you are in the habit of using confidence intervals to do tests with (a habit you would do well to get out of, as explained above), this interval only tells you that the P-value for your test is greater than 0.10 (because \(90% = 1 - 0.10\)). If your reader happens to want \(\alpha = 0.20\),¹ they will not know whether they should be rejecting or not.

5. (2 points) Do you have any concerns about the \(t\)-procedures (test and confidence interval) for these data? Explain briefly why or why not.

The assumption hiding behind the \(t\)-test is that the data values (MVV rates) have a normal distribution, or sufficiently close to it. The thing that guides the “sufficiently close” is the sample size as it applies to the Central Limit Theorem.

In this case, we have:

• data that are already close to normal (or only slightly non-normal)², from the graph you drew earlier
• a sample size of \(n = 20\), which is reasonably large in Central Limit Theorem terms.

So we should be happy to run the \(t\)-test and confidence interval, on the basis that we do have data that are sufficiently close to normal in shape. (Either the data are already normal enough, or they are close enough to normal that the sample size will take care of the short tails, if that’s what you saw.)

A third consideration here is:

• our decision to not reject the null hypothesis was very clear-cut, because the P-value was much bigger than 0.05. As a result, the normality would have to be off by a large amount to have any chance of changing our decision.

Nursing homes

A nursing home is a place where elderly or disabled people can live when they can no longer be cared for at home. The data set in http://datafiles.ritsokiguess.site/nursing-homes.csv contains information about 52 nursing homes in the US state of New Mexico. Our concern here is about nursing salaries, in the column NurseSalaries. Nursing homes can be in rural areas, labelled rural, or in cities or suburbs, labelled non-rural, in the column Location. We wish to determine whether there is evidence that nursing salaries are lower on average in the rural nursing homes. These nursing homes are a random sample of a larger number of nursing homes.

6. (1 point) Read in and display (some of) the data.

A distinct absence of surprises:

my_url <- "http://datafiles.ritsokiguess.site/nursing-homes.csv"
homes <- read_csv(my_url)

Rows: 52 Columns: 7
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): Location
dbl (6): Beds, InPatientDays, AllPatientDays, PatientRevenue, NurseSalaries,...

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

homes

nursing_homes would also be a good name.

Note for yourself that you correctly have 52 rows (nursing homes), and the columns you need: NurseSalaries and Location (you might have to scroll across to see them). NurseSalaries is quantitative, and Location is categorical.

Extra 1: these are the Nursing data from the Stat2Data package. In the original, what I called Location is a variable Rural coded as 1 (rural) and 0 (not rural). I like categorical variables to look like categorical variables, so I created a new column Location³ with text saying what the values represent, and got rid of Rural when I was done.

Extra 2: There are different ways of handling names that are multiple words. One way, as seen in NurseSalaries, is to glue the words together with uppercase letters at the start of words. This is called “camel case” because the uppercase letter(s) in the middle look like a camel’s hump(s), at least if you squint. A commoner way in R is to separate words using underscores, as in nursing_homes. This is called “snake case”, because you are making a long “snake” out of all the words. A third option can be used in things like file names, nursing-homes.csv, where the words are separated by minus signs. This is called “kebab case” because you can imagine the letters as pieces of food on a kebab stick, and the minus signs are where the kebab stick shows through.⁴

Another one you will sometimes see in older R code and data sets is a dot, like nursing.homes (or t.test) where the words are separated by dots. This works in R because a dot has no special meaning; in other languages, or in file names, the dot can be important, so you don’t see this much elsewhere.

Blind people (and others) use screen readers to read text on the screen out loud, because they can’t see it well enough to read it themselves. According to this, camel case works better than snake case for those people. Something to consider if screen reader users are part of your audience.

7. (2 points) Make a suitable graph of the variables of interest.

One quantitative variable NurseSalaries and one categorical variable Location, so a boxplot is the obvious thing:

ggplot(homes, aes(x = Location, y = NurseSalaries)) + geom_boxplot()

To note (for yourself): it does look as if the nurses in rural areas have lower salaries on average, but we don’t know yet whether that difference is real or just chance.

Another graph you might draw (thinking ahead to checking assumptions) is a normal quantile plot. This needs to have the rural and non-rural graphs separately (best, in facets) and, for the use we will be making of the normal quantile plot(s), it is better to use the same scale for each distribution:

ggplot(homes, aes(sample = NurseSalaries)) +
  stat_qq() + stat_qq_line() + facet_wrap(~ Location)

Extra: I needed this in order to write the next question:

homes %>% 
  count(Location)

8. (3 points) There are 34 rural nursing homes in the sample, and 18 non-rural ones. Using this information and your graph (from the previous question), justify a decision to run a \(t\)-test, and explain briefly which kind of \(t\)-test would be better.

1. To justify a \(t\)-test, we need sufficient normality in both groups, relative to the sample size. If you drew a boxplot:

• the salaries for the non-rural nurses are already close to normal (or, at least, symmetric) so no problem there.
• the salaries for the rural nurses are right-skewed, but there are 34 of them, so the Central Limit Theorem ought to be enough to overcome this. Hence, a \(t\)-test is fine.

If you drew a normal quantile plot, you get to the same place, but by a different route:

• you could say that both distributions are short-tailed, and therefore there will be no problem with the mean in either case.
• or you could say:
  • the non-rural distribution is close to normal, or is close enough to normal for a sample of size 18
  • the rural distribution is a bit further from normal, but the sample size of 34 should be large enough to overcome that.

Two points for appropriate discussion of normality for both groups based on your graph and the sample sizes.

(Aside: the boxplot says that the rural salaries have a right-skewed distribution, but the normal quantile plot says that this is not because the upper tail is long, but because the lower tail is short compared to a normal distribution.)

2. To justify the kind of \(t\)-test, decide whether the two groups have equal spread or not. On a boxplot, this is just comparing the heights of the boxes: the non-rural box is taller. On a normal quantile plot, the thing is to compare the slopes of the lines: the non-rural line is steeper. If you don’t see that, note that the non-rural salaries fill more of their box vertically than the rural ones do.

Either way, the two groups do not have equal spread, so we need to run a Welch \(t\)-test. One point for getting there by a reasonable route.

9. (3 points) Run your preferred \(t\)-test. What do you conclude from it, in the context of the data?

Start with the hypotheses. For a two-sample test, I find it easier to do these in words:

• null hypothesis: the rural and non-rural nurses have the same mean salary.
• alternative hypothesis: the rural nurses have a lower mean salary (or, the non-rural nurses have a higher mean salary).

If you want to do these in symbols, define your symbols carefully first. In particular, you need to be clear about which is rural and which is non-rural, such as this: let \(\mu_1\) be the population mean salary for rural nurses and \(\mu_2\) be the population mean salary for non-rural nurses. Then the null hypothesis is \(H_0: \mu_1 = \mu_2\) and the alternative is \(H_a: \mu_1 < \mu_2\). (If you define 1 and 2 the other way around, your alternative needs to be “greater than”.)

Now run the test. This is a one-sided alternative, so you need the right alternative = in your code. To see which group is “first” as far as t.test is concerned, one of these:

• on your graph, non-rural is first (on the left), which is the same as t.test wants
• alphabetically, non-rural is before rural

so your alternative has to say how non-rural compares to rural in that order in terms of what you are looking for, namely greater:

t.test(NurseSalaries ~ Location, data = homes, alternative = "greater")

    Welch Two Sample t-test

data:  NurseSalaries by Location
t = 3.3215, df = 26.11, p-value = 0.001326
alternative hypothesis: true difference in means between group non-rural and group rural is greater than 0
95 percent confidence interval:
 778.8986      Inf
sample estimates:
mean in group non-rural     mean in group rural 
               4859.500                3258.706 

The P-value, 0.0013, is less than 0.05, so we reject the null hypothesis in favour of the alternative, and conclude that rural nurses (in the population of all nurses, or all New Mexico nurses) do have a lower mean salary than non-rural nurses.

Extras (there are lots of them):

Extra 1: if you get your alternative the wrong way around, this happens:

t.test(NurseSalaries ~ Location, data = homes, alternative = "less")

    Welch Two Sample t-test

data:  NurseSalaries by Location
t = 3.3215, df = 26.11, p-value = 0.9987
alternative hypothesis: true difference in means between group non-rural and group rural is less than 0
95 percent confidence interval:
    -Inf 2422.69
sample estimates:
mean in group non-rural     mean in group rural 
               4859.500                3258.706 

and the very large P-value should make you stop and think, because that seems very inconsistent with the boxplot: this P-value, taken at face value, says that there is almost no evidence of a difference in means.

Extra 2: doing a one-sided test gets you a “one-sided confidence interval”, which is not the way we do them. If you want a confidence interval (to learn how far different the means could be in the population), re-do the test two-sided and ignore the P-value:

t.test(NurseSalaries ~ Location, data = homes)

    Welch Two Sample t-test

data:  NurseSalaries by Location
t = 3.3215, df = 26.11, p-value = 0.002652
alternative hypothesis: true difference in means between group non-rural and group rural is not equal to 0
95 percent confidence interval:
  610.3339 2591.2543
sample estimates:
mean in group non-rural     mean in group rural 
               4859.500                3258.706 

This 95% confidence interval⁵ says that the non-rural nurses have a mean salary that is between $610 and $2590 higher than the rural ones. (To judge by the means, this is per month.)

Extra 3: this time, we assessed the normality and equal spreads first by looking at the graph, so that we are happy that the assumptions for the test have been met. In the one-sample \(t\) question, we assessed the assumptions afterwards. For a two-sample question like this one, it is usually more natural to assess the assumptions at the same time as assessing equal spread (done using the same graph), but for a one-sample test you could do it either way. At this point in the course, we don’t have any alternative tests to do if we don’t like the normality, so the question is likely to be “justify the test you know how to do” (as here), especially if we are assessing assumptions first.

Extra 4: there are actually only 60 licenced nursing homes in New Mexico. It is likely that these 52 are the ones that completed the survey whose results are summarized here (and the other 8 might be in some way different), but even if these 52 are a random sample of the 60, we have sampled a large fraction of a population whose size we know, so to do this properly we should apply a finite population correction along the lines of STAC53. For the purposes of this question, however, we are pretending that our 52 nursing homes are a random sample from a population that is (much) larger than 52, hence the last sentence in the data description.

Footnotes

1. A rather insane choice, but you never know.

2. If you did this on r.datatools, your histogram will be something like my alternative version shown in question 2, and will not look so normal. Having said that, it is close to having a “flat top”, which is not as badly normal as something more obviously skewed (with a long tail at one end), or having outliers.

3. Using recode_values.

4. You cannot use kebab case in R variable names, however, since the - will be interpreted as a minus.

5. If you want something other than 95%, use conf.level as usual.