arXiv:2307.06024v2 [stat.

CO] 13 Jul 2023

balance - a Python package for balancing biased

data samples
Tal Sarig∗ , Tal Galili∗ , and Roee Eilat
Meta

June 2023

Abstract

Surveys are an important research tool, providing unique measure-

ments on subjective experiences such as sentiment and opinions that can-
not be measured by other means. However, because survey data is col-
lected from a self-selected group of participants, directly inferring insights
from it to a population of interest, or training ML models on such data,
can lead to erroneous estimates or under-performing models. In this pa-
per we present balance, an open-source Python package by Meta, offering
a simple workflow for analyzing and adjusting biased data samples with
respect to a population of interest.
The balance workflow includes three steps: understanding the ini-
tial bias in the data relative to a target we would like to infer, adjusting
the data to correct for the bias by producing weights for each unit in
the sample based on propensity scores, and evaluating the final biases
and the variance inflation after applying the fitted weights. The package
provides a simple API that can be used by researchers and data scientists
from a wide range of fields on a variety of data. The paper provides the
relevant context, methodological background, and presents the package’s
API.
∗ The first two authors contributed equally to this work

1
Contents
1 Introduction 3
2 Related Work 4

3 Methodological Background 6
3.1 The Total Survey Error framework . . . . . . . . . . . . . . . . . 6
3.2 Definitions and notations . . . . . . . . . . . . . . . . . . . . . . 8
3.3 Estimation of the survey weights . . . . . . . . . . . . . . . . . . 9
3.3.1 Post-stratification . . . . . . . . . . . . . . . . . . . . . . 10
3.3.2 Raking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
3.3.3 Inverse Propensity score Weighting (IPW) . . . . . . . . . 12
3.3.4 Covariate Balancing Propensity Score (CBPS) . . . . . . 14
3.4 Evaluation of survey weights . . . . . . . . . . . . . . . . . . . . 15
3.4.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.4.2 Visualizing Distributions . . . . . . . . . . . . . . . . . . . 15
3.4.3 Diagnostics for the covariates using ASMD . . . . . . . . 17
3.4.4 Diagnostics for the weights . . . . . . . . . . . . . . . . . 18
3.4.5 Diagnostics for the outcome . . . . . . . . . . . . . . . . . 20
4 The balance workflow 21
4.1 The workflow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
4.2 An end-to-end example . . . . . . . . . . . . . . . . . . . . . . . 22
4.2.1 Understanding the initial bias . . . . . . . . . . . . . . . . 22
4.2.2 Fitting survey weights . . . . . . . . . . . . . . . . . . . . 26
4.2.3 Evaluating the Results . . . . . . . . . . . . . . . . . . . . 26
4.3 How does balance implement the adjustment? . . . . . . . . . . 30
5 Future directions 31

References 33
Appendices 38
A Acknowledgments 38

B Limitations of the ASMD 38

C Kish’s design effect 39
D Estimating the variance of the weighted mean 42

2
1 Introduction

Surveys play an important role in the study of social phenomena across research
fields and industries. From their traditional usage by statistics bureaus in pro-
ducing population estimates, through the long history of public opinion surveys
in political science, to more recent applications like studying user experience
in online services and even playing part in epidemiological studies [1]. The
widespread use of surveys, and their unique role in providing measurements on
subjective indicators such as sentiment and opinions, makes the field abundant
with methodological research.

A central challenge in designing and analyzing survey data stems from bias
due to sampling limitations and non-response. Since the data is collected from
a self-selected group of participants, directly inferring insights or training ML
models on such data can result in erroneous estimates or under-performing mod-
els. An insightful theoretical framework for the sources of error present in sur-
vey data is given in the ”Total Survey Error” framework [2]. While the sources
might be different, similar manifestations of bias are often present in obser-
vational studies when comparing treatment groups, and in any data produced
through self-selection processes.

The field of survey statistics offers methods for mitigating bias in samples,
at least partially, by relying on auxiliary information (i.e., “covariates” or “fea-
tures”). When such information is available for all items in the sample as well as
for the population from which it was sampled, it can be used to create weights.
Under some assumptions on the relation between the auxiliary information, the
response mechanism, and the survey responses, applying the weights to the data
will produce less biased estimates or models. Different approaches were pro-
posed for the task, from simple post-stratification [3] to methods more suitable
for high dimensional covariates space such as raking [4, 5, 6], inverse propensity
weighting [7, 8, 9], covariate balancing methods [10], outcome regression based
approaches [11], and others. Weighting methods have been shown to be effective
in reducing bias of survey estimates [12].

Following methodological advancements in survey statistics, statistical soft-

ware packages were developed to allow researchers and practitioners to apply
these methodologies to survey data and observational data. Most software pack-
ages for this aim have R implementations, and other implementations in envi-
ronments such as SPSS, stata, SAS exist as well. In recent years a rich ecosystem
of data science software has been developed for Python, and its usage has be-
come prevalent among researchers and data scientists. This shift created a need
for a reliable Python package for working with survey data, and more generally
with biased data sets. Here we introduce balance - a Python package for bal-
ancing biased data samples. balance offers a simple easy-to-use framework for
weighting data and evaluating its biases. The package is designed to provide
best practices for weights fitting and offers several modeling approaches. The

3
methodology in balance can support ongoing automated survey data process-
ing, as well as ad-hoc analyses of survey data.

The main workflow API of balance includes three steps: (1) understanding
the initial bias in the data relative to a target population as observed by the
differences in covariate distribution (2) adjusting the data to correct for the bias
by producing weights for each unit in the sample based on propensity scores, and
(3) evaluating the final biases and the variance inflation after applying the fitted
weights. The adjustment step provides a few alternatives for the researcher to
choose from: Inverse propensity weighting using logistic regression model based
on LASSO (Least Absolute Shrinkage and Selection Operator [13]), Covariate
Balancing Propensity Scores [10], Raking, and post-stratification. The focus
is on providing a simple to use API, based on Pandas’s DataFrame structure,
which can be used by researchers from a wide spectrum of fields.

In this paper we describe the balance workflow in more detail and provide
guidance on how to implement it using the package. We include details on
methods, assumptions, and model choices made in the package. The method-
ological background part of the paper is an accessible review of the theoretical
frameworks, methods, and practices often used in survey statistics. We invite
readers new to the field to use it as a short and effective introduction.

The rest of this paper is structured as follows. We discuss related work in Sec-
tion 2, focusing on software packages available in the R and Python ecosystems
for survey data analysis and related use cases. In Section 3 we provide details
on the statistical background that guided the implementation of the package,
including theoretical frameworks, estimation methods, and diagnostic tools. In
Section 4 we present the balance workflow and provide an end-to-end walk
through using code snippets that are applied to simulated data. We conclude
with a discussion on future directions for the package in Section 5.

2 Related Work

The open-source ecosystem offers a variety of packages for weighting biased

data. This section gives a brief survey of prominent tools in this space and
describes some of their capabilities. We find the R ecosystem to be the most
developed in terms of packages for survey analysis. The Python ecosystem
has some packages for survey statistics. It also has several, well developed,
packages for casual inference, which employs similar models (e.g.: propensity
scores models, outcome models, etc.). While various R packages exist with
similar capabilities to what is available in balance, no Python package (that we
are aware of) gives a comprehensive end-to-end coherent solution for researchers.

The R ecosystem is exceptionally rich and diverse when it comes to survey

4
statistics. The most comprehensive review can be seen in the CRAN task view of
”Official Statistics & Survey Statistics” [14]. To date, it includes over 130 pack-
ages - ranging from the classical survey package [15] to more niche packages.
Similarly, the CRAN task view of ”Causal Inference” [16] also includes over 130
packages that offer related methods. A short review on the current state of R
packages can be found in the PSweight R package [17], which compares 9 R
packages that implement propensity score weighting with discrete treatments.
For survey weights diagnostics, the cobalt package [18] offers many options, in-
cluding balance tables and plots for covariates of multiple groups. This package
includes various capabilities that could inspire future development of balance.

For Python, the ipfn package [19] specializes in implementing a fast iterative
proportional fitting (raking). This package is utilized in balance and used
as the back-end for the raking implementation we rely on. The quantipy3
package [20] is designed to support data processing, analysis and reporting for
survey data using pandas and numpy. It supports native handling of special data
types like multiple choice variables, statistical analysis using case or observation
weights, DataFrame metadata and different data exports. quantipy3 seems to
be the most similar to what balance tries to achieve but lacks many of the
capabilities balance has in all stages of the workflow. The samplics package
offers a comprehensive solutions for dealing with complex sampling designs [21],
with various overlapping and non-overlapping capabilities between this package
and balance. The package offers tooling for random selection techniques used
to draw a sample from a population, and sample size calculations. It also
provides methods for weight adjustments using post stratification or calibration.
Additional capabilities in samplics include functions for estimation of statistics
and their variance (beyond just the Taylor linearization estimation in balance).
These include bootstrap, balanced repeated replication and Jackknife. Other
packages we found seem to be only lightly maintained, and do not provide
additional capabilities of relevance to our use-case. These include PySurvey
[22], Surveyweights [23], pscore match [24], pymatch [25], and causal nets [26].

Stepping aside from survey statistics, several Python packages offer tools for
casual inference that can be repurposed for adjusting biased samples. The
DoWhy package [27], developed by Microsoft, is a well maintained package with
the focus on causal inference. It models a given problem as a causal graph to
help explore assumptions clearly, estimate causal effects, and test assumptions’
validity for robustness. It offers a variety of methods for estimation includ-
ing Propensity-based Stratification, Propensity Score Matching, and Inverse
Propensity Weighting (similar to balance). It also offers outcome based mod-
els (currently not implemented in balance) using Linear Regression or General-
ized Linear Models, and supports other methods such as Instrumental Variable
methods. The package emphasizes graphical interpretation of causal inference.
It also gives various refutation methods (dummy outcome, simulated outcome,
etc.) and basic visualizations (e.g.: barplots of treatment and control). The
Empirical Calibration package [28], developed by Google, provides a method

