A simulated dataset in which source 2 recorded a count on every unit and
source 1 on a random 20% of units in each site. It is generated from the
model of binegbin_partialobs() with site effects on congruence and bias,
so every example can be run, and its estimates compared with known values.
The latent components behind each pair are in bicount_sim_latent.
Format
A data frame with 540 rows (three sites of 180 units) and 5 columns:
- site
Factor,
A,BorC.- unit
Integer, the unit within its site, 1 to 180.
- y1
Integer, the count recorded by source 1;
NAwhere source 1 did not record the unit.- y2
Integer, the count recorded by source 2, on every unit.
- y1_obs
Integer, 1 where
y1was recorded and 0 where it was not.
Details
Each unit's counts were constructed by trivariate reduction,
y1 = Ns + N1 and y2 = Ns + N2, from three independent negative-binomial
components with means mu = M f, lambda1 = M (1 - f)(1 + beta) and
lambda2 = M (1 - f)(1 - beta), and inverse dispersions phi = 1 / kappa^2.
Here M is the midpoint of the two expected counts, f the congruence, and
beta = tanh(delta) the bounded source bias. The generating values were:
| site | M | f | beta | kappas | kappax1 | kappax2 |
| A | 5 | 0.5 | 0 | 0.5 | 0.5 | 0.5 |
| B | 5 | 0.7 | 0.3 | 0.5 | 0.5 | 0.5 |
| C | 5 | 0.9 | -0.3 | 0.5 | 0.5 | 0.5 |
Source 1 was recorded on exactly 36 units per site, drawn at random within
site, so whether y1 was recorded does not depend on either count.
Generated by data-raw/bicount_sim.R in the package source, with
set.seed(20261002).
Before fitting with binegbin_partialobs(), replace the NA in y1 with a
placeholder such as 0: brms drops rows with a missing response, which would
discard the observed y2 on those rows. The likelihood does not read y1
where y1_obs is 0.
Examples
str(bicount_sim)
#> 'data.frame': 540 obs. of 5 variables:
#> $ site : Factor w/ 3 levels "A","B","C": 1 1 1 1 1 1 1 1 1 1 ...
#> $ unit : int 1 2 3 4 5 6 7 8 9 10 ...
#> $ y1 : int NA NA NA 5 NA 6 NA NA NA NA ...
#> $ y2 : int 6 6 7 3 4 8 5 4 12 4 ...
#> $ y1_obs: int 0 0 0 1 0 1 0 0 0 0 ...
table(bicount_sim$site, bicount_sim$y1_obs)
#>
#> 0 1
#> A 144 36
#> B 144 36
#> C 144 36
# Every recorded pair is the sum of its latent components
rec <- bicount_sim$y1_obs == 1
all(bicount_sim$y1[rec] ==
bicount_sim_latent$Ns[rec] + bicount_sim_latent$N1[rec])
#> [1] TRUE
if (FALSE) { # \dontrun{
# Fitting compiles a Stan model and samples, which takes minutes.
# The model is written in midpoint M, congruence f and source bias through
# non-linear formulas, with logM = log M, logitf = logit f,
# delta = artanh(beta), the three overdispersions kappa, and site effects on
# f and delta.
library(brms)
dat <- bicount_sim
dat$y1[dat$y1_obs == 0] <- 0L
fit <- brm(
bf(y1 | vint(y2, y1_obs) ~ 1, nl = TRUE) +
nlf(mu ~ logM + log_inv_logit(logitf)) +
nlf(lambdaone ~ log(2) + logM + log_inv_logit(-logitf) +
log_inv_logit(2 * delta)) +
nlf(lambdatwo ~ log(2) + logM + log_inv_logit(-logitf) +
log_inv_logit(-2 * delta)) +
nlf(shapes ~ log(1 / kappas^2)) +
nlf(shapexone ~ log(1 / kappaxone^2)) +
nlf(shapextwo ~ log(1 / kappaxtwo^2)) +
lf(logM ~ 1) + lf(logitf ~ 0 + site) + lf(delta ~ 0 + site) +
lf(kappas ~ 1) + lf(kappaxone ~ 1) + lf(kappaxtwo ~ 1),
data = dat,
family = binegbin_partialobs(),
stanvars = binegbin_partialobs_stanvars(),
prior = c(
set_prior("normal(0.5, 2)", class = "b", nlpar = "logM"),
set_prior("normal(0, 1.5)", class = "b", nlpar = "logitf"),
set_prior("normal(0, 0.5)", class = "b", nlpar = "delta"),
set_prior("normal(0, 1)", class = "b", nlpar = "kappas", lb = 0),
set_prior("normal(0, 1)", class = "b", nlpar = "kappaxone", lb = 0),
set_prior("normal(0, 1)", class = "b", nlpar = "kappaxtwo", lb = 0)
),
chains = 4, cores = 4, seed = 1
)
# Congruence by site, against 0.5, 0.7 and 0.9
cols <- c("Estimate", "Q2.5", "Q97.5")
plogis(fixef(fit)[c("con_siteA", "con_siteB", "con_siteC"), cols])
# Bias beta = tanh(delta) by site, against 0, 0.3 and -0.3
tanh(fixef(fit)[c("methd_siteA", "methd_siteB", "methd_siteC"), cols])
# The unrecorded source-1 counts, drawn conditionally on each unit's y2
y1_draws <- posterior_predict(fit)[, dat$y1_obs == 0]
} # }