Introduction to mediaPlanR

Units and the models of random duplication

Every function takes its inputs in explicit units. An audience is the number of people exposed by one insertion of a vehicle, and the population is the planning universe. The audiences of all the vehicles must measure the same opportunity of exposure and the same period: contacts are opportunities of exposure, and they can be read as advertising impressions only if the audiences measure that event (see the section “Notation and units” of ?mediaPlanR).

The simplest models need only the audience of each vehicle and the number of insertions:

result <- calc_sainsbury(
  audiences = c(300000, 180000, 120000),
  population = 1000000,
  insertions = c(4, 6, 10)
)
result
#> SAINSBURY MODEL
#> ===============
#> Description: random duplication and accumulation, with heterogeneous vehicles
#> 
#> HEADLINE METRICS:
#> -----------------
#> Total reach: 97.97% (979672 people)
#> 
#> EXPOSURE DISTRIBUTION:
#> ----------------------
#> (Percentage of the population receiving exactly N exposures)
#> 1 exposure: 8.93% (89343 people)
#> 2 exposures: 18.40% (184035 people)
#> 3 exposures: 23.63% (236318 people)
#> 4 exposures: 21.22% (212238 people)
#> 5 exposures: 14.18% (141765 people)
#> 6 exposures: 7.31% (73105 people)
#> 7 exposures: 2.98% (29815 people)
#> 8 exposures: 0.98% (9771 people)
#> 9 exposures: 0.26% (2600 people)
#> 10 exposures: 0.06% (565 people)
#> 11 exposures: 0.01% (100 people)
#> 12 exposures: 0.00% (15 people)
#> 13 exposures: 0.00% (2 people)
#> 14 exposures: 0.00% (0 people)
#> 15 exposures: 0.00% (0 people)
#> 16 exposures: 0.00% (0 people)
#> 17 exposures: 0.00% (0 people)
#> 18 exposures: 0.00% (0 people)
#> 19 exposures: 0.00% (0 people)
#> 20 exposures: 0.00% (0 people)
#> 
#> CUMULATIVE DISTRIBUTION:
#> -------------------------
#> (Percentage of the population receiving N or more exposures)
#> >= 1 exposure: 97.97% (979672 people)
#> >= 2 exposures: 89.03% (890329 people)
#> >= 3 exposures: 70.63% (706294 people)
#> >= 4 exposures: 47.00% (469976 people)
#> >= 5 exposures: 25.77% (257738 people)
#> >= 6 exposures: 11.60% (115973 people)
#> >= 7 exposures: 4.29% (42868 people)
#> >= 8 exposures: 1.31% (13053 people)
#> >= 9 exposures: 0.33% (3282 people)
#> >= 10 exposures: 0.07% (682 people)
#> >= 11 exposures: 0.01% (117 people)
#> >= 12 exposures: 0.00% (16 people)
#> >= 13 exposures: 0.00% (2 people)
#> >= 14 exposures: 0.00% (0 people)
#> >= 15 exposures: 0.00% (0 people)
#> >= 16 exposures: 0.00% (0 people)
#> >= 17 exposures: 0.00% (0 people)
#> >= 18 exposures: 0.00% (0 people)
#> >= 19 exposures: 0.00% (0 people)
#> >= 20 exposures: 0.00% (0 people)
#> 
#> SUMMARY STATISTICS:
#> --------------------
#> Average exposures per person reached: 3.55

The distribution has one value for each number of exposures from one to the total number of insertions, and its sum is the reach. The average frequency is the number of contacts divided by the people reached.

sum(result$distribution$percent)
#> [1] 97.96717
result$reach$percent
#> [1] 97.96717
sum(seq_along(result$distribution$people) * result$distribution$people) /
  result$reach$people
#> [1] 3.552211

calc_sainsbury() keeps each vehicle’s own audience (heterogeneous vehicles); calc_binomial() gives every insertion the insertion-weighted mean probability (homogeneous vehicles). Both assume random duplication and random accumulation: a repeat insertion in the same vehicle is treated as an independent trial, exactly like an insertion in a different vehicle. This is a simplification, and it is the hypothesis that the cited models state explicitly (Lee & Burkart, 1960; Chandon, 1976, as reviewed in Aldás Manzano, 1998, Sections 3.3.1.1 and 3.3.1.2, where it is cited as Chandon, 1985). It is also what allows them to run without any observed duplication data. With one insertion per vehicle, both reduce to their one-insertion-per-vehicle formulation (Aldás Manzano, 1998, Sections 3.2.2.1 and 3.2.2.2). The other models need observed duplications or the reach after two insertions, and take their own inputs, as the rest of this vignette shows.

The family of reach and exposure-distribution models

All these models answer the same question: what proportion of the population is reached, and how many times, by a given media plan? They differ in the shape of plan they were derived for, and in their assumptions about the duplication of audiences across vehicles and about the heterogeneity of individual exposure probabilities. Plan shape follows Aldás Manzano’s (1998) three-way split. The two Hofmans models carry their domain in the function name (calc_hofmans_accumulation(), calc_hofmans_duplication()) because they are two applications of one basis, Agostini’s formula, to different plan shapes, with different data, and they are not interchangeable (they coincide only when the vehicles are equal and every pair duplicates equally); calc_beta_binomial() is also an accumulation model, although its name, which is the standard one in the literature, does not say so. Accumulation models (Beta-Binomial, Hofmans accumulation) are derived for one vehicle with several insertions. Duplication models (Agostini, Hofmans duplication) are derived for several vehicles with one insertion each. Net cumulative audience models cover several vehicles with several insertions. Sainsbury and Binomial cover both simple cases because each is one model at two levels of generality (see the insertions argument): they accept one or several insertions in each vehicle.

