Is R's weighted sample without replacement function misleading?

Background

The 2023 article "Remarks on some misconceptions about unequal probability sampling without replacement" by Tillé suggests the sample function in R is misleading, when sampling without replacement and with unequal probabilities because the inputted vector of probability weights does not agree with the inclusion probabilities. It provides code to show this. It also recommends to use the R package sampling, co-authored by Tillé, and one of its many functions to sample without replacement and with unequal probabilities.

The documentation of sample states that

If replace is false, these probabilities are applied sequentially, that is the probability of choosing the next item is proportional to the weights amongst the remaining items. The number of nonzero weights must be at least size in this case*".

Tillé argues that some users assume, incorrectly, that this sequential procedure will lead to inclusion probabilities which agree with the inputted vector of probability weights.

Questions

• Does the base R sample function respect the inclusion probabilities when sampling without replacement and with unequal probabilities?
• If not, what are the inclusion probabilities for an inputted vector of probability weights?
• Are frequent R users aware of this discrepancy?
• In what situation would you actually prefer the method of the sample function? (the function seems faster than those from the package sampling, so it may be beneficial approximate inclusion probabilities are good enough).
• Can the inputted probabilities weights be transformed, prior to sampling, to achieve the expected inclusion probabilities?

There was work at the R Sprint earlier this year on adding an option to specify marginal sampling probabilities for unequal-probability sampling without replacement. It hasn't happened yet, but is likely to happen in the future.

One challenge is finding a method that doesn't break down when given larger sampling problems. There are many approaches that work at the scale needed for survey sampling but that don't scale well to larger sample and population sizes. For example, Tillé's method in UPtille accumulates rounding error for larger problems. Another issue is whether we need to use a method that produces non-zero pairwise sampling probabilities.

We may also add Poisson sampling as another option: it is fast and reproduces the input probabilities but doesn't fix the sample size.

The code below is adapted from the 2023 article "Remarks on some misconceptions about unequal probability sampling without replacement" by Tillé and generates multiple samples using the base R sample function and the UPbrewer function from the package sampling. The inclusion probabilities are determined after sampling and compared to the inputted vector of probability weights in the first plot. If the inclusion probabilities are correct, any residual discrepancies should be due to the finite sampling of this simulation and would follow a normal distribution, shown in the second graph.

I first sample with function UPbrewer from the package sampling. The function UPtille gives similar results but is over twice as slow.

library(sampling)
library(ggplot2)
set.seed(123)
N <- 20;
n <- 8;
ptm <- proc.time()
piks <- c()
pikests <- c()
distances <- c()
for(i in seq(100)){
pik <- runif(N, 0.2)
pik <- pik / sum(pik) * n
pikest <- rep(0, N)
SIM <- 12000
for(j in seq(SIM)){
s <- which(UPbrewer(pik) == 1)
pikest[s] <- pikest[s] + 1 / SIM
}
piks <- c(piks, pik)
pikests <- c(pikests, pikest)
distance <- (pikest - pik) / sqrt(pik * (1 - pik) / SIM)
distances <- c(distances, distance)
}
proc.time() - ptm # 105.66
ggplot() +
geom_point(aes(x = piks, y = pikests)) +
geom_smooth(aes(x = piks, y = pikests)) +
geom_abline(aes(intercept = 0, slope = 1), col = "red") +
expand_limits(x = c(0, 1), y = c(0, 1))
mean(distances) # -0.0009655566
sd(distances) # 0.9993592
ggplot() +
geom_density(aes(x = distances)) +
stat_function(fun = dnorm, col = "red")


I now sample with the function sample in base R.

library(sampling)
library(ggplot2)
set.seed(123)
N <- 20;
n <- 8;
ptm <- proc.time()
piks <- c()
pikests <- c()
distances <- c()
for(i in seq(100)){
pik <- runif(N, 0.2)
pik <- pik / sum(pik) * n
pikest <- rep(0, N)
SIM <- 12000
for(j in seq(SIM)){
s <- sample(x = seq(N), size = n, replace = FALSE, prob = pik)
pikest[s] <- pikest[s] + 1 / SIM
}
piks <- c(piks, pik)
pikests <- c(pikests, pikest)
distance <- (pikest - pik) / sqrt(pik * (1 - pik) / SIM)
distances <- c(distances, distance)
}
proc.time() - ptm # 26.78
ggplot() +
geom_point(aes(x = piks, y = pikests)) +
geom_smooth(aes(x = piks, y = pikests)) +
geom_abline(aes(intercept = 0, slope = 1), col = "red") +
expand_limits(x = c(0, 1), y = c(0, 1))
mean(distances) # 0.3682103
sd(distances) # 7.909408
ggplot() +
geom_density(aes(x = distances)) +
stat_function(fun = dnorm, col = "red")


The inclusion probabilities indeed do not agree with the inputted vector of probability weights when using the function sample. As mentioned in the article:

This method does not respect the inclusion probabilities. Moreover, the computation of the design and of the inclusion probabilities requires a sum over all the permutations of the sample, which quickly becomes impossible when the sample size is large. ...these authors all realized that this intuitive method is false. They all tried to find solutions, which marked the beginning of a prolific research in the field of the sampling methods.

I do not know if most users are aware that the sequential sampling does not respect the inclusion probabilities, but the R manual is not very explicit about it. I also do not know in what situations this type of sampling is beneficial.

Could the inputted probabilities weights be transformed, prior to sampling, to achieve the expected inclusion probabilities? Possibly. The code below fits a flexible polynomial through the previous output and uses this to transform the inputted vector. Following this transformation, the outputted inclusion probabilities do seem to agree with the original input. I do not know if this transformation can be obtained analytically, but it shows that, at least in principle, the correct inclusion probabilities can be recovered.

transformation.model <- lm(pik ~ poly(pikest, 4))
ptm <- proc.time()
piks <- c()
pikests <- c()
distances <- c()
for(i in seq(100)){
pik <- runif(N, 0.2)
pik <- pik / sum(pik) * n
transformed.pik <- predict(transformation.model,
newdata = data.frame("pikest" = pik))
pikest <- rep(0, N)
SIM <- 12000
for(j in seq(SIM)){
s <- sample(x = seq(N), size = n, replace = FALSE,
prob = transformed.pik)
pikest[s] <- pikest[s] + 1 / SIM
}
piks <- c(piks, pik)
pikests <- c(pikests, pikest)
distance <- (pikest - pik) / sqrt(pik * (1 - pik) / SIM)
distances <- c(distances, distance)
}
proc.time() - ptm # 29.97
ggplot() +
geom_point(aes(x = piks, y = pikests)) +
geom_smooth(aes(x = piks, y = pikests)) +
geom_abline(aes(intercept = 0, slope = 1), col = "red") +
expand_limits(x = c(0, 1), y = c(0, 1))
mean(distances) # 0.09198351
sd(distances) # 1.622353
ggplot() +
geom_density(aes(x = distances)) +
stat_function(fun = dnorm, col = "red")