H L Test (Hosmer-Lemeshow test ) - H0: the logistic regression fits well

http://www.stat.sc.edu/~hitchcock/516_sas_LOGISTIC_examp.txt



/* Data for the logistic regression example */

/* involving the cities' use of TIF */

/* that we studied in class */



/* Entering the data and naming the variables: */

/* Y is a binary response variable here. */

/* Y = 1 if a city used TIF. */

/* Y = 0 if the city did not use TIF. */



DATA citydata;

INPUT Y $ income;

CARDS;

0 9.2

0 12.9

0 9.2

1 9.6

0 9.3

1 10.1

0 9.4

1 10.3

0 9.5

1 10.9

0 9.5

1 10.9

0 9.5

1 11.1

0 9.6

1 11.1

0 9.7

1 11.1

0 9.7

1 11.5

0 9.8

1 11.8

0 9.8

1 11.9

0 9.9

1 12.1

0 10.5

1 12.2

0 10.5

1 12.5

0 10.9

1 12.6

0 11

1 12.6

0 11.2

1 12.6

0 11.2

1 12.9

0 11.5

1 12.9

0 11.7

1 12.9

0 11.8

1 12.9

0 12.1

1 13.1

0 12.3

1 13.2

0 12.5

1 13.5

;



run;



/* A simple plot of the data set */



goptions reset=all;

PROC GPLOT DATA=citydata;

PLOT Y*income;

run;



/* PROC LOGISTIC will fit a logistic regression model. */

/* The DESCENDING option is important! */

/* It defines mu as P(Y=1) as we did in class, rather than as P(Y=0). */



/* We create an output data set called NEW that contained the predicted probabilities. */

/* It also contains lower and upper bounds for a 95% CI for the true probability. */





PROC LOGISTIC DESCENDING DATA=citydata;

MODEL Y = INCOME / LACKFIT;

OUTPUT OUT=NEW P=PRED L=LOWER U=UPPER;

RUN;



/* Output: With the LACKFIT option, SAS provides a Hosmer-Lemeshow test for */

/* "H0: the logistic regression fits well". */

/* With a P-value of 0.4603, we cannot reject this null hypothesis. */

/* We have no reason to doubt the logistic model's fit. */

/* So it's fine to use the logistic regression model. */



/* The estimates beta_0-hat and beta_1-hat are -11.347 and 1.002. */

/* Estimated odds ratio = 2.723, and 95% CI for odds ratio is (1.526, 4.858). */



/* Output: SAS provides a likelihood-ratio test of H0: beta_1 = 0. Since the P-value */

/* is very small ( < .0001), we reject H0, conclude beta_1 is not zero, and conclude */

/* that income has a significant effect on the probability a city uses TIF. */



/* PLOTTING THE ESTIMATED LOGISTIC CURVE */



/* PROC PLOT gives a crude plot of the estimated logistic regression curve. */

/* The 95% CIs are also plotted. */



PROC PLOT DATA=NEW;

PLOT PRED*INCOME='P' LOWER*INCOME='L' UPPER*INCOME='U' / OVERLAY;

RUN;





/* Using PROC GPLOT gives a nicer-looking plot of the estimated logistic regression curve */

/* and the 95% CIs. */

/* This plot comes up in a separate window from the OUTPUT window. */

/* First the data must be sorted based on the X variable before using PROC GPLOT. */

/* (since we are going to connect the dots). */



/* We use a circle as the plotting symbol for the first curve (the predicted probabilities) */

/* and a star as the plotting symbol for the 2nd and 3rd curves (the CIs). */



PROC SORT DATA=NEW;

BY INCOME;

symbol1 i = join v=circle l=32 c = black;

symbol2 i = join v=star l=32 c = black;

symbol3 i = join v=star l=32 c = black;

PROC GPLOT DATA=NEW;

PLOT PRED*INCOME LOWER*INCOME UPPER*INCOME / OVERLAY;

RUN;

BASEL II and PD EAD LGD

Help banks to satisfy the three pillars of Basel II.




Pillar I: Minimum capital requirement:-

Through comprehensive end to end risk management system; banks would be able to allocate the required minimum capital for each risk



Pillar II: Supervisory review process:-

Real time monitoring of clients credit position and banks financial position through “Risk Analysis Engine”



Pillar III: Market discipline requirements-

Improved transparency, sound financial system, effective risk management and mitigation process, and disclosure of bank’s financial stability to public, through “Risk Analysis Dashboard”







PD - Probability of Default - The likelihood of a borrower to become default over a performance window after the observation point.

EAD - Exposure at Default - The outstanding balance at the time when the borrower becomes default.