What the models assume

Each model rests on assumptions about four things; the table below summarizes them for every model, and ?mediaPlanR gives the same explanation.

Individuals. Homogeneous: every person has the same probability of being exposed to an insertion of a vehicle. Heterogeneous: the probability differs between people, as between a loyal reader and an occasional one. With heterogeneity the same people tend to be exposed again and again, so reach grows more slowly with the number of insertions than under homogeneity.
Vehicles. Homogeneous: all vehicles have the plan’s average audience. Heterogeneous: each vehicle keeps its own audience.
Duplication between vehicles. Random: people choose vehicles independently, so the audience shared by two vehicles is the product of their audiences (two vehicles with audiences of 30% and 20% would share 6%). Observed: the shared audience is taken from data.
Accumulation in the same vehicle. Random: a repeat insertion in a vehicle is as independent as an insertion in another vehicle. Observed: the reach after two insertions (R2) shows how much a vehicle’s audience renews.

All the models also assume that exposure probabilities are stationary: they do not change over the period of the plan.

Model	Function	Plan shape	Duplication	Heterogeneity
Sainsbury	`calc_sainsbury()`	1 or more insertions, m vehicles	Random (independence and accumulation)	Vehicles only
Binomial	`calc_binomial()`	1 or more insertions, m vehicles	Random (independence and accumulation)	None (average vehicle)
Agostini (duplication)	`calc_agostini_duplication()`	1 insertion, m vehicles	Observed pairwise, one coefficient `k`	Not modeled (reach only)
Hofmans (duplication)	`calc_hofmans_duplication()`	1 insertion, m vehicles	Observed pairwise, pairwise coefficient	Not modeled (reach only)
Beta-Binomial	`calc_beta_binomial()`	1 vehicle, n insertions	Implicit (fitted to `R1`, `R2`)	Individuals (Beta-Binomial)
Hofmans (accumulation)	`calc_hofmans_accumulation()`	1 vehicle, n insertions	Implicit (fitted to `R1`, `R2`, optionally `R3`)	Not modeled (cumulative reach only)
Metheringham	`calc_metheringham()`	n insertions, m vehicles	Observed, averaged over insertion pairs	Individuals (Beta-Binomial)
CANEX	`calc_canex()`	n insertions, m vehicles	Observed, canonical correlations	Per vehicle (Beta-Binomial)
CBD (conditional Beta)	`calc_cbd()`	n insertions, m vehicles	Observed; canonical expansion on the zero/one grid of first-insertion exposure	Per vehicle (Beta-Binomial)
Beta-Binomial fitted to a reach	`fit_bbd_to_reach()`	n insertions, m vehicles	Reach supplied as a constraint	Individuals (Beta-Binomial)
CSD (canonical sequential aggregation)	`calc_csd()`	n insertions, m vehicles	Observed; canonical reach per sub-schedule	Per vehicle (Beta-Binomial)
MSAD (Morgensztern)	`calc_msad()`	n insertions, m vehicles	Observed; Morgensztern reach per sub-schedule	Per vehicle (Beta-Binomial)
MBD (multivariate Beta-Binomial)	`calc_mbd()`	n insertions, m vehicles	Observed pairs; imputed for 3+ vehicles (Waring’s theorem)	Per vehicle (Beta-Binomial)
Negative-Binomial approximation	`nbd_exposure_distribution()`	Continuous process, not a finite schedule	None across vehicles	Poisson-Gamma, unbounded

Example: three vehicles, one insertion each, different assumptions

Sainsbury and Binomial need only the audiences. Agostini and Hofmans also use the audience duplicated between each pair of vehicles, and Metheringham uses the average duplication.

data(ratings_example)
data(duplication_example)

reach_sainsbury <- do.call(calc_sainsbury, ratings_example)
reach_binomial <- do.call(calc_binomial, ratings_example)
reach_agostini <- do.call(calc_agostini_duplication, duplication_example)
reach_hofmans <- do.call(calc_hofmans_duplication, duplication_example)
reach_metheringham <- calc_metheringham(
  audiences = duplication_example$audiences,
  insertions = c(1, 1, 1),
  duplication_matrix = duplication_example$duplication_matrix,
  population = duplication_example$population
)

data.frame(
  model = c("Sainsbury", "Binomial", "Agostini (k = 1.125)", "Hofmans (duplication)",
            "Metheringham"),
  reach_pct = round(c(
    reach_sainsbury$reach$percent,
    reach_binomial$reach$percent,
    reach_agostini$reach$percent,
    reach_hofmans$reach$percent,
    reach_metheringham$reach$percent
  ), 2)
)
#>                   model reach_pct
#> 1             Sainsbury     66.40
#> 2              Binomial     65.70
#> 3  Agostini (k = 1.125)     65.45
#> 4 Hofmans (duplication)     64.17
#> 5          Metheringham     63.64

The observed duplication in duplication_example is higher than random duplication would produce for every pair of vehicles, so the models that use it estimate a lower reach than Sainsbury and Binomial. Agostini’s formula applies one empirical coefficient (k = 1.125, Agostini’s own estimate) to the total duplication, whereas Hofmans’ formula computes a coefficient from each pair’s own audiences, so no coefficient has to be fitted. Both are ad hoc formulas: they return reach only, and calc_agostini_duplication() and calc_hofmans_duplication() warn when the reach they compute leaves its logical range.

