STAC32 Assignment 5

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.

Packages

library(tidyverse)
library(smmr)

Stainless steel

The data in http://datafiles.ritsokiguess.site/stainless-steel.csv are from an investigation of the manufacture of stainless steel. In particular, the chromium content of the stainless steel manufactured by this process is supposed to be 19% by weight on average. (“Average” could be either mean or median.) Twelve specimens of steel were selected from a larger batch, and the chromium content of each specimen was measured (percent by weight). In the data file, the first column identifies the specimen and the second column contains the chromium content of that specimen.

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

A quick one-pointer for you:

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

Rows: 12 Columns: 2
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (2): sample, percent

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

steel

The specimen is identified in sample (which we will ignore) and the chromium content is in percent. (Given the column names, you could call the dataframe chromium if you wanted.)

Extra: this is the stainless data from the nsm3data package.

(2 points) Draw a suitable graph of the chromium content values, and explain briefly why a \(t\)-test for the mean would not be appropriate.

You have two choices here. The obvious one is a histogram, with not too many bins because there are only 12 observations:

ggplot(steel, aes(x = percent)) + geom_histogram(bins = 6)

This distribution is skewed to the right (or, there are two observations much higher than the others).

The second choice (by thinking ahead and realizing that normality or lack thereof is likely to be an issue here) is a normal quantile plot:

ggplot(steel, aes(sample = percent)) + stat_qq() + stat_qq_line()

This shows two high outliers rather more clearly, or you might read the whole plot as having a curved shape, in which case right skewness is what you find.

In order to run a \(t\)-test for the mean chromium content, we need the distribution to be normal enough given the sample size. The sample size here, 12, is not especially large, and the distribution is badly non-normal, so the Central Limit Theorem is not going to rescue us, and a \(t\)-test will be a bad idea.

Extra: the above is all I need, but you might also think of finding a bootstrap sampling distribution of the sample mean. The code is as before, but I am going to be lazy and use my package bstrap to do the heavy lifting:

library(bstrap)
plot_dist(steel$percent)

The upper end of the distribution is not far off the line, but the lower end definitely veers away from the line. With all the extreme points above the line, I would call this a curved shape and say that this distribution is still right-skewed; we would need a larger sample size to be happy running a \(t\)-test here.

You might also notice that this bootstrap sampling distribution has “waves”: it goes up less steeply than than the line for a while, then suddenly jumps up more steeply, then is less steep than the line again, then suddenly jumps more steeply again, and so on. You might remember¹ from the postage stamp question on the worksheet, where we were looking at a distribution with multiple modes, that a sudden steep bit on a normal quantile plot indicates a “hole” where there is less data than you’d expect, and a flat bit indicates a “peak” where there is more data than you’d expect. With that in mind, we’d expect holes in the bootstrap sampling distribution around 18, 18.5, 19.2, 20, 20.7, 21.4 and maybe another one around 22, and peaks in between those. Remember that the actual values are on the \(y\)-axis, so you want the \(y\)-coordinate of the steep parts.

I can actually verify this. bstrap will also let you save the bootstrap sampling distribution to do other things with:

d <- make_dist(steel$percent)
d

and now we can make a histogram of this and look for extra peaks. For this kind of thing, we need a lot of bins, but we have a lot of observations from the bootstrap sampling distribution of the sample mean, so:²

ggplot(d, aes(x = x_mean)) + geom_histogram(binwidth = 0.2)

I think I got the low points between the modes in about the right places before, and now you can see there are several peaks, at about 17.7, 18.3, 19, 19.6, 20.2, 20.8.³

So now, why? I think this is definitely something real rather than just chance, because I obtained 10,000 bootstrap samples, and the peaks are much higher than the surrounding troughs.

Let’s go back to the original distribution of chromium content values:

ggplot(steel, aes(sample = percent)) + stat_qq() + stat_qq_line()

All the observations were around 17.5, except for the two values 24.5 and 26 that were much bigger than the others. When we draw a bootstrap sample, what we get are the original data values, possibly with some missing ones and some repeats, so the mean of the bootstrap sample is likely going to be determined by which of the big ones we got, and how many of each. As a gross oversimplification, let’s pretend our observed sample has 10 “small” observations that are all 17.5, and two big ones that are 25.3, and see what happens to the bootstrap samples for varying numbers of the big ones.

if there are no big ones at all, the bootstrap sample mean will be 17.5.
if there is one big one, the sample mean will be

(17.5 * 11 + 25.3) / 12

[1] 18.15