LGD - Loss Given Default - The loss proportion of an exposure, e.g. EAD, after a recovery window from the time of default.



PD estimates are 1 year forward-looking probabilities of default (Default – fails to repay borrowings)

EAD estimates are a long-run default weighted average EAD; and

LGD estimates reflect economic downturn conditions

http://www.manyppt.com/05/The-New-Basel-Capital-Accord-and-Questions-for-Practice-and.html





Value at Risk (VAR) is used to measure the market risk .

VAR summarizes the likely loss in value of a portfolio over a given time period with specified probability.

Historical simulation, Model building approach, Montey Carlo simulation – are some of the VAR techniques.

Compute biserial, point biserial, and rank biserial correlations

Source -


http://support.sas.com/kb/24/991.html



%macro biserial(version, data= ,contin= ,binary= ,out=);




%if &version ne %then %put BISERIAL macro Version 2.2;



options nonotes;

* exclude observations with missing variables *;

data &out;

set &data;

where &contin>.;

if &binary>.;

run;



* compute the ranks for the continuous variable *;

proc rank data=&out out=&out ;

var &contin;

ranks r_contin;

run;



* compute proportion of binary, std of contin, and n *;

proc means data=&out noprint;

var &binary &contin;

output out=_temp_(keep=p stdy n) mean=p std=stdx stdy n=n;

run;



* sort by the binary variable *;

proc sort data=&out;

by descending &binary;

run;



* compute mean of contin and rank of contin var *;

proc means data=&out noprint;

by notsorted &binary;

var &contin r_contin;

output out=&out mean=my r_contin;

run;



* restructure the means computed in the step above *;

proc transpose data=&out out=&out(rename=(col1=my1 col2=my0));

var r_contin my;

run;



* combine the data needed to compute biserial correlation *;

data &out;

set &out(drop= _name_ _label_);

retain r1 r0 ;

if _n_=1 then do;

r1=my1;

r0=my0;

end;

else do;

set _temp_;

output;

end;

run;



* compute point biserial correlation *;

proc corr data=&data noprint outp=_temp_;

var &binary &contin;

run;







* extract the point biserial correlation from the matrix *;

data _temp_(keep=pntbisrl);

set _temp_(rename=(&contin=pntbisrl));

if _TYPE_='CORR' and &binary<>1 then output;



run;



options notes;

* compute biserial and rank biserial *;

data &out;

merge _temp_ &out;

if pntbisrl=1 then delete;

h=probit(1-p);

u=exp(-h*h/2)/sqrt(2*arcos(-1));

biserial=p*(1-p)*(my1-my0)/stdy/u;

rnkbisrl=2*(r1-r0)/n;



keep biserial pntbisrl rnkbisrl;

label biserial='Biserial Corr'

pntbisrl='Point Biserial Corr'

rnkbisrl='Rank Biserial Corr';

run;



%mend;







data k;

length x1 $ 1;

input x1 length;

event=(x1='y');

cards;

y 14.8

n 13.8

y 12.4

y 10.1

y 7.1

y 6.1

n 5.8

y 4.6

n 4.3

n 3.5

n 3.3

y 3.2

y 3.0

n 2.8

n 2.8

n 2.5

y 2.4

y 2.3

y 2.1

n 1.7

n 1.7

n 1.5

n 1.3

n 1.3

n 1.2

n 1.2

n 1.1

y 0.8

n 0.7

n 0.6

n 0.5

n 0.2

n 0.2

y 0.1

;



/* Define the BISERIAL macro */

%inc "";



%biserial(data=k, contin=LENGTH, binary=EVENT, out=out1);

*********************
data= SAS data set to be analyzed.
binary =  Name of dichotomous variable which must be numeric with values 0 and 1.
contin=  Name of continuous variable. Ranks of this variable will be computed to produce the rank biserial corr.
out= Output data set name.
*****************



proc print data=out1 label noobs;

title 'Point Biserial, Biserial and Rank Biserial Correlations';

run;

Logistic Regression - Assumptions

Residuals follow a binomial rather than a normal distribution. Normality of variables is not a stringent requirement.

A “nonlinear” (specifically an S-shaped, or sigmoidal, curve) relationship between IVs and the DV; however, this represents a linear relationship between the logit (natural log of the odds of the dependent occurring or not) and the set of IVs. See Addendum 2 for an illustration that compares probabilities, odds, and the logit.

Uses a maximum-likelihood rather than least-squares statistical model. In least squares, we select
regression coefficients that result in the smallest sum of squared differences between the observed
and the predicted values of the DV. In maximum-likelihood, the coefficients that make our observed results “most likely” are selected.

