The problem
BiasedMatch implements the two-urn biased sampling model
of Puza and Bonfrer (2018). Urn 1 contains N1 distinct
items and urn 2 contains N2 distinct items, with
N items common to both urns. A sample of size
m1 is drawn without replacement from urn 1 and a sample of
size m2 from urn 2. Each item common to both urns is given
a selection weight w, while every other item has weight 1,
so w > 1 makes common items more likely to be selected.
The quantity of interest is m, the number of items selected
in both samples (the number of matches).
The model has several natural readings: two retailers independently
choosing overlapping products to promote, two committees drawn from a
shared population of candidates, or classic capture-recapture sampling.
In each case the number of matches, and any bias toward the common
items, is what we want to understand.
The distribution of the number of matches
PROBM() returns the probability of a specific number of
matches, and DISTM() returns the full distribution together
with its mean and variance. When w = 1 the model reduces to
unbiased simple random sampling without replacement (a mixture of
hypergeometrics, Theorem 1 of the paper); when w > 1 it
uses the univariate and bivariate Wallenius’ non-central hypergeometric
distributions (Theorem 2).
The example below reproduces Table 1 of the paper, contrasting
w = 1 with w = 2 for N1 = 37,
N2 = 45, N = 16, m1 = 12,
m2 = 8.
d1 <- DISTM(N1 = 37, N2 = 45, N = 16, m1 = 12, m2 = 8, w = 1)
d2 <- DISTM(N1 = 37, N2 = 45, N = 16, m1 = 12, m2 = 8, w = 2)
round(
data.frame(
m = d1$mvec,
`w = 1` = d1$fmvec,
`w = 2` = d2$fmvec,
check.names = FALSE
),
4
)
#> m w = 1 w = 2
#> 1 0 0.3573 0.1182
#> 2 1 0.4137 0.3164
#> 3 2 0.1836 0.3309
#> 4 3 0.0404 0.1750
#> 5 4 0.0047 0.0507
#> 6 5 0.0003 0.0081
#> 7 6 0.0000 0.0007
#> 8 7 0.0000 0.0000
#> 9 8 0.0000 0.0000
Raising w from 1 to 2 shifts the distribution of matches
upward: the mean rises from 0.9225 to 1.7511, and the mode moves from 1
to 2.
plot(d1$mvec, d1$fmvec,
type = "b", pch = 16, ylim = c(0, 0.45),
xlab = "Number of matches (m)", ylab = "Probability"
)
lines(d2$mvec, d2$fmvec, type = "b", pch = 17, lty = 2)
legend("topright", c("w = 1", "w = 2"), pch = c(16, 17), lty = c(1, 2), bty = "n")
For the unbiased case w = 1, Vmf() provides
the closed-form variance of Theorem 1(b), which agrees with the exact
variance from DISTM():
c(closed_form = Vmf(N1 = 37, N2 = 45, N = 16, m1 = 12, m2 = 8), exact = d1$Vm)
#> closed_form exact
#> 0.7441474 0.7441474
Inference on the bias parameter
The inferential problem “inverts” the model: given the fixed urn
constants and a single observed number of matches m, what
can we say about the bias w? inferW() returns
the maximum likelihood estimate, the posterior mode under the
Jeffreys-type prior f(w) proportional to
1 / w, and a p-value for testing w = 1 (no
bias) against w > 1. likW() returns the
standardised likelihood and posterior curves for plotting.
Marketing application (Problem 3)
Masked promotion data from two competing electronics retailers give
N1 = 292, N2 = 350, N = 192,
m1 = 90, m2 = 59, with m = 13
matching promotions observed.
fit <- inferW(m = 13, N1 = 292, N2 = 350, N = 192, m1 = 90, m2 = 59)
fit
#> $mle
#> [1] 1.526549
#>
#> $posterior_mode
#> [1] 1.302926
#>
#> $p_value
#> [1] 0.1736506
#>
#> $alternative
#> [1] "greater"
#>
#> $m
#> [1] 13
curves <- likW(
m = 13, N1 = 292, N2 = 350, N = 192, m1 = 90, m2 = 59,
wgrid = seq(0.05, 8, by = 0.05)
)
plot(curves$w, curves$likelihood,
type = "l",
xlab = "w", ylab = "standardised value"
)
lines(curves$w, curves$posterior, lty = 2)
abline(v = c(fit$mle, fit$posterior_mode), col = "grey50", lty = 3)
legend("topright", c("likelihood", "posterior kernel"), lty = c(1, 2), bty = "n")
The MLE is about 1.53 and the posterior mode about 1.3, but the
p-value of 0.174 gives only weak evidence against w = 1. At
the 5% level we would not conclude that common items are preferentially
selected for promotion.
Committee gender bias (Problem 2)
Two committees of sizes C = 9 and K = 6 are
drawn from a common population of M = 22 men and
W = 30 women; we observe B = 2 men on both
committees. This maps onto the two-urn model with
N1 = N2 = M + W = 52, N = M = 22,
m1 = C, m2 = K, and w measuring
bias toward selecting men.
inferW(m = 2, N1 = 52, N2 = 52, N = 22, m1 = 9, m2 = 6)
#> $mle
#> [1] 16.9866
#>
#> $posterior_mode
#> [1] 2.291593
#>
#> $p_value
#> [1] 0.0609043
#>
#> $alternative
#> [1] "greater"
#>
#> $m
#> [1] 2
Sweeping over every possible value of B reproduces Table
2 of the paper. Note that once B is large enough the
likelihood is strictly increasing and the MLE becomes infinite, while
the posterior mode stays finite.
B_values <- 0:6
tab <- do.call(rbind, lapply(B_values, function(B) {
r <- inferW(m = B, N1 = 52, N2 = 52, N = 22, m1 = 9, m2 = 6)
data.frame(B = B, MLE = r$mle, posterior_mode = r$posterior_mode, p_value = r$p_value)
}))
round(tab, 4)
#> B MLE posterior_mode p_value
#> 1 0 0.0001 0.0001 1.0000
#> 2 1 2.5635 0.6447 0.3735
#> 3 2 16.9866 2.2916 0.0609
#> 4 3 Inf 5.1150 0.0047
#> 5 4 Inf 9.1058 0.0002
#> 6 5 Inf 13.4722 0.0000
#> 7 6 Inf 17.6357 0.0000
Hierarchical inference with MCMC
When several matching experiments are available, the weight can be
treated hierarchically and estimated by the Metropolis-Hastings samplers
MHALG() (single shared weight), RMMHALG()
(weight depends on covariates) and SMRMMHALG() (weights
shared within covariate-defined classes). restable()
summarises the draws and estprop() reports
serial-correlation and convergence diagnostics. Use many more iterations
than shown here for real inference.
set.seed(4172)
res <- MHALG(
B = 500, J = 5000,
N1vec = 37, N2vec = 45, N0vec = 16,
M1vec = 12, M2vec = 8, mvec = 3,
Cvec = 1, del = 0.1
)
restable(res, Bn = 500)
estprop(res, Bn = 500, Jn = 4000)