5
to compute empirical calibration weights using convex optimization. This ap-
proach balances out the marginal distribution of covariates directly while reduc-
ing the inflation of variance. This is similar to performing raking while trying
to keep the weights to be as equal as possible. It offers a bias correction solu-
tion that resembles the raking and CBPS methods that are implemented in the
balance package. The causalml package [29] provides a set of modeling and
causal inference methods for analyzing observational data using machine learn-
ing algorithms. It provides tool to estimate the Conditional Average Treatment
Effect (CATE) and the Individual Treatment Effect (ITE). This package offers
a wide variety of ML algorithms, including tree-based algorithms, meta-learner
algorithms, instrumental variables algorithms, and neural-network-based algo-
rithms. While these packages are comprehensive, there is still an overhead and
complexity for using them for balancing data for the workflow balance is opti-
mized in handling with the focus on surveys data.

3 Methodological Background

Before diving into the workflow and the implementation details of balance, we
introduce a brief description of the methodological background concerning the
representation error problem in surveys, weights estimation and tools to evaluate
survey weights.

3.1 The Total Survey Error framework

Figure 1: A flow diagram of ”Total Survey Error”, illustrating the different

components of surveys’ representation error.

The ”Total Survey” Error framework [2] provides a theoretical framework to

describe statistical properties of surveys. It is used as a conceptual tool for
researchers when designing and analyzing surveys to minimize estimation errors
and biases. While the research goal is to estimate a population parameter, such
as average or ratio, surveys only provides a glimpse on this parameter through
the survey responses and are subject to a range of sources for statistical errors,
as described by the ”Total Survey Error” concept.

6
 The ”Total Survey Error” has two main components: representation error and
measurement error [30]. Since neither can be overcome by increasing the sample
size, researchers should be aware of these as early as the survey design stage.
Measurement error deals with potential biases introduces to the estimation due
to the instrument of measurement. It includes questions about the validity of
the responses, the phrasing of the questions and how it is affecting what we
are trying to estimate, and similar questions related to whether we measure the
exact quantity we aim for. balance is focused on addressing and correcting the
representation errors in this framework, and hence for the rest of the section we
will focus on the representation part of the framework.

The representation error deals with how to infer from a subset of people to
the whole population on which we would like to learn, referred to as the target
population. The magnitude of the error depends on the group of respondents
to the survey and depends on how similar or different this group is from the
target population. Formally, Figure 1 shows the different sources of represen-
tation error and illustrates a breakdown of the difference between the group of
respondents and the target population.

The first error we consider is the coverage error. Its driver is the misalignment
between who can be sampled for the survey (the ”sampling frame”) and the
target population. In today’s world, where many, if not most, surveys are
conducted through the internet, a common sampling frame is people with access
to the internet. Since this sampling frame may not be representative of the whole
population, a caution should be taken when conducting survey over the web.

The most canonical example for a sampling frame that is not fully covering the
target population is the ”Literary Digest” 1936 poll [31]. During this year the
Literary Digest magazine ran a poll to predict the result of the U.S. election.
Franklin Delano Roosevelt was the Democrats candidate and the Republican
candidate was Governor Alfred Landon of Kansas. The magazine predicted a
decisive victory for Landon with a poll that was based on roughly 2.4 million
voters, but, as history tells us, Roosevelt won 62% of the votes. Even though
the poll sampled 10 million people, the sampling frame was skewed. The sample
included magazine readers, and people from phone lists and club memberships
lists. However, since people of lower socioeconomic status were disproportion-
ately not part of the magazine’s audience those days, the poll missed a significant
and unique portion of the U.S. voters population Due to setting a sampling frame
that ignores the target population definition, the magazine’s coverage error was
large and led to mis-prediction of the elections’ results.

Once the sampling frame is set, the researcher samples a certain amount
of people from the frame to ask to reply on the survey. This is the sample
population, or the group of people that have a ”real” opportunity to reply on
the survey. When doing so, the researcher reveal another gap where error can
occur due to sampling, which is the sampling error. This might be small if
we are able to sample either completely at random from the sampling frame

7
or by designing the sampling with the right sampling probabilities, but can be
significant given wrong assumptions on the structure of the sampling frame or a
complex mechanism of sampling. This error can be reduced as we increase the
sample size (and will be 0 if the sample population is the same is the sampling
frame), and is the one captured by the margin of error, often reported with a
survey results.

Once researcher has sent out the invitation for filling the survey to the sample
population, most often only a portion will choose to take part in the survey.
These are the respondents, or the observed sample. This self selection behavior
causes another error component which is the non-response error. This error can
be substantial depending on how the survey is conducted, the survey questions,
and other issues related to the instrument. The percent of non-response can
give us some intuition of how large this error is but the actual size of the bias
depends only on the properties of the people who chose to respond. In the case
of the Literary Digest poll the response rate was only 24%. In fact, research
suggests that the primary source of the error in the poll originated in the non-
response bias. Specifically, people who strongly disliked Roosevelt were more
willing to take the time to mail back their response [32, 33].

balance aims to correct for all types of representation errors at once (see
Figure 1). Using additional assumptions, as described in the next section, we
are able to make the group of respondents, i.e. the observed sample, similar
in properties to the target population and hence overcome some parts of the
representation error. However, it is important to note that there are cases where
it is impossible to fully correct the representation error. Such cases occur when
the assumptions on the missingness are not satisfied. The simplest example
of such case is when there is a substantial coverage error for which we cannot
overcome using auxiliary data. For example, if we want to learn about North
America’s population but survey people only from The United States. Even
given lots of auxiliary information we will likely not be able to adjust the sample
such that it correctly represents Canada’s population as well.

3.2 Definitions and notations

With the Total Survey Error framework in mind, we will now set definitions
to be used throughout the paper. Let S denote a sample of respondents to a
survey consisting of n respondents (sometimes referred to as the sample), and
let T represent a target population with N units.

Furthermore, we assume we have some auxiliary data on all units in sample

and target population, represented by a covariates (or features) vector attached,
Xi . Note that we assume that the same covariates are available for the sam-
ple and the target, otherwise we ignore the non-overlapping covariates. This
framework is applicable when we have the auxiliary data at the unit level for

8
the entire population or for a representative random sample from it. For ex-
ample, when we sample from a list of customers for which we have auxiliary
information available, or in cases when a reference survey is available (a survey
with better sampling properties to be used for correcting biases [34]). Another
common use-case is when census data of the population is available.

We define R to be an indicator for inclusion in the respondents group1 , i.e.

Ri = 1 if i ∈ S and Ri = 0 if i ∈ / S. Furthermore, we define y to be the
answer to one item of the survey. The answer can be numeric or discrete, and
is observed only for S. In our setup, we think about y as a constant and the
random variable, later considered for statistical properties, is R.

Our objective is to estimate a certain parameter of the population. The

simplest example is estimating the mean of one item of P the survey, i.e. the
mean of y. InPthis case, a natural estimate to ȳ = N1 i∈T yi is the sample
mean ȳS = n1 i∈S yi . However, due to the non-random sampling of S from the
population T , the proposed estimate will be biased, such that E [ȳ − ȳS ] ̸= 0.

3.3 Estimation of the survey weights

Weights are a common way to overcome survey error, and are essential when
estimating a population parameter [12]. This is generally done by incorporating
the weights, wi where i ∈ S, into the parameter estimation procedure2 . An
example is the case of estimating a parameter of the population using a weighted
mean of the sample: Pn
wi yi
ȳw = Pi=1
n (1)
i=1 wi
Further details about using the weights for estimations are described in section
3.4.5. One of the advantages of using weights, unlike alternative methods, is
the flexibility it gives in estimation. The weights depend only on the group
of respondents and not on the outcome itself, hence give the researcher the
flexibility to use the same set of weights for multiple outcomes, combine multiple
outcomes into one parameter, or consider the outcome in only in specific cuts
of the population.

We typically employ weights to adjust from the observed sample to match

the target population, so to overcome the representation bias of the sample. In
scenarios where the sampling procedure is set by design and therefore known, we
define the inverse of the sampling (or selection) probabilities as design weights 3 .
1 For estimation in balance, we think of the target population as a reference group, and

hence units of the target population are distinct from the units of the observed sample.
2 In the balance package, we choose to scale the weights to sum to the population size.

This way, each sample unit weight represents the number of corresponding units from the
population.
3 We may assume each unit i in S and/or T has a design weight, d . These are the
i

9
However, to overcome the full gap between the respondents and the target
population we need to estimate the weights according to the actual realization of
the observed sample. When estimated against the complete target population,
the weights can help address the non-response error, the by-design sampling
error, the ”unknown” sampling error and the coverage error.

A few assumptions are required to utilize the estimated weights for valid
estimations of the parameter of interest and mitigate the representation bias.
Under these assumptions, the estimation of the weights relies on the auxiliary
data Xi .

The first assumption is the Missing At Random assumption (MAR) [35]. The
MAR assumption states that the response mechanism is independent of the
survey responses conditional on the auxiliary data. In other words, Y ⊥ ⊥ R | X,
which means that given the covariates the likelihood of a person to respond to
the survey doesn’t depend on their answer. This assumption is also known as
the ignorability assumption or conditional unconfoundedness in causal inference
literature [36]. It is worth noting that recent research (such as [37]) proposes
alternative approaches to address missing values created by design in surveys
analysis.

The second assumption is positivity 0 < P (Ri = 1|Xi ) < 1 for all units in S
and T . 0 < P (Ri = 1|Xi ) means that given the auxiliary data, every unit of the
target population has non-zero probability to be include in the sample. In other
words, in a counterfactual world, any unit i ∈ T could have participated in the
survey given their covariates. Conversely, we also assume P (Ri = 1|Xi ) < 1,
implying that every unit in the observed sample is also in the target population.
The combination of the MAR assumption and the positivity assumption is often
known as the ”strong ignorability” assumption [36].

Given the assumptions we are now left with the question of how to esti-
mate the weights. The balance package currently supports 4 different methods
for estimating the weights: post-stratification, raking, Inverse Propensity score
Weights (IPW or IPSW) and Covariate Balancing Propensity Score (CBPS).
Next, we provide more background about the estimation process in each method
and describe advantages and limitations of each.

3.3.1 Post-stratification

Post-stratification [3] is one of the most common weighting approaches in sur-