From a reach to a distribution: Hofmans Beta-Binomial (HBBD)

Hofmans’ formula returns only the reach. The Hofmans Beta-Binomial (HBBD) procedure of Leckenby and Boyd (Aldás Manzano, 1998, Section 3.3.1.10; Kim, 2005, pp. 52-53) turns that reach into a frequency distribution in a second step: a Beta-Binomial whose mean is the mean exposure implied by the audiences and whose reach equals the Hofmans reach. The package offers the two steps as separate functions, calc_hofmans_duplication() and fit_bbd_to_reach(), and they are composed as follows (one insertion per vehicle, as in Hofmans’ duplication formula):

# Step 1: reach of the plan from the audiences and the pairwise duplications
hofmans <- do.call(calc_hofmans_duplication, duplication_example)

# Step 2: Beta-Binomial with that reach and the mean implied by the audiences
hbbd <- fit_bbd_to_reach(
  insertions = rep(1, length(duplication_example$audiences)),
  audiences = duplication_example$audiences,
  reach = hofmans$reach$people,
  population = duplication_example$population
)

hbbd$parameters$converged
#> [1] TRUE
round(hbbd$distribution, 4)
#>   contacts probability cumulative_probability
#> 1        0      0.3583                 1.0000
#> 2        1      0.4172                 0.6417
#> 3        2      0.1906                 0.2245
#> 4        3      0.0338                 0.0338

The distribution gives the share of the population exposed 0, 1, 2 and 3 times, from which the average frequency and the effective reach are read. Two cautions apply. Hofmans’ formula is empirical and can return a reach outside the logical range of the plan (with a warning), in which case fit_bbd_to_reach() stops rather than fitting an impossible reach. And the composition follows the published description of the procedure, but it has not been validated against a published numeric example (see the table in the Scope section).

Several insertions per vehicle: Metheringham

calc_metheringham() averages audiences and duplications over all pairs of insertions, both between vehicles and within a vehicle, to define an “average” vehicle, and fits a Beta-Binomial to the plan’s N = sum(insertions) insertions (Aldás Manzano, 1998, Section 3.3.1.5). The diagonal of the duplication matrix holds the audience duplicated between two insertions in the same vehicle.

data(metheringham_example)
metheringham <- do.call(calc_metheringham, metheringham_example)
metheringham$reach$percent
#> [1] 74.34177
length(metheringham$distribution$percent)   # zero to N = 12 exposures
#> [1] 12

# The mean of the distribution equals the plan's gross exposures per person
sum(seq_along(metheringham$distribution$percent) *
      metheringham$distribution$percent) / 100
#> [1] 1.44
with(metheringham_example, sum(audiences * insertions) / population)
#> [1] 1.44

One vehicle, several insertions: Beta-Binomial and Hofmans

The Beta-Binomial accumulation model needs the reach after one and two insertions. Hofmans’ accumulation model, with a constant coefficient, needs the same two figures; an observed third insertion (R3) lets the coefficient vary with the number of insertions.

beta_binomial <- calc_beta_binomial(R1 = 0.50, R2 = 0.55, insertions = 5,
                                    population = 1000000)
beta_binomial$reach$percent
#> [1] 60.33654

hofmans <- calc_hofmans_accumulation(R1 = 0.06, R2 = 0.103, insertions = 5)
hofmans$results
#>   insertions        RN
#> 1          1 0.0600000
#> 2          2 0.1030000
#> 3          3 0.1353285
#> 4          4 0.1605195
#> 5          5 0.1807018

Beta-Binomial versus Negative-Binomial: finite opportunities and count processes

Both calc_beta_binomial() and the Negative-Binomial distribution assume that each individual has a personal exposure probability (or rate) that varies across the population, but they mix different distributions. The Beta-Binomial mixes a Binomial with a fixed number of insertions with a Beta; the Negative-Binomial mixes a Poisson with a Gamma (Ehrenberg, 1959). The Negative-Binomial has unbounded support, so it suits page views, ad-server exposures and other processes with continuous opportunities. It should not be used to suggest that a finite-insertion schedule can register more exposures than insertions.

When individual counts are observed, fit_nbd_exposure() estimates the mean and the heterogeneity by maximum likelihood. For scenarios with external parameters, nbd_exposure_distribution() keeps the last cell as an open tail and never presents it as an exact count. Danaher (2007) develops the multivariate Negative-Binomial model that represents dependence across websites; the functions here are deliberately univariate.

counts <- c(rep(0, 40), rep(1, 25), rep(2, 15), rep(3, 8), 5, 7)
nbd_fit <- fit_nbd_exposure(counts)
nbd_fit$estimates
#> $mean_contacts
#> [1] 1.011111
#> 
#> $size
#> [1] 2.010449
#> 
#> $size_interval
#>     lower     upper 
#> 0.7361054 5.4909324 
#> 
#> $variance
#> [1] 1.519627

nbd_scenario <- nbd_exposure_distribution(
  mean_contacts = 2.5, size = 1.7, report_max = 9, opportunities = 9
)
nbd_scenario$diagnostics$probability_above_opportunities
#> [1] 0.01888247

The last diagnostic is the probability mass that the unbounded model assigns to more exposures than the finite opportunities allow.

CANEX: heterogeneity and jointly observed duplication

