# Determine sample size using reversed goodness of fit

I am asked to work on a specific problem in which I have to calculate certain expenditures for an industry, consisting of a population of about 400 companies. Although I already suggested to conduct a survey including all companies, I am asked to take a sample (preferable a small one, max N=30). Furthermore, the population can be divided in nine strata, and I should take this into account (so 30 observations in total, divided over 9 strata).

The variable I am collecting from the companies are expenditures on clean up operations, and are measured as continuous variable (in currency units). We want to calculate the abolute total expenditures for the entire industry based on the sample companies.

The sample I am asked to take seems too small to make a good estimate, so I am looking for a way to justify taking a small sample and still have a rather accurate estimate (best estimate under suboptimal project conditions, so to say). I am not very certain on how to proceed here, but I tried a few methods. I have some population data, one variable should correlate (at least, in theory) with the expenditures I am looking for. My question is: do these methods make sense? Is one better than the other? Are there better methods?

Method 1:

Use the formula N = (Z * SD / E)^2 for each strata (using SD of the variable that should correlate with expenditure) and tweak E so that N becomes about 30. This led me to have an E which I calculated as (0,2 * mean of stratum). In other words, in order to have a sample that fits the project's budget, we have to put up with a margin of error that is 20% of the stratum mean. Seems a bit steep.

Method 2:

I put a goodness of fit test in reverse, and used a couple of population variables to calculate the minimum proportion of the population I need in order to stay within the critical value limit. There are 2 x 5 strata, of which 1 is empty, so there are (2-1)(5-1) = 4 degrees of freedom. The critical value for 4 degrees of freedom (p = 0.05) is 9.49. In order words: I made a small model that calculates a 'hypothethical' chi2 value for several variables, assuming a certain proportion of the population is sampled. Lucky (I guess?) for me, N=30 falls within the critical value limit of 9.49. This takes the distribution of observations over the different samples into account: it calculates per sample what proportion of the population in each sample should be in order to stay within the critical value limit.

I don't like either of them, but I am in a bit of a pickle here due to project budgets and what not. Any help will be gladly received!

• Sorry. Voting to put this on hold as 'unclear what you're asking', pending clarifications: There are 9 strata. Are you allowed to sample 30 from each stratum (up to 270 out of 400)? Or a total of 30? // What info would you collect on each company sampled? // How is expenditure measured? Hi/Med/Lo or in dollars (or euros or pounds)? // What do you mean by 'critical value limit'? – BruceET Aug 2 at 20:27
• Thank you for your reply and question for clarification. I updated the question and added clarifications on all your points. I hope it is more clear now. If any questions remain, I will gladly answer them. – Eelco Aug 3 at 12:10

The problem is not so much to determine a sample size--that is specified as $$30$$--as it is to allocate the sampling budget to the strata. Let's solve this problem from first principles--the same ones used to develop the standard sample size formulas.

This answer, although developed for wide application, focuses throughout on the situation of sampling a small stratified population, where potentially the sample may include appreciable subsets of some strata; the usual asymptotic formulae might not be adequate; and we should expect to apply "finite population correction factors" to those formulae. These factors will show up in the form $$n_i-k_i$$ at places where we might otherwise expect to see $$n_i$$ alone.

The fundamental ideas that we need to employ come down to these:

1. Because the objective is to estimate the total industry cleanup expenditure (presumably during a well-defined period), it would be good to find a sampling procedure--among all the possible sample designs--that comes close to minimizing the variance of that estimate.

2. The variance depends both on the sample design and on the distributions of total expenditures within each sampling stratum.

To analyze and understand the implications, let's begin with some notation.

• The population is the disjoint union of strata $$\mathcal{S}_i,$$ $$i=1, 2, \ldots, m=9.$$

• Stratum $$i$$ includes $$n_i \ge 1$$ companies--and you know, with no uncertainties, what stratum each company is in. This enables us to designate companies by a pair of indices: $$i$$ for its stratum and a second index $$j$$ running from $$1$$ to $$n_i$$ for the company itself.

• The cleanup expenditure for company $$(i,j)$$ is $$z_{ij}$$ (and the company will assess and report this value accurately, with no error). Let $$T$$ be the total expenditure.

