Problem 1

Let \(N = 6\) and \(n = 3\). The population values are \(y_1 = 98\), \(y_2 = 102\), \(y_3 = 154\), \(y_4 = 133\), \(y_5 = 190\), \(y_6 = 175\). We compare two sampling plans for estimating \(\bar{y}_U\).

Setup

# Population values
y <- c(98, 102, 154, 133, 190, 175)
N <- length(y)

# Plan 1: eight samples, each with probability 1/8
plan1_samples <- list(c(1,3,5), c(1,3,6), c(1,4,5), c(1,4,6),
                      c(2,3,5), c(2,3,6), c(2,4,5), c(2,4,6))
plan1_probs   <- rep(1/8, 8)

# Plan 2: three samples with unequal probabilities
plan2_samples <- list(c(1,4,6), c(2,3,6), c(1,3,5))
plan2_probs   <- c(1/4, 1/2, 1/4)

(a) Population mean

\[\bar{y}_U = \frac{1}{N}\sum_{i=1}^{N} y_i\]

ybarU <- mean(y)
ybarU
[1] 142

The population mean is \(\bar{y}_U = 142\).

(b) Sampling distribution of \(\bar{y}\)

For each plan, we use:

  • \(E[\bar{y}] = \sum_{S} P(S)\,\bar{y}_S\)
  • \(V[\bar{y}] = \sum_{S} P(S)\,\big(\bar{y}_S - E[\bar{y}]\big)^2\)
  • \(\text{Bias}(\bar{y}) = E[\bar{y}] - \bar{y}_U\)
  • \(\text{MSE}(\bar{y}) = V[\bar{y}] + \big[\text{Bias}(\bar{y})\big]^2\)
# Table of each sample, its values, sample mean, and probability
sample_table <- function(samples, probs) {
  data.frame(
    Sample      = sapply(samples, function(s) paste0("{", paste(s, collapse = ","), "}")),
    Values      = sapply(samples, function(s) paste(y[s], collapse = ", ")),
    Sample_Mean = round(sapply(samples, function(s) mean(y[s])), 3),
    P_S         = as.character(MASS::fractions(probs))
  )
}

# E, V, Bias, MSE of the sample mean under a plan
eval_plan <- function(samples, probs) {
  ybar <- sapply(samples, function(s) mean(y[s]))
  E    <- sum(probs * ybar)
  V    <- sum(probs * (ybar - E)^2)
  Bias <- E - ybarU
  MSE  <- V + Bias^2
  c(E = E, V = V, Bias = Bias, MSE = MSE)
}

Plan 1

sample_table(plan1_samples, plan1_probs)
   Sample        Values Sample_Mean P_S
1 {1,3,5}  98, 154, 190     147.333 1/8
2 {1,3,6}  98, 154, 175     142.333 1/8
3 {1,4,5}  98, 133, 190     140.333 1/8
4 {1,4,6}  98, 133, 175     135.333 1/8
5 {2,3,5} 102, 154, 190     148.667 1/8
6 {2,3,6} 102, 154, 175     143.667 1/8
7 {2,4,5} 102, 133, 190     141.667 1/8
8 {2,4,6} 102, 133, 175     136.667 1/8
plan1 <- eval_plan(plan1_samples, plan1_probs)
round(plan1, 3)
      E       V    Bias     MSE 
142.000  18.944   0.000  18.944
  1. \(E[\bar{y}] = 142\)
  2. \(V[\bar{y}] = 18.944\)
  3. \(\text{Bias}(\bar{y}) = 0\)
  4. \(\text{MSE}(\bar{y}) = 18.944\)

Plan 2

sample_table(plan2_samples, plan2_probs)
   Sample        Values Sample_Mean P_S
1 {1,4,6}  98, 133, 175     135.333 1/4
2 {2,3,6} 102, 154, 175     143.667 1/2
3 {1,3,5}  98, 154, 190     147.333 1/4
plan2 <- eval_plan(plan2_samples, plan2_probs)
round(plan2, 3)
      E       V    Bias     MSE 
142.500  19.361   0.500  19.611
  1. \(E[\bar{y}] = 142.5\)
  2. \(V[\bar{y}] = 19.361\)
  3. \(\text{Bias}(\bar{y}) = 0.5\)
  4. \(\text{MSE}(\bar{y}) = 19.611\)

Summary

results <- round(rbind(`Plan 1` = plan1, `Plan 2` = plan2), 3)
results
           E      V Bias    MSE
Plan 1 142.0 18.944  0.0 18.944
Plan 2 142.5 19.361  0.5 19.611