if there are two big ones:

(17.5 * 10 + 25.3 * 2) / 12

[1] 18.8

if there are three big ones:

(17.5 * 9 + 25.3 * 3) / 12

[1] 19.45

if there are four big ones:

(17.5 * 8 + 25.3 * 4) / 12

[1] 20.1

and if there are five:

(17.5 * 7 + 25.3 * 5) / 12

[1] 20.75

That, despite the gross oversimplification involved, reproduces pretty well where the peaks are in the bootstrap sampling distribution of the sample mean, and explains why it came out the way it did. There is some variability around the peaks because the small observations are actually not all exactly 17.5 and the two big ones are actually slightly bigger and smaller than 25.3.

Hence, if you have outliers in your data, you can expect the bootstrap sampling distribution of the sample mean to have extra modes according to whether the bootstrap sample from which the mean came had zero, one, two, three,… copies of the outliers. This is as well as (most likely) also having an overall skewed shape.

(3 points) Is there any evidence that the median chromium content differs from the desired value? Run a suitable test and explain briefly.

The right test for the median is the sign test, and the desired “average” is 19 percent, so use this for your null median. We are looking for any evidence of a difference, so a two-sided test is called for.

Make sure smmr is loaded (having previously installed it), and run sign_test with the dataframe, quantitative column and null median in that order:

sign_test(steel, percent, 19)

$above_below
below above 
   10     2 

$p_values
  alternative    p_value
1       lower 0.01928711
2       upper 0.99682617
3   two-sided 0.03857422

The two-sided P-value is 0.039, which is less than 0.05, so there is evidence that the median chromium content is different from the desired value of 19 percent.

Extra: the table of aboves and belows tells you something about what is going on: 10 of the 12 sample values were below 19, and only 2 were above 19. With this much imbalance, it is not surprising that we reject the null hypothesis.

(2 points) Obtain a 80% confidence interval for the median chromium content.

This uses ci_median. Make sure to specify the confidence level, since it is not 95%:

ci_median(steel, percent, conf.level = 0.80)

[1] 17.50313 18.09961

17.50 to 18.10 percent.

Make sure to state your interval, rounded off suitably. The data values are given to one decimal, so state at the outside two decimals in your confidence interval. Your reader is not going to care about the other decimal places, and they probably don’t carry any meaning anyway, so don’t make your reader look at them. You’ll recall that in the confidence interval for the median, the end points are data values, so it is not surprising that the second decimal place is zero.

Extra 1: this interval is consistent with what we concluded from the test, in that 19 is outside the interval. Both the discussion in the above Extra and the result here point to the median being less than 19 percent rather than more.

Extra 2: in the problem about the postage stamps on the worksheet (specifically extra 1 to question 8), I talked about getting more accurate decimal places on the confidence interval. The default is two decimals, but if you want more, you use the input tol. Here that would go like this, for four decimals:

ci_median(steel, percent, conf.level = 0.80, tol = 1e-4)

[1] 17.50005 18.09996

and at this point you see that the endpoints of the interval are definitely 17.5 and 18.1. These must be two of the data values:

steel %>% 
  arrange(percent)

namely the fourth smallest and the fourth largest ones. There is actually a way of working out that it must be the fourth largest and smallest data values, but I like our trial and error procedure better when it comes to understanding what is going on.

Extra 3: I suppose I now need to show you how the exact way would go. The starting point is that P-values for the sign test come from the null distribution of the number of observations above (or below) the hypothesized median, which is (as we saw in lecture) binomial with \(n\) being the sample size (12 here) and \(p\) being 0.5.⁴

I’m going to make a table of probabilities from this binomial distribution. It will be most convenient to work with cumulative probabilities (that is, the probability of a certain value or less), which pbinom will give us:

tibble(x = 0:12) %>% 
  mutate(prob = pbinom(x, 12, 0.5))

We made an 80% confidence interval above. If we were using something like normal tables, we’d find the values that enclose the middle 80% of the distribution, leaving 10% each end. However, the binomial distribution is discrete (unlike the normal distribution that is continuous), so we won’t be able to get exactly 80%. The story is to get at least 80% in the middle, which means getting 10% or less in each end.

If you look at my table, you’ll see that \(x = 3\) would give you probability less than 10% at the low end, but \(x = 4\) gives you more (it gives you nearly 20% in the low end, and with the same at the high end, this would only be about a 60% interval). Hence we use \(x = 3\), which is the 4th row of this table (the first row being \(x = 0\)) and that will give us an interval of level

