Calibration-weighting a stratified simple random sample with SUDAAN

Raking and Raking to a Size Variable

One starts by downloading the data set (DAWN) in the Appendix and creating some additional variables in a new data set (R) with this SAS code:

DATA R; SET DAWN;

Q = FRAME_VISITS/1000;

W1 = W − 1;

QW1 = Q*W1;

PUBLICQ = PUBLIC*Q;

PUBLICQW1 = PUBLIC*Q*W1;

METROQ = METRO*Q;

METROQW1 = METRO*Q*W1;

One can estimate totals for the United States and within the four regions with this (SAS-callable SUDAAN) code for the expansion estimator:

PROC DESCRIPT DATA = R;

DESIGN = STRWOR; * the design is stratified simple random sampling without replacement;

NEST STRATUM;

WEIGHT W;

TOTCNT BIG_N; * these stratum population sizes are needed for finite-population correction;

CLASS REGION; * estimates will be generated by region;

VAR DR_VISITS; * the survey variable we are estimating;

OUTPUT TOTAL SETOTAL/FILENAME = OUT0 REPLACE;

RUN;

Rather than looking at the output, which the DESCRIPT procedure usually produces, we save it (i.e., the estimated regional and US totals and their standard errors) in OUT0. If one runs SAS STUDIO, then standard SUDAAN output will not appear. Instead, the user can output results into a SAS data set as is in the code above.

The strata in the sample are almost completely cross-classified by region, public/private, and metro/non-metro—but not quite. There is only one hospital in the sample in the East (REGION = 1) that is private (PUBLIC = 0) and not urban (METRO = 0), but it is not in its own stratum.

Raking is an iterative procedure to try to get the weighted estimates of the numbers of EDs in each category to equal the frame numbers in those categories. From the frame, we know the total number of public EDs (1642), metro EDs (856), and EDs in each of the four regions (489, 1636, 3124, and 1051). One rakes the weights to those totals using the WTADJUST procedure with the following code and outputs the results (estimated totals and their standard errors) into OUT1:

PROC WTADJUST DATA = R DESIGN = STRWOR;

ADJUST = POST; * calibration will be to totals produced outside the sample;

NEST STRATUM; WEIGHT W;

TOTCNT BIG_N;

CLASS REGION;

VAR DR_VISITS;

MODEL _ONE_ = PUBLIC METRO REGION/NOINT; *_ONE_ = 1;

POSTWGT 1642 856 489 1636 3124 1051;

OUTPUT TOTAL SE_TOTAL/FILENAME=OUT1 REPLACE;

RUN;

Note that here, SUDAAN treats PUBLIC and METRO as continuous variables whereas REGION is a categorical variable with four levels. Because it is a categorical variable, it appears in the CLASS statement, and WTADJUST will produce estimated means (drug-related visits per ED) and totals at the regional and US (all regions) levels.

WTADJUST treats raking like a regression with the raking variables as the explanatory variables and 1 (_ONE_) as the dependent variable. /NOINT is added to the model statement because there is no intercept.

The zeros in the sixth column in the output table (Table 1) tell us that calibration worked: The weighted numbers metro, public, and EDs in each region (the “sum[s] of the final adjusted weights over respondents”) match the numbers in the frame (the “control totals”).

Alternatively, the statement

OUTPUT TOTORIG TOTTRIM TOTFINAL CNTLTOTAL DIFFWT/ FILENAME DIFFWT REPLACE;

will put the second through sixth column from the table into an output data set DIFFWT.

Expressed mathematically, calibration weighting uses row vectors of calibration variables (PUBLIC, METRO, and REGION) denoted by xk for each k and the row vector Tx for the frame totals of components of xk to convert the design weights (dk) into calibrated weights (wk) by finding a column vector g—if one exists—that solves the calibration equation. Here, it is

where $w_k=d_k\exp (x_{k}g)$ are the calibrated weights for the k in S.

WTADJUST produces two other tables. One displays the components of g, which it calls “beta coefficients.” In this case, they are “estimates” of 0 and uninteresting in themselves. Another table produced by WTADJUST displays the estimates of the total $(\sum w_k{y}_k)$ and mean $(\sum w_k{y}_k/\sum w_k)$ of DR_VISITS by region and across all regions (the summations are over the sample in the relevant domain).

The hope is that by calibrating the weights to the numbers of public, metro, and regional EDs standard errors will decrease relative to the expansion estimator. It turns out, as we will see later in the section, they do not. Instead, when the goal is to estimate DR_VISITS, it makes more sense to find a g satisfying the following calibration equation:

where zk = xk ^. FRAME_VISITS (of ED k), Tz is the population total of zk, and, as before, $w_k=d_k\exp (x_{k}g)$. The number of components in zk and xk need to be equal.

Equation (2) forces the number of weighted estimates of public, metro, and regional frame visits to match their full-population targets. Satisfying the calibration equation in (2) is more sensible than the one in (1) because annual drug-related ED visits are more closely related to previous annual ED visits than the number of EDs. Observe that if Ty were exactly equal to Tz b* for some vector b*, ${t_y}^{\left(w\right)}={\sum }_S{w}_k{y}_k$ with the calibrated weights satisfying Equation (2) would estimate Ty perfectly.

The following code produces the calibrated weights and resulting estimates (in OUT2) from calibrating with Equation (2):

PROC WTADJX DATA = R DESIGN = STRWOR ADJUST = POST;

NEST STRATUM;

WEIGHT W;

TOTCNT BIG_N;

CLASS REGION; VAR DR_VISITS;

MODEL _ONE_ = PUBLIC METRO REGION/NOINT;

CALVARS PUBLICQ METROQ REGION*Q/NOINT;

POSTWGT 58000 44000 22000 43000 33000 36000;

OUTPUT TOTAL SE_TOTAL/FILENAME=OUT2 REPLACE;

RUN;

The MODEL statement is the same as in the previous WTADJUST, but a CALVARS statement has been added that contains the new calibration variables: PUBLICQ, METROQ, and REGION * Q (note that FRAME_VISITS has been divided by 1000 in each of these variables). The associated population totals, which have been computed using frame information, appear in the POSTWGT statement.

Solving for g, which the procedure does implicitly, is akin to a logistic regression without an intercept. The NOINT in both the MODEL and CALVARS statements indicates that there is no intercept.

A Short Divergence Into Theory

The $\exp (x_{k}g)$ in Equations (1) and (2) and special cases of the following weight-adjustment function are as follows:

where L = 0 and U = infinity (actually, it equals 10²⁰, which is close enough to infinity for all practical purposes). Those are the defaults for WTADJUST and WTADJX. The values LOWERBD (L) and UPPERBD (U) bound the size of the weight-adjustment function between L and U. LOWERBD cannot be negative, and UPPERBD must exceed LOWERBD. There is often a centering parameter, CENTER, in what SUDAAN calls “the general exponential model” or “GEM.” For our purposes, CENTER need only be defined when UPPERBD takes on its near-infinite default value and LOWERBD is not less than 1. GEM also allows the LOWEBD, UPPERBD, and CENTER to vary across the k, but that wrinkle does not concern us here.

Observe that with calibrated weights satisfying Equation (2),

where b* is the probability limit of

as the sample and population size grow arbitrarily large. We assume that b* exists.

Now Tzb* is a constant, and under mild conditions (that we assume to hold) b ≈ b* and wk ≈ dk, so the variance of ty^(w) is nearly equal to the variance of $t_e={\sum }_S{w}_k{e}_k$ where ek = yk − zkb; b is treated as if it were a constant. This is what WTADJX computes. So does WTADJUST when zk = xk. It appears that to keep standard errors low, one should find calibration variables that produce ek with small absolute values.

Two Examples of Pseudo-Optimal Calibration

Although b in Equation (5) looks like an estimated regression coefficient, it is not necessarily attached to a linear model. In fact, were the sample drawn using Poisson sampling, a pseudo-optimal version of b given the calibration vector zk sets xk = zk(dk − 1) (see, for example, Kott, 2011) and will usually produce smaller standard errors for estimated totals.

In the previous run of WTADJX the weight adjustment was a function of the model variables, PUBIC, METRO, and the four regions. Replacing each with a model variable of the form Variable * Q * W – 1 will usually result in smaller standard errors. This pseudo-optimal calibration is coded below with estimated totals and their estimated standard errors outputted into data set OUT3.

PROC WTADJX DATA = R DESIGN = STRWOR ADJUST = POST;

NEST STRATUM;

WEIGHT W;

TOTCNT BIG_N;

CLASS REGION; VAR DR_VISITS;

MODEL _ONE_ = PUBLICQW1 METROQW1 REGION*QW1/NOINT;

CALVARS PUBLICQ METROQ REGION*Q/NOINT;

POSTWGT 58000 44000 22000 43000 33000 36000;

OUTPUT TOTAL SE_TOTAL/FILENAME = OUT3 REPLACE;