(c) Which sampling plan is better?

Plan 1 is better. Under Plan 1, \(\bar{y}\) is unbiased (\(E[\bar{y}] = 142 = \bar{y}_U\)). Under Plan 2, \(\bar{y}\) has a bias of 0.5. Plan 1 also has a smaller variance (18.944 vs. 19.361) and a smaller MSE (18.944 vs. 19.611). On average, Plan 1 gives estimates that are both centered on the true mean and closer to it.

Problem 5

An SRS of size \(n = 30\) is taken from a population of size \(N = 100\). The sample values are in the data file srs30.csv.

Setup

if (file.exists("srs30.csv")) {
  srs30 <- read.csv("srs30.csv")
} else if (file.exists("srs30-2.csv")) {
  srs30 <- read.csv("srs30-2.csv")
} else {
  srs30 <- data.frame(y = c(8, 5, 2, 6, 6, 3, 8, 6, 10, 7, 15, 9, 15, 3, 5,
                            6, 7, 10, 14, 3, 4, 17, 10, 6, 14, 12, 7, 8, 12, 9))
}
y <- srs30$y

N <- 100           
n <- length(y)     
n
[1] 30
y
 [1]  8  5  2  6  6  3  8  6 10  7 15  9 15  3  5  6  7 10 14  3  4 17 10  6 14
[26] 12  7  8 12  9

(a) Sampling weights

\[w_i = \frac{1}{\pi_i} = \frac{N}{n}.\]

w <- rep(N / n, n)
unique(w)
[1] 3.333333

Each unit in the sample has a sampling weight of \(w_i = 100/30 = 3.3333\). Each sampled unit represents about 3.33 units in the population.

(b) Estimated population total

\[\hat{t} = \sum_{i \in S} w_i\, y_i = N\bar{y}\]

t_hat <- sum(w * y)
t_hat
[1] 823.3333
# Check: same as N * ybar
N * mean(y)
[1] 823.3333

The estimated population total is \(\hat{t} = 823.33\).

(c) 95% confidence interval for \(t\)

With the finite population correction (fpc), the standard error of \(\hat{t}\) is

\[SE(\hat{t}) = N\sqrt{\left(1 - \frac{n}{N}\right)\frac{s^2}{n}}.\]

The 95% CI is \(\hat{t} \pm t_{0.025,\,n-1}\, SE(\hat{t})\).

s2    <- var(y)                  
fpc   <- 1 - n / N               
tcrit <- qt(0.975, df = n - 1)   

# With fpc
se_fpc <- N * sqrt(fpc * s2 / n)
ci_fpc <- t_hat + c(-1, 1) * tcrit * se_fpc

# Without fpc
se_nofpc <- N * sqrt(s2 / n)
ci_nofpc <- t_hat + c(-1, 1) * tcrit * se_nofpc

results <- data.frame(
  Method = c("With fpc", "Without fpc"),
  SE     = round(c(se_fpc, se_nofpc), 2),
  Lower  = round(c(ci_fpc[1], ci_nofpc[1]), 2),
  Upper  = round(c(ci_fpc[2], ci_nofpc[2]), 2)
)
knitr::kable(results)
Method SE Lower Upper
With fpc 61.06 698.45 948.21
Without fpc 72.98 674.07 972.59

Using the fpc, a 95% confidence interval for the population total is (698.45, 948.21).

Does the fpc make a difference? Yes. The sampling fraction is \(n/N = 30/100 = 0.30\), so the fpc is \(1 - 0.30 = 0.70\). That shrinks the standard error by a factor of \(\sqrt{0.70} \approx 0.837\), from 72.98 to 61.06. The CI with the fpc is about 16.3% narrower than the CI without it. Because we sampled 30% of the population, the fpc should not be ignored here.

Problem 6

A university has \(N = 807\) faculty members. An SRS of \(n = 50\) faculty members was taken, and the number of refereed publications was recorded for each.

Setup

# Frequency table from the problem
pubs  <- 0:10
freq  <- c(28, 4, 3, 4, 4, 2, 1, 0, 2, 1, 1)

y <- rep(pubs, times = freq)

N <- 807
n <- length(y)
n
[1] 50
knitr::kable(data.frame(Publications = pubs, Faculty = freq))
Publications Faculty
0 28
1 4
2 3
3 4
4 4
5 2
6 1
7 0
8 2
9 1
10 1

(a) Histogram and shape