Does not assume homoscedasticity.

Assumes that there is little or no multicollinearity

Predicts the odds of an event occurring (see Addendum 1), which is based on the probability of

that event occurring. Precisely, the odds of an event occurring is: Odds=P/(1-P)
={prob. of event occurring}/[1-prob. of event occurring]
--------------------------------------------


Terminology for use with Logistic Regression -

Probability = P = probability of an event occurring (range of 0 - 1)


 Odds = P/(1-P) = ratio of the probability of an event occurring to the probability of the event not occurring (range of 0 – positive infinity)

Odds ratio = Odds1/ Odds2= ratio of two odds

Logit = ln(Odds) = predicted logged odds (range of neg. infinity - pos. infinity)


******************


References:


Hosmer, D. W., & Lemeshow, S. (2000). Applied logistic regression (2nd ed.). New York: John Wiley &

Sons, Inc.

Menard, S. (1995). Applied logistic regression analysis. Thousand Oaks, CA: Sage Publications.

Pampel, F. C. (2000). Logistic regression: A primer. Thousand Oaks, CA: Sage Publications.

Quantiles into new dataset

data tt;


do i=1 to 100;

x1=rannor(0);output;

x2=rannor(0);output;

end;

run;





ODS TRACE on / listing;



proc univariate data=tt;

var x1 x2;

run;



ODS TRACE off;







proc univariate data=tt;

var x1 x2;

ODS OUTPUT Quantiles=qntls;

run;

http://rfd.uoregon.edu/files/rfd/StatisticalResources/110_ods.txt

C - Stat Calculation Code for Validation and OOT

/* STEP - 1 Modeling for Development Dataset */



proc logistic data=DEV outmodel=logist_param descending;

model depvar=&var_r. /

outroc=roc

stb

scale=none

lackfit

selection=stepwise sle=.05 sls=.001;

output out=est_dev P=score;

run;



/* KS */

proc npar1way edf data=devscr;

class depvar;

var score;

run;





/* STEP - 2 Modeling for Validation Dataset */



proc logistic descending inmodel =logist_param;

score data = VAL out = est_val;

run;



/* KS */

proc npar1way data = est_val edf;

class depvar;

var P_1;

run;



/* STEP - 3 Modeling for OOT Validation Dataset */

proc logistic descending inmodel =logist_param;

score data = OOT out = est_val_OOT;

run;



/* KS */

ODS GRAPHICS on;

proc npar1way data = est_val_OOT edf;

class depvar;

var P_1;

run;





/* STEP - 4 C Statistic for Validation and Out of Time Sample DataSet */



%macro C_stat(input,var1,var2,output);



data d nd (rename=(&var1.=&var1._2 &var2.=&var2._2));

set &input.;

if &var1. then output d; else output nd;

run;



proc sql;

create table pairs as select a.&var1.,a.&var2.,b.&var1._2,b.&var2._2

from d as a, nd as b;

quit;



data class;

set pairs;

if &var2. gt &var2._2 then concordant =1; else concordant =0;

if &var2. = &var2._2 then tied =1; else tied =0;

if &var2. lt &var2._2 then discordant =1; else discordant =0;



proc summary data = class;

output out = sum

sum(concordant)=concordant

sum(tied)=tied;

run;



data &output.;set sum;

C = (concordant+tied)/_freq_;

run;

%mend;





%C_stat(est_val,depvar,p_1,C_for_val_data);

%C_stat(est_val_OOT,depvar,p_1,C_for_oot_data);

C - Stat Calculation Code

data test(drop=i);


do i=1 to 200;

y=round(ranuni(9),1);

x=ranuni(3);output;

end;

run;



proc logistic data=test outmodel=logist_param descending;

model y=x;

output out=est_dev P=score;

run;





%macro C_stat(input,var1,var2,output);



data d nd (rename=(&var1.=&var1._2 &var2.=&var2._2));

set &input.;

if &var1. then output d; else output nd;

run;



proc sql;

create table pairs as select a.&var1.,a.&var2.,b.&var1._2,b.&var2._2

from d as a, nd as b;

quit;



data class;

set pairs;

if &var2. gt &var2._2 then concordant =1; else concordant =0;

if &var2. = &var2._2 then tied =1; else tied =0;

if &var2. lt &var2._2 then discordant =1; else discordant =0;



proc summary data = class;

output out = sum

sum(concordant)=concordant

sum(tied)=tied;

run;



data &output.;set sum;

C = (concordant+tied)/_freq_;

run;

%mend;





%C_stat(est_dev,y,score,out_c);