Improved calibration estimators for the total cost of health programs and application to immunization in Brazil

Multi-stage/level sampling designs have been widely used by survey statisticians as a means of obtaining reliable and efficient estimates at a reasonable implementation cost. This method has been particularly useful in National country-wide surveys to assess the costs of delivering public health programs, which are generally originated in different levels of service management and delivery. Unbiased and efficient estimates of costs are essential to adequately allocate resources and inform policy and planning. In recent years, the global health community has become increasingly interested in estimating the costs of immunization programs. In such programs, part of the cost correspond to vaccines and it is in most countries procured at the central level, while the rest of the costs are incurred in states, municipalities and health facilities, respectively. As such, total program cost is a result of adding these costs, and its variance should account for the relation between the totals at the different levels. An additional challenge is the missing information at the various levels. A variety of methods have been developed to compensate for this missing data. Weighting adjustments are often used to make the estimates consistent with readily-available information. For estimation of total program costs this implies adjusting the estimates at each level to comply with the characteristics of the country. In 2014, A National study to estimate the costs of the Brazilian National Immunization Program was initiated, requested by the Ministry of Health and with the support of international partners. We formulate a quick and useful way to compute the variance and deal with missing values at the various levels. Our approach involves calibrating the weights at each level using additional readily-available information such as the total number of doses administered. Taking the Brazilian immunization costing study as an example, this approach results in substantial gains in both efficiency and precision of the cost estimate.


Estimation of total cost T
The scientific goal is the estimation of a total. This total is composed of totals at different levels. For example, for the Brazilian immunization study there are costs that are specific to the regions, some other costs are specific to the municipalities and some others are specific to the facilities. We introduce the following notation in order to develop the details for a three multi-stage design.
Let y Frmf be the value specific to health facility f , in municipality m, in r.
Let y M rm represents the value specific to municipality m, in region r.
Let y Rr represents the value specific to region r.
In addition, let U R be the set of all regions, U Mr be the set of all municipalities in region r and U Fmr be the set of all facilities in municipality m, region r. There are in total N R regions, N Mr municipalities in region r and N Fmr facilities in municipality m in region r. The total cost can be represented as (1) where T C represents the additional cost accrued centrally and it is assumed known.

The design
To estimate (1), a multi-stage simple sampling design is used. Let s R represent the sub sample of regions, s Mr represent the sub sample of municipalities in region r and s Fmr represent the sub sample facilities in municipality m in region r, respectively. The design is as follows. At first stage, a simple sample, s R , of n R regions is selected from U R with replacement. At second stage, for each region in s R , a sample of n Mr municipalities is selected from each U Mr . Finally, at third stage, within each municipality in s Mr , a sample of n Fmr facilities is selected from U Fmr . Each facility has a probability of being included in the final sample. This design can be characterize by its sampling probabilities. These probabilities are used in the estimation of the total T . This sampling probability for facility f from municipality m, from region r, can be written as These probabilities can be written in this form only if we can assume that the sampling designs within each region was planned before the regions were selected and that the sampling design was implemented as planned, then we have a principle called invariance. Otherwise, it is not guaranteed that the probabilities can be calculated. Note that where I F rmf , I Mrm and I Rr respectively represent the inclusion indicators for region r, municipality m in region r and faiclity f in municipality m in region r. The conditional inclusion probabilities are defined as π f |mr = E I F rmf =1 I Mrm=1 = n Fmr N Fmr ; π Rr = E (I Rr =1 ) = n R N R and π m|r = E ( I Mrm=1 | I Rr =1 ) = n Mr N Mr . Using this, the number of facilities sampled in municipality, region r can be written as n Frm =

Estimating T F
We now present the estimation of the total and its standard error in the next sections.
• Under the design above, the weighted estimator of T F is given by where each w Rrmf represents the sampling weight for facility f from municipality m, from region r.
• These weights are given by where π Rrmf is the inclusion probability for facility f in municipality m in region r.

Estimating T M
• Let w M rm = 1 π m|r π Rr ; that is, the weight for municipality m, which belongs to region r.
• The weighted estimator of T M is given by This estimator is unbiased for T M . The reason why we express (5) in terms of three sums is because it is used to calculate the variance.