hist(y,
     breaks = seq(-0.5, 10.5, by = 1),
     col    = "lightpink",
     border = "white",
     main   = "Refereed Publications per Faculty Member (n = 50)",
     xlab   = "Number of Refereed Publications",
     ylab   = "Number of Faculty Members",
     xaxt   = "n")
axis(1, at = 0:10)

The distribution is strongly skewed to the right. Most faculty in the sample (28 of 50, or 56%) have zero refereed publications. The counts drop off quickly, and a few faculty members have many publications (up to 10), which makes a long right tail.

(b) Estimated mean and SE

\[\bar{y} = \frac{1}{n}\sum_{i \in S} y_i, \qquad SE(\bar{y}) = \sqrt{\left(1 - \frac{n}{N}\right)\frac{s^2}{n}}\]

ybar <- mean(y)
s2   <- var(y)
fpc  <- 1 - n / N

se_ybar <- sqrt(fpc * s2 / n)

round(c(ybar = ybar, s2 = s2, fpc = fpc, SE = se_ybar), 4)
  ybar     s2    fpc     SE 
1.7800 7.1955 0.9380 0.3674

The estimated mean number of publications per faculty member is \(\bar{y} = 1.78\), with \(SE(\bar{y}) = 0.3674\).

(c) Is \(\bar{y}\) approximately normal?

# Sample skewness (a measure of how lopsided the data are)
skew <- mean((y - ybar)^3) / (mean((y - ybar)^2))^(3/2)
round(skew, 3)
[1] 1.545

Probably not, or only roughly. The Central Limit Theorem says \(\bar{y}\) becomes approximately normal as \(n\) grows. But how large \(n\) has to be depends on how skewed the population is. These data are heavily right-skewed (sample skewness \(\approx 1.54\)), with over half the values equal to 0. A sample of \(n = 50\) may not be large enough for the sampling distribution of \(\bar{y}\) to be close to normal. It will likely still be somewhat right-skewed. A common rule of thumb (Sugden et al., cited in Lohr) is that \(n\) should be greater than \(28 + 25 \times (\text{skewness})^2\):

n_needed <- 28 + 25 * skew^2
round(n_needed)
[1] 88

Because \(n = 50\) is less than the rule-of-thumb value of about 88, we should be cautious about assuming \(\bar{y}\) is normally distributed.

(d) Proportion with no publications and 95% CI