1 - 2 * 0.073

[1] 0.854

that is, at least the 80% we are claiming.

To complete the procedure: the actual interval is then the fourth smallest and fourth largest data values. I sorted the data values above (in extra 2); the 4th smallest data value is 17.5 and the fourth largest one is 18.1, so that’s the 80% confidence interval.⁵

According to this, this procedure is a simplified version of the Hettmansperger-Sheather interpolation method (just the first two steps), but I have also seen it elsewhere under other names; the ones I remember seeing first are Hahn and Meeker. It’s in their book “Statistical Intervals: A Guide for Practitioners”. I think the same procedure is the one given on page 82 there, but it’s given under the more general case of working out a CI for any percentile.

Babies

The data in http://datafiles.ritsokiguess.site/babies.csv are counts of the number of babies born in a certain year in 60 community hospitals (that are a random sample of “all community hospitals”).

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

As ever:

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

Rows: 60 Columns: 1
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (1): count

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

babies

I suppose “births” would be a better name, but something describing what you have is called for, one way or another.

Extra: these are the data from Exercise 2.3.6 in Biostatistics by Daniel. Where I got the values from, the one column was unhelpfully called V1 so I named it something more suitable for you to work with.

(2 points) Make a suitable histogram of the data, and from it explain briefly why the mean number of babies born would not be a sensible summary.

There are 60 observations, which would normally suggest about 8 bins for your histogram. There is a suggestion that the distribution is going to be somehow non-normal, though (otherwise the mean would be a sensible summary), and in that case, having a few more bins will show you the shape better. I went with 10 bins:

ggplot(babies, aes(x = count)) + geom_histogram(bins = 10)

This distribution is bimodal: it has two peaks, one around 30 and the other around 55, with a “hole” at around 38 between them.⁶

The mean will not be a sensible summary here because the distribution is clearly not normal: the mean is

babies %>% summarize(mean_babies = mean(count))

which is in between the peaks.

Extra 1: we are going to work with the median later, which I’m not really sure is any better. It looks as if the community hospitals we sampled from really divide into large ones that have about 55 births a year, and small ones that have about 30. A probably better way of analyzing data like these is to fit a “mixture model” in which an observation has a certain probability \(p\) of being from a normal distribution with mean \(\mu_1\) and SD \(\sigma_1\), and probability \(1-p\) of being from a normal distribution with mean \(\mu_2\) and SD \(\sigma_2\). To fit a mixture model, you have to estimate all five parameters (guessing that \(\mu_1\) will be about 30 and \(\mu_2\) about 55), but the problem is that you don’t know for sure whether a particular community hospital is a small one or a large one, especially if it had about 45 births: is it a small one that had a lot of births, or a large one that didn’t have very many? This is unlike a two-sample problem, where you know for sure which sample each observation came from.

In any case, mixture models are beyond our scope right now, so we will work with the median.

Extra 2: I didn’t ask you to draw a normal quantile plot, even though normality is definitely relevant here, because if you do, this is what you get:

ggplot(babies, aes(sample = count)) + stat_qq() + stat_qq_line()

This says “short tails”. I don’t know, looking at your histogram, whether you agree with that. The problem is that on a normal quantile plot, it’s very difficult to distinguish a short-tailed distribution like the uniform⁷ from a bimodal one.

Let me generate 60 random uniform observations⁸ and see whether we can tell the difference:

uniform <- tibble(x = runif(60, 0, 1)) 
ggplot(uniform, aes(sample = x)) + stat_qq() + stat_qq_line()

The uniform data more or less hugs the line in the middle and deviates only in the tails, while our bimodal data dip below the line then jump up above the line. A subtle difference indeed.

Like the stamp data on the worksheet, or the bootstrap sampling distribution of the sample mean of the stainless steel data on this assignment, we can look for where our normal quantile plot of babies born is steeper than the line (indicating less data than expected). This, reading on the \(y\) (data) axis, is between about 35 and 40, and that corresponds to the “hole” on our histogram, with a peak on each side. It is hard to judge where the normal quantile plot is flattest (which is where the peaks are), but if there is one hole, there is likely to be one peak on each side of it (that is to say, two modes).

In comparison, the normal quantile plot of the uniform data doesn’t really have anywhere (beyond chance) where it is clearly steeper than the line, and so doesn’t have a “hole”:

ggplot(uniform, aes(x = x)) + geom_histogram(bins = 8)

Those two modes are (I am asking you to believe) just chance, in that if we were to generate more uniform data, they would move or disappear. Try this yourself, and see whether the modes always appear and always in the same place.