RUN;

An intercept is added to the code below by inserting W1 into the model statement and _ONE_ into the CALVARS statement and the number of EDs on the frame (6300) into the POSTWGT statement, while retaining the NOINTs. The estimated totals and their estimates standard errors outputted into data set OUT4.

PROC WTADJX DATA = R DESIGN = STRWOR ADJUST = POST;

NEST STRATUM; WEIGHT W;

TOTCNT BIG_N;

CLASS REGION; VAR DR_VISITS;

MODEL _ONE_ = W1 PUBLICQW1 METROQW1 REGION*QW1/NOINT;

CALVARS _ONE_ PUBLICQ METROQ REGION*Q/NOINT;

POSTWGT 6300 58000 44000 22000 43000 33000 36000;

OUTPUT TOTAL SE_TOTAL/FILENAME = OUT4 REPLACE;

RUN;

We are now ready to compare the estimated totals and their coefficients of variation (CVs) computed in the five sets of alternative weights (i.e., original design, raked, raked to a size variable, quasi-optimal calibrated without an intercept, and quasi-optimal calibrated with an intercept).

DATA OUT0; SET OUT0;

DESCV = SETOTAL/TOTAL;

DESTOT = TOTAL; RUN;

DATA OUT1; SET OUT1;

RAK1CV = SE_TOTAL/TOTAL;

RAK1TOT = TOTAL; RUN;

DATA OUT2; SET OUT2;

RAK2CV = SE_TOTAL/TOTAL;

RAK2TOT = TOTAL; RUN;

DATA OUT3; SET OUT3;

QO1CV = SE_TOTAL/TOTAL;

QO1TOT = TOTAL; RUN;

DATA OUT4; SET OUT4;

QO2CV = SE_TOTAL/TOTAL;

QO2TOT = TOTAL; RUN;

DATA C; MERGE OUT0 OUT1 OUT2 OUT3 OUT4; BY VARIABLE REGION;

DESCV = ROUND(DESCV * 100, .01);

RAK1CV= ROUND(RAK1CV * 100, .01);

RAK2CV= ROUND(RAK2CV * 100, .01);

QO1CV = ROUND(QO1CV * 100, .01);

QO2CV = ROUND(QO2CV * 100, .01);

PROC PRINT; ID REGION; VAR DESTOTAL RAK1TOTAL RAK2TOTAL QO1TOTAL QO2TOTAL;

PROC PRINT; ID REGION; VAR DESCV RAK1CV RAK2CV QO1CV QO2CV; RUN;

The results are shown in Table 2.

The CVs from raking (second column of CV results) are no better than those of the expansion estimator (first column). Raking to a size variable (third column) decreases the CVs considerably, and quasi-optimal calibration (fourth and fifth columns) usually decreases CVs even further but not so dramatically.

Some More Theory

WTADJUST and WTADJX can be used to adjust weights to compensate for unit nonresponse. For example, setting L in Equation (3) at 1 and U at its near-infinite default (10²⁰) treats the unit response/nonresponse mechanism as virtually a logistic function of the components of xk with the signs on the regression coefficients reversed.

More generally, let rk = 1 when k is a unit respondent and 0 otherwise. Assume the probability of k responding when sampled is independent of whether any other k’ responds; $p_k=E\left(r_k\right)=1/\alpha (x_k\gamma ;L,U)$ can range from 1/U to 1/L.

The g satisfying the calibration equation:

is a consistent estimator for γ under mild conditions (i.e., its mean squared error tends to zero as the sample and population grow arbitrarily large). In Equation (7), $t_z={\sum }_S{d}_k{z}_k$is a vector of estimated totals based on the full sample before unit nonresponse. For simplicity, we are assuming no item nonresponse.

Equation (6) is called “calibration to the population” and is coded ADJUST = POST in the SUDAAN calibration procedures. Equation (7) is called “calibration to the full sample” as is coded ADJUST = NONRESPONSE. Either way, the estimated population total is

where ${\alpha }_k=\alpha (x_{k}g;L,U)$ is called the adjustment factor for ED k. Note that $w_{k \approx }{d}_k{r}_k/p_k.$

When calibrating to the full sample with Equation (7),

where, for technical reasons explained elsewhere in the literature (e.g., Kott & Liao, 2012), b* is the probability limit of

Unlike b, b* does not depend on which k are in the sample.