vey statistics. It originates from a stratified sample or probability sampling,
where the population is divided into sub-populations (strata) and a sample is
sampling weights that are based on the known sampling procedure from the sampling frame
to the sample population. These are often set to 1 for the respondents when the sampling
probabilities are unknown.

10
independently drawn from each. However, in post-stratification, the stratifica-
tion is done after the sample has been selected. This is done to overcome errors
originating in mechanisms outside the sampling design, such as non-response.

The idea behind post-stratification is straightforward. For each cell (strata)

in the population, calculate the percentage it represents of the total population.
Then fit weights so that they adjust each stratum in the sample so to have the
same proportions as each strata as in the population.

Let H be the group items which represent some stratum in the population,
and PH represent the proportion of this stratum in the target population T ,
i.e. ph = |H|
N . Let nH be the number of respondents from stratum H in the
observed sample. We also define the ”inflation factor” as I = N/n, i.e. the factor
indicating by how much we need to multiply the total observed sample size to
get to the total target population size. Consequently, the post-stratification
weight for each unit i from stratum H in the observed sample is:

n
wi = PH ∗I ∀i ∈ H (2)
nH

Note that the multiplication by I is a result of the arbitrary choice to scale

the weights to the population size, and could be omitted.

The goal of post-stratification is to have the sample exactly match the joint-
distribution of the auxiliary data of the target population. Hence it requires the
researcher to know the joint distribution of the covariates to weight on. This
level of resolution for the target population may not always be available. When
only marginal distributions of the covariates are available then raking might
serve as an alternative method to estimate the weights. Raking is described in
sub-section 3.3.2.

Another limitation of post-stratification is on the number of covariates that

can be used to correct the biases due to the limitations of having enough respon-
dents in each of the cells. Having a cell with very few respondents could easily
lead to a handful of respondents that receive very large weights - which leads
to inflated variance of estimation based on such weights. Furthermore, when
continuous variables are required for weighting, the researcher must decide on
the thresholds for bucketing the variables into. A more general approach is the
inverse propensity score weighting described in sub-section 3.3.3.

3.3.2 Raking

Raking [4, 5, 6], also known as Iterative Proportional Fitting procedure (IPF),
is a method that fits the sample data to a target population using only the
marginal distributions of the population’s covariates. Typically, we have access

11
to these marginal distributions but often not to their joint distribution. Since
raking weights do not represent the joint distribution, this can be thought of
as a type of regularized model. This approach helps to avoid over-fitting small
cells as in post-stratification and instead focuses only on the marginals [38].

Raking essentially applies post-stratification sequentially over all covariates

using only the marginal distributions. This is done repeatedly until a conver-
gence is achieved. If exist, the design weights of the sample are used as the
starting point of the algorithm. For example, we may have the marginal dis-
tribution of gender, age, and education. Raking would first adjust weights to
match the gender distribution and then take these weights as input to adjust
for age, and then for education. It would then adjust again to gender and then
again to age, and so forth until it converges. This process will repeat until one
of three stopping criteria are met: (1) we reached a pre-defined number of it-
erations,(2) the maximum difference in proportions between sample and target
marginal distribution on any covariate is smaller then a preset convergence rate,
or (3) the weights have converged and the change from one iteration to another
is smaller then a preset rate tolerance parameter.

The resulting weights will be close to the marginal distribution of the popu-
lation covariates. However, one cannot assume that the weighted sample joint
distribution is the same as the joint distribution in the target population. Hence,
if one wants to infer only for a sub-group of the population (such as young-adults
only), it is less recommended to use raking weights, and, if possible, one should
prefer a method that take into account the joint distribution of the covariates
if such data exists.

Similar to post-stratification, raking is limited by the number of covariates

that can be included in the model, due to the need to have enough respon-
dents in each margin cell. In addition, raking may be sensitive to the order in
which covariates are adjusted for, which may lead to under-correction of some
covariates.

3.3.3 Inverse Propensity score Weighting (IPW)

A natural expansion to post-stratification and raking is the inverse propensity

score weighting that can be viewed as a continuous extension of post stratifica-
tion.

The Propensity Score is defined as the conditional probability to be part of

the observed sample given the covariates:

p(X) := P r(R = 1|X) (3)

It was first suggested by Rosenbaum and Rubin [35, 36] as a method to perform
matching for causal effects estimations in observational studies, and was later

12
adopted to weighting survey data [7, 8, 9]. Rosenbaum and Rubin [36] showed
that the assumptions of ”strong ignorability” (uncounfoundness: Y ⊥ ⊥ R | X
and positivity: 0 < P (Ri = 1|Xi ) < 1) implies that Y ⊥ ⊥ R | p(X), and
that p(X) is the coarsest balancing score (a score B(X) that satisfies Y ⊥
⊥R|
B(X)). This means that the propensity score is an inexpensive way, in terms
of dimension, to estimate the selection probabilities. Hence, in the spirit of
”Horvitz–Thompson estimator” [39] of using the inverse selection probabilities
as weights, the inverse of the propensity score was suggested as a weighting
procedure to adjust for non-response bias [8].

The estimation of the propensity scores can be done in any standard tools for
classification, such as logistic regression, decision trees and random forests (such
as in [40]) or others. The choice of the model depends on the researcher assump-
tions regarding the parametric model of the non-response and the number and
types of features used. In balance, we chose to implement the estimation of
the propensity scores through a (regularized) logistic regression. The logistic
regression model assumes a linear relation between the covariates and the log
odds, of the form:
pi
log( ) = βT X (4)
1 − pi

Once the propensity scores are estimated, the weight of unit i is calculated by
wi = 1− pˆi
pˆi . This is because we define the target population as a reference group
and we don’t assume the target doesn’t include the observed sample (i.e. we
don’t exclude units from the target based on their appearance in the sample).
In this case, borrowing concepts from the causal inference literature [41], the
estimation we care about is only the estimation of the average treatment effect
for the ”control” (”untreated”) group (the target population), and hence we use
1−pˆi
pˆi as the weights.

One challenge when including many covariates is that the estimation of the
propensity scores (and hence, the weights), can have a a high variance, which
may lead to unnecessary inflation of the survey estimates4 [42]. In balance
we try to mitigate this by applying regularization to the logistic model using
LASSO (Least Absolute Shrinkage and Selection Operator) [43]. This either
excludes or reduces the magnitude of the covariates’ coefficients that are not
predictive for the response mechanism in the propensity model. This helps to
minimize the variance of the estimated weights, at the potential expense of some
consistent (hopefully small) bias in their estimated values. However, this process
doesn’t exclude covariates that are uncorrelated with the response itself. These
should be excluded by the researcher in order to avoid variance inflation [42].
Another protective measure against variance inflation and extreme weights is
4 Note that the variance estimator of the weighted mean presented in subsection 3.4.5 is a

closed form formula that assumes fixed weights, and hence the variability in the estimation of
the weights is not reflected in the formula. A more accurate estimator of the variance would
rely on bootstrap samples, which are more computationally expensive.

13
weight trimming. balance offers automatic trimming, for details see subsection
4.3.

Another weakness of inverse propensity score weighting is that it may strongly

depend on the specification of the model of the propensity scores, as shown in
a simulation study in [10]. Imai and Ratkovic [10] have suggested the method
of Covariate Balancing Propensity Score (CBPS) (described in the next sub-
section) to overcome this issue. Fitting tree-based methods for the propensity
scores have also shown to be a good alternative [44, 40].

3.3.4 Covariate Balancing Propensity Score (CBPS)

Covariate Balancing Propensity Score (CBPS), suggested by Imai and Ratkovic

[10], is a method to estimate the propensity score in a way that will also result
in maximizing the covariate balance. The method is preferable in the cases of
misspecification of the propensity score model, which may lead to a bias in the
estimated weights (and consequently, the estimated survey statistic). CBPS
is described in details in [10] and implemented in the R package [45]. We give
here a short summary of the method for completeness of the estimation methods
section.

The CBPS method is an expansion of the maximization problem of logistic

regression. The propensity score of the logistic regression model is modeled by:
exp(β T Xi )
pβ (Xi ) = ∀i ∈ S, T (5)
1 + exp(β T Xi )

By the maximum-likelihood approach, β is estimated by maximizing the log-

likelihood, which results in:
X X
β̂M LE = arg max log(pβ (Xi )) + log(1 − pβ (Xi ))
β
i∈S i∈T

At the maximum of the log-likelihood β satisfies first order condition:

" #
1 X p′β (Xi ) X −p′β (Xi )
+ =0
n pβ (Xi ) 1 − pβ (Xi )
i∈S i∈T

where p′β (Xi )

is the derivative of p by β T . This condition can be viewed as
the condition that balances a certain function of the covariates, in this case the
derivative of the propensity score p′β (Xi ).

Generally, one can expand the above to hold for any function f of the covari-
ates X (f (X)) and depends on the researcher’s goals and assumptions:
( )
X f (Xi ) X f (Xi )
E + =0 (6)
pβ (Xi ) 1 − pβ (Xi )
i∈S i∈T

14
CBPS method chooses f (X) = X as the balancing function f in order to bal-
ance the first moment of each covariate in addition to the to derivative of the
propensity score parametric model. The estimation of the propensity score is
then done by using Generalized Methods of Moments (GMM) [46] and is de-
scribed in [10].

3.4 Evaluation of survey weights

3.4.1 Overview

As mentioned, survey weights are essential to improve the accuracy of survey

estimates, but their reliability and validity hinge on several assumptions and
modeling decisions.

Survey weights are valuable when: (a) the non-response pattern is sufficiently
captured by the measurable covariates, (b) the covariates are accurately repre-
sented in the fitted propensity score model, ensuring that the weighted distribu-
tion of covariates in the sample closely resembles that in the target population,
and (c) the survey weights correlate with the outcome of interest to an extent
that justifies the increased variance resulting from the weights [47, 42].

To see the level to which survey data empirically complies with the above
criteria, diagnostics measures can be applied on each of the main elements:
covariates, weights, and outcome.

Both covariates and outcomes can be checked before and after applying the
weights, allowing for a comprehensive assessment of the weights influence on
the data. Such evaluations helps confirm whether the weights have successfully
enhanced the representativeness of the sample data in a way that also sub-
stantially influences the outcome of interest. Additionally, various diagnostics
can be performed on the weights themselves to understand if their are extreme
weights which dominate the sample, as well as the overall impact the weights
have on the effective sample size.

