LOGSELECT Procedure

(View the complete code for this example.)

If are independent binary (Bernoulli) random variables with common success probability , then their sum is a binomial random variable. In other words, a binomial random variable with parameters n and can be generated as the sum of n Bernoulli() random experiments. The LOGSELECT procedure uses a special syntax, the events/trials syntax, to express data in binomial form.

Consider the following data, taken from Cox and Snell (1989, pp. 10–11), of the number of ingots not ready for rolling (r) out of n tested for a number of combinations of heating time and soaking time. If each test is carried out independently and if for a particular combination of heating and soaking time there is a constant probability that the tested ingot is not ready for rolling, then the random variable r follows a binomial distribution, where the success probability is a function of heating and soaking time.

data Ingots;
   input Heat Soak r n @@;
   Obsnum= _n_;
   datalines;
7 1.0 0 10  14 1.0 0 31  27 1.0 1 56  51 1.0 3 13
7 1.7 0 17  14 1.7 0 43  27 1.7 4 44  51 1.7 0  1
7 2.2 0  7  14 2.2 2 33  27 2.2 0 21  51 2.2 0  1
7 2.8 0 12  14 2.8 0 31  27 2.8 1 22  51 4.0 0  1
7 4.0 0  9  14 4.0 0 19  27 4.0 1 16
;
data mylib.Ingots;
   set Ingots;
run;

The following statements show the use of the events/trials syntax to model the binomial response. The events variable in this situation is r (the number of ingots not ready for rolling), and the trials variable is n (the number of ingots that are tested). The dependency of the probability of not being ready for rolling is modeled as a function of heating time, soaking time, and their interaction. The ASSOCIATION option displays ordinal measures of association between the observed responses and the predicted probabilities. The CTABLE option saves statistics to the mylib.Roc data table that are used to evaluate the predictive power of the model. The OUTPUT statement stores the linear predictors and the predicted probabilities in the mylib.Out data table along with the other variables from the input data table.

proc logselect data=mylib.Ingots association ctable(out=mylib.Roc nocounts tpf fpf);
   model r/n = Heat Soak Heat*Soak;
   output out=mylib.Out xbeta predicted=Pred copyvars=(Heat Soak);
run;

The "Model Information" table shows that the data are modeled as binomially distributed with a logit link function (Output 18.2.1). This is the default link function in the LOGSELECT procedure for binary and binomial data. The procedure estimates the parameters of the model by a Newton-Raphson algorithm.

The second table in Output 18.2.1 shows that all 19 observations in the input data table were used in the analysis and the "Response Profile" table shows that the total number of events and nonevents equals 12 and 387, respectively.

Output 18.2.2 displays the convergence status table for this run. The LOGSELECT procedure satisfies the GCONV= convergence criterion.

Output 18.2.3 displays the "Dimensions" table for the model. The design matrix of the model (the matrix) has four columns, which correspond to the intercept, the Heat effect, the Soak effect, and the interaction of the Heat and Soak effects. The model is nonsingular because the rank of the crossproducts matrix equals the number of columns in . All parameters are estimable, and all participate in the optimization.

Output 18.2.4 displays the "Fit Statistics" table for this run. Evaluated at the converged estimates, –2 times the value of the log-likelihood function equals 27.9569. Further fit statistics are also displayed, all of them in "smaller is better" form. The AIC, AICC, and SBC criteria are used to compare non-nested models and to penalize the model fit for the number of observations and parameters. You can use the –2 log-likelihood value to compare nested models by way of a likelihood ratio test.

Output 18.2.5 shows the test of the global hypothesis that the effects jointly do not affect the probability of ingot readiness. You can obtain the chi-square test statistic by comparing the –2 log-likelihood value of the model with covariates to the value in the intercept-only model. The test is significant, with a p-value of 0.0082. One or more of the effects in the model have a significant impact on the probability of ingot readiness.

The "Parameter Estimates" table in Output 18.2.6 displays the estimates and standard errors of the model effects.

Output 18.2.7 displays the "Association Statistics" table, which is produced when you specify the ASSOCIATION option in the PROC LOGSELECT statement. The table contains four measures of association for assessing the predictive ability of a model. For more information, see the section Association Statistics.

You can construct the prediction equation of the model from the "Parameter Estimates" table. For example, an observation with Heat equal to 14 and Soak equal to 1.7 has the linear predictor

The probability that an ingot with these characteristics is not ready for rolling is

PROC LOGSELECT computes these linear predictors and probabilities and stores them in the mylib.Out data table. The following PROC PRINT statement produces Output 18.2.8, which shows the probability that an ingot with Heat equal to 14 and Soak equal to 1.7 is not ready for rolling:

proc print data=mylib.Out;
   where Heat=14 & Soak=1.7;
run;

The CTABLE option computes statistics for binary response models (on the basis of classifying observations according to whether their predicted probabilities exceed certain values) and stores the results in the mylib.Roc data table. For more information, see the section Classification Table and ROC Curves. You can use this data table to display the ROC curve by using the SGPLOT procedure as follows:

ods graphics on;
proc sgplot data=mylib.Roc aspect=1 noautolegend;
   title 'ROC Curve';
   xaxis values=(0 to 1 by 0.25) grid offsetmin=.05 offsetmax=.05;
   yaxis values=(0 to 1 by 0.25) grid offsetmin=.05 offsetmax=.05;
   lineparm x=0 y=0 slope=1 / lineattrs=(color=ligr);
   series x=FPF y=TPF;
   inset 'Area under the curve=0.7706' / position=bottomright;
run;

Output 18.2.9: Receiver Operating Characteristics Curve

data mylib.Ingots2;
   set Ingots;
   a = n - r;
run;

proc logselect data=mylib.Ingots2 association ctable(out=mylib.Roc nocounts tpf fpf);
   model r a = Heat Soak Heat*Soak;
   output out=mylib.Out xbeta predicted=Pred copyvars=(Heat Soak);
run;

data mylib.Ingots_binary;
   set Ingots;
   do i=1 to n;
     if i <= r then y=1; else y = 0;
     output;
   end;
run;

proc logselect data=mylib.Ingots_binary;
   model y(event='1') = Heat Soak Heat*Soak;
run;