calc_canex() implements Danaher’s (1991) canonical expansion model. It combines a Beta-Binomial distribution for each vehicle with a correlation term that adjusts the joint probability according to the duplication observed between each pair of vehicles.

data(canex_example)
canex_result <- do.call(calc_canex, canex_example)
canex_result
#> CANEX MODEL (Canonical Expansion)
#> =================================
#> Description: heterogeneous vehicles with observed duplication between vehicles
#> 
#> HEADLINE METRICS:
#> -----------------
#> Total reach: 53.96% (539628 people)
#> 
#> MODEL PARAMETERS:
#> -----------------
#> Probability of 0 exposures (%): 46.04
#> 
#> EXPOSURE DISTRIBUTION:
#> ----------------------
#> (Percentage of the population receiving exactly N exposures)
#> 0 exposures: 46.04% (460372 people)
#> 1 exposure: 19.64% (196407 people)
#> 2 exposures: 14.87% (148676 people)
#> 3 exposures: 14.65% (146530 people)
#> 4 exposures: 3.34% (33424 people)
#> 5 exposures: 1.46% (14591 people)
#> 
#> CUMULATIVE DISTRIBUTION:
#> -------------------------
#> (Percentage of the population receiving N or more exposures)
#> >= 0 exposures: 100.00% (1000000 people)
#> >= 1 exposure: 53.96% (539628 people)
#> >= 2 exposures: 34.32% (343221 people)
#> >= 3 exposures: 19.45% (194545 people)
#> >= 4 exposures: 4.80% (48015 people)
#> >= 5 exposures: 1.46% (14591 people)
#> 
#> SUMMARY STATISTICS:
#> --------------------
#> Average exposures per person reached: 2.11
canex_result$diagnostics
#> $negative_mass_truncated
#> [1] 0
#> 
#> $mass_before_renormalization
#> [1] 1
#> 
#> $correlation_min_eigenvalue
#> [1] 0.9059873
#> 
#> $mean_exposures_expected
#> [1] 1.14
#> 
#> $mean_exposures_result
#> [1] 1.14

The number of exposure combinations grows exponentially with the number of vehicles and insertions, so the model is intended for small and medium schedules. For very large schedules calc_canex() stops with an informative error instead of exhausting memory. Because the truncated expansion (pairwise terms only; see the section “Order of the expansion” of ?calc_canex) can produce negative probabilities, calc_canex() sets them to zero and renormalizes; the diagnostics report how much mass was truncated, and the mean number of exposures that the expansion preserves and the one it returns.

CSD and MSAD: two reach estimators, one sequential architecture

calc_csd() implements the Canonical Sequential Aggregation Distribution of Kim (2005): Beta-Binomial marginals, reach from the pairwise canonical expansion, and non-random convolution that preserves the margins, and with them the mean number of exposures. csd_example is an original three-vehicle scenario.

data(csd_example)
csd_result <- do.call(calc_csd, csd_example)
csd_result$reach
#> $probability
#> [1] 0.6809823
#> 
#> $percent
#> [1] 68.09823
#> 
#> $people
#> [1] 0.6809823
csd_result$distribution
#>    contacts  probability     percent       people cumulative_probability
#> 1         0 0.3190176501 31.90176501 0.3190176501           1.0000000000
#> 2         1 0.1660488500 16.60488500 0.1660488500           0.6809823499
#> 3         2 0.1690685122 16.90685122 0.1690685122           0.5149334999
#> 4         3 0.1935314665 19.35314665 0.1935314665           0.3458649876
#> 5         4 0.0780020953  7.80020953 0.0780020953           0.1523335211
#> 6         5 0.0470451034  4.70451034 0.0470451034           0.0743314258
#> 7         6 0.0164951066  1.64951066 0.0164951066           0.0272863224
#> 8         7 0.0080081268  0.80081268 0.0080081268           0.0107912158
#> 9         8 0.0020895005  0.20895005 0.0020895005           0.0027830890
#> 10        9 0.0006935885  0.06935885 0.0006935885           0.0006935885
csd_result$steps
#>   step added_vehicle target_reach zero_probability random_zero_probability
#> 1    1             2    0.6044618        0.3955382               0.3888701
#> 2    2             3    0.6809823        0.3190177               0.3084787
#>   expansion_adjustment duplicated_reach row_margin_error column_margin_error
#> 1           0.01714722        0.1238669     0.000000e+00        0.000000e+00
#> 2           0.03416422        0.1302102     1.665335e-16        1.387779e-17

The function reproduces the worked example of Kim (2005) to within a few hundredths of a percentage point in the maintainer’s local tests (the figures of that example are not distributed with the package). The small differences arise because the dissertation rounds parameters and intermediate matrices to four decimals, whereas the function keeps the precision of its inputs.

MSAD: Morgensztern’s sequential aggregation

calc_msad() implements the structure described by Leckenby and Rice (1986) and Kim (2005): a Beta-Binomial marginal for each vehicle, the Morgensztern reach for each sub-schedule, and a non-random convolution that preserves the marginals and sets the zero-exposure probability to 1 - reach.

data(msad_example)
msad_result <- do.call(calc_msad, msad_example)
msad_result$reach
#> $probability
#> [1] 0.6383175
#> 
#> $percent
#> [1] 63.83175
#> 
#> $people
#> [1] 0.6383175
msad_result$steps
#>   step added_vehicle target_reach duplicated_reach row_margin_error
#> 1    1             2    0.6383175       0.09001128     1.110223e-16
#>   column_margin_error
#> 1                   0
msad_result$diagnostics
#> $probability_sum
#> [1] 1
#> 
#> $minimum_probability
#> [1] 0.02012174
#> 
#> $gross_mean_contacts
#> [1] 1.41
#> 
#> $distribution_mean_contacts
#> [1] 1.41
#> 
#> $mean_error
#> [1] 4.440892e-16

