Blending of Probability and Non-Probability Samples: Applications to a Survey of Military Caregivers

TL;DR

This paper develops novel statistical methods to combine probability and non-probability samples, improving inference accuracy for small subpopulations, demonstrated through a survey of military caregivers with enhanced insights and reduced bias.

Contribution

The paper introduces new weighting techniques for blending probability and non-probability samples, addressing bias and variance issues in small subpopulation surveys.

Findings

Blended samples yield stronger conclusions about caregiver hardships.

Proposed methods improve estimate precision while controlling bias.

Simulation studies validate the effectiveness of the blending techniques.

Abstract

Probability samples are the preferred method for providing inferences that are generalizable to a larger population. However, when a small (or rare) subpopulation is the group of interest, this approach is unlikely to yield a sample size large enough to produce precise inferences. Non-probability (or convenience) sampling often provides the necessary sample size to yield efficient estimates, but selection bias may compromise the generalizability of results to the broader population. Motivating the exposition is a survey of military caregivers; our interest is focused on unpaid caregivers of wounded, ill, or injured servicemembers and veterans who served in the US armed forces following September 11, 2001. An extensive probability sampling effort yielded only 72 caregivers from this subpopulation. Therefore, we consider supplementing the probability sample with a convenience sample from the same subpopulation, and we develop novel methods of statistical weighting that may be used to combine (or blend) the samples. Our analyses show that the subpopulation of interest endures greater hardships than caregivers of veterans with earlier dates of service, and these conclusions are discernably stronger when blended samples with the proposed weighting schemes are used. We conclude with simulation studies that illustrate the efficacy of the proposed techniques, examine the bias-variance trade-off encountered when using inadequately blended data, and show that the gain in precision provided by the convenience sample is lower in circumstances where the outcome is strongly related to the auxiliary variables used for blending.

Tables

Table 1. Table 1: The number of respondents by domain and panel.

Domain	KP	WWP
Veterans	251	–
Post-9/11 Caregiver	72	281
Pre-9/11 Caregiver	522	3
Civilian Caregiver	1828	–
Non-Caregiver	1163	–

Table 2. Table 2: Means and standard errors of selected survey variables for post-9/11 caregivers calculated using various sources of data.

	KP Only		WWP Only		Blended: Unweighted
Variable Description	Mean	s.e.	Mean	s.e.	Mean	s.e.
Caregiver lives with veteran	0.453	0.085	0.858	0.021	0.755	0.032
Veteran deployed to war zone	0.577	0.086	0.929	0.015	0.839	0.029
Vet. has a service-related disability rating	0.563	0.095	0.883	0.016	0.801	0.031
Vet.’s disability rating is 70+%	0.280	0.086	0.715	0.023	0.605	0.033
Vet. has serv.-rel. traumatic brain injury	0.164	0.069	0.687	0.027	0.554	0.032
Vet. has serv.-rel. mental health problems	0.450	0.097	0.897	0.021	0.783	0.033
Caregiver is female	0.476	0.082	0.940	0.014	0.822	0.032
Caregiver depression level†	7.071	0.808	9.485	0.446	8.871	0.356
Caregiver anxiety level†	38.409	5.221	50.896	1.727	47.722	3.887
Caregiving caused quitting work entirely	0.204	0.065	0.502	0.028	0.426	0.031
Caregiving disturbs sleep	0.474	0.090	0.808	0.028	0.723	0.030
Caregiving causes financial strain	0.588	0.085	0.762	0.021	0.718	0.028
Caregiver feels overwhelmed by caregiving	0.398	0.085	0.815	0.024	0.709	0.032

Table 3. Table 3: Benchmarks and weighted averages among post-9/11 caregivers for selected auxiliary variables ( 𝐱 i subscript 𝐱 𝑖 {\bf x}_{i} ). Three schemes for calculating blending weights are shown: simultaneous propensity scores (SPS), disjoint propensity scores (DPS), and simultaneous calibration (SC).

	Bench-	Blended: Weighted
Variable Description	marks	SPS	DPS	SC
Caregiver lives with veteran	0.453	0.467	0.573	0.489
Veteran deployed to war zone‡	0.559	0.592	0.631	0.591
Veteran has a service-related disability rating	0.563	0.577	0.648	0.591
Veteran’s disability rating is 70+%	0.280	0.275	0.320	0.308
Vet. has service-related traumatic brain injury	0.164	0.133	0.166	0.189
Vet. has service-related mental health problems‡	0.510	0.491	0.529	0.549
Caregiver is female‡	0.555	0.509	0.557	0.592

Table 4. Table 4: Means and standard errors of selected outcome variables for post-9/11 caregivers from blended data with various weighting schemes. The p 𝑝 p -values are derived from a test of the hypotheses in ( 11 ). The weighting schemes shown are simultaneous propensity scores (SPS), disjoint propensity scores (DPS), and simultaneous calibration (SC).

	Blended: SPS			Blended: DPS			Blended: SC
Variable Description	Mean	s.e.	$p$-val.	Mean	s.e.	$p$-val.	Mean	s.e.	$p$-val.
Caregiver depression level†	7.661	0.575	0.109	7.574	0.895	0.513	7.801	0.668	0.112
Caregiver anxiety level†	41.551	3.887	0.084	43.657	4.827	0.344	41.706	3.933	0.340
Caregiving caused quitting work	0.206	0.072	0.009	0.237	0.066	0.489	0.237	0.069	0.114
Caregiving disturbs sleep	0.499	0.069	0.005	0.536	0.089	0.298	0.554	0.068	0.061
Caregiving causes financial strain	0.593	0.064	0.174	0.627	0.089	0.557	0.575	0.067	0.327
Feels overwhelmed by caregiving	0.443	0.104	0.007	0.480	0.078	0.291	0.480	0.087	0.010

Table 5. Table 5: Results for the models in ( 13 ) for various weighting schemes. See the caption to Table 2 for a description of the acronyms.

		(Intercept)	ERA	AGE	SEX	INC	EMP
KP Only	Estimate	11.4897	1.5102	-0.7542	-0.7812	-1.4843	-0.6085
	Std. Error	1.8073	0.9375	0.3330	0.6134	0.4848	0.7070
	$p$-value	0.0000	0.1078	0.0239	0.2034	0.0023	0.3898
Blended: Unweighted	Estimate	12.2617	2.3185	-0.8561	-1.6111	-1.6121	-0.1945
	Std. Errors	1.3262	0.5632	0.2727	0.4990	0.3540	0.5224
	$p$-value	0.0000	0.0000	0.0017	0.0013	0.0000	0.7098
Blended: SPS	Estimate	11.1934	1.6835	-0.7676	-1.1130	-1.3655	-0.2015
	Std. Error	2.6438	0.9258	0.4008	0.7327	0.8287	1.1329
	$p$-value	0.0002	0.0778	0.0639	0.1380	0.1086	0.8599
Blended: DPS	Estimate	11.5365	1.9340	-0.7633	-1.2863	-1.4557	-0.3199
	Std. Error	2.6331	0.6637	0.3003	0.4554	0.6961	1.1046
	$p$-value	0.0001	0.0063	0.0158	0.0079	0.0440	0.7739
Blended: SC	Estimate	11.6872	2.1431	-0.8125	-1.3930	-1.4348	-0.3934
	Std. Error	1.8824	0.7785	0.3195	0.5647	0.7694	0.6950
	$p$-value	0.0000	0.0094	0.0157	0.0188	0.0709	0.5751

Equations

$$d_i=P(S_1|{\mathbf{x}}_i),q_i=P(S_2|{\mathbf{x}}_i),\text{mbox}and{p}_i=P(S|{\mathbf{x}}_i),$$

$$d_i^{\ast }=P(S_1^{\ast }|{\mathbf{x}}_i),$$

$$r_i=P(S_1|S_1^{\ast },{\mathbf{x}}_i),$$

$${\gamma }_i=P(S_2|S_1\cup S_2,{\mathbf{x}}_i).$$

$$P(S_2|{\mathbf{x}}_i,y_i)=P(S_2|{\mathbf{x}}_i).$$

$$f(y_i|{\mathbf{x}}_i,S_2)=f(y_i|{\mathbf{x}}_i),$$

$${\gamma }_i=\frac{P(S_2|{\mathbf{x}}_i)}{P(S_1\cup S_2|{\mathbf{x}}_i)}=\frac{P(S_2|{\mathbf{x}}_i)}{P(S_1|{\mathbf{x}}_i)+P(S_2|{\mathbf{x}}_i)}=\frac{q_i}{d_i+q_i},$$

$$q_i=\frac{d_i{\gamma }_i}{1-{\gamma }_i}=\frac{r_i{d}_i^{\ast }{\gamma }_i}{1-{\gamma }_i}.$$

$$\text{mbox}E\left[\frac{{\sum }_{i\in S_1}{d}_i^{-1}{y}_i}{{\sum }_{i\in S_1}{d}_i^{-1}}\right]\approx \text{mbox}E\left[\frac{{\sum }_{i\in S_2}{q}_i^{-1}{y}_i}{{\sum }_{i\in S_2}{q}_i^{-1}}\right]\approx \frac{{\sum }_{i\in \Omega }{y}_i}{N},$$

$$\kappa =\frac{({\sum }_{i\in S_1}{d}_i^{-1})({\sum }_{i\in S_2}{q}_i^{-2})}{({\sum }_{i\in S_1}{d}_i^{-1})({\sum }_{i\in S_2}{q}_i^{-2})+({\sum }_{i\in S_1}{d}_i^{-2})({\sum }_{i\in S_2}{q}_i^{-1})}.$$