In estimating the variance of ty^(w), one cannot treat the square-bracketed term in Equation (8) as a constant because rk is a random variable (when calibrating to the population, rkαk(yk − zkb*) replaces the term in the squared brackets). This can cause difficulty when finite-population correction matters, and it often does in establishment surveys based on stratified simple random samples. The SUDAAN calibration weighting procedures handle this by adding the statement VARNONADJ. This adds ∑S dkrkαk²(1 − 1/αk)(yk − zkb)²(n/nr) to the variance estimator, where n is the original sample size, and nr is the respondent sample size (n/nr is an ad hoc adjustment for replacing b* with b). The addition removes the impact of the original sample finite-population-correction factors (1 − 1/dk, which is a constant within each stratum) on that part of the original sample that did not respond. Note that 1 − 1/αk estimates the probability that sampled ED k is a nonrespondent. It is analogous to the second-stage finite-population-correction factors in a two-stage sample.

Fitting a Logistic Response Model

The SAS data set R contains a variable RESPONDENT that is equal to 1 if the ED is a unit respondent to the survey and 0 otherwise. Assuming the probability that an ED responds to the survey is a logistic function of its frame visits (FRAME_VISITS) and that the estimated number of frame visits can be determined for the full sample before nonresponse, there are (at least) two ways that SAS-callable SUDAAN can estimate the annual number of drug-related ED visits in the United States (and within census regions) in a nearly unbiased fashion. The first way uses PROC RLOGIST, and the second uses PROC WTADJUST. The former fits a weighted maximum-likelihood equation (∑S dk (rk − 1/[1 + exp(-zkg)])zk^T = 0) rather than the calibration equation. The respective codes follow.

PROC RLOGIST DATA = R;

DESIGN = STRWOR VARNONADJ;

NEST STRATUM;

WEIGHT W;

TOTCNT BIG_N;

CLASS REGION;

VAR DR_VISITS;

MODEL RESPONDENT = LOG_FRAME;

RUN;

PROC WTADJUST DATA = R ADJUST = NONRESPONSE;

DESIGN = STRWOR VARNONADJ;

NEST STRATUM;

WEIGHT W;

TOTCNT BIG_N;

CLASS REGION;

VAR DR_VISITS;

LOWERBD 1; CENTER 2;

/* UPPERBD is set at its near-infinite default; dding CENTER 2 assumes virtually the same response function fit by RLOGIST */

MODEL RESPONDENT = LOG_FRAME;

RUN;

Tables 3 and 4 show the key RLOGIST results. A 1 percent increase in ED k’s frame visits results in an estimated 0.27 percent increase (the value in the LOG_FRAME row and Beta Coeff. column) in its odds of response, pk/(1-pk).

Tables 5 and 6 show the analogous WTADJUST results.

WTADJUST models the adjustment factor (αk), which is the inverse of the estimated response probability. The WTADJUST code above estimates the percent decrease in the odds of an ED responding caused by a one percent increase in its frame visits as −0.30, which is not exactly −0.27. Nevertheless, both values are statistically consistent estimates of the same parameter value. The estimated means and totals from the two procedures are likewise not the same. In all cases, using RLOGIST appears to produce smaller standard errors.

Despite this, an advantage of using WTADJUST over RLOGIST is that one can directly output the nonresponse-adjusted weights by adding a statement to WTADJUST like the following:

OUTPUT IDVAR WTFINAL ADJFACTOR/FILENAME=OUT REPLACE;

WTFINAL are the adjusted weights (wk), ADJFACTOR the adjustment factors (rkα(xk^Tg)), and OUT is the data set containing both and other variables listed on a separate IDVAR statement. There is no parallel statement in RLOGIST. Moreover, one cannot bound the probabilities of response by 1/U and 1/L with RLOGIST like one can with WTADJUST.

A third method of estimating drug-related ED visits, under the same logistic response model employs WTADJX. It again features LOG_FRAME as the sole non-intercept in the model statement but adds the statement CALVARS FRAME_VISITS, which means WTADJX attempts to calibrate on frame visits rather than on the log of frame visits (LOG_FRAME) even though response is assumed to be a logistic function of the latter. By calibrating on FRAME_VISITS rather than LOG_FRAME, one is attempting to reduce the size of the terms within the squared brackets of Equation (8) but not necessarily the precision of the regression coefficient.

The code for the WTADJX procedure described above is as follows:

PROC WTADJX DATA = R DESIGN = STRWOR ADJUST = NONRESPONSE VARNONADJ;

NEST STRATUM; WEIGHT W; TOTCNT BIG_N;

CLASS REGION;