By default the vehicles with the highest audience are aggregated first, as Kim indicates. An explicit permutation tests the sensitivity to the order. The function stops if the duplications or the Morgensztern reach violate probabilistic bounds.

Kim specifies MSAD but refers to Lee (1988) for the complete numerical example. calc_msad() reproduces the three-vehicle numerical example of Lee (1988, pp. 81-93) with a largest absolute difference of 0.001077973 in probability in the cell of three exposures, where Lee prints three decimals; that comparison is part of the maintainer’s local tests, because the example is not distributed with the package (see inst/DATA-PROVENANCE.md).

Validation against observed data

Observed data are always supplied by the analyst. For one schedule they must be a table with the columns contacts and observed, explicitly including zero exposures, and the scale of the values (counts, probabilities or percentages) must be declared.

# Invented counts, for illustration only: they are not real data
observed_counts <- data.frame(
  contacts = 0:9,
  observed = c(520, 150, 130, 80, 50, 40, 20, 5, 3, 2)
)

csd_evaluation <- evaluate_exposure_model(
  observed = observed_counts,
  predicted = csd_result,
  observed_scale = "count"
)

csd_evaluation$summary
#> $schedules
#> [1] 1
#> 
#> $kim_aer
#> [1] 0.4187132
#> 
#> $kim_ape
#> [1] 0.4425541
#> 
#> $mean_total_variation
#> [1] 0.2067042
#> 
#> $mean_cell_mae
#> [1] 0.04134083
#> 
#> $mean_reach_absolute_error
#> [1] 0.2009823
#> 
#> $mean_contact_bias
#> [1] 0.563
csd_evaluation$by_schedule
#>   schedule observed_reach predicted_reach reach_signed_error
#> 1        1           0.48       0.6809823          0.2009823
#>   reach_absolute_error kim_relative_reach_error kim_distribution_error
#> 1            0.2009823                0.4187132              0.4425541
#>   total_variation   cell_mae observed_mean_contacts predicted_mean_contacts
#> 1       0.2067042 0.04134083                  1.247                    1.81

If the prediction is also a table, its columns must be named contacts and predicted, and predicted_scale must be declared. To evaluate several schedules at once, add an identifier column to both tables and name it in schedule_col.

The evaluator requires identical supports. It does not interpret an open tail of the Negative-Binomial approximation as an exact frequency: observed and predicted tails must first be collapsed at the same threshold. Besides Kim’s average percentage error in reach (AER) and in the exposure distribution (APE), it returns absolute reach error, total variation, cell mean absolute error and mean-exposure bias.

MBD: the multivariate Beta-Binomial model

calc_mbd() implements Cheong’s (2007) model. Pairwise co-exposure probabilities are observed; those of three or more vehicles are imputed from a Beta-Binomial fitted to the subset’s own mean audience and duplication. They are combined into the joint zero/one grid by Waring’s (1792) inclusion-exclusion theorem. Vehicles are then peeled off one at a time, in reverse aggregation order, expanding each remaining pattern with a conditional Beta-Binomial allocation.

data(mbd_example)
mbd_result <- do.call(calc_mbd, mbd_example)
mbd_result$reach
#> $probability
#> [1] 0.5922363
#> 
#> $percent
#> [1] 59.22363
#> 
#> $people
#> [1] 0.5922363
mbd_result$distribution
#>   contacts probability   percent     people cumulative_probability
#> 1        0  0.40776371 40.776371 0.40776371             1.00000000
#> 2        1  0.18055145 18.055145 0.18055145             0.59223629
#> 3        2  0.18470533 18.470533 0.18470533             0.41168484
#> 4        3  0.17493820 17.493820 0.17493820             0.22697951
#> 5        4  0.03498325  3.498325 0.03498325             0.05204131
#> 6        5  0.01705806  1.705806 0.01705806             0.01705806

Against the three-vehicle example of Cheong (2007, p. 75), also published as Table 4 of Cheong, Leckenby and Eakin (2011), the function reproduces the published distribution (in the maintainer’s local tests, because that example is not distributed with the package) with a largest absolute difference of 0.0020348568 in probability (0.20348568 percentage points, in the cell of three exposures), larger than the rounding of one printed cell (0.0005). The regression test accepts differences below 0.0025, a threshold above the real maximum. The printed intermediate tables contain arithmetic inconsistencies, which are a plausible explanation for part of the differences; their full causal attribution has not been established, and this is an approximate replication of a printed example, not an independent validation.

Unlike CSD and MSAD, the model can produce negative probabilities when it combines four or more vehicles, because the co-exposure of three or more vehicles is imputed rather than observed. Cheong’s four-vehicle example is corrected by a first-order consistency check. The five-vehicle example is not: Cheong notes that a higher-order check would be needed, which was never implemented, and proposes iterative proportional fitting (Deming, 1943; the method is due to Deming and Stephan, 1940) as future work. calc_mbd() applies exactly the checks Cheong specifies (first order only) and then, as a package safety net, sets any cell that remains negative in the final distribution to zero and redistributes that mass proportionally. This is not Cheong’s MBD-ADJ, which applies the correction to the table that Cheong calls the UD (univariate distribution) before it is expanded and collapsed, so the published evidence for MBD-ADJ does not transfer to it. The diagnostics report it:

mbd_result$diagnostics$negative_mass_adjusted
#> [1] 0
mbd_result$diagnostics$cells_adjusted
#> [1] 0

With four or more vehicles calc_mbd() warns about this limitation. The aggregation order can also alter the result (Cheong did not investigate this systematically), and the model is limited to 12 vehicles because of the exponential cost of the exposure grid, the same practical limit Cheong applied in the dissertation’s own evaluation.

The safety net is not the only way in which the result can depart from the inputs. When a vehicle is expanded, the conditional probability that splits each row of the joint table between the two Beta-Binomial expansions of that vehicle is bounded to the interval that the marginal of the vehicle allows (Cheong, 2007, pp. 65-66). If the value implied by the observed co-exposures falls outside the interval, replacing it by the nearest end changes the marginals and the mean number of exposures of the result, even when no probability is negative. The following data are compatible (a joint distribution of the numbers of exposures of the two vehicles exists, and the duplication, 0.15, is above the random one, 0.06) and imply one exposure per person, but the distribution returned has a different mean, with no cell adjusted:

vehicles <- data.frame(insertions = c(2, 2), R1 = c(0.30, 0.20),
                       R2 = c(0.40, 0.30))
duplications <- matrix(c(0, 0.15, 0.15, 0), nrow = 2)
mbd_mean <- calc_mbd(vehicles, duplications, aggregation_order = "given")
unlist(mbd_mean$diagnostics[c("negative_mass_adjusted", "gross_mean_contacts",
                              "distribution_mean_contacts", "mean_error")])
#>     negative_mass_adjusted        gross_mean_contacts 
#>                 0.00000000                 1.00000000 
#> distribution_mean_contacts                 mean_error 
#>                 1.05378788                 0.05378788

Read mean_error every time, whether or not the safety net was engaged. The average frequency of the result is the mean of its own distribution divided by its reach, so it cannot be combined with the gross rating points of the inputs (100 here, against about 105 implied by the distribution): the identity between rating points, reach and average frequency requires the same distribution. CSD and MSAD conform each step to the margins that they combine and conserve the mean; CANEX, CBD and MBD can alter it.

CBD: the conditional Beta distribution

calc_cbd() implements the Conditional Beta Distribution of Kim (1994), reviewed in Kim (2005, pp. 59-65). It is a multivariate model of conditional expansion, not a sequential-aggregation model like CSD and MSAD. In Kim’s (2005) evaluation CBD is among the three most accurate models, together with CSD and CANEX; in Cheong’s (2007) it is second only to CANEX for the complete exposure distribution.

The model has three steps. Danaher’s (1991) canonical expansion truncated to its pairwise terms, the mechanism used by calc_canex() and calc_csd(), gives the joint distribution of exposure to one insertion of each vehicle (the “(0,1) grid”: every combination of exposed, 1, or not exposed, 0, to the first insertion of each vehicle). Conditionally on that state, the remaining insertions of each vehicle follow a Beta-Binomial distribution whose Beta parameters are updated by the first insertion, and the vehicles are independent given the grid (Danaher, 1992a). Finally, the zero cell of the grid is conformed to the zero probability of the schedule with all the insertions and its remaining levels are rescaled so that the cell keeps its probability. In its general formulation CBD lets any existing exposure model supply that zero probability (Kim, 1994, p. 124; Kim, 2005, p. 65, gives CANEX as an example); the version implemented here, like Kim’s worked example, takes it from the canonical expansion of all the insertions. The result does not depend on the aggregation order.

cbd_result <- do.call(calc_cbd, csd_example)
cbd_result$reach
#> $probability
#> [1] 0.6809823
#> 
#> $percent
#> [1] 68.09823
#> 
#> $people
#> [1] 0.6809823
cbd_result$distribution
#>    contacts  probability     percent       people cumulative_probability
#> 1         0 0.3190176501 31.90176501 0.3190176501           1.0000000000
#> 2         1 0.1707909669 17.07909669 0.1707909669           0.6809823499
#> 3         2 0.1684272430 16.84272430 0.1684272430           0.5101913830
#> 4         3 0.1898293130 18.98293130 0.1898293130           0.3417641399
#> 5         4 0.0758276599  7.58276599 0.0758276599           0.1519348270
#> 6         5 0.0474020446  4.74020446 0.0474020446           0.0761071671
#> 7         6 0.0168683896  1.68683896 0.0168683896           0.0287051225
#> 8         7 0.0085833116  0.85833116 0.0085833116           0.0118367329
#> 9         8 0.0023937018  0.23937018 0.0023937018           0.0032534213
#> 10        9 0.0008597195  0.08597195 0.0008597195           0.0008597195

The implementation reproduces the CBD columns of the two-vehicle comparisons of Kim (2005, Appendix B) and Hong (1998, Appendix E) within the rounding of those sources (0.01 percentage points), taking the duplication that each source implies (the sources do not print it: it was reconstructed from another model’s published column, so this checks consistency between published columns, not predictive accuracy), and the three-vehicle worked example of Kim (1994, pp. 125-139) except for its zero cell. In that example the zero of the canonical expansion is printed as 0.5066, whereas the formula applied to the printed inputs gives about 0.510; calc_cbd() evaluates the formula, so its zero and one-contact probabilities differ from the printed ones by 0.36 and 0.29 percentage points there. The comparisons with those published cases are made only in the maintainer’s local runs, because their input tables are not distributed; the public tests check that the (0,1) grid recovers each vehicle’s own R1 as its marginal and that, with zero correlation, the model reduces to the exact convolution of the vehicles’ own Beta-Binomial marginals.