$$\bar{\theta }=\bar{\kappa }{\hat{\theta }}_1+(1-\bar{\kappa }){\hat{\theta }}_2$$

$$\bar{\kappa }=\frac{\text{mbox}Var({\hat{\theta }}_2)-\text{mbox}Cov({\hat{\theta }}_1,{\hat{\theta }}_2)}{\text{mbox}Var({\hat{\theta }}_1)+\text{mbox}Var({\hat{\theta }}_2)-2\text{mbox}Cov({\hat{\theta }}_1,{\hat{\theta }}_2)}.$$

$$p_i=\frac{d_i}{1-{\gamma }_i}=\frac{r_i{d}_i^{\ast }}{1-{\gamma }_i}$$

$$\text{mbox}E\left[\frac{{\sum }_{i\in S_1}{w}_i{y}_i}{{\sum }_{i\in S_1}{w}_i}\right]\ne \frac{{\sum }_{i\in \Omega }{y}_i}{N}$$

$${\mathbf{t}}_x=\sum _{i\in S}{v}_i{\mathbf{x}}_i.$$

$${\mathbf{t}}_x=\sum _{i\in S_1}{v}_{1i}^{\ast }{\mathbf{x}}_i=\sum _{i\in S_2}{v}_{2i}^{\ast }{\mathbf{x}}_i.$$

$$\kappa =\frac{({\sum }_{i\in S_1}{v}_{1i}^{\ast })({\sum }_{i\in S_2}(v_{2i}^{\ast }{)}^2)}{({\sum }_{i\in S_1}(v_{1i}^{\ast }{)}^2)({\sum }_{i\in S_2}{v}_{2i}^{\ast })+({\sum }_{i\in S_1}{v}_{1i}^{\ast })({\sum }_{i\in S_2}(v_{2i}^{\ast }{)}^2)}$$

$$y_i=\mu +\delta 1_{\{i\in S_2\}}+{ϵ}_i,$$

$${\mathcal{H}}_0:\delta =0\text{mbox}against{\mathcal{H}}_1:\delta \ne 0.$$

$$\log \left(\frac{{\gamma }_i}{1-{\gamma }_i}\right)={\zeta }_0+{\mathbf{\zeta }}^{\prime }{\mathbf{x}}_i,$$

$$\text{mbox}{DEP}_i={\eta }_0+{\eta }_1\text{mbox}{ERA}_i+{\eta }_2\text{mbox}{AGE}_i+{\eta }_3\text{mbox}{SEX}_i+{\eta }_4\text{mbox}{INC}_i+{\eta }_5\text{mbox}{EMP}_i+{ϵ}_i.$$

$$\log \left(\frac{{\rho }_i}{1-{\rho }_i}\right)=b_0+{\mathbf{b}}_1^{\prime }{\mathbf{v}}_i,$$

$$\text{mbox}{bias}_j=\frac{1}{K}\sum _{k=1}^K\text{mbox}{e}_j^{[k]}\text{mbox}and\text{mbox}{rMSE}_j=\sqrt{\frac{1}{K}\sum _{k=1}^K\text{mbox}{SRE}_j^{[k]}}.$$

$${\hat{p}}_j=\frac{1}{K}\sum _{k=1}^K{1}_{\{p_j^{[k]}\le \alpha \}}$$

Full text

Keywords: Inverse probability weighting, Calibration, Propensity scores.

1 Introduction

Probability (or representative) sampling, in which each member of a population is sampled with a known and nonzero probability, has long served as the gold-standard for producing inferences that are generalizable to the whole population. However, researchers are increasingly interested in deriving estimates on small segments of a population, and it is often infeasible to use probability sampling to produce samples of sufficient size from such segments due to the lack of prevalence of individuals therein.

One small segment of the population for which there is great interest is service members and veterans who served in US military operations in Iraq and Afghanistan since 2001. Policymakers are increasingly interested in identifying the health, educational, labor, and other social needs of this population so that they can form policies or create programs to better address these needs. In 2014, there were 2.6 million post-9/11 veterans, representing 0.8% of the total US population. Although the number of veterans is expected to increase to 3.6 million in 2019, it is still projected to only comprise 1.0% of the population (NCVAS, 2015). Oftentimes, researchers are interested in even further narrowly defined segments of veterans. Specifically, this article is motivated by efforts to draw inferences regarding the wellbeing of unpaid caregivers of wounded, ill, and injured service members veterans of the US armed forces who served after September 11, 2001 (post-9/11) and caregivers of those who did not (pre-9/11). Sampling of veterans of the Armed Forces is aided by databases that are compiled by organizations that provide services to such individuals (e.g., Iraq and Afghanistan Veterans of America, Disabled Veterans of America, American Legion). Utilization of these services is voluntary; therefore, these databases do not yield sampling frames that are representative of the veteran population, and therefore samples drawn using these lists are considered non-probability samples.

Non-probability (or convenience) samples are becoming an increasingly efficient and cost-effective alternative to probability samples. Baker et al. (2013) provide an outline of types of convenience samples (e.g., mall intercepts, volunteer panels, river samples and respondent driven samples). Of particular interest in the current literature are volunteer (opt-in) web panels (e.g., Baker et al., 2010; Biffignandi and Bethlehem, 2012; Tourangeau et al., 2013). Convenience samples are usually not representative of the complete population from which they were drawn, and inferences extracted from them do not necessarily generalize to that population (i.e., such samples may be fraught with selection bias). Therefore, a wealth of literature (e.g., Schonlau et al., 2004; Duffy et al., 2005; Schonlau et al., 2007; Chang and Krosnick, 2009; Blasius and Brandt, 2010; Yeager et al., 2011) has been devoted to evaluating the quality of inferences drawn from convenience samples.

It is possible to use non-representative samples to produce estimates that generalize to a population even in circumstances where micro-level data that are representative of that population are not available. For instance, if the non-representative sample is differentiated from the population solely across characteristics for which aggregate measures are available for the entire population, a host of statistical weighting techniques, wherein certain observations are provided greater emphasis than others in estimator calculations, may be applicable (see Wang et al., 2015, for example). However, it is often the case that no comprehensive summary quantities are available from the entire segment of interest or that the population and convenience sample are differentiated on the basis of subtle characteristics that are not commonly reported in census-type summaries. These disadvantages can be overcome by using a (potentially small) probability sample in tandem with non-representative data.

The idea of combining data from two or more surveys has a lengthy history within extant literature (e.g., Zieschang, 1986; Renssen and Nieuwenbroek, 1997; Merkouris, 2004; Rao and Wu, 2010). This discussion has spurned the development of methods for combining (or blending) data from probability and convenience samples so as to yield generalizable inferences. Elliott and Haviland (2007) and Ghosh-Dastidar et al. (2009) discuss composite estimators—i.e., estimation involving outcomes is performed on each sample individually and the results are combined in some manner. Here, we focus on weighting techniques that enable the two datasets to be analyzed as a single dataset (following weighting). We separate these methods into two general classes of procedures: 1) propensity score methods (Lee, 2006; Lee and Valliant, 2009; Schonlau et al., 2009; Valliant and Dever, 2011), which model the mechanism that differentiates the two samples, and 2) calibration weighting (Lee and Valliant, 2009; DiSogra et al., 2011), which determines weights that satisfy the constraint that totals calculated across a set of variables using a weighted sample equal corresponding totals calculated using the entire population. Each of these weighting procedures is designed to match the convenience sample to the population across a set of auxiliary variables; congeniality between the sample and the remainder of the population across other variables (e.g., outcomes) is assumed to hold following weighting adjustments.

Our first methodological contribution is to introduce a novel and easily applicable manner of using propensity scores to calculate weights by directly approximating inclusion probabilities. We rigorize the framework for combining samples by developing and distinguishing between methods that weight the samples separately (wherein each sample is weighted to be individually representative of the population) and simultaneously (in which the samples are only representative when used in tandem), and we illustrate the comparative advantages of these two types of weights. Specifically, weights based on disjoint blending enable evaluation of the representativeness of the blended sample through statistical testing of the adequacy of blending (i.e., are there latent characteristics that differentiate the two samples?). In contrast, weights based on simultaneous blending appear to yield smaller variances than their disjoint counterparts and do not require that the convenience sample covers the population.

For the application discussed herein, we utilize data from a survey administered as part of an extensive study of military caregivers; the survey and a wealth of findings from the study are described in detail in Ramchand et al. (2014). The caregivers survey was administered to a probability-based panel of approximately 45,000 American adults (GfK KnowledgePanel); however, this yielded only 72 post-9/11 caregivers. A convenience sample that provided an additional 281 post-9/11 military caregivers was derived from a database provided by the Wounded Warrior Project (WWP). Response patterns from caregivers sampled from the WWP are markedly different from the those given by respondents in the probability sample. Nonetheless, the proposed methods for weighting are shown to be effective at combining data from the two sources. Using only the probability sample, we estimate that post-9/11 caregivers score 1.5 points higher average than pre-9/11 caregivers on a 24-point measure of depression levels. However, a 95% confidence interval on this estimate ranges from -0.36 to 3.38 points. Using the blending methods outlined here, we see a marked improvement in precision of this estimate (see Section 3 for details).

The weighting methodology (and techniques for variance estimation) is detailed in Section 2. In Section 3, we relay results from the application of these methods to the survey of military caregivers. Section 4 presents simulation studies performed using data from the caregiver survey and purely synthetic data, and Section 5 provides concluding points of discussion. Supplementary materials consider extensions to settings involving multiple convenience samples, provide details not given in the main article, and tie together loose ends.