Distributions can be compared using summary statistics and plots, with calcu-
lations incorporating the fitted survey weights. The following sections describes
various methods for that purpose.

3.4.2 Visualizing Distributions

Distribution plots are effective tools for visualizing the covariates and outcomes
in the data, offering insights that extend beyond basic summary statistics. For
numerical variables there are Kernel Density Estimator (KDE) plots (see an
example in Fig 2), histograms, and quantile-quantile (QQ) plots. For categor-

15
ical variables it is common to use bar-plots (see an example in Fig 2). These
distribution plots enable users to observe the differences between the observed
sample and the target population, as well as the influence of the applied weights.

(a) Bar plot for gender

(b) KDE for income

Un-weighted sample Weighted sample True value

Figure 2: Examples (from simulated data) of diagnostic plots for covariates

The advantage of visualizations lies in their ability to reveal unexpected pat-

terns in the complete range of data, as opposed to looking on summary statistics
only. However, scaling these visualizations can be challenging. For instance,
while examining KDE plots for each covariate comparing the sample and tar-
get population is informative, it is often more efficient for the researcher to
have summary statistics that can quickly convey the extent of bias in different
features. This is particularly useful when evaluating multiple weighting solu-
tions. The following sections discuss particular summary statistics that helps
in addressing this need.

16
3.4.3 Diagnostics for the covariates using ASMD

A fundamental statistic for comparing distributions is the first moment, i.e. the
mean, of each covariate for the target and the sample (weighted or unweighted).
For each covariate, it is insightful to observe how much closer the application
of weights brings us to the target mean. The Absolute Standardized Mean
Deviation (ASMD) can be used to summarize this effect.

The Absolute Standardized Mean Deviation (ASMD) is a statistical measure

employed to compare the means of two groups (in our case, the sample and the
target). It is computed as follows:

X̄Sample − X̄T arget

ASM D = (7)
SD

where X̄Sample and X̄T arget are the means of the sample and target. The
SD can be either the pooled standard deviation of sample and target, or the
standard deviation of the target population. In balance we use the standard
deviation of the target population.

The concept of ASMD is derived from the standardized mean difference, which
is a measure of effect size used to compare the means of two groups, expressed
in the standard deviation units. This is often referred to as Cohen’s d [48], a
standardized measure of the magnitude of the difference between two means.
ASMD values range from 0 to infinity, with larger values indicating greater
differences between the means of the two groups.5

The ASMD can be calculated using the unweighted and weighted mean of
the Sample, and these two quantities can be compared. If applying the weights
lead to an ASMD value that is closer to 0 than the ASMD of the unadjusted
sample, then it is an indication the weights help to reduce the bias. The level
of adjustment can be measured by taking the difference of these two ASMD
values.

ASM Ddif f = ASM Dunadjusted − ASM Dweighted (8)

The more the adjusted ASMD (the ASMD that is based on the weighted mean
of the sample) is smaller than the unadjused ASMD (based on the unweighted
mean) - i.e.: the closer the diff is to 0 - the stronger the indication we have
of the potential benefit of the weights for adjusting a bias in the covariates.
The magnitude of the difference of the two ASMD values is a measure of the
5 A value of 0 signifies no difference between the means, while a value of 1, for example, indi-

cates that the difference between the means is equal to one standard deviation. ASMD is most
easily conceptualized when the distributions being compared are unimodal and symmetric.

17
impact of the weights. If ASM Ddif f is positive then it means the weights have
helped reduce the bias, while a negative value indicates that the weights have
potentially increased the bias.

Since we often wish to adjust over many covariates, then the ASMD differ-
ence from each covariate can be summarized by taking the average ASMD (or
ASM Ddif f ) over all covariates. This gives a single summary statistic to mea-
sure the level of impact the weights had on reducing the bias of the sample in
the covariates.

For categorical variables, one possible behavior for ASMD calculation is to

use dummy variables and calculates the ASMD for each of them. The ASMD
per dummy variable approach could lead to over-weighting categorical variables
with many categories when calculating the mean ASMD. A possible solution is
to aggregate these ASMD per covariate. I.e.: calculates a single mean ASMD
value for each categorical variable, and then the general mean ASMD will give
each variable the same weight in the final calculation.

For some limitations of using ASMD, see the appendix section B.

3.4.4 Diagnostics for the weights

Kish’s design effect

One important aspect to consider when using survey weights is the potential
increase in variance of some estimate of interest (e.g.: the mean) due to the
variability in the weights. This is measured by a quantity known as design
effect, which is generally defined as the variance of the weighted estimator to
the variance expected with simple random sample (SRS) without replacement
[49, 50]. It assesses the potential impact that weights might have on the variance
of estimating the weighted mean.

Kish’s design effect [51] is a widely known and commonly used design effect
measure for the potential impact that weights might have on the variance of
estimating the population mean using the weighted mean. Kish’s design ef-
fect assumes that there is no correlation between the weights and the outcome
variable, also known as ”haphazard weights.”, its formula is:

Pn 1
Pn 2
n w2 i=1 wi w2
Def f = Pn i=1 i2 = n

2 = (9)
w2
Pn
( i=1 wi ) 1
n i=1 wi

The effective sample size proportion (ESSP) indicates what is the effective
proportion of sample size we’ll keep after applying the weights. It’s simply the
inverse of Def f :

18
 1
ESSP = (10)
Def f

The effective sample size is a related measure that takes into account both
the design effect and the actual sample size. It can be used to approximate the
number of independent observations that would yield the same variance as the
weighted sample. The effective sample size is calculated as follows (where n is
the sample size):

n
nef f = ESS = (11)
Def f

The effective sample size provides a useful way to gauge the impact of the
weights on the precision of the estimates. A smaller effective sample size in-
dicates that the weights have introduced greater variability in the estimates,
potentially requiring a larger actual sample size to achieve the desired preci-
sion.

Further details on assumptions and proofs are available in appendix C.

Summary Statistics for the Distribution of Weights

While Kish’s design effect can be used to estimate an effective sample size as
a summary measure for the impact of using weights, it may also be beneficial
to examine the distribution of weights using other summary statistics. For in-
stance, extremely large or small weights could indicate potential issues with the
weighting process or the presence of outliers in the data used for estimating the
weights. Furthermore, the distribution of weights can help determine whether
the weights are concentrated on a small number of observations or more evenly
distributed across the sample. These observations are often easier to infer from
summary statistics than from distribution plots of the weights. Understanding
the distribution of the weights can also help to better understand Kish’s design
effect (and effective sample size) value, which may indicate whether follow-up
manipulation of the weights is necessary (such as using an alternative weighting
model or weight trimming).

For diagnostic purposes, it is often more convenient to examine the weights af-
ter they have been normalized so that their sum equals the sample size. i.e.: by
dividing each weight in the sample by the average of the weights (wi∗ = wi /w̄).
When weights are normalized to sum to the sample size, they have the appealing
property of being more or less informative as they deviate from 1. A weight
smaller than 1 for an observation indicates that the weighting procedure consid-
ers this observation less informative than the average observation. Conversely,
a weight larger than 1 suggests that this observation is more informative, on
average, than other observations.

19
 For instance, if we have weights based on gender and find that males have
weights smaller than 1 while females have weights larger than 1, we can infer
that our sample has an over-representation of males and an under-representation
of females - an imbalance that the weights attempt to rectify.

It is helpful to look at the distribution of the weights. Looking at the KDE

plot can help detect multimodal distribution (which might indicate clusters of
users of higher/lower representativeness of the population). It is also helpful to
look at basic summary statistics, such as the main quartiles (25%, 50%, and
75%) as well as the proportion of weights above and below certain values (e.g.,
over 2 and under 0.5, along with other similar quantities). This can help identify
which proportions of the responses might be over/under weighted. Such insights
could lead to followup changes to the final weighting model. For example, if we
find out a handful of users have weights that are extremely large we might decide
to look at the skewed features. We might find a need to bucket some classes in a
covariate together, remove some features from the weighting model, use weight
trimming, or some other post-processing manipulation to the weights.

3.4.5 Diagnostics for the outcome

The entire procedure of fitting weights and diagnostics is geared towards an

impactful change in the outcome (or outcomes) of interest towards reducing the
estimation bias. A common population parameter of interest is the mean. The
statistics used to review it are the sample weighted mean, the variance of the
weighted mean, as well as asymptotic confidence intervals.

The formula for the weighted mean, using a Horvitz–Thompson estimator

[39], is simply:

Pn
wi yi
ȳw = Pi=1
n (12)
i=1 wi

The variance of the weighted mean is based on the π-estimator for the ratio-
mean:[52]

n
1 X
V\
(ȳw ) = Pn w2 (yi − ȳw )2 (13)
( i=1 wi )2 i=1 i

This estimator works for cases when the probability of selection for each yi
are not identical, treating the yi values themselves as fixed.6 See section D for
more details.
6 The formula presented for the variance of the weighted mean assumes that the weights

are known and fixed quantities. Hence, this formula does not account for the uncertainty that
is introduced from the estimation of the weights. If measuring this uncertainty is of interest,

20
 The confidence intervals (CI) available uses the above formula and are the
standard approximate CI based on the central limit theorem:

q
CI(µ) : ȳw ± zα/2 V\
(ȳw ) (14)

In an applied setting, it is advisable to calculate the weighted mean and

their CI after applying the weights, and also without weights, and compare the
quantities to each other. The difference of the weighted and unweighted mean
could be thought of as an estimator of the potential bias reduced by using the
weights (assuming the general trend of the ASMD calculations on the covariates
indicate a positive improvement in their imbalance). This estimated bias can be
compared to the effective sample size to allow a rough decision if the increase
in variance due to the weights is adequately compensated by the reduction in
bias.

4 The balance workflow

4.1 The workflow

Survey data weighting using balance is achieved with the following three main
steps:

1. Understanding the initial bias in the data relative to a tar-