Example datasets

Models with complex inputs come with datasets whose elements match the function’s arguments, so they can be used with do.call():

data(canex_example)
str(canex_example)
#> List of 3
#>  $ vehicles_data:'data.frame':   2 obs. of  3 variables:
#>   ..$ insertions: num [1:2] 3 2
#>   ..$ R1        : num [1:2] 0.3 0.12
#>   ..$ R2        : num [1:2] 0.4 0.18
#>  $ duplications : num [1:2, 1:2] NA 0.05 0.05 NA
#>  $ population   : num 1e+06
do.call(calc_canex, canex_example)$reach$percent
#> [1] 53.96284

Every dataset is an original illustrative scenario, meant for quick experimentation; none reproduces inputs published by a third party, whose redistribution basis is not established (the comparisons with published examples are local-only tests of the maintainer). ratings_example serves calc_sainsbury() and calc_binomial(), and duplication_example serves calc_agostini_duplication() and calc_hofmans_duplication(), because each pair shares the same input shape. No dataset is named binomial, so that stats::binomial is never masked.

The package has two Hofmans models, derived from the same basis (Agostini’s formula) but for different plan shapes and with different data, so both function names carry their domain: calc_hofmans_accumulation() (one vehicle, several insertions, with hofmans_accumulation_example) and calc_hofmans_duplication() (several vehicles, one insertion each, with duplication_example).

Domain of validity of the multivariate models

In plain words: these models take the overlaps between vehicles from your data and build the whole distribution of exposures from them. For that the overlaps must be possible, and plausible together. The functions react in two ways, and each reaction asks for a different action:

The function	It means	What to do
Stops with an error naming a pair or a triple of vehicles	The overlaps are impossible: a pair overlaps more than its smaller audience, or less than the two audiences force; or the overlaps of three vehicles cannot hold together	Review the data of the vehicles named: units (proportions, not percentages), and whether every overlap is measured on the same population
Stops with an error about `R2`	`R2` is outside the range allowed for its `R1` (between `R1` and `2 * R1 - R1^2`, the range of the Beta-Binomial family: a larger `R2`, up to `min(1, 2 * R1)`, can be possible data that this family cannot represent)	Check that `R2` is the cumulative reach after two insertions, as a proportion; if it is, review how the exposures were measured and whether the model suits the vehicle, without lowering `R2` just to make it run
Warns about a low duplication	A pair overlaps clearly less than chance would (less than the product of their audiences). That can be valid data, but the studies that evaluated the models mostly used overlaps at least as large as chance, so the approximations can differ a lot from one model to another	Check how the duplication of the pair named was defined (a zero is often a measurement problem, but not always), read the diagnostics, and compare the result with another model before using it; do not replace valid observations automatically
Warns about truncated probabilities (`calc_canex()`)	The approximation produced negative probabilities that were set to zero and the distribution was normalized again; the mean number of exposures may have changed (the warning gives the expected and the returned mean, and `diagnostics` the truncated mass)	Compare `mean_exposures_expected` with `mean_exposures_result` in `diagnostics`; treat the distribution with caution

The technical detail follows.

calc_canex(), calc_csd(), calc_msad(), calc_cbd() and calc_mbd() take the observed one-insertion duplications as given. They stop when a duplication violates the Fréchet bounds implied by the audiences, or when some triple of vehicles has no non-negative joint exposure probability compatible with its pairwise duplications (a necessary condition that is not sufficient for four or more vehicles). They warn when the canonical correlation of a pair is below -0.1. Hong (1998, pp. 238-239) formulates, for the models and data of his study, the hypothesis of a duplication not below the random one and treats null duplications as measurement problems; it is not a logical condition (the limits of a duplication are the Fréchet bounds, and a lower duplication can be valid data), but the approximations can diverge considerably in that range. The threshold is a decision of this package, and the published reference plans stay above it. calc_canex() also warns when truncating the negative probabilities of its expansion is material, and reports in diagnostics the mean number of exposures that the expansion preserves and the one it returns. ?calc_canex documents the domain in full.

Scope

The package adds a model only when its primary specification and a numeric reference that allows an independent regression test are available.

Model	Numeric evidence	Status
CSD	Complete example, Kim (2005), Tables 4.2.2.1-4.2.2.10; 40 plans of Appendix B	Implemented and validated
MSAD	Algorithm given by Kim; detailed example in Lee (1988, pp. 81-93); 40 plans of Kim (2005, Appendix B)	Implemented and validated against Lee’s example and Kim’s plans (plan 19 differs by 0.099 percentage points)
CANEX	Danaher’s (1991) formula, the CSD example that uses its reach, and the CANX column of Kim (2005, Appendix B)	Implemented; reports approximation diagnostics
Beta-Binomial and Binomial	Closed-form formulas	Implemented and tested at their limits
MBD	Complete specification and three-vehicle example, Cheong (2007) and Cheong, Leckenby and Eakin (2011); four- and five-vehicle examples show the model’s limits	Implemented, with a package safety net (not Cheong’s MBD-ADJ) and warnings; approximately reproduces the published example (largest absolute difference 0.0020348568 in probability)
CBD	Specification and three-vehicle example in Kim (1994, pp. 115-139); two-vehicle outputs in Kim (2005, Appendix B) and Hong (1998, Appendix E)	Implemented; reproduces those outputs (a consistency check between published columns); zero cell of Kim’s example differs by about 0.4 percentage points
HBBD (Hofmans Beta-Binomial)	Kim (2005, pp. 52-53) gives the Hofmans reach formula followed by iterative mean-zero Beta-Binomial fitting	Available as a two-step composition, `calc_hofmans_duplication()` followed by `fit_bbd_to_reach()` (worked example in this vignette and in `?fit_bbd_to_reach`); there is no single function, and it is not validated against a published numeric example in this release
DMD (Dirichlet multinomial)	Kim (2005, pp. 54-56) specifies the distribution and refers to detailed steps in Boyd (1985) and Kim (1994, pp. 89-98)	Not implemented or validated in this release
HBD	Approximation of MBD described by Cheong (2007), without a verified numeric example	Not implemented
Multivariate Negative-Binomial	Requires Danaher’s (2007) formulation for online processes	Not implemented