2 Weighting Schemes

Prior to proposing weighting methods for blending probability and convenience samples, we introduce some fundamental concepts and notation that will be used throughout. To simplify the exposition, we assume unstratified, unclustered designs (however, stratification and clustering are discussed in Section 5). We let $\Omega$ denote the set of units (individuals) within a population, whereas $S_{1}\subset\Omega$ and $S_{2}\subset\Omega$ denote the set of units within the probability and convenience samples, respectively. (In the supplemental materials, we extend theory outlined here to allow for multiple convenience samples.) Under the assumption that the two samples are disjoint (i.e., $S_{1}\cap S_{2}=\emptyset$), we notate the set of respondents within the combined (or blended) sample as $S=S_{1}\cup S_{2}$. We assign a weight $w_{i}$ to each unit $i\in S$. The problem of interest involves inferences for a variable $\{y_{i}\}_{i\in S}$, whereas a set of auxiliary variables (that does not include $y_{i}$), which are denoted by the vector ${\bf x}_{i}$ for unit $i$, are used to capture discrepancies between the two samples. Although the theory presented here is applicable for multivariate outcomes, we assume univariate $y_{i}$ to simplify the exposition. The notation $\{\cdot\}_{i\in A}$ indicates a set of data values over a subset $A\subseteq\Omega$; when referencing all sampled units, the subscript $i\in S$ is excluded. Both $y_{i}$ and ${\bf x}_{i}$ are observed for all $i\in S$. For the remaining units of the population, neither case-level data values nor population totals are assumed to be known.

We also assume that $S_{1}$ has been subjected to an estimable non-response mechanism. That is, $S_{1}\subseteq S_{1}^{*}$, where $S_{1}^{*}$ is the set of individuals selected for participation within probability sampling. Because the analyst has little control over the composition of the convenience sample, we do not consider non-response for $S_{2}$ (i.e., we embed non-response within the mechanism that generates $S_{2}$).

Note that in lieu of variances, we often discuss design effects. The design effect, as proposed by Kish (1965), quantifies the inflation in variance that is a consequence of a weighting scheme. Specifically, it represents the ratio of the variance of a weighted estimator with the variance of an equivalent unweighted estimator (i.e., under a simple random sampling design).

2.1 Weights Based on Propensity Scores

If the probability of inclusion into a sample is known for each sampled unit and positive for all units in the population, one can use the reciprocal of such probabilities as sampling weights; these weights may be used to provide unbiased estimators of population quantities within sampling schemes that are not simple random sampling (Kish, 1965). Generally, $P(i\in A)$ denotes the probability that individual $i$ is in set $A$; this expression is written equivalently as $P(A)$ for simplicity. Our first set of procedures are based on the quantities,

[TABLE]

which represent the probability of inclusion in the probability, convenience, and blended samples, respectively, for individual $i$ given the individual’s set of auxiliary characteristics ${\bf x}_{i}$.

In the presence of a non-response mechanism, $\{d_{i}\}$ is derived as follows. The practitioner assigns a sampling probability, denoted

[TABLE]

to all $i\in\Omega$. Next, $S_{1}$ is generated by the response mechanism

[TABLE]

and as a result $d_{i}=d_{i}^{*}r_{i}$. For convenience, we write all probabilities listed above as being conditioned upon the same set of auxiliary variables ${\bf x}_{i}$. Of course, this is not required.

Since $\{q_{i}\}$ and $\{p_{i}\}$ cannot be directly determined, and to help approximate them, we use estimated values of the probability that a sampled unit was included into the convenience sample, which is given by

[TABLE]

The quantity $\gamma_{i}$ is analogous to a propensity score as commonly used in treatment/control studies (Rosenbaum and Rubin, 1983) and to differentiate nonrespondents from respondents (Little and Rubin, 2002; Bethlehem et al., 2011)—therefore, weights that incorporate this quantity are referred to as propensity score-based weights. Our approaches based on these weights mandate the following assumptions:

We assume the selection probabilities for the probability sample are based only on a set of design variables, ${\bf x}_{i}^{*}$, so that $d^{*}_{i}=P(S_{1}^{*}|{\bf x}_{i}^{*})$ where ${\bf x}_{i}^{*}\subseteq{\bf x}_{i}$ and ${\bf x}_{i}^{*}$ is observed for all $i\in\Omega$. Also, the response mechanism is dependent only on a set of variables, $\widetilde{\bf x}_{i}$, so that $r_{i}=P(S_{1}|S_{1}^{*},\widetilde{\bf x}_{i})$ where $\widetilde{\bf x}_{i}\subseteq{\bf x}_{i}$ and ${\bf x}_{i}^{*}$ is observed for all $i\in S^{*}_{1}$.

Assumption 1**.**

Design of $S_{1}^{*}$: The selection probabilities for the probability sample depend only on a set of design variables ${\bf x}_{i}^{*}$, so that $d^{*}_{i}=P(S_{1}^{*}|{\bf x}_{i})=P(S_{1}^{*}|{\bf x}_{i}^{*})$ where ${\bf x}_{i}^{*}\subseteq{\bf x}_{i}$ and ${\bf x}_{i}^{*}$ is observed for all $i\in S_{1}^{*}\cup S_{2}$ with $d_{i}^{*}>0$ for all $i\in\Omega$. Furthermore, $d_{i}^{*}$ is known for all $i\in S_{1}^{*}\cup S_{2}$.

Assumption 2**.**

MAR of $S_{1}^{*}$: The response mechanism that yields $S_{1}$ is dependent only on a set of variables $\widetilde{\bf x}_{i}$, where $\widetilde{\bf x}_{i}\subseteq{\bf x}_{i}$ and $\widetilde{\bf x}_{i}$ is observed for all $i\in S^{*}_{1}\cup S_{2}$. That is, $r_{i}=P(S_{1}|S_{1}^{*},{\bf x}_{i})=P(S_{1}|S_{1}^{*},\widetilde{\bf x}_{i})$. Furthermore, $r_{i}$ is expressed as an estimable function of $\widetilde{\bf x}_{i}$ with $r_{i}>0$ for each $i\in\Omega$.

It follows from Assumptions 1 and 2 that $d_{i}>0$ for all $i\in\Omega$, and consequentially the probability sample covers the population.

Assumption 2 establishes that non-response in $S_{1}^{*}$ is missing at random, in the nomenclature of Little and Rubin (2002), with respect to $\widetilde{\bf x}_{i}$. That is, $P(S_{1}|S_{1}^{*},{\bf x}_{i},z_{i})=P(S_{1}|S_{1}^{*},{\bf x}_{i})$, where $z_{i}$ is exogenous to $\widetilde{\bf x}_{i}$. This observation and Assumption 1 establish that the probability sample is conditionally independent the outcome of interest in that $P(S_{1}|{\bf x}_{i},y_{i})=P(S_{1}|{\bf x}_{i})$. We next impose a similar expression for the convenience sample.

Assumption 3**.**

Ignorability:* It holds that*

[TABLE]

Letting $f(\>\cdot\>|\>\cdot\>)$ denote a conditional density, (2) can be equivalently stated as

[TABLE]

which is the exchangeability (or ignorability) assumption that underpins traditional propensity score techniques (Rosenbaum and Rubin, 1983).

Note that the theory that forms the methodology presented here is valid in the event that the outcome $y_{i}$ is included in the set of auxiliary variables used for blending. However, as illustrated in Section 4.2, it is not beneficial to do so (specifically, no precision is gained by utilization of the convenience sample in such circumstances). To ensure validity when $y_{i}\notin{\bf x}_{i}$, we impose Assumption 3. To assess whether this assumption holds in practice, we suggest the test for the adequacy of blending as described in Section 2.3.

Assumption 4**.**

Model Adequacy:* The models used to estimate $r_{i}$ and $\{\gamma_{i}\}$ are each correctly specified.*

It is common practice in missing data analysis and in propensity score approaches to assume that the models used (i.e., non-response and propensity score models) are not misspecified—certainly, our procedures require the same. In that vein, if the mechanism that generates the convenience sample is thought to be independent of some elements of ${\bf x}_{i}^{*}$ or $\widetilde{\bf x}_{i}$, the model for $\gamma_{i}$ should be adjusted accordingly.

Assumption 5**.**

Positivity:* For all $i\in\Omega$, $q_{i}>0$.*

Traditional propensity score methods require a positivity assumption (i.e., the propensity scores are non-zero). Here, positivity is implied by assuming that the convenience sample covers the population. However, as is clarified shortly, this assumption is not needed for all methods we present.

Robustness to these assumptions is assessed in Section 4. Given the assumptions, we present two options for weighting via inclusion probabilities based on propensity scores; the first weights the samples separately, whereas the second jointly weights the samples to obtain representativeness.

2.1.1 Disjoint Weighting

Our first propensity score weighting scheme involves estimation of $\{q_{i}\}$, the probability of inclusion into the convenience sample. Note that

[TABLE]

where we have made use of the fact that $S_{1}\cap S_{2}=\emptyset$. Solving the above expression for $q_{i}$ gives

[TABLE]

Note that (3) does not indicate that $q_{i}$ is dependent upon $d_{i}$ (i.e., any alterations in $d_{i}$ while leaving $q_{i}$ unchanged will be offset by changes in the value of $\gamma_{i}$).

It follows from Assumptions 1 and 2 that if weights calculated using $d_{i}^{-1}$ are applied to the probability sample, the probability sample will be representative of the population. Likewise, Assumption 5 implies that weights calculated using $q_{i}^{-1}$, when applied to the convenience sample, will ensure that the convenience sample is representative. That is,