A stratified (simple random) sampling plan will independently sample $$k_i$$ (a number between $$1$$ and $$n_i$$ inclusive) of the companies in stratum $$\mathcal{S}_i$$ without replacement. Let $$X_{ij}$$ be the binary random variable indicating whether company $$(i,j)$$ is in the sample. The best (least-variance) estimator of the total expenditure in each stratum is the average expenditure observed in the sample, multiplied by the stratum size $$n_i:$$

$$\widehat T = \sum_{i=1}^m n_i\left( \frac{1}{k_i} \sum_{j=1}^{n_i} X_{ij} z_{ij} \right) = \sum_{i=1}^m\sum_{j=1}^{n_i} \frac{n_i z_{ij}}{k_i} X_{ij}.$$

Idea $$(1)$$ above requires us to develop an expression for the variance of this estimator. The independent choices of samples without replacement in each stratum determine this variance because (a) independence allows us to compute the variance as a sum of contributions within each stratum and (b) sampling without replacement implies

$$E[X_{ij}] =\frac{k_i}{n_i}$$

(each company has a $$k_i/n_i$$ chance of being in the sample) and when $$k\ne j,$$

\eqalign{ E[X_{ij} X_{ik}] &= \Pr(j\text{ in the sample and } k\text{ in the sample})\\ &= \Pr(j\text{ in the sample} \mid k\text{ in the sample})\Pr(k\text{ in the sample})\\ &= \frac{k_i-1}{n_i-1}\frac{k_i}{n_i}. }

Consequently for $$k\ne j$$

$$\operatorname{Cov}(X_{ij}, X_{ik}) = E[X_{ij} X_{ik}] - E[X_{ij}]E[X_{ik}] = -\frac{k_i(n_i-k_i)}{n_i^2(n_i-1)}.$$

When $$k=j,$$

$$\operatorname{Cov}(X_{ij}, X_{ik}) = \operatorname{Var}(X_{ij}) = E[X_{ij}^2] - E[X_{ij}]^2 = \frac{k_i(n_i-k_i)}{n_i^2}.$$

A useful way to express all covariances at once employs the Kronecker delta $$\delta_{jk}$$ which is equal to $$1$$ when $$j=k$$ and otherwise is zero. This enables the two different expressions $$-1/(n_i-1)$$ and $$1$$ to be written