get population: First, the survey data is loaded for both respondents
and the target population. A pandas DataFrame can be created using
[Link] csv() and converted into a balance Sample class object
with [Link] frame. A similar step is repeated for the target pop-
ulation’s data, and then the two Sample objects can be combined by as-
signing the target object as the target of the sample object. Once the
data is loaded, we can conduct a diagnostic evaluation of the sample-
vs-target covariates’ distributions to determine if weighting is necessary.
These include ASMD and distribution plots such as bar-charts and kernel-
density-estimation plots.
2. Adjusting the sample to the target: next, we generate weights for
the sample to more accurately represent the target population’s distribu-
tions. Currently, the package implements the following methods: Inverse
then it is possible to perform an end to end bootstrap simulation which includes re-sampling
from the sample, calculating the weighted mean estimation, and then repeating the process a
few times, and using the bootstrap estimations of the mean to estimate the variance.

21
 Probability Weighting (IPW) using LASSO regression, Covariate Balanc-
ing Propensity Score (CBPS), Post-stratification, and raking. These are
all available through the adjust method in the Sample class.
3. Results evaluation: Once the weights are estimated, their effect is eval-
uated on the covariates imbalance (again, using ASMD and plots), the
effective sample size, and the change in weighted mean of the outcome as
well as their confidence intervals.

The next section gives a detailed example for applying this workflow.

4.2 An end-to-end example

4.2.1 Understanding the initial bias

Loading simulated data

This section presents an example of simulated data extracted from the balance
tutorial page [53]. The data set is comprised of two pandas DataFrames:
one for the target population and the other for the sample population. Both
DataFrames contain an identifier column (id), three covariate columns (gender,
age group, and income), and an outcome variable (happiness)7

In this particular simulation, we intentionally designed the outcome to be

associated with all covariates, ensuring that this relationship remains consistent
for both the target and sample populations. It is important to note that in
real-world data sets, we generally don’t observe the outcome for the target pop-
ulation. However, in this simulated data set we have included it for illustrative
purposes. This setup allows us to later demonstrate how weighting methods can
mitigate bias and approximate population-level parameters more accurately.

In real-world use-cases the data is often loaded using [Link] csv().

Here, we use pre-made DataFrames that can be loaded (and inspected) using
the following Python code:
1 from balance import load_data
2 # INFO (2023 -05 -14 09:00:15 ,410) [ __init__ / < module > ( line 52) ]:
Using balance version 0.9.0
3
4 target_df , sample_df = load_data ()
5
6 print ( " sample_df : \ n " , sample_df . head () )

sample_df:
7 Code for creating the distributions is available here: [Link]
facebookresearch/balance/blob/main/balance/datasets/__init__.py#L17.

22
 id gender age_group income happiness
0 0 Male 25-34 6.428659 26.043029
1 1 Female 18-24 9.940280 66.885485
2 2 Male 18-24 2.673623 37.091922
3 3 NaN 18-24 10.550308 49.394050
4 4 NaN 18-24 2.689994 72.304208

# The target_df DataFrame looks similarly to sample_df.

Creating instances of the Sample class with the DataFrames

The main class for our analyses is the Sample class from the balance package.
The following illustrates how we incorporate the DataFrames into this class:
1 from balance import Sample
2
3 sample = Sample . from_frame ( sample_df , ou tc om e _c ol u mn s =[ " happiness "
])
4
5 target = Sample . from_frame ( target_df , ou tc om e _c ol u mn s =[ " happiness "
])
6 # Usually the code will be simply :
7 # target = Sample . from_frame ( target_df )
8 # This is since most times we do not have the outcome for the
target . In the example in this paper we have added it just to
validate later that the weights indeed help us reduce the bias
of the outcome .
9
10
11 # Following this , we associate the Sample object instance of sample
with that of the target object , enabling us to adjust the
sample to match the target .
12 s a m p l e _ w i t h _ t a r g e t = sample . set_target ( target )

The Sample class provides a wide range of attributes, methods, and properties.
For instance, the df property can reveal the DataFrame encapsulated within the
instance of the Sample class (e.g.: sample with [Link]):

Invoking the Sample object directly provides a concise summary of its at-
tributes:
1 sample_with_target

(balance.sample_class.Sample)

balance Sample object with target set

1000 observations x 3 variables: gender,age_group,income
id_column: id, weight_column: weight,
outcome_columns: happiness

23
 target:

balance Sample object

10000 observations x 3 variables: gender,age_group,income
id_column: id, weight_column: weight,
outcome_columns: happiness

3 common variables: gender,age_group,income

Exploring the imbalances in covariates

We can use methods such as .covars() with .plot(), .mean(), and .asmd()
to get some diagnostics about the imbalance.

We can use the .plot() method to look at the distributions of covariates in

the sample versus the target data.
1 s a m p l e _ w i t h _ t a r g e t . covars () . plot ()

The output in Figure 3 helps to easily identify imbalance. For example, we

can see the sample has many more males than females, as opposed to a 50%-50%
split in the target population. And for age group we can see how the sample is
skewed towards younger respondents, as compared to the target population.

The package leverages plotly [54] (as the default) to create interactive visu-
alizations, but it also supports static figures using the seaborn package [55] for
added flexibility.

The default asmd method uses ASMD to compare sample (which is un-
weighted) with the target using dummy variables for categorical variables, and
calculates the ASMD for each of them. The aggregate ASMD per covariate
can be achieved using the aggregate by main covar = True argument, as de-
scribed in section 3.4.3.
1 print ( s a m p l e _ w i t h _ t a r g e t . covars () . asmd ( a g g r e g a t e _ b y _ m a i n _ c o v a r =
True ) . T . round (2) )

source self
age_group 0.23
gender 0.25
income 0.49
mean(asmd) 0.33

The ASMD helps quantify the levels of imbalance in each covariate.

24
 (a) Bar plot for gender

(b) Bar plot for age group

(c) KDE for income

Un-weighted sample Weighted sample

Figure 3: Examples (from simulated data) of diagnostic plots for covariates

(unweighted sample vs target)

25
 4.2.2 Fitting survey weights

In order to estimate weights for the sample the .adjust() method as used on
the Sample object. The default is ipw, and other methods could be invoked
using the method argument.
1 # Using ipw to fit survey weights
2 adjusted = s a m p l e _ w i t h _ t a r g e t . adjust ()

4.2.3 Evaluating the Results

Covariates

We can get a basic summary of the results using the .summary() method:
1 print ( adjusted . summary () )

Covar ASMD reduction: 59.7%, design effect: 1.897

Covar ASMD (7 variables): 0.327 -> 0.132
Model performance: Model proportion deviance explained: 0.172

It shows that the weights led to an improvement of around 60% reduction

in the mean ASMD (from 0.327 to 0.132), and that the price we paid for it is
an increasing the variance of the estimator by 1.897 in comparison to a random
sample (as calculated using Kish’s design effect, if assuming haphazard weights).

The same tools used to evaluate the bias before adjustment can be used for
evaluating the effect of the weights on the balance after adjustment.
1 adjusted . covars () . plot ()

The output in Figure 4 shows how the weights help mitigate some (though
not all) of the bias, for all three covariates (gender, age and income).

We can also see the improvement per caovariate (averaged across category)
using the .asmd() method:
1 print ( adjusted . covars () . asmd ( a g g r e g a t e _ b y _ m a i n _ c o v a r = True ) . T .
round (2) )

source self unadjusted unadjusted - self

age_group 0.06 0.23 0.18
gender 0.10 0.25 0.16
income 0.24 0.49 0.25
mean(asmd) 0.13 0.33 0.20

26
 (a) Bar plot for gender

(b) Bar plot for age group

(c) KDE for income

Un-weighted sample Weighted sample True value

Figure 4: Examples (from simulated data) of diagnostic plots for covariates

(unweighted and weighted sample vs target)

27
 We can see that while we got improvements in all covariates, there is still
some imbalance that remained, especially in the income variable.

Weights

Next, we wish look at the diagnostics of the weights to identify if there are
any extreme weights or signs of issue that requires further investigation. This
can be done by using the summary method on the .weights() method of the
adjusted object.
1 print ( adjusted . weights () . summary () . round (2) )

var val
0 design_effect 1.90
1 effective_sample_proportion 0.53
2 effective_sample_size 527.04
...
7 describe_min 0.31
11 describe_max 11.65
16 prop(w < 1) 0.65
21 prop(w >= 10) 0.00

We can see a design effect of 1.9 which corresponds with an effective sample
size proportion of 53%. Since the size of the sample was 1000, it means that the
effective sample size is 527. We can also see that 65% of the weights are below
1, meaning that we down-sized 65% of our sample. The minimal weight is 0.31
and the max weight is 11.65, with almost no weights above 10. A conclusion
here is that the weights are not too extreme and we get some sense of the cost
that using the weights would incur on the precision of our estimates.

Outcome

The summary method on the outcomes method gives us the weighted means
and confidence intervals for the adjusted sample, the target, and the unadjusted
sample data.

From the results below we can see that the real population level mean of
happiness in the simulation was 56.2. In our (unweighted/unadjusted) sample
it was 48.5. Meaning, the bias was roughly 7.7 points. After applying the
weights, we got a value of 53.3, reducing the bias to roughly only 2.9 points.
Note that this comparison is only possible in a simulated environment and is
given here for a proof of concept of the effect of the weights. We can also see
that the CI of the self and unadjusted show very different ranges of bands,
indicating how the weights clearly got us a significant change in the estimated
mean8 . While the model improved the bias, we know it didn’t fix it completely.
8 Comparing the CI of the data with and without the weights is a good approximation

28
 This is because the model also did not perfectly fix the covariate imbalance,
since it used some regularization in the process.
1 print ( adjusted . outcomes () . summary () )
2 adjusted . outcomes () . plot ()

Source Happiness
Self 53.389
Target 56.278
Unadjusted 48.559
Self CI (52.183, 54.595)
Target CI (55.961, 56.595)
Unadjusted CI (47.669, 49.449)

The output of .plot is in Figure 5. It shows that we got a relatively sym-

metrical uni-modal distribution (before and after applying the weights). So we
don’t observe and strong irregular behavior of the outcome. Note that we are
able to compare the outcome with and without the weights in the sample to the
real outcome distribution in the target population only because this is simulated
data. In real-world cases, we are not expected to have access to the outcome
distribution of the target population. Also, it is relatively common to get out-
come responses in binary or likert scales, and not a continuous variable. The
.plot would work with these just as well.

Un-weighted sample Weighted sample True value