Estimating T R
• Let w Rr = 1 π Rr , that is, the weight corresponding to region r.
• The weighted estimator of T R is given by This estimator is unbiased for T R .

The estimator of the total
Usually, IPW estimators only involve data at the last stage of the design, such as T F . In this case, the variance of the estimator is well known and straight forward to calculate following standard properties of multi stage designs . Our estimator of interest, however, involves information at every stage of the design. In order to find an expression for the variance of T , we define the variables: Then, the estimator T be written as and the total can be expressed as Note that E T u = T.

Variance of T
The variance of (9) can be calculated as where the first term represents the variation due to the sampling design and can be written as where ∆ Rrs = π Rrs − π Rr π Rs , π Rrs = P (I Rr = 1, I Rs = 1) which represents the pairwise inclusion probabilities at the regional level,ť Rr = w Rr The term V Rr is calculated in the same way as (12) . For a three multi-stage design, we can express (12) as This expression is found using conditional expectations properties, where V R represents the variation due to the sampling of regions, V M represents the variation due to the sampling of municipalities V F is the variation due to the sampling of health facilities (see, e.g.  ). Under the underlying design, they have the forms where and The second term of (11) represents the probabilistic variation that would have been obtained had we had complete data. This can be expressed as 3 Estimation of the variance of T In practice, only the final sample of facilities is known, as so u rmf is only known for the sample. This implies that the variance (13) cannot be calculated. However, we can estimate it using different approaches. The first approach, is to analytically find an unbiased estimator. A second approach is to use bootstrap re sampling, which return a consistent estimator.

Analytic estimation of E V ar T u
Under the underlying design, an unbiased estimator of the term (12) can be obtained analytically. An unbiased estimator of this is given by where and

Analytic estimation of V ar E T u
Under the underlying design, an unbiased estimator of the term (20) can be obtained analytically. An consistent estimator of this is given by Or an unbiased estimator is given by

Post-design methods: Calibrated weights
The main purpose of post-designs methods is to improve the estimates by reducing the variance, i. e to produce more efficient estimators. Post-design methods include techniques such as post stratification, estimation of weights and calibration of weights. The latter is known to yield more efficient estimates by adjusting the weights using information readily-available at first phase. In Brazil, for example, the total number of doses administered is known and this can be use to adjust the weights. Breslow et al. [2009a] and Breslow et al. [2009b] proposed using calibration to adjust the sampling weights in case-cohort designs. Let v F , v M and v R denote readily available variables at first phase and at the different levels. The former is a variable available for all units at the final stage, the second is available for all second stage units and the latter is available for all first-stage units, e.g regions. The main goal is to adjust weights such that estimates for the totals T v F , T v M and T v R equal the actual totals T v F , T v M ) and T v R ), but forcing the new weights ( w) to be as close as possible to the design weights. For example, consider IPW estimates T v F and T v M based on the same weights used to estimate T with T , i.e.
Since we know T v F , T v M and T v R we can modify the weights to force (30) Intuitively, the modified weights ( w F rmf , w Mrm ) contain information about the entire population of facilities and we can use the new weights to find a new estimate of T .
Mathematically, a distance function d(w, w) measuring the distance from original weights(w) to new weights( w) is minimized subject to the constrains (30). There are several options for choosing this distance. Substantial gains in accuracy depend on how correlated the calibration variables are with the quantity of interest [Deville and Sarndal, 1992]. In order to find the asymptotic variance of T cal , we follow [Deville and Sarndal, 1992]. Choosing one the distance functions suggested by them ensures that we can write where T F,reg , T M,reg and T R,reg are the so-called regression estimators shown in [Deville and Sarndal, 1992] given by and Using (31) we can write: and

Variance estimation
In order to estimate the variance in practice we use the expression where, using Taylor expansion as in , we obtain that the approximate variance is th same as the variance of w Rr E Rr , with We now define u rmf = e F rmf + w Mrm w F rmf n Fmr e Mrm + w Rr w F rmf n Fr e Rr , where e F rmf , e Mrm and e Rr are estimates for E F rmf , E Mrm and E Rr . They are given by  (40) is found in a similar way as the variance in section 2.