\[\hat{p} = \frac{\#\{y_i = 0\}}{n}, \qquad SE(\hat{p}) = \sqrt{\left(1 - \frac{n}{N}\right)\frac{\hat{p}(1-\hat{p})}{n-1}}\]

p_hat <- mean(y == 0)
se_p  <- sqrt(fpc * p_hat * (1 - p_hat) / (n - 1))
ci_p  <- p_hat + c(-1, 1) * qnorm(0.975) * se_p

round(c(p_hat = p_hat, SE = se_p, Lower = ci_p[1], Upper = ci_p[2]), 4)
 p_hat     SE  Lower  Upper 
0.5600 0.0687 0.4254 0.6946

The estimated proportion of faculty with no refereed publications is \(\hat{p} = 0.56\). A 95% confidence interval is (0.425, 0.695). We are 95% confident that between about 42.5% and 69.5% of the university’s faculty members have no refereed publications.

Problem 10

Which SRS design gives the most precision for estimating a population mean, assuming every population has the same population variance \(S^2\)?

  1. An SRS of size 400 from a population of size 4,000
  2. An SRS of size 30 from a population of size 300
  3. An SRS of size 3,000 from a population of size 300,000,000

Approach

For an SRS, the variance of the sample mean is

\[V(\bar{y}) = \left(1 - \frac{n}{N}\right)\frac{S^2}{n}.\]

Calculations

designs <- data.frame(
  Design = 1:3,
  n      = c(400, 30, 3000),
  N      = c(4000, 300, 300000000)
)

designs$sampling_fraction <- designs$n / designs$N
designs$fpc               <- 1 - designs$sampling_fraction
designs$var_over_S2       <- designs$fpc / designs$n        
designs$se_over_S         <- sqrt(designs$var_over_S2)      

knitr::kable(designs,
             digits    = 6, format.args = list(big.mark = ","),
             col.names = c("Design", "n", "N", "n/N", "fpc",
                           "V(ybar) / S^2", "SE(ybar) / S"))
Design n N n/N fpc V(ybar) / S^2 SE(ybar) / S
1 400 4,000 0.10000 0.90000 0.002250 0.047434
2 30 300 0.10000 0.90000 0.030000 0.173205
3 3,000 300,000,000 0.00001 0.99999 0.000333 0.018257 
best <- designs$Design[which.min(designs$var_over_S2)]
best
[1] 3

Relative precision

# How many times larger each design's variance is compared to the best design
designs$relative_variance <- designs$var_over_S2 / min(designs$var_over_S2)
knitr::kable(designs[, c("Design", "n", "N", "relative_variance")],
             digits    = 2, format.args = list(big.mark = ","),
             col.names = c("Design", "n", "N", "Variance relative to best"))
Design n N Variance relative to best
1 400 4,000 6.75
2 30 300 90.00
3 3,000 300,000,000 1.00

Conclusion

Design 3 (an SRS of 3,000 from a population of 300,000,000) gives the most precision.

Designs 1 and 2 both sample 10% of their populations, so they have the same fpc (0.90). Design 3 samples only a tiny fraction of its population (fpc \(\approx 1\)), but its sample size is much larger. Precision depends mainly on the sample size \(n\), not on the population size \(N\) or the sampling fraction. The fpc only matters when \(n/N\) is large, and here it shrinks the variance by at most 10%. Since Design 3 has by far the largest \(n\), its variance is the smallest. It is about 6.8 times smaller than Design 1 and about 90 times smaller than Design 2.

Problem 12

We want to estimate the percentage of people immunized against measles in each of 5 Arizona cities. Each estimate should have a margin of error of \(e = 0.04\) (4 percentage points) at 95% confidence.

Approach

We don’t know \(p\) ahead of time, so we use \(p = 0.5\). This gives the largest possible value of \(p(1-p)\), so it’s the most conservative choice. Without the fpc, the required sample size is

\[n_0 = \frac{z_{\alpha/2}^2\, p(1-p)}{e^2}.\]

With the finite population correction, the required sample size is

\[n = \frac{n_0}{1 + n_0/N}.\]

Calculations

cities <- data.frame(
  City       = c("Casa Grande", "Gila Bend", "Jerome", "Phoenix", "Tempe"),
  Population = c(48571, 1922, 444, 1445632, 161719)
)

z <- qnorm(0.975)   
e <- 0.04          
p <- 0.5            

n0 <- z^2 * p * (1 - p) / e^2
n0
[1] 600.2279
cities$n_without_fpc <- ceiling(n0)
cities$n_with_fpc    <- ceiling(n0 / (1 + n0 / cities$Population))
cities$difference    <- cities$n_without_fpc - cities$n_with_fpc
cities$pct_sampled   <- round(100 * cities$n_with_fpc / cities$Population, 2)

knitr::kable(cities,
             format.args = list(big.mark = ","),
             col.names = c("City", "Population", "n (no fpc)", "n (with fpc)",
                           "Difference", "% of city sampled"))
City Population n (no fpc) n (with fpc) Difference % of city sampled
Casa Grande 48,571 601 593 8 1.22
Gila Bend 1,922 601 458 143 23.83
Jerome 444 601 256 345 57.66
Phoenix 1,445,632 601 600 1 0.04
Tempe 161,719 601 599 2 0.37

Conclusion

Without the fpc, every city needs a sample of about 601 people. The fpc reduces the required sample size by \(n_0/(1 + n_0/N)\), which only matters when \(n_0\) is a meaningful fraction of \(N\).

  • The fpc makes a big difference for Gila Bend and Jerome. These are small towns. Gila Bend only needs about 458 people instead of 601, and Jerome only needs about 256 of its 444 residents.
  • The fpc makes only a small difference for Casa Grande (8 fewer people) and almost none for Tempe and Phoenix. Their populations are so large compared to the sample that the fpc is close to 1.

The sample size needed for a given precision depends very little on the population size unless the population is small.

Problem 13: Objects in a Bin

The bin contains \(N = 20{,}000\) objects (15,000 squares and 5,000 circles). Each object has a shape, a color (black or gray), and an area.

Setup

# Read the population.
if (file.exists("shapespop.csv")) {
  shapespop <- read.csv("shapespop.csv")
} else {
  if (!requireNamespace("SDAResources", quietly = TRUE)) {
    install.packages("SDAResources", repos = "https://cloud.r-project.org")
  }
  data(shapespop, package = "SDAResources")
}
str(shapespop)
Classes 'tbl_df', 'tbl' and 'data.frame':   20000 obs. of  5 variables:
 $ ID   : num  1 2 3 4 5 6 7 8 9 10 ...
  ..- attr(*, "format.sas")= chr "BEST"
 $ shape: chr  "square" "square" "square" "square" ...
  ..- attr(*, "format.sas")= chr "$"
 $ color: chr  "black" "black" "black" "black" ...
  ..- attr(*, "format.sas")= chr "$"
 $ area : num  41 37 19 21 14 58 49 38 58 10 ...
  ..- attr(*, "format.sas")= chr "BEST"
 $ conv : num  1 1 0 0 0 0 0 1 0 1 ...
  ..- attr(*, "format.sas")= chr "BEST"
 - attr(*, "label")= chr "SHAPESPOP                       "
N <- nrow(shapespop)
N
[1] 20000

(a) Select an SRS of size 200 and find the sampling weights

n <- 200

set.seed(3003)   
index  <- sample(1:N, size = n, replace = FALSE)
shapes <- shapespop[index, ]

shapes$weight <- N / n
unique(shapes$weight)
[1] 100
write.csv(shapes, "shapes_sample.csv", row.names = FALSE)

head(shapes)
         ID  shape color area conv weight
15403 15403 circle black   10    1    100
10002 10002 square  gray   23    0    100
10787 10787 square  gray   20    1    100
18994 18994 circle black   15    1    100
14049 14049 square  gray   23    1    100
3340   3340 square black    6    0    100

Each sampled object has a sampling weight of \(w_i = N/n = 20{,}000/200 = 100\). Each object in the sample represents 100 objects in the bin. The sample is saved as shapes_sample.csv.

(b) Histogram of the areas

hist(shapes$area,
     col    = "lavender",
     border = "white",
     main   = "Areas of Objects in the SRS (n = 200)",
     xlab   = "Area",
     ylab   = "Number of Objects")

summary(shapes$area)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   3.00   16.00   24.00   25.38   32.00   67.00

(c) Average area, 95% CI, and total area

\[\bar{y} = \frac{1}{n}\sum_{i\in S} y_i, \qquad SE(\bar{y}) = \sqrt{\left(1-\frac{n}{N}\right)\frac{s^2}{n}}, \qquad \hat{t} = N\bar{y}\]

fpc   <- 1 - n / N
tcrit <- qt(0.975, df = n - 1)

ybar    <- mean(shapes$area)
se_ybar <- sqrt(fpc * var(shapes$area) / n)
ci_ybar <- ybar + c(-1, 1) * tcrit * se_ybar

t_area    <- N * ybar
se_t_area <- N * se_ybar

round(c(Mean = ybar, SE = se_ybar, Lower = ci_ybar[1], Upper = ci_ybar[2]), 3)
  Mean     SE  Lower  Upper 
25.380  0.977 23.453 27.307 
round(c(Total_Area = t_area, SE_Total = se_t_area), 1)
Total_Area   SE_Total 
  507600.0    19545.4

The estimated average area is \(\bar{y} = 25.38\). A 95% CI is (23.45, 27.31). The estimated total area covered by all objects in the bin is \(\hat{t} = N\bar{y} = 507,600\).

(d) Total number of gray objects and 95% CI

\[SE(\hat{t}) = N\sqrt{\left(1-\frac{n}{N}\right)\frac{\hat{p}(1-\hat{p})}{n-1}}.\]

p_gray    <- mean(shapes$color == "gray")
t_gray    <- N * p_gray
se_t_gray <- N * sqrt(fpc * p_gray * (1 - p_gray) / (n - 1))
ci_gray   <- t_gray + c(-1, 1) * tcrit * se_t_gray

round(c(p_hat = p_gray, Total = t_gray, SE = se_t_gray,
        Lower = ci_gray[1], Upper = ci_gray[2]), 3)
   p_hat    Total       SE    Lower    Upper 
   0.350 7000.000  672.840 5673.189 8326.811

The estimated number of gray objects in the bin is 7,000. A 95% CI is (5,673, 8,327).

(e) Total number of circles and 95% CI

p_circ    <- mean(shapes$shape == "circle")
t_circ    <- N * p_circ
se_t_circ <- N * sqrt(fpc * p_circ * (1 - p_circ) / (n - 1))
ci_circ   <- t_circ + c(-1, 1) * tcrit * se_t_circ

round(c(p_hat = p_circ, Total = t_circ, SE = se_t_circ,
        Lower = ci_circ[1], Upper = ci_circ[2]), 3)
   p_hat    Total       SE    Lower    Upper 
   0.295 5900.000  643.319 4631.402 7168.598 
# Does the CI contain the true value of 5,000?
covers <- ci_circ[1] <= 5000 & 5000 <= ci_circ[2]
covers
[1] TRUE

The estimated number of circles in the bin is 5,900. A 95% CI is (4,631, 7,169). Yes, the CI includes the true population value of 5,000 circles.