Department of Statistics, Virginia Tech

Recently, the transit system increased their bus fare by $0.10 in an effort to fix their budget without significantly affecting the monthly expense of their passengers. Since most of their passengers take the Metro train, they did not increase the rail fare. Now, the administration wants to determine whether the median amount of fare spent in a month has significantly increased to a value higher than $200. A sample of regular customers revealed the following observed monthly fares.

`fares <- c(205, 188, 192, 220, 203, 197, 206, 209, 184, 200, 214, 217)`

- From this data, is there evidence to conclude that there has indeed been an increase in revenue?

One way to settle this is via a hypothesis test where \(\mathcal{H}_0\): the 50% quantile is less than 200; versus \(\mathcal{H}_1\): the 50% quantile is greater than 200. I.e., a lower-tailed test.

- In that case the \(p\)-value is \(\mathbb{P}(T_1 \leq t_1)\), where \(T_i \sim \mathrm{Bin}(n, p^*)\).

So in R we have

```
pstar <- 0.5
xstar <- 200
n <- length(fares)
t1 <- sum(fares <= xstar)
```

and then \(p\)-value is thus

```
phi <- pbinom(t1, n, pstar)
phi
```

`## [1] 0.387207`

- I.e., not enough evidence to reject \(\mathcal{H}_0\) at the 5% level.
- We conclude that the median amount of fare spent in a month has
*not*significantly increased to a value higher than $200.

We should not be surprised to find that this is the same as what we get out of the library function `binom.test`

.

`binom.test(t1, n, pstar, "less")`

```
##
## Exact binomial test
##
## data: t1 and n
## number of successes = 5, number of trials = 12, p-value = 0.3872
## alternative hypothesis: true probability of success is less than 0.5
## 95 percent confidence interval:
## 0.0000000 0.6847622
## sample estimates:
## probability of success
## 0.4166667
```

- But we have to be careful, because
`binom.test`

does not recognize the*definition*of a quantile, i.e., it does not differentiate between \(T_1\) and \(T_2\). We got the right answer here because only \(T_1\) was involved. We have to be more careful when performing upper-tailed and two-tailed tests. - We even get a 1-sided confidence interval, but it is not at all meaningful in our quantile context. It is for \(p\), and we want one for \(x_{p^*}\).

Entering college, freshmen must take an exam, and the upper quartile (75th quantile) is established at a score of 193. A particular high school sends 112 of its graduates to college. Suppose 100 students have a score less than or equal to 193, and 98 students have a score strictly less than 193. Conduct the quantile test and determine whether this established quantity holds for this high school.

The hypotheses are \(\mathcal{H}_0\): The 75th quantile is 193, versus \(\mathcal{H}_1\): the 75th quantile is not 193.

- Since this is a two-tailed test, the \(p\)-value is the minimum of two CDF calculations: \[ \varphi = 2* \min \left\{\mathbb{P}(T_1 \leq t_1), \mathbb{P}(T_2 \geq t_2) \right\}, \quad \mbox{ where } T_1, T_2 \sim \mathrm{Bin}(n, p^*). \]

So in R we have

```
pstar <- 0.75
t1 <- 100
t2 <- 98
n <- 112
```

and then the \(p\)-value is thus

```
phi <- 2*min(pbinom(t1, n, pstar), 1-pbinom(t2-1, n, pstar))
phi
```

`## [1] 0.001709746`

- An easy rejection of \(\mathcal{H}_0\) at the typical 5% level.
- This high school is exceptional (on the low end) for this exam.

What happens when we use `binom.test`

?

Sixteen transistors are selected at random from a large batch of transistors and are tested. The number of hours until failure is recorded for each one. We wish to find an approximately 90% confidence interval for the upper quartile. The data follows.

`fail <- c(67.1, 63.7, 46.9, 56.5, 63.2, 59.9, 47.2, 73.3, 63.3, 49.1, 78.5, 56.8, 63.4, 59.2, 64.1, 67.7)`

- To calculate the \(r\) and \(s\) indexes into the relevant order statistics we need to set \(\alpha\), \(n\) and \(p^*\).

```
alpha <- 0.1
n <- length(fail)
pstar <- 0.75
```

- Then the indexes are determined by the \(\alpha/2\) and \(1-\alpha/2\) quantiles of a \(\mathrm{Bin}(n, p^*)\) random variable.

```
r <- qbinom(alpha/2, n, pstar)
s <- qbinom(1-alpha/2, n, pstar)
c(r,s)
```

`## [1] 9 15`

- Then, we sort the observed data values and extract out the \(r\)th and \(s\)th ones.

```
fail.sorted <- sort(fail)
c(fail.sorted[r], fail.sorted[s])
```

`## [1] 63.3 73.3`

So our CI is [63.3, 73.3], but it doesn’t have cover exactly 90%. To figure out the exact confidence level requires inverting the quantile calculations.

```
alpha1 <- pbinom(r-1, n, pstar)
alpha2 <- 1-pbinom(s-1, n, pstar)
1 - alpha1 - alpha2
```

`## [1] 0.9093936`

- So to we actually got a 91% CI (rather than “at least” 91% since the data can reasonably be assumed to be continuous).

Lets revisit one of our earlier examples (we can’t do the Metro one because we don’t have the raw observations) and find a 95% confidence interval for the median.

```
alpha <- 0.05
n <- length(fares)
pstar <- 0.5
r <- qbinom(alpha/2, n, pstar)
s <- qbinom(1-alpha/2, n, pstar)
fares.sorted <- sort(fares)
c(fares.sorted[r], fares.sorted[s])
```

`## [1] 192 209`

- Clearly the $200 mark is right in there, so we’re not going to reject that \(\mathcal{H}_0\).

Lets calculate the actual confidence level.

```
alpha1 <- pbinom(r-1, n, pstar)
alpha2 <- 1-pbinom(s-1, n, pstar)
1 - alpha1 - alpha2
```

`## [1] 0.9077148`

- a bit narrower than we asked for, but that’s ok.

It easy to cut-and-paste to construct a new CI for a new data set, but it would be even easier if we didn’t have to. Lets make our own quantile CI function, with defaults for the median and 95% level.

```
quantile.confint <- function(x, pstar=0.5, level=0.95)
{
alpha <- 1-level
n <- length(x)
xs <- sort(x)
r <- qbinom(alpha/2, n, pstar)
s <- qbinom(1-alpha/2, n, pstar)
CI <- c(xs[r], xs[s])
alpha1 <- pbinom(r-1, n, pstar)
alpha2 <- 1-pbinom(s-1, n, pstar)
level.act <- 1 - alpha1 - alpha2
return(list(CI=CI, level.act=level.act))
}
```

Now lets try it out.

`quantile.confint(fail, 0.75, 0.9)`

```
## $CI
## [1] 63.3 73.3
##
## $level.act
## [1] 0.9093936
```

`quantile.confint(fares)`

```
## $CI
## [1] 192 209
##
## $level.act
## [1] 0.9077148
```

- There we go!