SAS-ohjelmia
Logistinen regressio/ kevät 2001

libname LogReg  'c:\kurssit\LogReg01\DataPA';
data LogReg.penalty;
infile 'c:\kurssit\LogReg01\DataPA\penalty.txt';
input death blackd whitvic serious culp serious2;
run;

PROC REG DATA=LogReg.penalty;
  MODEL death=blackd whitvic serious;
RUN;

PROC FREQ DATA=LogReg.penalty;
  TABLES death*blackd / NOCOL NOROW NOPERCENT;
RUN;

PROC LOGISTIC DATA=LogReg.penalty DESCENDING;
  MODEL death=blackd whitvic serious;
RUN;

PROC GENMOD DATA=LogReg.penalty;
  MODEL death=blackd whitvic serious / DIST=BINOMIAL;
RUN;

PROC GENMOD DATA=LogReg.penalty;
  CLASS culp;
  MODEL death = blackd whitvic culp / D=B;
RUN;

PROC GENMOD DATA=LogReg.penalty;
  MODEL death = blackd whitvic culp blackd*whitvic / D=B;
RUN;

PROC GENMOD DATA=LogReg.penalty;
  MODEL death = blackd whitvic culp blackd*whitvic / D=B TYPE3;
RUN;

PROC LOGISTIC DATA=LogReg.penalty DES;
  MODEL death = blackd whitvic culp / CLPARM=BOTH;
RUN;

PROC LOGISTIC DATA=LogReg.penalty DES;
  MODEL death = blackd whitvic culp / ALPHA=.10 CLPARM=BOTH;
RUN;

PROC LOGISTIC DATA=LogReg.penalty DES;
  MODEL death = blackd whitvic culp / CLODDS=WALD CLODDS=PL;
RUN;

PROC LOGISTIC DATA=LogReg.penalty DES;
  UNITS culp=2 / default=1;
  MODEL death = blackd whitvic culp / CLODDS=WALD CLODDS=PL;
RUN;

PROC GENMOD DATA=LogReg.penalty;
  MODEL death = blackd whitvic culp / D=B WALDCI LRCI;
RUN;

DATA compsep;
  INPUT x y;
  DATALINES;
1 0
2 0
3 0
4 1
5 1
6 1
;
PROC LOGISTIC DES;
  MODEL y = x;
RUN;

PROC GENMOD DATA=compsep;
 MODEL y = x / D=B;
RUN;

DATA quasisep;
  INPUT x y;
  DATALINES;
1 0
2 0
3 0
4 0
4 1
5 1
6 1
;
PROC LOGISTIC DES;
  MODEL y = x;
RUN;
PROC GENMOD;
  MODEL y = x /D=B;
RUN;

PROC GENMOD DATA=LogReg.penalty;
  WHERE blackd=0;
  CLASS culp;
  MODEL death = culp serious / D=B TYPE3;
RUN;

PROC FREQ DATA=LogReg.penalty;
  WHERE blackd=0;
  TABLES culp*death / NOCOL NOROW NOPERCENT;
RUN;

DATA;
  SET LogReg.penalty;
  IF culp=1 THEN culp=2;
PROC GENMOD;
  WHERE blackd=0;
  CLASS culp;
  MODEL death=culp serious /D=B;
RUN;

PROC GENMOD DATA=LogReg.penalty;
  WHERE blackd=0 AND culp > 1;
  CLASS culp;
  MODEL death=culp serious /D=B;
RUN;

PROC GENMOD;
  WHERE blackd=0;
  CLASS culp;
  MODEL death=culp serious /D=B TYPE3;
RUN;

PROC GENMOD DATA=LogReg.penalty;
  WHERE blackd=0;
  CLASS culp;
  MODEL death=culp serious /D=B;
  CONTRAST 'culp1 v. culp5' CULP 1 0 0 0 -1;
RUN;

PROC LOGISTIC DES DATA=LogReg.penalty;
  MODEL death = blackd whitvic serious serious2;
RUN;

PROC REG DATA=LogReg.penalty;
  MODEL death = blackd whitvic serious serious2 / TOL VIF; 
RUN;

PROC LOGISTIC DES DATA=LogReg.penalty;
  MODEL death = blackd whitvic serious serious2;
  OUTPUT OUT=a PRED=phat;
DATA b;
  SET a;
  w = phat*(1-phat);
PROC REG DATA=b;
  WEIGHT w;
  MODEL death = blackd whitvic serious serious2 / TOL VIF; 
RUN;

PROC LOGISTIC DES DATA=LogReg.penalty;
  MODEL death = blackd whitvic culp / AGGREGATE SCALE=NONE;
RUN;

PROC GENMOD DATA=LogReg.penalty;
CLASS culp;
MODEL death = blackd whitvic culp blackd*whitvic blackd*culp
              whitvic*culp blackd*whitvic*culp / D=B;
RUN;

PROC GENMOD DATA=LogReg.penalty;
MODEL death = blackd whitvic culp / D=B;
RUN;

PROC LOGISTIC DES DATA=LogReg.penalty;
  MODEL death = blackd whitvic culp / LACKFIT;
RUN;

PROC LOGISTIC DES DATA=LogReg.penalty;
  MODEL death = culp whitvic blackd / RSQ;
RUN;

PROC GENMOD DATA=LogReg.penalty;
  MODEL death = culp blackd whitvic / D=B OBSTATS RESIDUALS;
  MAKE 'OBSTATS' OUT=a;
RUN;

PROC LOGISTIC DES DATA=LogReg.penalty;
  MODEL death = culp whitvic blackd / INFLUENCE;
RUN;