[TABLE]

Note that the first and second expectations in (4) are taken over the random variables used to generate $S_{1}$ and $S_{2}$, respectively. The above result follows from the unbiasedness of the Horvitz-Thompson estimator and illustrates that representativeness for both panels can be obtained without combining the two samples. Nonetheless, fusion of the data sources is preferred since blending is expected to yield estimators with lower variance than those which are constructed using the samples individually.

We next propose so-called disjoint propensity score weights—the phrase “disjoint” blending is used since such procedures involve combination of two samples that are individually representative. Analysis of multiple surveys using disjoint combing methods has been considered previously (e.g., Kott and Vogel, 1995; Renssen and Nieuwenbroek, 1997; Merkouris, 2004; Lee, 2006). Although the two samples are essentially weighted separately in disjoint blending, the weights are concatenated to yield a single set of weights that enables the two datasets to be analyzed in tandem. Note that this mandates that the analyst decide how much emphasis to give one sample over the other in the concatenated weights. Therefore, our first set of blending weights (notated $\{w^{*}_{i}\}$ for $i\in S$), is calculated by setting $w^{*}_{i}=\kappa d_{i}^{-1}$ for $i\in S_{1}$ and $w^{*}_{i}=(1-\kappa)q_{i}^{-1}$ for $i\in S_{2}$ for some constant $\kappa\in[0,1]$.

As illustrated in the supplemental materials to this article, $\{w^{*}_{i}\}$ may be used to obtain unbiased estimators of population quantities for any $\kappa\in[0,1]$.

In order to minimize the Kish approximation of the design effect (i.e., $\mbox{deff}=n\sum w_{i}^{2}/(\sum w_{i})^{2}$), we suggest

[TABLE]

We do not recommend selection of $\kappa$ so as to minimize the true design effect (or variance for that matter) since doing so would require a different value of $\kappa$ for each estimate produced. However, such approaches have been considered previously (Hartley, 1974; Rivers et al., 2003; Lohr, 2011). Those authors suggest post hoc blending of estimators taken from multiple frames. Therein, the two samples are combined after analyses involving outcomes begins. (In our proposed methods, blending occurs prior to such analysis.) Consider the estimators $\hat{\theta}_{1}$ and $\hat{\theta}_{2}$ as described in the previous paragraph. The post hoc blended estimator of $\theta$ is

[TABLE]

where $\bar{\kappa}$ is selected to minimize the mean squared error of $\bar{\theta}$. That is, set

[TABLE]

This calculation incorporates the assumed unbiasedness of $\hat{\theta}_{1}$ and $\hat{\theta}_{2}$ but does not assume $\theta_{1}$ and $\theta_{2}$ are independent (as $\theta_{2}$ is calculated using weights that invoke $S_{1}$).

Blending via disjoint weights as proposed with the use of the Kish approximation should yield a nearly minimum MSE estimator, and it does not require analyzation of the samples separately, nor does it require a different choice of $\kappa$ for each estimate.

2.1.2 Simultaneous Weighting

Weights may also be calculated using the probability of inclusion into the blended sample. Note that $p_{i}=P(S_{1}\cup S_{2}|{\bf x}_{i})=P(S_{1}|{\bf x}_{i})+P(S_{2}|{\bf x}_{i})=d_{i}+q_{i}$. Hence,

[TABLE]

for all $i\in S$. Therefore, our second proposed choice of weights for blending probability and convenience samples is determined by setting $w_{i}=p_{i}^{-1}$ for all $i\in S$. We refer $\{w_{i}\}$ as simultaneous propensity score weights. Note that Assumption 5 is not required for simultaneous weighting. That is, if $q_{i}=0$, it follows that $\gamma_{i}=0$ and thus $p_{i}=d_{i}>0$ (and $S$ covers $\Omega$).

Simultaneous blending, which has been described previously in Lee and Valliant (2009) in the context of calibration weighting, is designed to make the two samples individually representative of the population when used in tandem. However, when the samples are weighted in this manner, they are not individually representative. That is, in general,

[TABLE]

where $N$ is the cardinality of $\Omega$ (i.e., the population size). Equivalence does not approximately hold in general in this expression. Therefore, it is not prudent to use $\{w_{i}\}$ in conjunction with the observed data for the purpose of testing for the presence of discrepancies between the probability and convenience samples that are unexplained (which would determine whether or not $\{{\bf x}_{i}\}$ is sufficient for modeling inclusion probabilities).

Note that both disjoint and simultaneous propensity score methods require a known or estimated value of $d_{i}$ for each $i\in S_{1}\cup S_{2}$. In practice, these probabilities are easily determined. From Assumption 1, $d_{i}^{*}$ is known for all $i\in\Omega$. (E.g., if $n$ individuals from a stratum of size $N$ are randomly selected asked to participate in the study within probability sampling, one may use $d_{i}^{*}=n/N$ for all $i$ in the stratum’s population.) It is also assumed that the response mechanism, $r_{i}$, that yields $S_{1}$ is estimated as a function of predictor variables—for example, $\hat{r}_{i}=1/(1+\exp\{-\hat{\alpha}_{0}-\widehat{\bm{\alpha}}^{\prime}\widetilde{\bf x}_{i}\})$. Since $\hat{r}_{i}$ can be calculated for all $i\in S$, we approximate $d_{i}$ for $i\in S$ using $d_{i}^{*}\hat{r}_{i}$.

2.2 Calibration Weighting

If population totals for a set of auxiliary variables ${\bf x}_{i}$, that is ${\bf t}_{x}=\sum_{i\in\Omega}{\bf x}_{i}$, are known, calibration weighting (Deville and Särndal, 1992; Särndal, 2007) is useful for deriving weights that can be used to calculate generalizable estimators from non-representative data. As such, it has a lengthy history as a tool for non-response adjustments (e.g. Kott, 2006) and has a similar utility for blending probability and convenience samples. Further, it should give results that are similar to those found using propensity score weighting (albeit through different computational procedures).

Calibration involves the determination of weights $\{v_{i}\}$ that satisfy

[TABLE]

Under a linear representation of $y_{i}$ as a function of ${\bf x}_{i}$ (e.g., $\mbox{E}[y_{i}]=\bm{\beta}^{\prime}{\bf x}_{i}$), it follows that $\hat{t}_{y}=\sum_{i\in S}v_{i}{y}_{i}$ is a unbiased estimator of $t_{y}=\sum_{i\in\Omega}{y}_{i}$. Similarly, unbiased estimates of $\mbox{E}[y_{i}]$ can be calculated using the weighted data.

Similar to the framework we established for weighting via propensity scores, we let weights $\{v_{t}\}$, which satisfy (8) for $S=S_{1}\cup S_{2}$, denote so-called simultaneous calibration weights; these weights jointly blend the probability and convenience samples to produce a blended sample that is representative of the population. However, it is also possible to individually weight the two samples so that each one is separately representative of the population. That is, we may use calibration to find weights $\{v^{*}_{1i}\}_{i\in S_{1}}$ and $\{v^{*}_{2i}\}_{i\in S_{2}}$ that satisfy

[TABLE]

If we define weights $\{v^{*}_{i}\}_{i\in S}$ via $v^{*}_{i}=\kappa v^{*}_{1i}$ for $i\in S_{1}$ and $v^{*}_{i}=(1-\kappa)v^{*}_{2i}$ for $i\in S_{2}$ for some $\kappa\in[0,1]$, then $\{v^{*}_{i}\}$ is conducive for calculation of unbiased estimates of $\mbox{E}[y_{i}]$ and $\mbox{E}[y_{i}|{\bf x}_{i}]$. We refer to $\{v^{*}_{i}\}$ as disjoint calibration weights. Similar to the disjoint propensity score-based weights, we set

[TABLE]

in order to minimize the Kish approximation of the design effect. As is the case with $\{w^{*}_{i}\}$ (the disjoint propensity score weights), $\{v^{*}_{i}\}$ can be used to test for post-weighting discrepancies between the two samples discussed in Section 2.3 below; the simultaneous calibration weights $\{v_{i}\}$ should not be used for this purpose. Unlike the weights based on selection probabilities, calibrated weights may adjust for non-response within the probability sample (if the non-response mechanism is explained by the auxiliary variables ${\bf x}_{i}$). Further, disjoint calibration weights can be used with post hoc blending estimation in the vein of (5).

Given initial design weights $\{\omega_{i}\}$, the calibration weighting procedure, denoted generalized raking by Deville and Särndal (1992), calculates the set of calibrated weights $\{v_{t}\}$ which minimize the distance between the initial and calibrated weights subject to a distance metric $G(\cdot)$ and the constraints imposed by (8). That is, the quantity $\sum_{i\in S}\omega_{i}G(v_{i}/\omega_{i})$ is minimized subject to (8). Calculations show that $\{\omega_{i}\}$ satisfies the expression $\sum_{i\in S}\omega_{i}F({\bf x}_{i}^{\prime}\bm{\xi}){\bf x}_{i}={\bf t}_{x}$, where $F(\cdot)$ is the inverse of the derivative of $G(\cdot)$ and where $\bm{\xi}$ is a vector of Lagrange multipliers (which may be extracted via Newton’s method). See Deville et al. (1993) for details. We use a truncated linear distance metric bounded below at zero (i.e., $G(x)=(1/2)(x-1)^{2}$ for $x\geq 0$ and $G(x)=\infty$ otherwise), which has been shown to help minimize the design effect (Robbins et al., 2017) while ensuring that $w_{i}>0$ for all $i\in S$. The propensity score-based weights $\{w_{i}\}$ or $\{w^{*}_{i}\}$ can be used as initial values within calculation of their respective calibration weights $\{v_{i}\}$ or $\{v^{*}_{i}\}$. However, it is also feasible to initialize the algorithm by giving each respondent equal weight.