Figure 5: Examples (from simulated data) of diagnostic plots for outcome (un-
weighted and weighted sample vs target)

Downloading data
for the impact of the weights, but is not statistically precise. Future work is planned for
introducing more formal confidence intervals of the impact of the weights by using paired
t-test style analysis. See the discussion and future work section for more.

29
 Once we are settled with the weights we got, we can download them as csv,
as follows:
1 adjusted . to_download () # Will create a download link in jupyter
2 # We can also prepare the data to be exported as csv
3 # The following code showes the first 500 characters for simplicity
:
4 adjusted . to_csv ()

4.3 How does balance implement the adjustment?

Pre-processing Before applying any of the adjustment methods, balance per-

forms a pre-processing step to improve models’ results. The pre-processing step
includes a few best practiced that makes the use of balance easy and automatic
for a default usage.

Transformations. balance applies the following default behaviours:

1. Handling missing values: balance handles missing values automatically

by adding a special indicator column to any variable that contains miss-
ing values. The advantage of this is that these are then considered as a
separate category for the adjustment.
2. Feature engineering: by default, balance applies feature engineering to be
able to fit the covariate distribution better, and not only the first moment.
Specifically, each continuous variable is bucketed into 10 quantiles buckets,
and rare categories variables are grouped together so to avoid overfitting9 .

Model matrix. The model matrix of the covariates used in balance for the
logistic regression in ipw and for CBPS propensity score is constructed before
the fitting is done using the transformed variables and one-hot encoding for
discrete variables. The default behaviour is an additive model including all
joint covariates of the target and the observed sample. However, thorough the
argument formula, one can input a formula for a specified relation between the
variables. The formulas adopt the notation from the patsy Python package [56],
facilitating a range of operations like addition, multiplication (for interaction
effects), and power transformations. A detailed example is available in the
”balance: transformations and formulas” tutorial [57].

Adjustment through ipw

ipw is implemented using LASSO regularized logistic regression. To avoid

non-balanced categories in the logistic regression, balance scales the prevalence
of the target population to be similar to the observed sample.
9 The user has also an option to change these default behaviours, through setting different

values to the transformations argument of [Link].

30
 The penalty factor λ of the LASSO is chosen through cross-validation. Two
methods for choosing the parameter are suggested:

1. Unbounded Design Effect: If one doesn’t want to bound the design effect
of the resulted weights (the default behaviour with max de=None), the
penalty if chosen using lambda 1se, which is the largest value of λ such
that the cross-validated error is within one standard error of the minimum
value.
2. Bounded Design Effect: If one chooses to bound the design effect (e.g. by
using max de=2), a grid search over 10 of the values of λ that brings the
largest design effect within the bound is done, where the λ is chosen to be
the one that brings the largest ASMD reduction.

In addition, a penalty factor argument can be also used to indicate how

much the model should focus to adjust each term of the formula. Larger penalty
factors means that the covariate is more likely to be regularized by the LASSO
penalty and as a result the adjustment of this covariate will be smaller, i.e. will
end in a less balanced covariate. This feature can be particularly useful when
certain components are believed to be more or less responsible for bias in the
data, or when the user wants to explore different adjustment scenarios.

Post-processing

Weights in balance are scaled to the population size after estimated, and
can be interpret as the number of people from the target the sample unit rep-
resent. After the adjustment and scaling is done, weights trimming from above
is performed. This is done in order avoid over fitting of the model and unnec-
essary variance inflation. The weights are trimmed and scales in a way that
keeps the mean and sum of the weights the same as before trimming, such that
the interpretation how many units in the target this unit represent holds after
trimming.

5 Future directions

The balance package offers benefits for researchers interested in analyzing data
with non-response bias in the Python environment by being easy to use, provid-
ing an end-to-end workflow, and released as open-source. While comprehensive,
there is still room for improvement and expansion. This section highlights sev-
eral possible areas for future development in the balance package.

1. Better Diagnostic Tools for Covariates: The current metric of ASMD

has limitations, especially when applied to a wide range of distributions

31
 and for categorical variables. Future versions could include more robust
measures like the Common Language Effect Size [58] and better methods
for handling categorical variables, such as Kullback-Leibler divergence.
There is also room for adding statistical hypothesis tests for the evalua-
tions, as well as more plots. The cobalt R package [18] is a good source
of inspiration.
2. Expanded Estimation and Diagnostic Tools for Outcomes: Cur-
rently, the package primarily provides the weighted mean and its confi-
dence intervals. A helpful improvement would be to directly measure the
estimated bias reduction caused by the weights, including a confidence
interval for this estimate. Also, current implementation focuses on the
weighted average and the linearization (Taylor) estimator for the vari-
ance. Other possible statistics, and estimations of variance exists. The
samplics package already implements some of these and would be a good
source of inspiration [21].
3. Diagnostics for the bias-variance trade-offs when using weights:
At present, the user is provided with a set of weights but with no easy
way to check the bias-variance tradeoffs for alternative levels of trimming
or tuning other parameters. A future version of the package could include
more diagnostics tools and allow automated functions for weight trimming,
such as based on empirical-MSE estimation for a given outcome over a
range of potential weight trimming values. This could lead to a better
balance between the variance induced by the weights and the bias they
reduce and save researcher’s time in manual tweaking.
4. Built-in Model Comparison for Multiple Weights: Our ultimate
goal is to allow the most flexibility to the user by conducting easy com-
parisons of multiple models and adjustments to the weights in order to
choose the model that best fits his/hers data.
5. Feature Selection for Propensity Score Models: When given several
potential models, it can be challenging to choose the best one. This choice
could depend on various factors, such as the balance between reduced
bias and incurred variance or the impact of different models on different
outcomes. Further development in this area could provide useful tools for
sensitivity analysis and decision making.
6. Expansion Beyond Propensity Score Models: The next step for the
package could be to include outcome models and doubly robust models.
Thus making the package more versatile and comprehensive.

These possible improvements represent exciting opportunities for the future of

the balance package, aiming to provide a more robust and user-friendly tool for
researchers in the Python environment. We welcome any feedback, suggestions,
and opportunities for collaborations.

32
References
[1] C. M. Astley, G. Tuli, K. A. Mc Cord, E. L. Cohn, B. Rader, T. J. Var-
relman, S. L. Chiu, X. Deng, K. Stewart, T. H. Farag, et al., “Global
monitoring of the impact of the covid-19 pandemic through online surveys
sampled from the facebook user base,” Proceedings of the National Academy
of Sciences, vol. 118, no. 51, p. e2111455118, 2021.

[2] R. M. Groves and L. Lyberg, “Total survey error: Past, present, and fu-
ture,” Public opinion quarterly, vol. 74, no. 5, pp. 849–879, 2010.
[3] R. J. Little, “Post-stratification: a modeler’s perspective,” Journal of the
American Statistical Association, vol. 88, no. 423, pp. 1001–1012, 1993.

[4] J.-C. Deville and C.-E. Särndal, “Calibration estimators in survey sam-
pling,” Journal of the American statistical Association, vol. 87, no. 418,
pp. 376–382, 1992.
[5] G. Kalton, Compensating for missing survey data. Inst for Social Research
the Univ, 1983.

[6] A. Mercer, A. Lau, and C. Kennedy, “For weighting online opt-in samples,
what matters most?,” 2018.
[7] M. David, R. J. Little, M. Samuhel, and R. Triest, “Nonrandom nonre-
sponse models based on the propensity to respond,” in Proceedings of the
Business and Economic Statistics Section, American Statistical Associa-
tion, pp. 168–173, 1983.
[8] R. J. Little, “Survey nonresponse adjustments for estimates of means,” In-
ternational Statistical Review/Revue Internationale de Statistique, pp. 139–
157, 1986.
[9] A. Ekholm and S. Laaksonen, “Weighting via response modeling in the
finnish household budget survey,” Journal of Official statistics, vol. 7, no. 3,
pp. 325–337, 1991.
[10] K. Imai and M. Ratkovic, “Covariate balancing propensity score,” Journal
of the Royal Statistical Society: Series B: Statistical Methodology, pp. 243–
263, 2014.

[11] A. Gelman and J. Hill, Data analysis using regression and multilevel/hier-
archical models. Cambridge university press, 2006.
[12] G. Solon, S. J. Haider, and J. M. Wooldridge, “What are we weighting
for?,” Journal of Human resources, vol. 50, no. 2, pp. 301–316, 2015.

[13] T. Robert, “Regression shrinkage and selection via the lasso,” Journal of
the Royal Statistical Society Series B, vol. 58, p. 267, 1994.

33
[14] M. Templ, A. Kowarik, and T. Schoch, “CRAN Task View: Official Statis-
tics & Survey Statistics.” Version 2023-02-19, 2023.
[15] T. Lumley, “Analysis of complex survey samples,” Journal of Statistical
Software, vol. 9, no. 1, pp. 1–19, 2004. R package verson 2.2.
[16] I. Mayer, P. Zhao, N. Greifer, N. Huntington-Klein, and J. Josse, “CRAN
Task View: Causal Inference.” Version 2022-12-07, 2022.
[17] T. Zhou, G. Tong, F. Li, L. E. Thomas, and F. Li, “The r journal: Psweight:
An r package for propensity score weighting analysis,” The R Journal,
vol. 14, pp. 282–300, 2022. [Link]
[18] N. Greifer, “Covariate balance tables and plots: a guide to the cobalt
package,” Accessed March, vol. 10, p. 2020, 2020. R package version 4.5.1.
[19] D. Forthomme, “ipfn: Iterative proportional fitting for python with n di-
mensions, 2021.” GitHub repository. [Link]
[20] G. Freysson, “Quantipy3.” GitHub repository. [Link]
Quantipy/quantipy3.
[21] M. S. Diallo, “samplics: a python package for selecting, weighting and
analyzing data from complex sampling designs.,” Journal of Open Source
Software, vol. 6, no. 68, p. 3376, 2021.
[22] J. Friedman, “Pysurvey.” PyPI, 2013. [Link]
PySurvey/.
[23] P. Wildeford, “Surveyweights.” GitHub repository. [Link]
rethinkpriorities/surveyweights.
[24] K. Ottoboni, “A python package for propensity score matching.” GitHub
repository. [Link]
[25] B. Miroglio, “A python package for propensity score matching.” GitHub
repository. [Link]
[26] P. Milica, “A python package for propensity score matching.” GitHub
repository. [Link]
[27] A. Sharma and E. Kiciman, “Dowhy: An end-to-end library for causal
inference,” arXiv preprint arXiv:2011.04216, 2020.
[28] X. Wang, J. Miao, and Y. Sun, “A python library for empirical calibration,”
arXiv preprint arXiv:1906.11920, 2019. [[Link]].
[29] H. Chen, T. Harinen, J.-Y. Lee, M. Yung, and Z. Zhao, “Causalml: Python
package for causal machine learning,” 2020.
[30] M. J. Salganik, Bit by bit: Social research in the digital age. Princeton
University Press, 2019.