$$\frac{\delta_{jk} n_i - 1}{n_i-1} = \left\{\matrix{-\frac{1}{n_i-1} & j\ne k \\ 1 & j=k.}\right.$$

The covariance of any two terms that form the expression for $$\hat T$$ is zero when they come from different strata; and in a common stratum $$\mathcal{S}_i,$$ the foregoing immediately gives

\eqalign{ \operatorname{Cov}\left(\frac{n_i z_{ij}}{k_i} X_{ij}, \frac{n_i z_{ik}}{k_i} X_{ik}\right) &= \frac{n_i z_{ij}}{k_i} \frac{n_i z_{ik}}{k_i} \frac{k_i(n_i-k_i)}{n_i^2} \frac{\delta_{jk} n_i - 1}{n_i-1} \\ &= \frac{n_i-k_i}{k_i(n_i-1)} z_{ij}z_{ik}(\delta_{jk} n_i - 1) . }

These enable the variance of $$\widehat T$$ (a linear combination of the $$X_{ij}$$) to be computed as

$$\operatorname{Var}(\widehat T) = \sum_{i=1}^{m} \frac{n_i-k_i}{k_i(n_i-1)} \left(\sum_{j=1}^{n_i} n_iz_{ij}^2 - \sum_{j,k}^{n_i}z_{ij}z_{ik} \right)= \sum_{i=1}^{m} \frac{n_i-k_i}{k_i}\frac{n_i^2}{n_i-1}\sigma_i^2,\tag{*}$$

where

$$\sigma_i^2 = \frac{1}{n_i}\sum_{j=1}^{n_i} z_{ij}^2 - \left(\frac{1}{n_i}\sum_{j=1}^{n_i} z_{ij}\right)^2$$

refers to the variance of that stratum. (In any stratum with $$n_i=1,$$ necessarily $$k_i=1$$ which assures that stratum contributes zero to the variance. In effect, a stratum of size $$1$$ designates a company that will be in the sample no matter what. We may as well ignore all such strata in our analyses and thereby assume $$n_i\ge 2$$ from now on.)

Since the strata are determined, their sizes $$n_i$$ and their variances are determined, too (although the latter is not known!), leaving the sample sizes $$k_i$$ as the only quantities left to be determined: you will choose them to minimize $$(*)$$ so that their total is no greater than your sampling budget $$K=30.$$

It seems we're at an impasse: the optimal sample sizes depend on properties of the strata, but the purpose of sampling is precisely to estimate properties of the population! The key to progress lies in recognizing that these are different properties: given assumptions you make about the strata variances, you can design a conditionally optimal sampling program for estimating the population total.

Let's put off the discussion of those assumptions and suppose you have a sense of the strata variances. In that case, ignoring momentarily the requirement that $$k_i\ge 1$$ and $$k_i$$ is integral, a simple application of Lagrange Multipliers shows that when no optimal $$k_i$$ is less than $$1,$$ the variance of $$\widehat T$$ is minimized when sample sizes are nearly proportional to the stratum sizes and standard deviations. Specifically, the squared sample sizes are

$$k_i^2\ \propto\ \frac{n_i}{n_i-1} n_i^2\sigma_i^2 = \omega_i^2\tag{1}$$

(this defines the weights $$\omega_i^2$$).

Thus, you can meet your sampling budget of no more than $$K$$ samples by setting $$k_i = \frac{\omega_i}{\sum_{i=1}^m \omega_i} K.\tag{2}$$

If it turns out some of these are less than $$1,$$ set them to $$1$$ and re-apply $$(1)$$ and $$(2)$$ to the remaining $$k_i.$$ Iterate (less than $$m-1$$ times) to a solution.

Observations:

1. $$(2)$$ only depends on the ratios of the stratum variances rather than the variances themselves. Thus, you don't need to know the actual stratum variances--you only need to have a reasonable guess about their relative values.

2. The solution given in $$(1)$$ and $$(2)$$ is unlikely to be entirely integral. We need the best integral solution, not the best solution in real numbers. People usually just round the real optimum, but there's no guarantee this will be best. Strictly speaking you need special techniques to find the integral optimum efficiently, but as a practical matter rounding may be good enough. (The code example below goes a little further by exploring all nearby integral solutions, but that still provides no guarantee of finding the very best integral solution.)

3. Provided you don't round down to zero, you can assure every stratum is sampled. (When the solution suggests several $$k_i$$ ought to be close to zero and less than $$1,$$ you might consider combining their strata if there's hope those strata have similar average values of $$z_{ij}$$ and then re-solving the problem with the smaller number of strata.)

4. If you err in estimating the stratum variances, you will wind up with a suboptimal sample--but it will still be useful. It merely won't achieve the best possible expected variance in the estimate of the total. The example below illustrates this point by comparing solutions based on very different stratum variance estimates.

5. For strata that are not small, $$(1)$$ says to pick the sample sizes $$k_i$$ in proportion to $$n_i \sigma_i.$$ This provides rules of thumb for such special cases. In particular,

• When you believe the strata variances $$\sigma_i^2$$ are all roughly the same, the optimum solution is proportional allocation where $$k_i \propto n_i.$$

• When the strata are all roughly the same size, the optimal solution is to devote samples to strata in proportion to their standard deviations, $$k_i \propto \sigma_i.$$

The initial term $$n_i/(n_i-1)$$ in $$(1)$$ indicates that very small strata deserve slightly more samples than otherwise indicated by these rules of thumb. This is due to the rapidly decreasing uncertainty in the stratum total as the sample comes close to exhausting a stratum.

To illustrate this analysis and the observations, consider an instance of the original problem in which

1. The stratum sizes are $$(n_i) = (4, 10, 19, 29, 40, 53, 67, 87, 96).$$ They sum to $$N=400.$$

2. The stratum variances are $$(\sigma_i^2) = (140,100,74,53,38,27,19,14,10).$$

This models a typical situation in which a few small strata exhibit relatively high variances. These strata might be comprised of, say, the largest companies, or of similar companies with the most uncertain liabilities based on some other study.