Assumptions 1-3 are necessary when using calibration. Note that unlike propensity scores, calibration may be applied without knowing $d_{i}$ for $i\in S_{2}$. Since we do not use a propensity score model, Assumption 4 is not needed. Analogously, validity of calibration estimators is frequently illustrated on outcomes that obey $y_{i}=\beta_{0}+\bm{\beta}^{\prime}{\bf x}_{i}+\epsilon_{i}$ (see Deville and Särndal, 1992, for example). However, it commonplace to apply calibration estimation to outcomes that do not obey this formulation (e.g., binary data). Assumption 5 is needed for disjoint (but not simultaneous) calibration.

Calibration weighting for the purpose of blending probability and convenience samples has been proposed by DiSogra et al. (2011); therein, a technique akin to the simultaneous calibration method presented here is suggested. Lee and Valliant (2009) propose a method akin to the disjoint calibration procedure discussed herein.

2.2.1 Estimation of Benchmarks

The vector of population totals ${\bf t}_{x}$ is said to be composed of benchmark values of the auxiliary variables. It is unlikely that this vector will be known in practice, especially if the population of interest is a rare segment of a larger population. In the event that no elements of ${\bf t}_{x}$ are known, one can approximate it using $\hat{\bf t}_{x}=\sum_{i\in S_{1}}d_{i}^{-1}{\bf x}_{i}$. where $\{d_{i}^{-1}\}$ are the design weights. Calibration weights calculated using $\hat{\bf t}_{x}$ offer no perceived advantage over the propensity score-based weights. Nonetheless, calibration weighting can be used to incorporate limited aggregate information regarding the population of interest. As an example: If ${\bf x}_{i}$ can be decomposed into two vectors ${\bf x}_{i}^{(1)}$ and ${\bf x}_{i}^{(2)}$, satisfying ${\bf x}_{i}=(({\bf x}_{i}^{(1)})^{\prime},({\bf x}_{i}^{(2)})^{\prime})^{\prime}$, where ${\bf t}_{x}^{(1)}=\sum_{i\in\Omega}{\bf x}_{i}^{(1)}$ is known but ${\bf t}_{x}^{(2)}=\sum_{i\in\Omega}{\bf x}_{i}^{(2)}$ is not, one can apply calibration to find weights $\{v_{i}^{(1)}\}$ that obey $\sum_{i\in S_{1}}v_{i}^{(1)}{\bf x}_{i}^{(1)}={\bf t}^{(1)}_{x}$. Setting $\tilde{\bf t}_{x}=\sum_{i\in S_{1}}v_{i}^{(1)}{\bf x}_{i}$, the vector $\tilde{\bf t}_{x}$ yields benchmark values that can be used in calculation of blending weights $\{v_{i}\}$, $\{v_{1i}^{*}\}$ and $\{v_{2i}^{*}\}$. Calibration weighting in this manner may account for non-response among the probability sample (so long as the non-response mechanism is influenced by only variables within ${\bf x}_{i}^{(1)}$).

2.3 A Test for the Adequacy of Blending

We are interested in knowing whether or not the blended sample is representative of the population (following application of the weighting methods described above). If the probability of inclusion in the convenience sample has not been appropriately modeled or if an insufficient scope of auxiliary variables has been used within calibration weighting, the blended sample is not assured to be representative. Therefore, we develop a method for statistically testing the adequacy of blending.

Let ${\theta}$ denote a population parameter of interest. Further, let $\hat{\theta}_{1}$ and $\hat{\theta}_{2}$ denote estimates of $\theta$ found using only $i\in S_{1}$ and only $i\in S_{2}$, respectively. We will test to see whether or not $\hat{\theta}_{1}$ and $\hat{\theta}_{2}$ are statistically equivalent. Lack of equivalence may indicate the presence of a latent characteristic that exists outside the set of auxiliary variables $\{{\bf x}_{i}\}$ which differentiates $S_{1}$ and $S_{2}$ (in which case (2) is violated). If $\hat{\theta}_{1}$ and $\hat{\theta}_{2}$ are calculated using simultaneous blending weights, we expect a non-negligible discrepancy between $\hat{\theta}_{1}$ and $\hat{\theta}_{2}$ when blending is adequate (even though $\hat{\theta}$, an estimator of $\theta$ determined using data from both samples with simultaneous blending weights, is expected to be unbiased). To observe $\hat{\theta}_{1}\approx\hat{\theta}_{2}$ under an adequate selection of ${\bf x}_{i}$, disjoint blending weights (determined using propensity scores or calibration) are needed.

We propose a test for the adequacy of blending that is based off of the mean value of an outcome of interest $\{y_{i}\}$. Letting $\mu=\sum_{i\in\Omega}y_{i}/N$, we evaluate differences between $\hat{\mu}_{1}$ and $\hat{\mu}_{2}$, which are estimates of $\mu$ based off of the probability and convenience samples, respectively, and which are calculated using a disjoint weighting scheme. We approach this problem by fitting the model

[TABLE]

for each $i\in S$ using weighted least squares. We let $\{\epsilon_{i}\}$ denote mean zero errors and, generally, $1_{\{A\}}$ represents the indicator of event $A$. We test

[TABLE]

This test may be performed using the statistic $z^{*}=\hat{\delta}/\widehat{\mbox{Var}}(\hat{\delta})^{1/2}$.

We formulate the test using the linear model in (10) as opposed to a basic two-sample formulation so as to enable easy extension to analogous tests that involve non-linear response (e.g., logistic regression) or that involve regression models with additional covariates. One may also develop an omnibus version of the above which can incorporate multiple outcomes (although this is not considered further here).

2.4 Variance Estimation

We briefly discuss techniques for variance estimation in the presence of weights calculated via the methodology outlined herein. We first describe (Taylor series) linearization and conclude with a discussion of resampling methods (with a focus on the jackknife). Further details on both techniques are given in the supplementary materials to this article.

Linearization provides algebraic approximations of variances of estimators found using complex survey data. The procedure invokes Taylor series expansions to translate an estimator into terms for which the variance can be estimated using straightforward procedures designed for complex surveys (e.g., the Horvitz-Thompson estimator Horvitz and Thompson, 1952). We implement linearization through the survey package in R (Lumley, 2004, 2011).

By using resampling techniques, one can incorporate aspects of an estimation process into variance calculations that are not easily captured algebraically. We focus on the jackknife (Quenouille, 1949, 1956), wherein data are segmented into replication groups, and the estimation process is applied to each replication group separately. Here, we consider a delete-a-group jackknife (Kott, 2001), wherein samples are segmented into $G$ mutually exclusive and exhaustive groups, and the weights are recalculated for each group.

3 A Survey of Military Caregivers

The early 21st century wars in Afghanistan and Iraq (Operation Enduring Freedom, Operation Iraqi Freedom, and Operation New Dawn) have produced a wave of returning wounded veterans and servicemembers. Many of those who have engaged in these modern conflicts have suffered unique traumas (e.g., traumatic brain injury associated with improvised explosive devices); fortunately, they have also benefited from advances in battlefield medicine and rehabilitative services (thereby reducing death rates overall but resulting in an increased proportion of veterans who have endured disabling injuries). As a result, America’s modern-day wounded warriors are often in need of a variety of caregiving services (e.g., bathing, dressing, and eating, etc.) upon their return to the homefront. However, these caregiving burdens often fall upon a range of family members, friends and other acquaintances (i.e., military caregivers) who do not receive wages for their services. As recently as a decade following the beginning of the conflicts in Iraq and Afghanistan, little was known about the average American military caregiver (e.g., their prevalence across the American population, the specific obstacles they face, etc.); this lack of information acted as an impediment for policymakers aiming to establish services to assist caregivers.

Ramchand et al. (2014) report findings from a large study of military caregivers performed from August to October 2013 by the RAND Corporation and sponsored by the Elizabeth Dole Foundation—our empirical application here utilizes data from this study. Substantive definitions utilized within Ramchand et al. (2014) are reviewed here. Generally, a caregiver is someone (a family member, friend, neighbor, etc.) who provides a broad range of care for an individual with a disabling (physical or mental) wound, injury or illness. A key characteristic of caregivers in our context is that they do not receive wages for caregiving. A military caregiver provides these services to a current or former member of the U.S. Armed Services. Further, a military caregiver is denoted as post-9/11 if the care recipient served in the armed forces at any point following September 11, 2001 (regardless of prior military service).

A crucial aspect of the RAND survey of military caregivers is that it was probability-based. A key policy question investigated by the military caregivers study involved the comparison of caregivers of servicemembers who served during the Afghanistan and Iraq conflicts (i.e., post-9/11 caregivers) to caregivers of servicemembers who served in prior eras (i.e., pre-9/11 caregivers). It was understood that probability-based sampling would not yield a sufficient number of post-9/11 caregivers; therefore, convenience samples from this segment were also drawn. Here, we apply the methodology described in Section 2 in order to blend data from the probability-based sample of caregivers with data from a convenience sample of post-9/11 caregivers and thus enable inferences that are generalizable to this sub-population.

Ramchand et al. (2014) describe a multitude of analyses that were performed using the caregiver data. Primarily, marginal means for respondents within the domains described above are calculated for several outcomes of interest. A primary outcome discussed here is the depression score, which is gauged using the eight-question Patient Health Questionnaire (Kroenke et al., 2009). A related outcome measure discussed here is anxiety levels, which are measured using the Mental Health Inventory anxiety subscale. On these two scales, higher scores indicate higher levels of depression and anxiety.