References

Agostini, J. M. (1961). How to estimate unduplicated audiences. Journal of Advertising Research, 1(3), 11-14. https://doi.org/10.1080/00218499.1961.12519620
Aldás Manzano, J. (1998). Modelos de determinación de la cobertura y la distribución de contactos en la planificación de medios publicitarios impresos (Models for determining reach and exposure distribution in print media planning). Doctoral dissertation, Universidad de Valencia, Spain.
Boyd, M. M. (1985). The Dirichlet multivariate multinomial distribution as a magazine exposure model. Unpublished doctoral dissertation, University of Illinois at Urbana-Champaign.
Chandon, J.-L. (1976). A comparative study of media exposure models. Unpublished doctoral dissertation, Northwestern University, Evanston, IL.
Cheong, Y. (2007). Multivariate Beta Binomial Distribution Model as a Web Media Exposure Model. Doctoral dissertation, The University of Texas at Austin. Handle 2152/3215 (University of Texas at Austin repository).
Cheong, Y., Leckenby, J. D., & Eakin, T. (2011). Evaluating the multivariate beta binomial distribution for estimating magazine and Internet exposure frequency distributions. Journal of Advertising, 40(1), 7-23. https://doi.org/10.2753/JOA0091-3367400101
Danaher, P. J. (1991). A canonical expansion model for multivariate media exposure distributions: A generalization of the “duplication of viewing law”. Journal of Marketing Research, 28(3), 361-367. https://doi.org/10.1177/002224379102800311
Danaher, P. J. (1992a). A Markov-chain model for multivariate magazine-exposure distributions. Journal of Business & Economic Statistics, 10(4), 401-407. https://doi.org/10.1080/07350015.1992.10509915
Danaher, P. J. (2007). Modeling page views across multiple websites with an application to Internet reach and frequency prediction. Marketing Science, 26(3), 422-437. https://doi.org/10.1287/mksc.1060.0226
Deming, W. E. (1943). Statistical adjustment of data. New York: Wiley (reprinted by Dover in 1964; Cheong, 2007, cites it with the Dover imprint).
Deming, W. E., & Stephan, F. F. (1940). On a least squares adjustment of a sampled frequency table when the expected marginal totals are known. Annals of Mathematical Statistics, 11(4), 427-444.
Ehrenberg, A. S. C. (1959). The pattern of consumer purchases. Applied Statistics, 8(1), 26-41. https://doi.org/10.2307/2985810
Hofmans, P. (1966). Measuring the cumulative net coverage of any combination of media. Journal of Marketing Research, 3(3), 269-278. https://doi.org/10.1177/002224376600300307
Hong, J. (1998). Advertising media models for Internet reach/frequency estimation. Unpublished doctoral dissertation, The University of Texas at Austin.
Kim, H. (1994). A conditional beta distribution model for advertising media reach/frequency estimation. Unpublished doctoral dissertation, The University of Texas at Austin.
Kim, H. G. (2005). A Canonical Sequential Aggregation Media Model. Doctoral dissertation, The University of Texas at Austin. Handle 2152/1590 (University of Texas at Austin repository).
Leckenby, J. D., & Rice, M. D. (1986). The declining reach phenomenon in exposure distribution models. Journal of Advertising, 15(3), 13-20. https://doi.org/10.1080/00913367.1986.10673014
Lee, A. M., & Burkart, A. J. (1960). Some optimization problems in advertising media planning. Operational Research Quarterly, 11(3), 113-122.
Lee, H.-K. (1988). Sequential aggregation advertising media models. Unpublished doctoral dissertation, The University of Texas at Austin.
Metheringham, R. A. (1964). Measuring the net cumulative coverage of a print campaign. Journal of Advertising Research, 4(4), 23-28. https://doi.org/10.1080/00218499.1964.12519751
Waring, E. (1792). On the principles of translating algebraic quantities into probable relations and annuities. Cambridge. (The work cited for the theorem by Cheong, 2007, and Aldás Manzano, 1998; not Waring’s Meditationes Algebraicae, 1770.)

Morgensztern (1970), whose reach formula underlies MSAD, is known to this package through Aldás Manzano (1998) and Kim (2005). Lee and Burkart (1960) and Lee (1988) are cited as they appear in Kim (2005). The bibliographic details of Kim (1994), Lee (1988) and Hong (1998) are those of the title pages of the dissertations; those of Chandon (1976) and Boyd (1985) are given in the reference list of Kim (2005).