The best non-integral solution to the constrained minimization of the variance of $$\widehat T$$ is $$k^{*} = (1, 1.59, 2.53, 3.24, 3.76, 4.19, 4.44, 4.65, 4.6).$$ Its first component was set by the constraint $$k_1 \ge 1$$ and, as intended, the values sum to $$K=30.$$ Among the $$256$$ nearby integral values, the best solution among those that sum to $$30$$ is $$k = (1,2,3,3,4,4,4,5,4).$$ Notice how $$k^{*}_9 = 4.6$$ was rounded down to $$4$$ (for otherwise the rounded values would sum to $$31,$$ which is too large.) As a basis for comparison, suppose you cannot estimate the variances and so you just take them all to be equal (which obviously is a gross error). The best integral solution is $$k_e = (1,1,1,2,3,4,5,6,7).$$ Compared to $$k,$$ this solution oversamples the high-variance strata and undersamples the low-variance ones.

To compare these solutions, consider the variance of $$\widehat T$$ in each case: for $$k$$ is is $$120857$$ while for $$k_e$$ it is $$140618,$$ about $$16\%$$ greater. That means the incorrect estimation of the stratum variances has made its $$30$$ samples about as effective as $$25$$ optimally-deployed samples. Considering how wrong the stratum variance estimates were, that's not much of a loss: it indicates both that (1) there is something to be gained from this sample design process but (2) the gains are relatively insensitive to the variance assumptions you need to make to carry it out.

### Reference

Thompson, Steven K. (1992), "Sampling" (Part 1), Wiley.

### Code

This R code created the example. It includes functions that are generally useful for finding optimal stratified samples of small populations. It is practical for up to about $$m\approx 20$$ strata; beyond that, the brute-force search for integral solutions is too expensive and would have to be replaced by a more sophisticated (branch-and-bound) algorithm or by sampling.

#
# Find k*, the optimal (non-integral) solution.
#
ss <- function(K, n, sigma.2) {
i <- n >= 1
repeat {
K0 <- K - sum(!i)
omega <- sqrt(n^3 * sigma.2 / (n-1))
omega <- omega / sum(omega[i])
j <- i & (omega * K0 >= 1)
if (all(j==i)) break
i <- j
}
k <- rep(1, length(n))
k[i] <- K0 * omega[i]
return(k)
}
#
# Evaluate the variance of the estimator of the population sum.
#
s.var <- function(k, n, sigma.2) {
sum((n-k)/k * n^2 / (n-1) * sigma.2)
}
#
# Create an example.
#
K <- 30                                           # Desired sample size
n <- round((1:9)^(3/2) * 400/sum((1:9)^(3/2)))
n[which.max(n)] <- n[which.max(n)] + 400 - sum(n) # Stratum sizes
sigma.2 <- signif(200*exp(-(1:9)/3), 2)           # Estimated stratum variances
#
# Find the (non-integral) solution.
#
k <- ss(K, n, sigma.2)
#
# Find the best among the closest integral solutions (brute force).
#
x <- mapply(function(x,y) list(c(x,y)), floor(k), ceiling(k))
X <- as.matrix(expand.grid(x[k > 1]))
X <- X[rowSums(X) + sum(k<=1) ==K,]
Y <- matrix(1, nrow(X), length(k))
Y[, k>1] <- X
SS <- apply(Y, 1, function(k0) s.var(k0, n, sigma.2))
#
# Report and compare the results.
#
k.even <- round(ss(K, n, 1)) # Assume all sample stratum are equal (to 1)

#-- The various solutions
rbind(Rounded=Y[which.min(SS), ],
Original=round(k, 2),
Equal variance=k.even)

#-- The associated variances of the population total estimator
logsig <- function(x, d=2) 1 + signif(x-1, d) # Formats values close to 1
var.t <- s.var(k, n, sigma.2)
var.t.even <- s.var(k.even, n, sigma.2)
rbind(Rounded=c(Value=min(SS), Ratio=logsig(min(SS)/var.t)),
Original=c(var.t, 1),
Equal variance=c(var.t.even, logsig(var.t.even/var.t)))

• Thanks, it took me some time to process, but it is a very detailed explanation of how to tackle the problem. Mathematics is such a beautiful language. However, I am not quite sure how to read the infinity sign in equation (1). I should brush up on my math, but maybe you like to explain? – Eelco Aug 6 at 11:44
• That's not an infinity sign: it's a standard symbol for "is proportional to." – whuber Aug 6 at 12:03
• Thanks for your help! – Eelco Aug 6 at 12:04
• Thank you for your comments. I made some editorial changes to clarify that point and to correct a few typographical errors in key formulas. – whuber Aug 6 at 12:17