34
[31] D. Freedman, R. Pisani, and R. Purves, Statistics Fourth Edition. WH
Norton & Company New York, 2007.
[32] P. Squire, “Why the 1936 literary digest poll failed,” Public Opinion Quar-
terly, vol. 52, no. 1, p. 125, 1988.
[33] D. Lusinchi, ““president” landon and the 1936 literary digest poll: Were au-
tomobile and telephone owners to blame?,” Social Science History, vol. 36,
no. 1, pp. 23–54, 2016.
[34] S. Lee and R. Valliant, “Estimation for volunteer panel web surveys us-
ing propensity score adjustment and calibration adjustment,” Sociological
Methods & Research, vol. 37, no. 3, pp. 319–343, 2009.
[35] D. B. Rubin, “Inference and missing data,” Biometrika, vol. 63, no. 3,
pp. 581–592, 1976.
[36] P. R. Rosenbaum and D. B. Rubin, “The central role of the propensity
score in observational studies for causal effects,” Biometrika, vol. 70, no. 1,
pp. 41–55, 1983.
[37] R. J. Little, “Missing data assumptions,” Annual Review of Statistics and
Its Application, vol. 8, pp. 89–107, 2021.
[38] M. P. Battaglia, D. C. Hoaglin, and M. R. Frankel, “Practical considera-
tions in raking survey data,” Survey Practice, vol. 2, no. 5, 2009.
[39] D. G. Horvitz and D. J. Thompson, “A generalization of sampling without
replacement from a finite universe,” Journal of the American statistical
Association, vol. 47, no. 260, pp. 663–685, 1952.
[40] S. Watkins, M. Jonsson-Funk, M. A. Brookhart, S. A. Rosenberg, T. M.
O’Shea, and J. Daniels, “An empirical comparison of tree-based methods
for propensity score estimation,” Health services research, vol. 48, no. 5,
pp. 1798–1817, 2013.
[41] F. Li, K. L. Morgan, and A. M. Zaslavsky, “Balancing covariates via
propensity score weighting,” Journal of the American Statistical Associ-
ation, vol. 113, no. 521, pp. 390–400, 2018.
[42] R. J. Little and S. Vartivarian, “Does weighting for nonresponse increase
the variance of survey means?,” Survey Methodology, vol. 31, no. 2, p. 161,
2005.
[43] R. Tibshirani, “Regression shrinkage and selection via the lasso,” Journal
of the Royal Statistical Society: Series B (Methodological), vol. 58, no. 1,
pp. 267–288, 1996.
[44] B. K. Lee, J. Lessler, and E. A. Stuart, “Improving propensity score weight-
ing using machine learning,” Statistics in medicine, vol. 29, no. 3, pp. 337–
346, 2010.

35
[45] C. Fong, M. Ratkovic, and K. Imai, CBPS: Covariate Balancing Propensity
Score, 2019. R package version 0.21.
[46] L. P. Hansen, “Large sample properties of generalized method of moments
estimators,” Econometrica: Journal of the econometric society, pp. 1029–
1054, 1982.

[47] L. Kish and J. O. Stat., “Weighting for unequal π,” Journal of Official
Statistics, vol. 8, pp. 183–200, 1992.
[48] J. Cohen, Statistical power analysis for the behavioral sciences. Academic
press, 2013.

[49] R. Arnab, “Chapter 19 - complex surveys: Categorical data analysis,” in

Survey Sampling Theory and Applications (R. Arnab, ed.), pp. 645–671,
Academic Press, 2017.
[50] K. A. Henry and R. Valliant, “A design effect measure for calibration
weighting in single-stage samples,” Survey Methodology, vol. 41, no. 2,
pp. 315–331, 2015.
[51] L. Kish, Survey Sampling. New York: John Wiley & Sons, Inc., 1965.
[52] C.-E. Sarndal, B. Swensson, and J. Wretman, Model Assisted Survey Sam-
pling. New York, NY: Springer, 1992.

[53] T. Galili and T. Sarig, “balance quickstart: Analyzing and adjusting the
bias on a simulated toy dataset.” [Link]
tutorials/quickstart/, 2023. [Online; accessed 14-May-2023].
[54] P. T. Inc., “Collaborative data science,” 2015.

[55] M. Waskom, O. Botvinnik, D. O’Kane, P. Hobson, S. Lukauskas, D. C.

Gemperline, T. Augspurger, Y. Halchenko, J. B. Cole, J. Warmenhoven,
J. de Ruiter, C. Pye, S. Hoyer, J. Vanderplas, S. Villalba, G. Kunter,
E. Quintero, P. Bachant, M. Martin, K. Meyer, A. Miles, Y. Ram,
T. Yarkoni, M. L. Williams, C. Evans, C. Fitzgerald, Brian, C. Fonnesbeck,
A. Lee, and A. Qalieh, “mwaskom/seaborn: v0.8.1 (september 2017),”
Sept. 2017.
[56] M. Wardrop, “Patsy: A python library for describing statistical models.”
[Link] 2023. Python library version (0.5.3).
[57] T. Galili and T. Sarig, “balance: transformations and formu-
las.” [Link]
transformations_and_formulas/, 2023. [Online; accessed 14-May-
2023].
[58] K. O. McGraw and S. P. Wong, “A common language effect size statistic,”
Psychological Bulletin, vol. 111, no. 2, pp. 361–365, 1992.

36
[59] M. S. Ben-Shachar, I. Patil, R. Thériault, B. M. Wiernik, and D. Lüdecke,
“Phi, fei, fo, fum: Effect sizes for categorical data that use the chi-squared
statistic,” Mathematics, vol. 11, no. 9, 2023.
[60] Wikipedia, “Design effect — Wikipedia, the free encyclopedia.”
[Link]
oldid=1137559092, 2023. [Online; accessed 27-April-2023].
[61] S. Gabler, S. Häder, and P. Lahiri, “A model based justification of kish’s
formula for design effects for weighting and clustering,” Survey Methodol-
ogy, vol. 25, pp. 105–106, 1999.

[62] Wikipedia, “Weighted arithmetic mean — Wikipedia, the free en-

cyclopedia.” [Link]
20arithmetic%20mean&oldid=1147722211, 2023. [Online; accessed 27-
April-2023].
[63] W. G. Cochran, Sampling Techniques. Nashville, TN: John Wiley & Sons,
3 ed., 1977.

37
Appendices
A Acknowledgments

The balance package was (and is) developed by many people, including: Roee
Eilat, Tal Galili, Daniel Haimovich, Kevin Liou, Steve Mandala, Adam Obeng
(author of the initial internal Meta version), Tal Sarig, Luke Sonnet, Sean Tay-
lor, Barak Yair Reif, and others.

The balance package was open-sourced by Tal Sarig, Tal Galili and Steve
Mandala, from Central Applied Science at Meta, in late 2022.

Branding created by Dana Beaty, from the Meta AI Design and Marketing
Team.

B Limitations of the ASMD

It is also worth noting some of the disadvantages in ASMD:

1. Sensitivity to extreme values: Since ASMD is based on the first moment,

ASMD can be sensitive to outliers or extreme values in the data, which
can lead to a distorted representation of the differences between the two
groups. This could be mitigated by turning to robust measure, but these
are currently not implemented in balance.
2. Inability to detect distributional differences: ASMD focuses solely on the
mean difference between the two groups, and does not account for differ-
ences in other distributional characteristics, such as variance, skewness or
number of modes. This means that two groups with similar means but
different variances or shapes may have a low ASMD value, which could be
misleading. This can be addressed by looking at distribution plots. The
next section discusses methods that are available in balance.

3. The need for context: ASMD values are unitless and can be difficult to in-
terpret without context. Does an ASMD value below 0.1 indicate an effect
size which is small or large? The interpretation of ASMD is often com-
parative within a specific research context. For example, ASMD changes
could be compared across different covariates and alternative weights, so
to identify which set of weight effects the bias of which covariate.
4. Limited applicability to categorical variables: ASMD is primarily designed
for comparing continuous variables, and its applicability to categorical

38
 variables is more limited. In such cases the covariate can be turned into
several dummy variables using one hot encoding and the ASMD can be
calculated on these values of zeros and ones. Alternative measures that
directly compare categorical distributions [59] are currently not imple-
mented in balance.

Despite these limitations, the ASMD can be a useful measure for comparing
the effect of the weights on the covariates.

C Kish’s design effect

C.1 Design effect in general

A design effect [60]10 is a measure of the increase in variance of an estimate due

to the use of survey weights compared to an equal probability sample of the
same size. Theoretically, it is calculated as follows:

V ar(θ̂weighted )
Def f = (15)
V ar(θ̂un−weighted )

where V arweighted and V arunweighted are the variances of the weighted and
unweighted estimates, respectively.

A design effect greater than 1 indicates that the variance of the weighted
estimate is larger than that of an unweighted estimate, while a design effect
less than 1 suggests that the variance of the weighted estimate is smaller (for
example, when using design weights based on stratified sampling). A design
effect of 1 implies that the use of weights does not affect the variance of the
estimate (which happens only if all weights are equal to the same, non 0, value).

A design effect that is larger than 1 does not necessarily imply that the weights
are undesirable, as they can still improve the accuracy and representativeness
of the estimates. Put differently, it may be that the bias corrected by applying
the weights is substantially larger than the variance added due to using them.

Kish’s design effect is a specific measure for a specific parameter (the popu-
lation mean), with specific assumptions. The following sections discusses these
some of these assumptions.
10 The text in this section is a modified version of the text we wrote for the Wikipedia article

