18 Sampling a Site or Region with Spatial Units (Ch. 18)
The realities of sampling are such that we often sample spatially, and as a result in practice deal with cluster samples more often than simple random samples.
18.1 Estimating Population Proportions
Calculating the standard error of a proportion from a cluster sample features a more complicated equation than simply calculating the standard error of a proportion, but needn’t in practice be difficult. That equation is \[SE=\sqrt{(\frac{1}{n})(\frac{\Sigma{(\frac{x}{y}-P)^2(\frac{yn}{Y})^2}}{n-1})(1-\frac{n}{N})}\] where SE is the standard error of the proportion, n is the sample size (the number of clusters sampled), N is the population size (the number of clusters in the population that could have been sampled), x is the quantity of object x in a given cluster, y is the quantity of object y in that same cluster, P is the estimate of the proportion \(\frac{x}{y}\) for the population, and Y is \(\Sigma{y}\) (the total quantity of object y in the population). Table 18.1 provides counts of sherds from ten randomly sampled excavation units, which Drennan uses to explore how to calculate the standard error of a proportion from a cluster sample.
sherds <- data.frame(Unit = c(7, 18, 29, 31,37, 56, 72, 83, 87, 91),
CMsherds = c(10, 13, 16, 19, 17, 21, 18, 30, 19, 12),
TotalSherds = c(32, 27, 38, 73, 55, 41, 63, 81, 56, 34))
kable_classic(kable(sherds), full_width = F)| Unit | CMsherds | TotalSherds |
|---|---|---|
| 7 | 10 | 32 |
| 18 | 13 | 27 |
| 29 | 16 | 38 |
| 31 | 19 | 73 |
| 37 | 17 | 55 |
| 56 | 21 | 41 |
| 72 | 18 | 63 |
| 83 | 30 | 81 |
| 87 | 19 | 56 |
| 91 | 12 | 34 |
Table 18.2, which contains the elements of the standard error calculation, we can produce through vectorized operations. Note that because we are defining some of these columns with reference to others, they have to be calculated in the sequence written (e.g., sherds$yn2 cannot be calculated until sherds$yn has been calculated). We could avoid this by writing more complicated code.
sherds$xy.prop <- round(sherds$CMsherds / sherds$TotalSherds, 3)
sherds$xy.prop.dev <- round((sherds$CMsherds / sherds$TotalSherds) -
(sum(sherds$CMsherds) / sum(sherds$TotalSherds)), 3)
sherds$xy.prop.dev2 <- round(sherds$xy.prop.dev^2, 6)
sherds$yn <- (sherds$TotalSherds*length(sherds$Unit)) / sum(sherds$TotalSherds)
sherds$yn2 <- (sherds$yn)^2
sherds$weights <- round(sherds$xy.prop.dev2 * sherds$yn2, 6)
head(sherds)## Unit CMsherds TotalSherds xy.prop xy.prop.dev xy.prop.dev2 yn yn2
## 1 7 10 32 0.312 -0.037 0.001369 0.64 0.4096
## 2 18 13 27 0.481 0.131 0.017161 0.54 0.2916
## 3 29 16 38 0.421 0.071 0.005041 0.76 0.5776
## 4 31 19 73 0.260 -0.090 0.008100 1.46 2.1316
## 5 37 17 55 0.309 -0.041 0.001681 1.10 1.2100
## 6 56 21 41 0.512 0.162 0.026244 0.82 0.6724
## weights
## 1 0.000561
## 2 0.005004
## 3 0.002912
## 4 0.017266
## 5 0.002034
## 6 0.017646Summing the last column provides the numerator for the middle term of the equation for calculating SE from cluster samples.
## [1] 0.053132That quantity - 0.053132 - can be plugged into the equation for calculating SE, along with values for n and N.
## [1] 0.02305038Calculating this way has the advantage of being easy to follow (and troubleshoot), but a function is much more efficient.
SECluster <- function(x, y, N, ...) {
P <- sum(x)/sum(y)
Y <- sum(y)
n <- length(x)
return(sqrt((1/n)*(sum(((x/y - P)^2*(y*n/Y)^2), na.rm=T)/(n - 1))*(1 - n/N)))
}
SECluster(sherds$CMsherds, sherds$TotalSherds, 100)## [1] 0.0230664318.2 Estimating Population Means
Estimating population means from cluster samples involves a similar procedure. The equation for standard error of the mean is \[SE=\sqrt{(\frac{1}{n})(\frac{\Sigma{(\bar{x}-\bar{X})^2(\frac{yn}{Y})^2}}{n-1})(1-\frac{n}{N})}\] where SE is the standard error of the mean, n is the number of clusters in the sample, N is the number of clusters in the population, \(\bar{x}\) is the mean of x in a cluster, \(\bar{X}\) is the estimated population mean of x (the overall sample mean), y is the quantity of objects in a cluster for which x has been measured, and Y is \(\Sigma{y}\) (the total quantity of y in all the clusters in the sample).
Table 18.3 provides some sample data, in a format that is a bit complicated: the second column (projectile point lengths) is itself a vector, with between 1 and 5 elements. We can produce this in R using a list, which is a useful thing to be familiar with.
ProjPointLengths <- data.frame(Unit = sherds$Unit, Lengths =
I(list(c(15, 19, 23), 17, c(18, 23),
c(18, 18, 27), c(18, 19), 24,
c(20, 21, 26, 28, 29), 16, 28, c(25, 26))))
ProjPointLengths## Unit Lengths
## 1 7 15, 19, 23
## 2 18 17
## 3 29 18, 23
## 4 31 18, 18, 27
## 5 37 18, 19
## 6 56 24
## 7 72 20, 21, ....
## 8 83 16
## 9 87 28
## 10 91 25, 26You can use this data frame to build Table 18.4, but we’ll skip straight to calculating standard error - with the exception of a quick note about the complications introduced by having a list within a data frame. We can no longer rely on simple vectorization (try mean(ProjPointLengths$Lengths)), because that column is a list rather than a simple vector; instead we have to use sapply() or something like it.
ProjPointLengths$Y <- sapply(ProjPointLengths$Lengths, length)
ProjPointLengths$xbar <- sapply(ProjPointLengths$Lengths, mean)
head(ProjPointLengths)## Unit Lengths Y xbar
## 1 7 15, 19, 23 3 19.0
## 2 18 17 1 17.0
## 3 29 18, 23 2 20.5
## 4 31 18, 18, 27 3 21.0
## 5 37 18, 19 2 18.5
## 6 56 24 1 24.0Rather than produce the remaining columns of Table 18.4 (and please note that the second column in Table 18.4 on Drennan p248 should be y, not Y), we’ll opt for a function that will take care of that for us.
SEmeanCluster <- function(x, N) {
y <- sapply(x, length)
xbar <- sapply(x, mean)
Xbar <- mean(unlist(x))
n <- length(x)
Y <- length(unlist(x))
return(round(sqrt((1/n)*(sum(((xbar - Xbar)^2*(y*n/Y)^2)/(n - 1))*(1 - n/N))),2))
}
ProjPointLengths_SE <- SEmeanCluster(ProjPointLengths$Lengths, 100); ProjPointLengths_SE## [1] 1.07We can use that standard error to calculate the error range in our estimate of the mean at a given confidence level.
## [1] 2.420508We can be 95% confident that the mean length of all the projectile points at the site is 21.8mm ± 2.4mm.
18.3 Densities
Note that whether a sample can be treated as a simple random sample or a cluster sample depends on the unit of analysis. As Drennan notes, for the example above a measure like density of sherds/unit could be treated as a simple random sample, because density is a property of the unit (e.g., the grid square), and since those units were selected randomly they can be treated as a simple random sample. You can follow Drennan’s logic (and math) on p249 using the methods covered in Ch.9.