3.1 The Probability Sample: KnowledgePanel (KP)

The GfK KnowledgePanel (KP; GfK, 2013), which is a nationally representative panel of approximately 45,000 American adults (during August 2013), was used as a probability sample for the military caregivers survey. Each member of the KnowledgePanel was given a screener and was consequentially placed into one of the following domains: 1) military veteran receiving unpaid care, 2) military caregiver, 3) civilian caregiver, and 4) other (i.e., non-veteran/non-caregiver). All respondents who fell into the first two domains listed were given a full survey. Caregivers were asked to report on themselves and their care recipient; likewise, veterans were asked to report on themselves and their caregiver. Random samples of members who fell into the final two domains were selected and asked to complete a full survey. Of the 41,163 panelists invited to complete the screener, 28,164 (68%) did so.

To facilitate blending of data from supplemental sources, military caregivers who completed the full survey are then further segmented on the basis of the criteria mentioned earlier. Table 1 lists the number of respondents from each domain that were sampled from the KP (as well as the respondents from the convenience samples, which are described later).

Pre-stratification weights (i.e., screener weights) for all screener respondents from the KP were calculated by using calibration with a set of standard demographic variables (e.g., age, gender, income, etc.) as auxiliary variables (see ${\bf x}_{i}$ in (8)) where the benchmark values (${\bf t}_{x}$) are determined using information (or the non-institutionalized U.S. adult population) from the 2010 Census and American Community Survey (see Ramchand et al., 2014, for futher details). These weights account for non-response at the screener stage. In the ensuing, military caregivers are segmented further by era of service (pre- vs. post-9/11), and we focus on the domain of post-9/11 caregivers.

3.2 The Convenience Sample: Wounded Warrior Project (WWP)

As indicated by Table 1, there are only 72 post-9/11 caregivers from the probability sample, which yields concern that efforts to compare this domain to other domains will be under-powered. To further supplement our data, we used a database from the Wounded Warrior Project (WWP). The WWP is a nonprofit veterans service organization that offers a variety of programs to support servicemembers who were wounded during the military actions following the events of September 11, 2001. The WWP maintains a database of individuals who have registered with the organization as caregivers of wounded, ill, or injured veterans of post-9/11 military actions—this database yielded a sampling frame for additional post-9/11 caregivers. The rate of response among individuals sampled from the WWP was 20%, which yielded 281 post-9/11 caregivers who completed the survey.

We compare the KP and WWP across a wide range of variables. Tables 2 and 4 show (among other information) means and standard errors of several variables from the caregiver survey for respondents from the KP and the WWP. We note that Table 2 includes auxiliary variables used in weighting, whereas Table 4 includes a selection of outcome variables (although this discrepancy is not yet germane to the discussion). Since the WWP contains primarily post-9/11 caregivers, the information in these tables is calculated using only data from post-9/11 caregivers. The table also lists unweighted means and standard errors of these variables when calculated using the KP and WWP jointly. Note that the standard errors presented in these tables (and in all subsequent analyses) are calculated using Taylor series linearization. As indicated by Tables 2, there are substantial differences across the KP and WWP panels. Specifically, the WWP caregivers are predominantly female (94%), whereas around half of the caregivers from the KP are female. Likewise, around 69% of the caregivers from the WWP are caring for a veteran with traumatic brain injury (TBI), whereas only 16% of care recipients from the KP have TBI. Furthermore, WWP caregivers report higher depression and anxiety levels. Note that all differences in means seen in Table 2 are statistically significant at the 5% level.

3.3 The Blending Weights

Since the unweighted respondents from the WWP report markedly different characteristics than those given by respondents from the probability sample (KP) and since policymakers need nationally representative estimates, we apply the weighting techniques of Section 2 for the purpose of blending KP and WWP samples.

3.3.1 Calculation of Propensity Scores

As outlined in Section 2, calculation of propensity score-based weights requires two sets of probabilities as input: $\{d_{i}\}$ and $\{\gamma_{i}\}$. Recall that it is required that $d_{i}$, which gives the likelihood of unit $i$ being selected into the probability sample, be known for each $i\in S$. Recall also that a pre-stratification weight was assigned to each post-9/11 caregiver within the KP and that these weights incorporate adjustments for non-response among the probability sample (i.e., non-compliance with the initial KP screener) and any residual inconsistencies between the KP and the full population; therefore, we set $\{d_{i}\}_{i\in S_{1}}$ equal to the inverse of the initial pre-stratification weights. Due to the proprietary nature of the KP, we lack data on screener non-respondents, and therefore cannot fully model screener non-response so as to calculate an informed value of $\{d_{i}\}_{i\in S_{2}}$. Hence, we set $d_{i}=n_{1}/\sum_{j\in S_{1}}d_{j}^{-1}$ for each $i\in S_{2}$, where $n_{1}$ is the sample size of $S_{1}$—this effectively assumes that all cases in $S_{2}$ had equal probability of inclusion into $S_{1}$. For non-response adjustments that are more rigorous, we rely on calibration weights.

Next, we calculate $\gamma_{i}$, the probability unit $i$ sampled unit being selected within the convenience sample for each $i\in S$ (i.e., the propensity score). These probabilities are unknown but can be easily estimated. Specifically, we use logistic regression:

[TABLE]

where ${\bf x}_{i}$ is a set of auxiliary variables (the same set is used within calibration weighting) and $(\zeta_{0},\bm{\zeta}^{\prime})^{\prime}$ is a vector of regression parameters. As a substitute for $\{\gamma_{i}\}$, we use $\{\hat{\gamma}_{i}\}$, which contains the predicted values (derived via the above regression) of the elements of $\{\gamma_{i}\}$. Non- and semi-parametric models for binary response could be used (e.g., Hahn, 1998; Hirano et al., 2003; Frölich, 2006), and boosted regression models (Friedman, 2001, 2002; Ridgeway et al., 2014) have become increasingly popular for propensity score estimation; however, such techniques are not given further consideration here.

3.3.2 Calibration Weighting

When applying calibration here and in Section 4, we use the function calibrate() within the survey package in R (Lumley, 2004, 2011). Specifically, we set calfun = ’linear’ and bounds = c(0,Inf). As initial weights in the simultaneous calibration method, weights from the propensity score-based scheme were used. Details regarding calculation of benchmarks for calibration in this application are provided in the supplemental materials.

3.3.3 Results from Blending

We calculated three sets of weights to blend the samples of post-9/11 caregivers: 1) simultaneous propensity scores (SPS), 2) disjoint propensity scores (DPS) and 3) simultaneous calibration (SC). Weights based on disjoint calibration could not be calculated since the algorithm to do so failed (this is likely due to KP and WWP respondents being substantially different and thereby disabling the feasibility of a solution to the calibration equations in (9)). Standard errors reported in this section were calculated using a delete-a-group jackknife with $G=40$ and with reweighting for each replicate group.

A subset of auxiliary variables (${\bf x}_{i}$) used in calculation of all weights includes all variables listed in Table 3. These variables include five questions that gauge the status of the caregiver as being an “early adopter” of technology. Such questions have proven to be useful for weighting respondents of online surveys (DiSogra et al., 2011) by helping to quantify their familiarity with and use of electronic media. Table 3 also lists the estimated mean of each variable when calculated using each of the three sets of blending weights. Standard errors of these auxiliary variables are not provided. Note that weight trimming (Lee et al., 2011; Potter and Zheng, 2015) was used to reduce the influence of outlying weights. Specifically, as the final step in weighting, the highest and lowest 1% of weights were truncated to match the corresponding percentile, and weights within these bounds were adjusted so that the sum of the weights remains unchanged. As a consequence, the mean of the auxiliary variables when calculated using calibration weighting differs slightly from the benchmarks values. Additionally, Table 3 indicates that weights based on propensity scores yield means of calibration variables that closely approximate the benchmark values. The full set of auxiliary variables for blending is shown in Table A.2 of the supplemental materials along with information analogous to that which is provided in Tables 2 and 3.

Our set of auxiliary variables chosen not only because the KP and WWP caregivers differ on the basis of these characteristics but also because the selection of auxiliary variables was limited to demographic characteristics of the caregiver and descriptors of the care recipient. Descriptors of caregiver well-being are outcome variables for this study. Table 4 lists selected outcome variables when calculated using only the KP caregivers, only the WWP caregivers, and while using the three sets of blended weights. Table A.3 provides this information for a larger set of outcome variables. For each method of blending (e.g., SPS, DPS, SC) and for each outcome variable, a $p$-value of the hypotheses in (11) is provided. Corresponding $p$-values calculated using unweighted data (not shown for brevity) quantify the discrepancies between the two samples. Specifically, of the 31 outcomes listed in Table A.3, 23 have an unweighted $p$-value less than 0.05 (16 have a $p$-value less than 0.01).