VAR DR_VISITS;

LOWERBD 1; CENTER 2;

MODEL RESPONDENT = LOG_FRAME;

CALVARS FRAME_VISITS;

RUN;

Unfortunately, the output table with the “Final Weight Sum Minus Control Totals” column reveals that calibration fails (not shown). When the numbers in that column are not all zeros (or very close to it), then the estimated totals and means are specious.

Simply replacing FRAME_VISITS by Q (which, recall, is FRAME_VISITS/1000) in the CALVARS statement fixes things. That is why we created Q.

Tables 7 and 8 show the key results after the correction:

Observe that the estimated standard error for the US-level estimated total (and mean) is the least across the three methods. At the same time, the p-value of the estimated LOG_FRAME coefficient is the highest across the methods, suggesting that although it is the most efficient method among the three in estimating means and totals, it is least efficient in estimating the response model.

Calibration Weighting When Nonresponse Is a Function of the Survey Variable

One can also use WTADJX, but not RLOGIST, when nonresponse is assumed to be a function of DR_VISITS itself—that is, when nonresponse is not assumed to be missing at random (Kott & Liao, 2017). In the code below, response is assumed to be a logistic function of an intercept and the log of DR_VISITS, which is denoted LOG_DR:

PROC WTADJX DATA = R DESIGN = STRWOR ADJUST = NONRESPONSE VARNONADJ;

NEST STRATUM; WEIGHT W; TOTCNT BIG_N;

CLASS REGION;

VAR DR_VISITS;

LOWERBD 1; CENTER 2;

MODEL RESPONDENT = LOG_DR;

CALVARS Q;

OUTPUT TOTAL SE_TOTAL/FILENAME=OUT0 REPLACE;

RUN;

Rather than printing the total here, the totals and their standard errors under the model have been outputted into OUT0 for future comparison.

As in the previous section, standard errors are likely reduced by calibrating to population totals (for the population size and FRAME_VISITS) rather than full-sample-estimated totals:

PROC WTADJX DATA = R DESIGN = STRWOR ADJUST = POST VARNONADJ;

NEST STRATUM; WEIGHT W; TOTCNT BIG_N;

CLASS REGION; VAR DR_VISITS;

LOWERBD 1; CENTER 2;

MODEL RESPONDENT = LOG_DR;

CALVARS Q;

POSTWGT 6300 134000;

OUTPUT TOTAL SE_TOTAL/FILENAME=OUT1 REPLACE;

RUN;

These totals and their standard errors under the model have been outputted into OUT1.

One can likely reduce the standard errors further by adding more population targets akin to raking to a size variable. At the same time, this adds potential dummies for public, metro, and region to the assumed response model (even when they are not needed in response modeling, which appears to be the case here):

PROC WTADJX DATA = R DESIGN = STRWOR ADJUST = POST VARNONADJ;

NEST STRATUM; WEIGHT W; TOTCNT BIG_N;

CLASS REGION; VAR DR_VISITS;

LOWERBD 1; CENTER 2;

MODEL RESPONDENT = LOG_DR PUBLIC METRO REGION;

CALVARS Q PUBLICQ METROQ REGION*Q;

POSTWGT 6300 134000 58000 44000 22000 43000 33000 36000;

OUTPUT TOTAL SE_TOTAL/FILENAME=OUT2 REPLACE;

OUTPUT ADJFACTOR/FILENAME = AFNR REPLACE;

These totals and their standard errors under the model have been outputted into OUT2.

The addition of the line “OUTPUT ADJFACTOR/FILENAME = AFNR REPLACE;” in the code above creates data set AFNR, with the adjustment factors. A PROC UNIVARIATE of the adjustment factors (not shown) suggests we can set the upper bound at 2.75 for the adjustment factors in the code below with the hope or containing the variability of the weights and thus standard errors:

PROC WTADJX DATA = R DESIGN = STRWOR ADJUST = POST VARNONADJ;

NEST STRATUM; WEIGHT W; TOTCNT BIG_N;

CLASS REGION; VAR DR_VISITS;

LOWERBD 1; CENTER 2; UPPERBD 2.75;

MODEL RESPONDENT = LOG_DR PUBLIC METRO REGION;

CALVARS Q PUBLICQ METROQ REGION*Q;

POSTWGT 6300 134000 58000 44000 22000 43000 33000 36000;

OUTPUT TOTAL SE_TOTAL/FILENAME=OUT3 REPLACE;

/* When running SAS STUDIO, we make sure calibration is successful with the following two optional lines */