(4 points) A hospital administrator tells you that the median number of births in community hospitals last year was 37, and asks you to assess the evidence that the median number of births has increased. What do you conclude, and why?

This means to do a suitable test for the population median, that is to say, a sign test, using last year’s value as the null median, and testing whether the population median has increased. Dataframe, quantitative column, null median in that order:

sign_test(babies, count, 37)

$above_below
below above 
   21    38 

$p_values
  alternative    p_value
1       lower 0.99083146
2       upper 0.01817159
3   two-sided 0.03634318

The null hypothesis is that the population median is 37, and the alternative is that it is greater than 37. The correct one-sided P-value is 0.018, which is smaller than 0.05, so we reject the null hypothesis and conclude that the population median is now greater than 37.

(2 points) The hospital administrator who asked you to assess the evidence vaguely remembers the sign test from their time at university. They say to you “only about 35% of the data values are below 37, which is not far from a 50-50 split. I don’t think you should be rejecting a median of 37.” How do you respond to this?

The “about 35%” is correct:

21 / (38 + 21)

[1] 0.3559322

What the test is saying is that if 37 is the correct population median, then it is actually unlikely to get as few as 21 values below the population median; we would expect to see more than that: that is, closer still to a 50-50 split than we observed. The reason that we can reject with a 65-35 split is that the sample size is on the large side (\(n = 60\)); the larger the sample size, the closer to a 50-50 split we would expect to see, if the null hypothesis is true.

Make the point that if the null is true, what we observed actually is unusual, and that’s why we were able to reject the null hypothesis.

(2 points) Obtain a 90% confidence interval for the population median.

Get the interval with the right level of confidence and round it off, bearing in mind that the numbers of babies born have to be whole numbers:

ci_median(babies, count, conf.level = 0.90)

[1] 40.00281 49.99585

so the interval is from 40 to 50: with 90% confidence, the median number of babies born in a community hospital is between 40 and 50.

Cite one decimal if you want (40.0 to 50.0),⁹ but not more than that; the decimals beyond that are literally noise caused by the workings of the trial and error.

Extra: let me see if I can get this exactly, in the same manner as for the stainless steel question (in Extra 3 there). This time the sample size is 60, so we need a binomial distribution with \(n = 60\) and \(p = 0.5\):

tibble(x = 0:60) %>% 
  mutate(prob = pbinom(x, 60, 0.5))

We just found a 90% interval, which means 90% (at least) in the middle and 5% (or less) in each end. 0.05 is \(5 \times 10^{-2}\) in scientific notation, and \(x = 23\) has cumulative probability just less than that, with the probability for \(x = 24\) being bigger. The table starts at zero, so we want the 24th smallest and largest data values.¹⁰

We don’t have the sorted data values, but we can get them easily enough:

babies %>% arrange(count)

The 24th smallest value is 40, and the 24th biggest value (go to the end and count backwards) is 50, reproducing our confidence interval. The actual confidence level is, working out the middle as the whole thing minus the two ends:

1 - 0.04623 * 2

[1] 0.90754

which is a little more than 90%. (Exercise for you: if you wanted a 95% confidence interval, show that you would need to use the 22nd smallest and largest values and that using the 23rd smallest and largest would not quite work.)

Footnotes

Hi! I’m Troy McClure. You might remember me from such films as…↩︎
I used binwidth rather than bins on this one because I saw before that the “holes” were about 0.7 apart and if I have a binwidth that is something like a quarter of that, I should be able to see both the holes and the peaks.↩︎
The locations of the peaks are easier to read off a histogram than to infer from a normal quantile plot, because on the latter, you have to judge where the trace of points is flattest, which is harder work than judging where it is steepest.↩︎
Because if the null hypothesis is true, each observation has probability 0.5 of being above or below the hypothesized median.↩︎
Which would also be the 85% confidence interval.↩︎
The exact place you think the peaks and/or hole are will depend on the number of bins you choose, but if you have enough bins you will see the two peaks with a hole in between.↩︎
Which you saw on assignment 4.↩︎
Because our counts of babies born also had 60 observations↩︎
This would confirm that the ends of the interval are actually data values, since the data values are whole-number counts.↩︎
There was a lot in that paragraph, and you might need to take it slowly. Make sure you get each thing before moving on to the next. When I wrote that paragraph, I had just worked out the procedure for the stainless steel data and was happy that I understood it, so that paragraph just came tumbling out.↩︎