As seen by comparing of Tables 2 and 4 (and by examining Table A.3, the means of the outcome variables when found using weighted (blended) data tend to be close to the means found using only the KP sample; uniformly, the blended means are well within the error bound for the values found from the KP. Additionally, the standard errors of the blended values imply that an improvement in precision has been obtained by using the convenience sample. However, we see that simultaneous blending offers smaller standard errors (and therefore a greater increase in precision) than disjoint blending. This observations yield the conclusion that the simultaneous blending weights are likely to provide smaller design effects—this claim is investigated further in Section 4. In some circumstances, disjoint blending yields larger standard errors than were observed when only the KP was used. This is a consequence of the blended data having a higher mean than was observed in only the probability sample, and is not evidence of a loss of precision.

The $p$-values reported in Table 4 (and Table A.3) are often close to zero when simultaneous weights are used—this is expected and is not a testament to the quality of weighting. For the disjoint blending method, the $p$-values are larger—these $p$-values give a more accurate assessment of the adequacy of weighting (with respect to the selection of auxiliary variables). None of the 31 outcomes in Table A.3 have a $p$-value less than 0.05 under disjoint blending—this provides evidence that the set of auxiliary variables in Table 2 sufficiently account for differences between KP and WWP post-9/11 caregivers.

When disjoint weights were used, we also calculated the post hoc blending estimator described in (5); however, full results are omitted for brevity. We note that the point estimates yielded by the post hoc method are consistently similar to those seen for the DPS procedure in Table 4. For example, in the depression and anxiety outcomes, the post hoc estimators are 7.30 and 42.97, respectively (compared to respective values of 7.57 and 43.66 for DPS). The post hoc blending estimator is examined further in Section 4.

As a final data analysis, we discuss the influence of weighting on inferences drawn from the following regression model of a military caregiver’s depression level:

[TABLE]

In the above, $\mbox{DEP}_{i}$ is the depression score for the $i^{\rm th}$ caregiver, $\mbox{ERA}_{i}$ is an indicator variable which is unity if the caregiver is post-9/11, and $\mbox{AGE}_{i}$, $\mbox{SEX}_{i}$, $\mbox{INC}_{i}$ and $\mbox{EMP}_{i}$ are the age (continuous), gender (binary), income (continuous) and employment status (binary), respectively, of the caregiver. Using all military caregivers (pre- and post-9/11) from the KP and WWP, we estimate the model in (13). The various blending weights are used for post-9/11 caregivers, and the pre-stratification weights are applied to the pre-9/11 caregivers; results are seen in Table 5.

Of primary interest is whether or not $\eta_{1}$ differs from zero; if $\eta_{1}>0$, then post-9/11 caregivers have higher depression levels (conditional on covariates) than pre-9/11 caregivers. Using weighted least squares (WLS) with the initial pre-stratification weights, the regression model is estimated while incorporating data from only military caregivers from the KP sample. When only military caregivers from the probability sample are used, the results illustrate that post-9/11 caregivers have higher levels of depression, but the estimate is not statistically different from zero at the 5% significance level. The lack of statistical significance is likely a consequence of sample size. All methods of calculating blending weights implicate that era of service is statistically significant (at the 5% level) with post-9/11 caregivers having higher levels of depression when controlling for other characteristics. However, the era of service coefficient is less statistically significant when calculated using disjoint blending. Further, the discrepancy between levels of depression for pre- and post-9/11 caregivers is smaller when any method of blending is used than when the unweighted WWP data are used. The emotional hardships endured by post-9/11 military caregivers are thought to be a consequence (in part) of the types of injuries suffered by soldiers that served in the recent campaigns and the fact that these caregivers are often spouses and family members tending to a young veteran with a debilitating injury.

The sample means and standard errors reported in this section are calculated using the svymean() function from the survey package in R. Regression results and $p$-values for comparing samples are derived using the svyglm() function.

4 Simulations

We perform simulation studies to examine efficacy of the weighting methodology in greater detail. First, we provide simulations performed using data from the caregiver study that are designed to assess bias and mean-squared error in situations where model assumptions are and are not met. We also perform simulations with purely synthetic data that are designed to illustrate the efficacy of variance estimators.

4.1 Simulations with caregiver data

Our first simulation study uses data from the caregiver survey. We form a pseudo-population of post-9/11 caregivers using the observed data and then repeatedly draw samples from this population. Prior to describing construction of the pseudo-population, we note that an opt-in sample military caregivers was drawn from a volunteer internet-based panel—this sample yielded 171 post-9/11 caregivers (see Table 1). The opt-in sample was found to contain data anomalies (e.g., clearly erroneous response patterns) and was disregarded from the analyses provided by Ramchand et al. (2014). However, it is used here to help build the pseudo-caregiver population. Specifically, the pseudo-population of post-9/11 caregivers consists of all observed military caregivers from the KP and opt-in samples. The WWP is excluded from this phase of our study because it represents a highly differentiated segment of the caregiver population (e.g., the WWP is nearly entirely post-9/11 and is heavily burdened). Our pseudo-population has 940 cases.

Approximately 150 caregiver are selected for inclusion in the probability sample by drawing units from the pseudo-population at random with probability $d_{i}^{*}=0.16$. Selected caregivers are considered respondents with a probability $r_{i}$ that satisfies $\log(r_{i}/(1-r_{i}))=(1/3)\mbox{F}_{i}-(2/3)\mbox{A}_{i}$, where $\mbox{F}_{i}$ is an indicator of whether or not the caregiver is female and $\mbox{A}_{i}$ is the caregiver’s age categorized in a 5-point Likert scale (both variables are centered). A convenience sample of is drawn from the pseudo-population by using the following logistic model:

[TABLE]

where $\rho_{i}$ is the probability that case $i$ is assigned to the convenience sample. Further, ${\bf v}_{i}$ is a set of variables upon which the probability and convenience samples may be differentiaed and $(b_{0},\bm{b}_{1}^{\prime})^{\prime}$ is a vector of coefficients. When applying (14), we use standardized versions of all variables so as to ensure that the coefficients $\bm{b}_{1}$ have interpretability. If a unit is selected into both the probability and convenience samples, it is assigned to the probability sample. Note that the selection mechanism in (14) is treated as being unknown when our methods are applied; therefore, a sample drawn in accordance with (14) is a convenience sample.

We draw convenience samples under five different settings for the purpose of exploring how the inclusion mechanism influences final inferences. The set of variables ${\bf v}_{i}$ and the values of the coefficients used in (14) are listed in Table 6 for each of the settings.

Across the settings, caregiver depression represents the outcome variable of interest and caregiver anxiety denotes a (sometimes) latent outcome variable that is strongly correlated with caregiver depression ($\rho\approx 0.65$). The remaining variables are used as auxiliary variables in all settings. Within Setting 1, neither depression nor anxiety influence selection into the convenience sample (although the remaining variables have a strong influence on the probability of inclusion). Setting 2 is the same as the first, except anxiety (which does not influence the selection mechanism) is used as an auxiliary variable (this setting is considered to help assess the need for parsimony in the selection of auxiliary variables). Setting 3 introduces circumstances where a latent outcome (anxiety) that is not the outcome of interest influences selection (note that anxiety is not used as an auxiliary variable for blending in this setting). Setting 4 presents circumstances where the outcome of interest directly influences inclusion probabilities. Setting 5 is identical to the fourth setting with the exception that anxiety is added to the list of auxiliary variables used for blending. This final setting is designed to investigate whether or not one can lessen the effect of bias by including covariates that are correlated with the outcome (but do not necessarily have a direct influence on the sampling mechanisms) as auxiliary variables within the weighting schemes. Within Settings 3, 4, and 5, the degree to which blending is insufficient (e.g., the samples are directly differentiated by depression levels) is controlled by the tuning parameter $\tau$. Although we use $\tau=1/2$ here, simulations are provided in the supplementary materials that consider various values of this coefficient.

A visualization of the discrepancies between the pseudo-population, the probability sample, and the convenience sample is provided in Table A.1 of the supplemental materials to this article.

Each setting of this simulation study involves $K=10,000$ independent iterations. In each iteration, probability and convenience samples are drawn from the pseudo-population of post-9/11 military caregivers. Each of these samples are drawn independently of one another within each iteration (although a unit that is selected for both the convenience and probability samples is assigned to the probability sample). The four weighting methods described in Section 2 are applied within each of the $K$ iterations of the simulation. To review, the methods include simultaneous propensity scores (SPS), disjoint propensity scores (DPS), simultaneous calibration (SC) and disjoint calibration (DC). The set of auxiliary variables (${\bf x}_{i}$) used in calculation of weights for each setting is listed in Table 6.

Estimated response probabilities, $\hat{r}_{i}$, for the probability sample are calculated using a logistic model, and the estimated probability of inclusion into the convenience sample, $d_{i}$, is then given by $\hat{d}_{i}=d_{i}^{*}\hat{r}_{i}$. When the propensity score-based methods are applied, the propensity scores are calculated using the logistic model in (12). However, $\{\gamma_{i}\}$ does not necessarily obey the stated logistic function under the schemes used to draw the probability and convenience samples here. When calibration is applied, totals are estimated using the probability sample (weighted with $\hat{d}_{i}^{-1}$). However, each data unit is given the same initial weight value (which is set to the number of units in the population divided by the number of unit sampled) in the calculation of calibrated weights.

To explain the bookkeeping of bias and error in parameter estimators, let $\mu$ denote the population parameter of interest, which in this case is the mean depression level of post-9/11 caregivers. Let $\hat{\mu}$ denote a value of $\mu$ calculated using the complete pseudo-populations ($\hat{\mu}$ represents a benchmark value of $\mu$). Further, let $\hat{\mu}^{[k]}_{j}$ denote the estimated value of $\mu$ calculated using the $k^{\rm th}$ replication (for $k\in(1,\ldots,K$) and the $j^{\rm th}$ weighting scheme (for $j\in(1,\ldots,4)$). The relative error (reported as a percent change from the benchmark) inherent in $\hat{\mu}^{[k]}_{j}$ is calculated via $\mbox{e}^{[k]}_{j}=100(\hat{\mu}^{[k]}_{j}-\hat{\mu})/\hat{\mu}$ whereas the squared relative error is given by $\mbox{SRE}^{[k]}_{j}=(\mbox{e}^{[k]}_{j})^{2}.$ To aggregate the findings across all iterations, we approximate the bias and root-mean squared error via

[TABLE]

Let $p^{[k]}_{j}$ denote the $p$-value of a test of the hypotheses in (11) (or related hypotheses if a more exhaustive regression model is used) for the $k^{\rm th}$ replication with the $j^{\rm th}$ weighting scheme. We report

[TABLE]

which denotes the rate at which the null hypothesis of adequate blending is rejected. Recall that $1_{\{A\}}$ denotes the indicator of event $A$. We use a significance level of $\alpha=0.05$. The rejection rate $\hat{p}_{j}$ approximates type I error in the event that ${\cal H}_{0}$ is true (e.g., Settings 1 and 2) and the power of the test when ${\cal H}_{1}$ is true (e.g., Settings 3–5). We also calculate the design effect that is seen in $\hat{\mu}$ and we calculate the standard error of this quantity; to avoid redundancy, we only report the design effect (the comparative results for the standard errors are similar). Lastly, we calculate the post hoc blending estimator from (5) when using the two sets of disjoint weights. To approximate the variances and covariances in (6) for the post hoc estimator, we use a jackknife. Only the bias and rMSE are reported for the post hoc estimator of $\mu$. Standard errors and the design effect corresponding to $\hat{\mu}$ are calculated using Taylor series linearization. To ensure computational feasibility of these simulations, we do not use resampling approaches.

Findings are reported in Table 7 for $\tau=1/2$. Within each of the settings, weighted blending (regardless of method) is preferable to unweighted blending. As expected, if the discrepancies between the probability and convenience samples have been appropriately modeled (i.e., Setting 1 and 2), weighted blending produces more accurate (lower bias) and more precise (lower rMSE) estimators than those which are found using only the probability sample; otherwise (i.e., Settings 3–5), it may be preferable to use only the probability sample. It appears that calibration and propensity scores produce similar results. However, we note that the propensity score-based methods tend to yield lower design effects, whereas the calibration methods often observe lower rMSE. Further, we see that simultaneous weighting yields lower design effects and consequentially lower rMSE than disjoint weighting methods. We also see that when blending is based on a sufficient model (i.e., Setting 1), the rejection rate of the hypotheses in (11) is close to its nominal level so long as disjoint blending is used; also, our expectation is validated that simultaneous weights are inappropriate for testing the sufficiency of blending.

From Setting 2, we learn that parsimony in the choice of auxiliary variables is beneficial since the second setting observes higher rMSE than the first (likewise, the test for the adequacy of blending is conservative in the Setting 2). However, the remaining settings illustrate that it is necessary to have a robust set of auxiliary variables. That is, even though anxiety is not the outcome of interest, we see in Setting 3 that allowing the probability of selection into the convenience sample to depend upon anxiety induces bias into estimators found by the weighting schemes. This bias would be smaller if depression and anxiety were not highly correlated. Similarly, Setting 5 (when compared to Setting 3) illustrates that if the probability of selection into the convenience sample depends on the outcome of interest (depression), bias can be reduced by using additional variables that are correlated with the outcome as auxiliary variables in the calculation of weights.

From Table 7, we also see that the post hoc blending estimator performs comparably to the corresponding estimators founding using the corresponding disjoint blending weights in all settings. Therefore, we conclude that there is little to no loss of efficiency that stems from the use of the approximated design effect in our disjoint blending technique in comparison to more rigorous methods of variance minimization.

The supplemental materials present further simulations in this setup that study the effect of the parameter $\tau$ on the performance of the methods.

4.2 Simulations with synthetic data

Here, we compare methods for variance estimation with blended data using synthetic data. In these simulations, we vary the degree to which auxiliary variables used for blending are correlated with the outcome of interest. For a population of size $N=10,000$, we generate a 2-dimensional vector of auxiliary variables via ${\bf X}=(X_{1},X_{2},X_{3})\sim N({\bf 0},{\bf I}_{3})$, where ${\bf I}_{k}$ is the identity matrix of dimension $k$. The outcome is simulated using $Y=\beta(X_{1}+X_{2})+\epsilon$, where $\epsilon\sim N(0,\sigma^{2}_{e})$. We choose $\beta$ and $\sigma_{e}^{2}$ so as to ensure the coefficient of determination, $R^{2}=2\beta^{2}/(2\beta^{2}+\sigma_{e}^{2})$, takes on a desired value while maintaining $\mbox{Var}(Y)=1$. That is, for a given value of $R^{2}$, we set $\beta=\sqrt{R^{2}/2}$ and $\sigma_{e}^{2}=1-R^{2}$. Next, we select $n_{1}=200$ case for the probability sample at random from the simulated population. Those cases are considered respondents (and consequentially elements of the probability sample) with probability $r_{i}$ that satisfies $r_{i}=1/\{1+\exp(-0.15X_{3})\}$. An element of the simulated population is assigned to the convenience sample with a probability determined by $p=1/\{1+\exp[-4.2-0.5(X_{1}+X_{2})]\}$. As before, data elements selected into both samples are assigned to the probability sample. Blending weights are then calculated using simultaneous propensity scores. Logistic modeling is used to estimate response probabilities and propensity scores as functions of $(X_{1},X_{2})$ and $X_{3}$, respectively.

Letting $\eta=\mbox{E}[Y]$, we calculate $\hat{\eta}$, the estimate of $\eta$, when found using only the probability sample and when found using the (weighted) blended sample. We also calculate the standard error of $\hat{\eta}$ when only the probability sample is used, when the blended weights are used with Taylor series linearization for variance estimation, and when a delete-a-group jackknife (as outlined in Section 2.4) with $G=40$ is used. Note that weights are recalculated for each replicate group in the jackknife. Using the estimated standard error (and the corresponding value of $\hat{\eta}$), we can calculate the upper and lower bounds of a $95\%$ confidence interval. For each confidence interval, we track the percent of iterations in which the true mean, $\eta=0$, falls in the interval—this is referred to as coverage.

Results are shown in Figure 1. We see that coverage remains at the nominal level of 95% when variance estimation is performed with a jackknife; however, coverage decreases with increasing $R^{2}$ when linearization is used. Therefore, one should not estimate variance with linearization when the auxiliary variables are moderately to strongly related with the outcome. Further, in this study, standard error remains constant across all values of $R^{2}$ when only the probability sample is used and when blending is performed with variance estimated using linearization. However, standard error increases with $R^{2}$ when the variance of the blending estimator is calculated with a jackknife. In fact, when $R^{2}=0$, the jackknife yields the same efficiency as linearization; however, when $R^{2}=1$, the jackknife yields a standard error equivalent to that provided by only the probability sample. Therefore, if one were to use the outcome as an auxiliary variable for blending, one would have little to no gain in precision from blending samples over what would be achieved with only the probability sample. These observations imply that the design effect and rejection rates reported for Settings 2 and 5 in Table 7 are likely understated since linearization was used (in Settings 2 and 5, $R^{2}=0.42$; otherwise, $R^{2}=0.14$).

5 Discussion

We introduced four methods for calculating weights that blend probability and convenience samples (using combinations of disjoint vs. simultaneous blending and propensity scores vs. calibration). Simultaneous methods appear to yield lower design effects (and variances). In addition, simultaneous blending does not require coverage of the convenience sample (i.e., Assumption 5). However, disjoint blending is needed to assess the adequacy of the set of auxiliary variables. Calibration and propensity score procedures perform comparably, with propensity scores yielding slightly smaller design effects in our empirical illustrations. Note that calibration requires a feasible solution to the calibration equations (which was not satisfied for disjoint calibration in our application). For variance estimation, a jackknife was shown to outperform linearization in simulations.

We revisit the assumptions mandated by our approaches and discuss the degree to which they were satisfied in our empirical examples. With respect to Assumptions 1 and 2 in the caregiver data, the probability sample appears to cover the population; however, as noted in Section 3.3.1, the proprietary nature of the KP prevented us from rigorously adjusting for non-response when estimating $d_{i}$ for the convenience sample. Our solution to this issue was to assume that KP non-respondents are not substantially different from respondents, which appears reasonable observationally. Further, $d_{i}$ is not needed for all $i\in S$ for calibration methods, and propensity score methods provided results that were similar to calibration in the military caregivers application. Assumption 3 appears upheld in our data application since the test for the adequacy of blending did not reject for any outcomes. Assumption 4 (adequacy of the propensity score model) may not have been upheld. In fact, in our simulation studies, it was not correct (since a logistic model was used to draw the convenience sample, the propensity score model will have a different form). The simulations showed that despite this, the propensity score-based techniques perform well. Lastly, although the WWP was strongly differentiated from the KP, we saw no evidence to suggest that Assumption 5 (positivity of the convenience sample) was not satisfied.

Although we made efforts within our simulation study to provide guidance regarding how one should select a choice of auxiliary variables that are to be used in the methods outlined here, we feel that more work in this vein is needed. For instance, Brookhart et al. (2006) provide a thorough assessment of model selection within traditional propensity score procedures. They suggest that variables that are not related to treatment status but related to the outcome should always be included. However, this guidance does not translate to our work: on account of the findings of Section 4.2, building propensity scores for blending while unnecessarily using auxiliary variables that are strongly related to the outcome will lead to a reduction in the gain in precision provided by the convenience sample. Therefore, we recommend being parsimonious when selecting auxiliary variables.