OUTPUT DIFFWTZ/FILENAME = DIFF REPLACE;

PROC PRINT DATA = DIFF; RUN; * All DIFFWTZ should be 0;

The new totals and their standard errors under the model have been outputted into OUT3.

The following code compares the various estimates of the totals. Results are shown in Table 9.

DATA OUT0; SET OUT0;

NONCV = SE_TOTAL/TOTAL;

NONTOTAL = TOTAL;

RUN;

DATA OUT1; SET OUT1;

POSTCV = SE_TOTAL/TOTAL;

POSTTOTAL = TOTAL;

RUN;

DATA OUT2; SET OUT2;

RSCV = SE_TOTAL/TOTAL;

RSTOTAL = TOTAL;

RUN;

DATA OUT3; SET OUT3; ;

RS2CV = SE_TOTAL/TOTAL;

RS2TOTAL = TOTAL;

RUN;

DATA C; MERGE OUT0 OUT1 OUT2 OUT3; BY VARIABLE REGION;

NONCV = ROUND(NONCV * 100, .01);

POSTZCV= ROUND(POSTCV * 100, .01);

RSCV= ROUND(RSCV * 100, .01);

RS2CV = ROUND(RS2CV * 100, .01);

PROC PRINT; ID REGION; VAR NONTOTAL POSTTOTAL RSTOTAL RS2TOTAL;

PROC PRINT; ID REGION; VAR NONCV POSTCV RSCV RS2CV; RUN;

where

NONRxx has been computed with ADJUST = RESPONSE and CALVARS Q,

POSTxx with ADJUST = POST and CALVARS Q,

PRSxx with ADJUST = POST and raking to size with default upper bound, and

PRS2xx with ADJUST = POST and raking to size with an upper bound of 2.75.

As expected, calibrating to the population tends to reduce standard errors compared with calibrating to the sample, although not in every region. Adding some raking-to-size variables reduces standard errors even more. Constraining the weight adjustments to be no greater than 2.75 decreases the estimated US-level standard error by around 10 percent but increases the estimated standard error in Region 1 by roughly 35 percent. There are small estimated standard error increases in two of the three other regions as well. This shows that smaller upper bounds (and presumably less-variable calibrated weights) do not necessarily translate into smaller standard errors.

One can check whether the response model with and without the upper bound of 2.75 produce significantly different totals in the following manner. Make two copies of each record in newly created data set R2. Place one copy into DOMAIN = 1 and the other into DOMAIN = 2. DOMAIN = 1 has no upper bound (i.e., the default is used), whereas DOMAIN = 2 has an upper bound of 2.75. (So long as it is finite, the value assigned to CENTER only affects the intercept.) DOMAIN replaces the previously missing intercept in the CALVARS statement, so there is now a NOINT at the end of both the MODEL and CALVARS statements. When comparing estimates, it is common to treat the design as with replacement:

DATA R2; SET R;

PSU + 1;

U = .; DOMAIN = 1; OUTPUT;

U = 2.75; DOMAIN = 2; OUTPUT;

PROC WTADJX DATA = R2 DESIGN = WR;

ADJUST = POST; * One can also use ADJUST= NONRESPONSE for this for test;

NEST STRATUM PSU; WEIGHT W;

CLASS DOMAIN REGION;

VDIFFVAR DOMAIN =(1 2); * This tells SUDAAN to compare estimates between domains;

VAR DR_VISITS;

LOWERBD 1; CENTER 2; UPPERBD U;

MODEL RESPONDENT = DOMAIN LOG_DR*DOMAIN PUBLIC*DOMAIN;

METRO*DOMAIN REGION*DOMAIN/NOINT;

CALVARS DOMAIN Q*DOMAIN PUBLICQ*DOMAIN;

METROQ*DOMAIN REGION*DOMAIN*Q/NOINT;

POSTWGT 6300 6300 134000 134000 58000 58000;

44000 44000 22000 22000 43000 43000 33000 33000 36000 36000;

OUTPUT TOTAL SE_TOTAL /FILENAME=TEST REPLACE;

We compute t-values for the differences, but they are not significant. In fact, all the t-values are well below 1 in absolute value. Similar results (Table 10) obtain when the ADJUST = NONRESPONSE option is used:

DATA TEST; SET TEST;

T_TOTAL = ROUND(TOTAL/SE_TOTAL, .01);

PROC PRINT; ID REGION; VAR TOTAL SE_TOTAL T_TOTAL;

RUN;