PROC LOGISTIC DES DATA=LogReg.penalty;
  MODEL death = culp blackd whitvic;
  OUTPUT OUT=a(KEEP= death pred up lo chi dev) P=pred U=up L=lo 
     RESCHI=chi RESDEV=dev;
  OUTPUT OUT=b(KEEP= death pred dculp dblack dwhitvic dev) P=pred 
     DFBETAS=int dculp dblack dwhitvic DIFDEV=dev;
RUN;

PROC GPLOT DATA=b;
  PLOT dev*pred / VAXIS=AXIS1 HAXIS=AXIS1;
  SYMBOL V=DOT HEIGHT=.35;
  AXIS1 MINOR=NONE WIDTH=2 MAJOR=(WIDTH=2);
RUN;

PROC LOGISTIC DES DATA=LogReg.penalty;
  MODEL death = culp blackd whitvic / LINK=PROBIT;
RUN;

PROC LOGISTIC DES DATA=LogReg.penalty;
  MODEL death = culp blackd whitvic / LINK=CLOGLOG;
RUN;

/*	******************	*/
/* 4. Luku	*/
/* 4.2 2x2 taulukkojen logitanalyysi */

PROC GENMOD DATA=LogReg.penalty;
  MODEL death = blackd / D=B;
RUN;


DATA tab4_1a;
  INPUT f blackd death;
  DATALINES;
22 0 1
28 1 1
52 0 0
45 1 0
;
PROC GENMOD DATA=tab4_1a;
  FREQ f;
  MODEL death = blackd / D=B;
RUN;



DATA tab4_1b;
  INPUT death total blackd;
  DATALINES;
22 74 0
28 73 1
;
PROC GENMOD DATA=tab4_1b;
  MODEL death/total = blackd / D=B;
RUN;


PROC LOGISTIC DATA=tab4_1b;
  MODEL death/total = blackd / AGGREGATE SCALE=N;
RUN;


PROC FREQ DATA=tab4_1a;
  WEIGHT f;
  TABLES blackd*death / CHISQ CMH;
RUN;

/* 4.3 3-ulotteiset taulukot  */

DATA interco;
  INPUT white male yes no;
  total=yes+no;
  DATALINES;
1 1 43 134
1 0 26 149
0 1 29 23
0 0 22 36
;
PROC GENMOD DATA=interco;
  MODEL yes/total=white male / D=B;
RUN;

DATA;
  CHI=1-PROBCHI(.0583,1);
  PUT CHI;
RUN;

PROC GENMOD DATA=interco;
  MODEL yes/total=white male white*male/ D=B;
RUN;

/* 4.3 4-ulotteiset taulukot  */

DATA working;
  INPUT france manual famanual total working;
  DATALINES;
1   1   1  107  85
1   1   0   65  44
1   0   1   66  24
1   0   0  171  17
0   1   1   87  24
0   1   0   65  22
0   0   1   85   1
0   0   0  148   6
;
PROC GENMOD DATA=working;
  MODEL working/total = france manual famanual / D=B;
RUN;

PROC GENMOD DATA=working;
  MODEL working/total = france manual famanual 
    france*manual france*famanual manual*famanual / D=B;
RUN;

DATA wisc;
  INPUT iq parent ses total coll;
  DATALINES;    
1   1   1   353     4
1   1   2   234     2
1   1   3   174     8
1   1   4    52     4
1   2   1    77    13
1   2   2   111    27
1   2   3   138    47
1   2   4    96    39
2   1   1   216     9
2   1   2   208     7
2   1   3   126     6
2   1   4    52     5
2   2   1   105    33
2   2   2   159    64
2   2   3   184    74
2   2   4   213   123
3   1   1   138    12
3   1   2   127    12
3   1   3   109    17
3   1   4    50     9
3   2   1    92    38
3   2   2   185    93
3   2   3   248   148
3   2   4   289   224
4   1   1    77    10
4   1   2    96    17
4   1   3    48     6
4   1   4    25     8
4   2   1    92    49
4   2   2   178   119
4   2   3   271   198
4   2   4   468   414
; 
PROC GENMOD DATA=wisc;
  CLASS iq ses;
  MODEL coll/total=iq ses parent / D=B TYPE3;
RUN;

PROC GENMOD DATA=wisc;
  CLASS iq ses;
  MODEL coll/total=iq|ses|parent @2/ D=B TYPE3;
RUN;

PROC GENMOD DATA=wisc;
  CLASS iq ses;
  MODEL coll/total=iq ses parent / D=B TYPE3;
  CONTRAST 'Melko alh. vs. melko kork.' ses 0 1 -1 0;
RUN;

/* 4.6 Ylihajonta */

DATA LogReg.nihdoc;
  INPUT nih docs pdoc;
  DATALINES; 
.5 8 1
.5 9 3
.835 16 1
.998 13 6
1.027 8 2
2.036 9 2
2.106 29 10
2.329 5 2
2.523 7 5
2.524 8 4
2.874 7 4
3.898 7 5
4.118 10 4
4.130 5 1
4.145 6 3
4.242 7 2
4.280 9 4
4.524 6 1
4.858 5 2
4.893 7 2
4.944 5 4
5.279 5 1
5.548 6 3
5.974 5 4
6.733 6 5
7 12 5
9.115 6 2
9.684 5 3
12.154 8 5
13.059 5 3
13.111 10 8
13.197 7 4
13.433 86 33
13.749 12 7
14.367 29 21
14.698 19 5
15.440 10 6
17.417 10 8
18.635 14 9
21.524 18 16
;
PROC GENMOD DATA=LogReg.nihdoc;
  MODEL pdoc/docs=nih / D=B;
RUN;