Introduction
This situation is also called oversampling, retrospective sampling, biased sampling, or choice-based sampling. In the biological sciences, studies using this kind of sampling are known as case-control studies.
Parameter and odds ratio estimates of the covariates (and their confidence limits) are unaffected by this type of stratified sampling, so no weighting is needed. However, the intercept estimate is affected by the sampling, so any computation that is based on the full set of parameter estimates is incorrect, such as the predicted event probabilities, differences or ratios of event probabilities (the ratio is called the relative risk), and false positive and negative rates. If you know the probabilities of events and nonevents in the population, then you can adjust the intercept either by weighting or by using an offset.
Adjusting the Intercept
To adjust by weighting, add a variable to your data set that takes the value p1/r1 in event observations, and the value (1-p1)/(1-r1) in nonevent observations, where p1 is the probability of an event in the population and r1 is the proportion of events in your data set. Specify this variable in the WEIGHT statement in PROC LOGISTIC.
Or, to adjust by using an offset, add a variable to your data set defined as log[(r1*(1-p1)) / ((1-r1)*p1)], where log represents the natural logarithm. Specify this variable in the OFFSET= option of the MODEL statement in PROC LOGISTIC.
These two methods are not equivalent. The weighting method seems to perform better when the model is not correctly specified. For details, see Scott and Wild (1986), "Fitting Logistic Models under Case-Control or Choice-Based Sampling," Journal of the Royal Statistical Society B, 48.
Adjusting the Predicted Probabilities
The weighting method also has the advantage of directly providing properly adjusted predicted probabilities via the PREDICTED= or PREDPROBS= option in the OUTPUT statement. You can use the PRIOR= or PRIOREVENT= option in the SCORE statement to input the population event and nonevent probabilities and obtain predicted probabilities that are equivalent to those produced by the offset method. These methods are illustrated below.
For more information about analyzing data from stratified sampling, see Hosmer and Lemeshow (2000) or Collett (2003).
Example
The following example illustrates obtaining predicted probabilities adjusted for oversampling. Data set FULL is created containing a binary response, Y (with event=1 and nonevent=0), and predictor, X. The true model from which the data is generated is logit(p) = log(p/(1-p)) = -3.35 + 2*X, resulting in approximately a 0.1 overall proportion of events. Note that the LOGISTIC function computes the inverse of logit(p): p = logistic(logit(p)) = 1/(1+exp(-logit(p))). PROC FREQ shows the proportions of events and nonevents.
data full;
do i=1 to 1000;
x=rannor(12342);
p=logistic(-3.35+2*x);
y=ranbin(98435,1,p);
drop i;
output;
end;
run;
proc freq data=full;
table y / out=fullpct(where=(y=1)
rename=(percent=fullpct));
title "Response counts in FULL data set";
run;
A subset data set, SUB, of data set FULL is obtained by oversampling. That is, all Y=1 observations and one-ninth of the Y=0 observations are retained resulting in a sample with approximately equal numbers of events and nonevents.
data sub;
set full;
if y=1 or (y=0 and ranuni(75302)<1/9) then output;
run;
proc freq data=sub;
table y / out=subpct(where=(y=1)
rename=(percent=subpct));
title "Response counts in oversampled, subset data set";
run;
A subset sample created by oversampling can also be produced using PROC HPSAMPLE. The following statements generate a similar, but not identical, subset. Option EVENT="1" specifies that the value 1 in the TARGET variable, Y, is considered the event level. The SAMPPCTEVT=100 option requests that 100 percent of the Y=1 observations be selected. EVENTPROP=0.5 causes selection of enough nonevents to result in the final sample containing 50% events – that is, equal numbers of events and nonevents.
proc hpsample data=full out=sub2 seed=13579
samppctevt=100 eventprop=0.5 event="1";
var x; class y; target y;
run;
Using subset data, SUB, the offset (OFF) and weights (W) are computed as discussed above using the known population and sample event probabilities. Unadjusted, offset-adjusted, and weight-adjusted logistic models are then fit yielding corrected intercepts. Predicted probabilities from each of the methods are computed as discussed above. All results are accumulated in data set OUT.
data sub;
set sub;
if _n_=1 then set fullpct(keep=fullpct);
if _n_=1 then set subpct(keep=subpct);
p1=fullpct/100; r1=subpct/100;
w=p1/r1; if y=0 then w=(1-p1)/(1-r1);
off=log( (r1*(1-p1)) / ((1-r1)*p1) );
run;
This LOGISTIC step fits the unadjusted logistic model to the subset data and saves the predicted probabilities. This model does not correct the intercept for the oversampling.
proc logistic data=sub;
model y(event="1")=x;
output out=out p=pnowt;
title "True Parameters: -3.35 (intercept), 2 (X)";
title2 "Unadjusted Model";
run;
This next LOGISTIC step fits the logistic model using the weight adjustment method and saves the adjusted predicted probabilities.
proc logistic data=out;
model y(event="1")=x; weight w;
output out=out p=pwt;
title2 "Weight-adjusted Model";
run;
This LOGISTIC step uses the offset adjustment method to fit the logistic model and saves the adjusted predicted probabilities.
proc logistic data=out;
model y(event="1")=x / offset=off;
output out=out xbeta=xboff;
title2 "Offset-adjusted Model";
run;
data out;
set out;
poff=logistic(xboff-off);
run;
These statements refit the unadjusted model and uses the proportions in the full data set as prior probabilities to compute and save the prior-adjusted probabilities.
proc freq data=full noprint;
table y / out=priors(drop=percent
rename=(count=_prior_));
run;
proc logistic data=out;
model y(event="1")=x;
score data=sub prior=priors out=out2;
title2 "Unadjusted Model; Prior-adjusted probabilities";
run;
data out;
merge out out2;
drop _lev:;
run;
Finally, these statements produce a graph comparing the true probabilities and all of the adjusted probabilities using PROC SGPLOT. NOTE
proc sort data=out;
by x;
run;
proc sgplot data=out;
title "Adjustments to Probabilities";
series y=p x=x / name="p" legendlabel="True"
lineattrs=(color=red pattern=solid);
series y=pnowt x=x / name="pnowt" legendlabel="Unadjusted"
lineattrs=(color=black pattern=solid);
series y=pwt x=x / name="pwt" legendlabel="Weighted"
lineattrs=(color=green pattern=solid);
series y=poff x=x / name="poff" legendlabel="Offset"
lineattrs=(color=brown pattern=solid);
series y=p_1 x=x / name="p_1" legendlabel="Priors"
lineattrs=(color=yellow pattern=solid);
keylegend "p" "pnowt" "pwt" "poff" "p_1" /
position=topleft across=3 down=3 border
title="Adjustments" location=inside;
run;
Note that the offset-adjusted probabilities match those from using the PRIOR= option in the SCORE statement. Both the offset and the weight adjustments come closer to the true probabilities than the unadjusted probabilities which are incorrect because of the unadjusted intercept. The offset adjustment seems slightly better than the weight adjustment in this case, possibly because the correct model is specified. With a larger initial data set (say, 10000) and sample (1000 events and nonevents), both adjustment methods yield very similar probabilities very close to the true probabilities.
__________
NOTE: These statements produce the same graph using a custom graph template.
proc template;
define statgraph LOGAdj;
begingraph;
entrytitle "Adjustments To Probabilities";
notes "plot of various adjustments for oversampling";
layout overlay;
seriesplot y=p x=x / name="p"
legendlabel="True" lineattrs=(color=red);
seriesplot y=pnowt x=x / name="pnowt"
legendlabel="Unadjusted" lineattrs=(color=black);
seriesplot y=pwt x=x / name="pwt"
legendlabel="Weighted" lineattrs=(color=green);
seriesplot y=poff x=x / name="poff"
legendlabel="Offset" lineattrs=(color=brown);
seriesplot y=p_1 x=x / name="p_1"
legendlabel="Priors" lineattrs=(color=yellow);
discretelegend "p" "pnowt" "pwt" "poff" "p_1" /
halign=left valign=top across=3 down=3 border=true
title="Adjustments" titleborder=false;
endlayout;
endgraph;
end;
run;
proc sgrender data=out template=LOGAdj;
run;