on Kish’s design effect [60].

39
C.2 Assumptions and derivation

The formula of Kish’s design effect computes the increase in variance of the
weighted mean due to ”haphazard” weights, which occur when y consists of
observations selected using unequal probabilities, without within-cluster cor-
relation or any relationship to the expected value or variance of the outcome
measurement. From a model-based perspective [61], the formula holds when all
n observations (y1 , ..., yn ) are (at least approximately) uncorrelated and have
the same variance for some response variable of interest (y). The formula also
assumes that the weights are not random variables but rather known constants.
These, for example, can be the inverse of the selection probability for a pre-
determined and known sampling design.

The conditions on y are trivially satisfied if the y observations are independent

and identically distributed (i.i.d.) with the same expectation and variance. It
is important to note that if y1 , ..., yn do not have the same expectations, the
estimated variance of the estimator cannot be used for calculating the variance
using a simple weighted variance formula, as the estimation assumes that all
yi s have the same expectation. Specifically, if there is a correlation between
the weights and the outcome variable y, the expectation of y is not the same
for all observations but rather depends on the specific weight value for each
observation. In such cases, while the design effect formula might still be accurate
(assuming other conditions are met), a different estimator for the variance of
the weighted mean may be needed, such as a weighted variance estimator.

If different yi ’s have distinct variances, the weighted variance might capture

the correct population-level variance, but Kish’s formula for the design effect
may no longer be valid. A similar issue occurs if there is a correlation structure
in the samples, such as when using cluster sampling.

Kish’s formula estimates the increase in the variance of the weighted mean
based on ”haphazard” weights. Let y be observations selected using unequal
selection probabilities (with no within-cluster correlation, and no relationship
to the expectancy or variance of the outcome measurement);[47] and let y ′ be
the observations we would have had if we got them from simple random sample,
then Kish’s formula for Deff is:
 ! !
n
P n
var wi yi /
P
wi
var(ȳw ) i=1 i=1
Def f (kish) = var(ȳ ′ ) = 
n

yi′/n
P
var
i=1

From a model based perspective[61], this formula holds when all n observa-
tions (y1 , ..., yn ) are (at least approximately) uncorrelated (∀(i ̸= j) : cor(yi , yj ) =
0), with the same variance (σ 2 ) in the response variable of interest (y). It also as-
sumes the weights themselves are not a random variable but rather some known
constants (E.g.: the inverse of probability of selection, for some pre-determined
and known sampling design).

40
 The conditions on y are trivially held if the y observations are i.i.d with
the same expectation and variance. In such case we have y = y ′ , and we can
estimate var (ȳw ) by using var (ȳw ) = var (ȳ) × Def f [47].

C.3 Proof

We present here a simplified proof to Kish’s formula: Def f := var(ȳ w) w 2

var(ȳ ′ ) = w̄2 for
the case when there are no clusters (i.e.: no intraclass correlation between the
elements of the sample), so that each strata includes only one observation. The
proof is shown in full in [61] .

n
P 
wi yi n
! n
1  2 3
 i=1 X X
var (ȳw ) = var 
 P n
 = var
 wi′ yi = var (wi′ yi )
wi i=1 i=1
i=1
n
wi2
P
n n n
4 5 6 7
X X X i=1
= wi′2 var (yi ) = wi′2 σ 2 = σ 2
wi′2 = 2
σ  2
n
i=1 i=1 i=1
P
wi
i=1
n n n
(16)
wi2
P
wi2 wi2
P P
2 i=1
8 9 10 σ
= σ 2  i=1 2 = σ 2  ni=1 2 = n
2
n
n P
P n
P n wi wi
wi n  i=1  n2  i=1 
i=1 n n

σ 2 w2 12
11
= = var (ȳ ′ ) Def f
n w̄2
var (ȳw )
=⇒ Def f (kish) =
var (ȳ ′ )

Transitions:

1. from definition of the weighted mean.

2. using normalized (convex) weights definition (weights that sum to 1): wi′ =
wi
Pn .
wi
i=1

3. sum of uncorrelated random variables.

4. If the weights are constants (from the basic properties of the variance).
Another way to say it is that the weights are known upfront for each
observation i. Namely that we are actually calculating var (ȳw |w)

41
 5. assume all observations have the same variance (σ 2 ).

D Estimating the variance of the weighted mean

D.1 Formulation

This section discusses the derivation of the formula presented in the paper for
the variance of the weighted mean, also known as π-estimator for ratio-mean.11

We are interested in estimating the variance of the weighted mean when the
various yi are not assumed to be i.i.d random variables. An alternative perspec-
tive for this problem is that of some arbitrary sampling design of the data in
which units are selected with unequal probabilities (with replacement) [63].

Unlike classical ”model based” approaches, in which the randomness is de-

scribed by the randomness of the y value, here we consider the value of yi as
constant, where the variability comes from the selection procedure. We let Ri
be the Bernoulli indicator that is equal to 1 if observation i is in the observed
sample, and 0 if not. The probability of a unit to be sampled given a sample S of
size n is denoted by πi := P (Ri = 1 | S). Furthermore, we denote the one-draw
probability of selection by pi := P (Ri = 1|one sample draw) ≈ πni . For the
following derivation we’ll assume that the probability of selecting each element
is fully represented by these probabilities [52], i.e. selecting some element will
not influence the probability of drawing another element (this doesn’t apply for
things such as cluster sampling design).

Since each outcome yi is fixed, and the randomness comes from unit i being
included in the sample or not (Ri ), we often talk about the multiplication of
the two, which is a random variable. To avoid confusion in what to follow, we
define: yi′ = yi · Ri . This satisfies: E[yi′ ] = yi E[Ri ] = yi πi and V[yi′ ] = yi2 V[Ri ] =
yi2 πi (1 − πi ).

In this ”design based” perspective, the weights are obtained by taking the
inverse of the selection probability (i.e.: the inflation factor), i.e. wi = π1i ≈
1
n×pi . The weights in this setting are considered fixed and known.

We assume
Pn that the target population size N is unknown, and is estimated
by N̂ = i=1 wi . Our parameter of interest is the weighted mean, that can be
written as a ratio: PN yi PN
i=1 πi wi yi
Ȳ = PN 1 = Pi=1 N
(17)
i=1 πi i=1 wi

11 The text in this section is a modified version of the text we wrote for the Wikipedia article

on the weighted mean [62].

42
This ratio is estimated by the observed sample using:

Pn yi Pn
ˆ i=1 πi wi yi′
Ȳ = Pn 1 = Pni=1 (18)
i=1 πi i=1 wi Ri

This is called a Ratio estimator and it is approximately unbiased for Ȳ [52,

p.182]

In this case, the variability of the ratio depends on the variability of the
random variables both in the numerator and the denominator - as well as their
correlation. Since there is no closed analytical form to compute this variance,
various methods are used for approximate estimation, primarily Taylor series
first-order linearization, asymptotics, and bootstrap/jackknife.[52, p. 172] The
Taylor linearization method could lead to under-estimation of the variance for
small sample sizes in general, but that depends on the complexity of the statistic.
For the weighted mean, the approximate variance is supposed to be relatively
accurate even for medium sample sizes [52]. For when the sampling has a random
sample size, such as in Poisson sampling, it is as follows: [52]

n
1 X
V\
(ȳw ) = Pn w2 (yi − ȳw )2 (19)
( i=1 wi )2 i=1 i

We note that if πi ≈ pi n, then either using wi = π1i or wi = p1i would give

the same estimator, since multiplying wi by some factor would lead to the same
estimator. It also means that if we scale the sum of weights to be equal to a
known-from-before population size N , the variance calculation would look the
same. When all weights are equal to one another, this formula is reduced to the
standard unbiased variance estimator.

Note that for the trivial case in which all the weights are equal to 1, the
above formula is just like the maximum-likelihood formula for the variance of
the mean (but not that it is not the unbiased variance, i.e. dividing it by n
instead of (n-1)).

D.2 Proof

We show here a short proof for the variance formula presented above:
n
1 X
V\
(ȳw ) = Pn w2 (yi − ȳw )2 (20)
( i=1 wi )2 i=1 i

The Taylor linearization states that for a general ratio estimator, Q, of two
sums, Y and Z, can be expressed by: [52, p.178]

43
 Pn ′ n 
1 X yi′ zi′

Ŷ i=1 wi yi
Q̂ = = n ′ ≈Q+ Z −Q (21)
P
Ẑ i=1 wi zi i=1
πi πi

And the variance can be approximated by: [52, p.178]

n n
!  
\ 1 X X ˇ yi − Q̂zi yj − Q̂zj 1 \ \
V (Q̂) = ∆ij = V (Ŷ ) + Q̂V (Ẑ) − 2Q̂Ĉ(Ŷ , Ẑ)
Ẑ 2 i=1 j=1 πi πj Ẑ 2
(22)

where Ĉ(Ŷ , Ẑ) is the estimated covariance between the Y and Z, and ∆ij =
C(Ri , Rj ).

Since Ĉ is the covariance of two sums of random variables, it would include

many combinations of covariances that will depend on the indicator variables.
If the selection probability are uncorrelated (i.e.: ∀i ̸= j : ∆ij = C(Ri , Rj ) = 0),
this term would include only the summation of n covariances for each element
i between yi′ = Ri · yi and zi′ = Ri · zi . This helps illustrate that this formula
incorporates the effect of correlation between y and z on the variance of the
ratio estimators.

When defining zi = 1 the above becomes: [52, p.182]

n n  
\ 1 X X ˇ yi − ȳw yj − ȳw
V (Q̂) = V\
(ȳw ) = ∆ij (23)
N̂ 2 i=1 j=1 πi πj

∆
ˇ ij = ij . If the selection probability are uncorrelated (i.e.: ∀i ̸= j :
where ∆ πij
∆ij = C(Ri , Rj ) = 0), and when assuming the probability of each element is
very small (i.e.: (1 − πi ) ≈ 1), then the above reduced to the following:

n  2 n
1 X yi − ȳw 1 X
\
V (ȳw ) = (1 − πi ) = Pn w2 (yi − ȳw )2 . (24)
N̂ 2 i=1 πi ( i=1 wi )2 i=1 i

44