/* Esimerkkiohjelmia - SAS */
/* 2-ulotteisen taulukon rakenne */

data Avadex;
	input kasvain $ altistus $ lkm;
	cards;
	on on 4  
	on ei 5 
	ei on 12
	ei ei 74
	;
proc freq order=data;
	weight  lkm;
	tables kasvain*altistus/nopercent norow nocol;
run;
proc freq order=data;
	weight  lkm;
	tables kasvain*altistus/nopercent norow;
run;

/*tables kasvain*altistus/chiq expected measures;*/

libname Fraadata  'c:\kurssit\Fraa\Fraa01\Datat';
data Fraadata.Avadex;
infile 'c:\kurssit\Fraa\Fraa01\Datat\Avadex.txt';
input kasvain $ altistus $ lkm;
run;
proc freq order=data;
	weight  lkm;
	tables kasvain*altistus/nopercent norow nocol;
run;

/* 2.2.1 Framinghamin sydäntutkimus */

data Framing;
input Vp $ Tauti $ lkm @@;  
cards;
_116    On  3 _116    Ei 153
117_126 On 17 117_126 Ei 235
127_136 On 12 127_136 Ei 272
137_146 On 16 137_146 Ei 255
147_156 On 12 147_156 Ei 127
157_166 On  8 157_166 Ei 77
167_186 On 16 167_186 Ei 83
187_    On  8 187_    Ei 35
;
proc freq order=data;
	weight  lkm;
	tables Tauti*Vp/nopercent norow nocol; 
	/* tables Tauti*Vp/nofreq nopercent nocol; */
run;

/* 2.2.3-2.2.4 Suhteellinen riski ja vedonlyöntisuhde */
data Avadex;
	input kasvain $ altistus $ lkm;
	cards;
	on on 4  
	on ei 5 
	ei on 12
	ei ei 74
	;
proc freq order=data;
	weight  lkm;
	tables altistus*kasvain/measures nopercent nocol;
run;

/* 2.3.1 Riippuvuus järjestysasteikolla: Työtyytyväisyys

data jobsat;
	input tulot $ tyytyv $ lkm;
	cards;
	i<6	ETyytym	20
	i<6	HTyytym	24
	i<6	MTyytyv	80
	i<6	ETyytyv 82 
	i6-14	ETyytym	22
	i6-14	HTyytym	38
	i6-14	MTyytyv	104	
	i6-14	ETyytyv 125 
	i15-24	ETyytym	13
	i15-24	HTyytym	28
	i15-24	MTyytyv	81
	i15-24	ETyytyv	113
	i25-	ETyytym	7
	i25-	HTyytym	18
	i25-	MTyytyv	54
	i25-	ETyytyv 92
	;
proc freq data=jobsat order=data;
	weight  lkm;
	tables tulot*tyytyv/measures nopercent norow nocol;
run;

/* 3.  Yhteensopivuuden testaus */
data vari;
	input hiukset $  lkm @@;
	cards;
	vaalea 30 ruskea 47 punainen 23
	;
proc freq order=data;
	weight  lkm;
	tables hiukset/  testf=(30 60 10);
run;
data varit;
	input hiukset $ silmat $ lkm @@;
	cards;
	vaalea sininen 15 vaalea ruskea 10 vaalea vihrea 5
 	ruskea sininen 18 ruskea ruskea 27 ruskea vihrea 2  
	punain sininen 10 punain ruskea  7 punain vihrea 6
	;
proc freq order=data;
	weight  lkm;
	tables hiukset*silmat/ nocol norow testf=(12 15 3 22 34 4 6 1 3);
run;
Chi=1-probchi(288.268,8);
Put chi;
run;
/* 3.3 Riippumattomuuden testaus */

data puoluek;
        input sukupuoli $ puolue $ lkm @@;
        cards;
        naiset demokraatti  279 naiset riippumaton 73 naiset republikaani 225
        miehet demokraatti  165 miehet riippumaton 47 miehet republikaani 191
        ;
proc freq data=puoluek order=data;
        weight  lkm;
        tables sukupuoli*puolue / chisq expected nopercent nocol norow;
run;

/* 3.3.3 Residuaalit */
proc genmod data=puoluek;
	class sukupuoli puolue;
	model lkm = sukupuoli puolue / dist=poi residuals;
run;
/* model lkm = sukupuoli puolue / dist=poi obstats; */

/* 3.4.2 Raskaudenaikainen alkoholin käyttö ja sikiön epämuodostumat */

/* Esimerkki. Äidin alkoholin käyttö ja vauvojen epämuodostumat. */
data infants;
input alkoholi $ epamuod $ lkm @@;
cards;
C=0 Ei 17066 C<1 Ei 14464 C1-2 Ei 788  C3-5 Ei 126 C>5 Ei 37
C=0 On 48    C<1 On 38    C1-2 On 5    C3-5 On 1   C>5 On 1 
;
proc freq order=data; weight lkm;
   tables alkoholi*epamuod / chisq cmh1 nopercent nocol;
run;

/* Seuraavassa korrelaatio ja M^2 laskettu luokkakeskipisteitä käyttäen */                                   
data infants;
input alkoholi epamuod $ lkm @@;
cards;
0 Ei 17066  0.5 Ei 14464  1.5 Ei 788  4.0 Ei 126  7.0 Ei 37
0 On  48    0.5 On 38     1.5 On 5    4.0 On 1    7.0 On 1 
;
proc freq; weight lkm;
   tables alkoholi* epamuod / chisq cmh1 noprint;
run;

/* Lasketaan korrelaatio keskijärjestyksiä käyttäen */                                   

proc freq; weight lkm;
   tables alkoholi* epamuod / cmh1 scores=ridit noprint;
run;

/* 3.6.3 Fisherin tarkka testi */ 
/* Fisherin teekoe */
/* =========================================================== */
 Data Fisher;
      input k_ensin $ arvaus $ lkm @@;
      cards;
      maito maito 3 maito tee 1
      tee   maito 1 tee   tee 3
      ;
proc freq order=data; weight lkm;
   tables k_ensin*arvaus / chisq nocol norow nopct;
run;

/* Tavallinen chisq-optio tables-komennossa tulostaa myös Fisherin */
/* tarkan testin p-arvot */

/* Voit kokeilla myös seuraavaa erillistä exact-komentoa */

proc freq order=data; weight lkm;
   tables k_ensin*arvaus / nocol norow nopct;
   exact chisq; 
run;

/* 4. Luku */
/* 3-ulotteiset taulukot */

/*  Kuolemanrangaistus uhrin ja syytetyn rodun mukaan */
Data Tau4_2;
input syytetty $ kuolema $ lkm @@;
cards;
Valkoinen On 19  Valkoinen Ei 141
Musta     On 17  Musta     Ei 149
;
proc freq order=data; weight lkm;
   tables syytetty*kuolema / measures nopercent nocol;
run;

Data Tau4_3; /* Uhri valkoinen */
input syytetty $ kuolema $ lkm @@;
cards;
Valkoinen On 19  Valkoinen Ei 132
Musta     On 11      Musta Ei  52
;
proc freq order=data; weight lkm;
   tables syytetty*kuolema / measures nopercent nocol;
run;

Data Tau4_4; /* Uhri musta */
input syytetty $ kuolema $ lkm @@;
cards;
Valkoinen On 0.5   Valkoinen Ei  9.5
    Musta On 6.5       Musta Ei 97.5
;
proc freq order=data; weight lkm;
   tables syytetty*kuolema / measures nopercent nocol;
run;

Data penalty;
input uhri $ syytetty $ kuolema $ lkm @@;
cards;
Valkoinen Valkoinen On 19 Valkoinen Valkoinen Ei 132
Valkoinen     Musta On 11 Valkoinen     Musta Ei 52
Musta     Valkoinen On  0     Musta Valkoinen Ei  9
Musta         Musta On  6     Musta     Musta Ei 97
;
proc freq order=data; weight lkm;
   tables uhri*syytetty*kuolema / cmh measures nopercent nocol;
run;

/* Taulukon 4.5 aineisto. Hoidon ja vasteen ehdollinen riippumattomuus */

data Tau4_5;
input hoito $ tulos $ lkm @@;
cards;
h1 onnist 20   h1  epaonn 20
h2 onnist 20   h2  epaonn 40
;
proc freq order=data; weight lkm;
   tables hoito*tulos / measures nopercent nocol;
run;

data Tau4_5;
input sairaala $ hoito $ tulos $ lkm @@;
cards;
s1 h1 onnist 18  s1  h1  epaonn 12
s1 h2 onnist 12  s1  h2  epaonn  8
s2 h1 onnist  2  s2  h1  epaonn  8
s2 h2 onnist  8  s2  h2  epaonn 32  
;
proc freq order=data; weight lkm;
   tables sairaala*hoito*tulos / chisq nopercent nocol;
run;
proc freq order=data; weight lkm;
   tables sairaala*hoito*tulos / cmh nopercent nocol;
run;


/* Taulukon 4.6 aineisto. Kiinalainen tutkimus tupakoinnin ja keuhkosyövän */    
/* välisestä yhteydestä */

/* =============== */
data cmh;
input center $ smokes $ cancer $ count @@;
cards;
Beijing Yes Yes 126  Beijing Yes No 100  Beijing No Yes 35 Beijing No No  61
Shanghai Yes Yes 908 Shanghai Yes No 688 Shanghai No Yes 497 Shanghai No No 807
Shenyang Yes Yes 913 Shenyang Yes No 747 Shenyang No Yes 336  Shenyang No No 598
Nanjing Yes Yes 235  Nanjing Yes No 172   Nanjing No Yes  58  Nanjing No No 121
Harbin Yes Yes 402  Harbin Yes No 308   Harbin No Yes 121  Harbin No No 215
Zhengzho Yes Yes 182 Zhengzho Yes No 156 Zhengzho No Yes 72  Zhengzho No No  98
Taiyuan Yes Yes  60  Taiyuan Yes No  99   Taiyuan No Yes  11  Taiyuan No No  43
Nanchang Yes Yes 104  Nanchang Yes No 89  Nanchang No Yes 21 Nanchang No No  36
;
proc freq order=data; weight count;
   tables center*smokes*cancer / cmh nopercent nocol;
run;
/* ************* */

/* ********************************* */
/* Esimerkki: Kuorsaus ja sydäntauti */
/* ********************************* */

Data kuorsaus;
input kuorsaus tauti lkm;
cards;
0 1  24  
0 0  1355
2  1  35   
2  0  603
4 1  21
4 0  192
5  1  30   
5  0  224
;
proc freq order=data;
 weight lkm;	 
   tables kuorsaus*tauti/ chisq nocol nopct;
run;

proc genmod data=kuorsaus descending;
freq lkm; 
model tauti = kuorsaus /  dist=bin  link=identity;/* Linearinen tn-malli */
run;

proc genmod data=kuorsaus descending;
freq lkm; 
model tauti = kuorsaus /  dist=bin  link=logit; /* Logistinen regressiomalli */
run;

Data kuorsaus;
input kuorsaus tauti ei_tauti;
kaikki=tauti+ei_tauti;
cards;
0  24  1355
2  35   603
4  21   192
5  30   224
;
proc genmod; 
model tauti/kaikki = kuorsaus /  dist=bin  link=identity predicted;/* Linearinen tn-malli */ 
run;
proc genmod; 
model tauti/kaikki = kuorsaus /  dist=bin  link=logit predicted; /* Logistinen regressiomalli */

proc genmod; model tauti/kaikki = kuorsaus /  dist=bin  link=probit predicted; /* Probit-malli */
run;

/* ******************* */
/* Poissonin regressio */
/* ******************* */

Kopio itsellesi SAS-aineisto crabs.sd2 hakemistosta
http://mtl.uta.fi/tilasto/frekvenssiaineistot/Aineistot/

Jos kopioit sen esimerkiksi hakemistoon 'c:\tilap', niin seuraavalla 
libname-komennolla voit luoda SAS-kirjaston 'fraadata'.

libname fraadata 'c:\tilap';
proc genmod data = fraadata.crabs;
model satell = width / dist=poi link=log; /* Logaritmisesti lineaarinen Poissonin malli */
run;
model satell = / dist=poi link=identity; /* Lineaarinen Poissonin malli */  
run;
  
/* Ikääntyneiden autoilijoiden (65-84) onnettomuusaste */

data edrivers;
input sex  years accident;
lyears=log(years);
cards;
1  21.4 320
0  17.3 175
;
proc genmod data=edrivers order=data;
model accident = sex  /  dist=poisson  link=log offset=lyears; /* Poissonin regressio */
run; /* Sex on dummy-muuttuja */

data edrivers;
input sex $ years accident;
lyears=log(years);
cards;
Men   21.4 320
Women 17.3 175
;
proc genmod data=edrivers;
class sex;
model accident = sex  /  dist=poisson  link=log offset=lyears; /* Sex luokitteluasteikollinen */
run;

/* Vakuutusaineisto */

data vakuutus;
input autot  korv koko $;
lautot=log(autot);
cards;
900  143 pieni  
1700 110 kkoko  
400   15 suuri 
;
proc genmod data=vakuutus order=data;
class koko;
model korv = koko/  dist=poisson  link=log offset=lautot;
run;

data vakuutus;
input autot  korv koko;
lautot=log(autot);
cards;
900  143 1  
1700 110 2  
400   15 3 
;
proc genmod data=vakuutus;
model korv = koko/  dist=poisson  link=log offset=lautot;
run;

/* Seuralaisten lkm pituuden eri arvoilla */

libname Fraadata 'c:\kurssit\Fraa\Fraa01\Datat';
data fraadata.crabs_w;
infile 'c:\kurssit\Fraa\Fraa01\Datat\crab_w.txt';
input  width cases tsatell;
l_cases=log(cases);
proc genmod;
   model tsatell = width / dist=poi link=log offset=l_cases;
run;
/* ******************************  */
/* 5.4.1 Testaus: Wald ja us-testi */

libname Fraadata 'c:\kurssit\Fraa\Fraa01\Datat';
proc genmod data = fraadata.crabs;
model satell = / dist=poi link=log;
run;
proc genmod;
   model satell = width / dist=poi link=log type3;
run;

/* ******************************  */
/* 5.4.2 Poissonin mallin sopivuus aineistoon */

/* Estimoidut frekvenssit */
libname Fraadata 'c:\kurssit\Fraa\Fraa01\Datat';
proc genmod data = fraadata.crabs;
   model satell = width / dist=poi link=log predicted;
run;

Data crabstab;
input lkm width satell var;
cards;
14 22.69 1.00 2.77
14 23.84 1.43 8.88
28 24.78 2.39 6.54
39 25.84 2.69 11.38
22 26.79 2.86 6.89
24 27.74 3.88 8.81
18 28.67 3.94 16.88
14 30.41 5.14 8.29
;
proc genmod;
   freq lkm;
   model satell = width / dist=poi link=log;
run;

/* residuaalit */
libname Fraadata 'c:\kurssit\Fraa\Fraa01\Datat';
proc genmod data = fraadata.crabs;
   model satell = width / dist=poi link=log residuals;
run;
/* obstats: kaikki havaintosuureet */
proc genmod data = fraadata.crabs;
   model satell = width / dist=poi link=log obstats;
run;

/* *********************    */
/*     6. Luku 		    */
/*     Logistinen regressio */

/* Logist. reg. mallin parametrit chd-datasta */

libname Fraadata 'c:\kurssit\Fraa\Fraa01\Datat';
data Fraadata.chddata;
infile 'c:\kurssit\Fraa\Fraa01\Datat\chddata.txt';
input id age chd;
run;

PROC GENMOD DATA=Fraadata.chddata DESCENDING;;
  MODEL chd=age / DIST=BINOMIAL;
RUN;

PROC LOGISTIC DATA=Fraadata.chddata DES;
  MODEL chd=age;
RUN;

/* Luokitellusta aineistosta sairastuneiden osuudet */

libname Fraadata 'c:\kurssit\Fraa\Fraa01\Datat';
data chd1;
set fraadata.chddata;
agec = 0;
if age ge 20 and age le 29 then agec=1;
if age ge 30 and age le 34 then agec=2;
if age ge 35 and age le 39 then agec=3;
if age ge 40 and age le 44 then agec=4;
if age ge 45 and age le 49 then agec=5;
if age ge 50 and age le 54 then agec=6;
if age ge 55 and age le 59 then agec=7;
if age ge 60 and age le 99 then agec=8;
run;
proc tabulate data=chd1;
class agec;
var age chd;
table agec all,
	N (age chd)*mean;
run;

/* 6.2.2 Luottamusvälit */

libname Fraadata 'c:\kurssit\Fraa\Fraa01\Datat';
PROC LOGISTIC DATA=Fraadata.chddata DES;
  MODEL chd=age / CLPARM=WALD CLPARM=PL;
RUN;

PROC GENMOD DATA=Fraadata.chddata DESCENDING;;
  MODEL chd=age / DIST=BINOMIAL WALDCI LRCI;
RUN;

PROC LOGISTIC DATA=Fraadata.chddata DES;
  MODEL chd=age / CLODDS=BOTH;
RUN;

libname Fraadata 'c:\kurssit\Fraa\Fraa01\Datat';
PROC LOGISTIC DATA=LogReg.chddata DES;
  UNITS age=5;
  MODEL chd = age / CLODDS=WALD CLODDS=PL;
RUN;

/* 6.3 Mallin tarkistukset */

libname Fraadata 'c:\kurssit\Fraa\Fraa01\Datat';

PROC LOGISTIC DATA=Fraadata.chddata DES;
  MODEL chd=age / aggregate scale=none;
RUN;

PROC GENMOD DATA=Fraadata.chddata DESCENDING;;
  MODEL chd=age / DIST=BINOMIAL;
RUN;

/* Luokiteltu aineisto */
data chdclass;
input total desease meanage;
cards;
10  1 25.40
15  2 32.00
12  3 36.92
15  5 42.33
13  6 47.23
 8  5 51.88
17 13 56.88
10  8 63.00

PROC GENMOD DATA=chdclass DESCENDING;
  MODEL desease/total=meanage / DIST=BINOMIAL;
RUN;
PROC LOGISTIC DATA=chdclass DES;
  MODEL desease/total=meanage / aggregate scale=none;
RUN;


/* 7 Logaritmisesti lineaariset mallit  */
/* 7.1 2-suuntaisten taulukoiden loglin mallit */ 

data Dpenalty;
input Trotu $ RKuo $ frekv @@;
cards;
Valk  Kylla 19 Valk  Ei 141
Musta Kylla 17 Musta Ei 149
;
proc freq data=Dpenalty order=data;
        weight  frekv;
        tables Trotu*RKuo / chisq expected nopercent nocol norow;
run;
proc genmod order=data; class Trotu  RKUO;
   model  frekv = Trotu RKuo / dist=poi link=log pred;

proc catmod order=data; weight frekv;
        model Trotu*RKuo= _response_ / pred=freq noiter;
        loglin Trotu RKuo;
run;

/* ***************** */
/* 3-suuntaiset taulukot  */

data aids;
input rotu $ azt $ oireet $ lkm @@; 
cards; 
valkoinen  kylla kylla 14 valkoinen  kylla ei  93  valkoinen  ei kylla 32 valkoinen  ei ei  81 
musta      kylla kylla 11 musta      kylla ei  52  musta      ei kylla 12 musta      ei ei  43
; 
proc genmod   order=data; class rotu  azt oireet; 
   model lkm = azt oireet rotu rotu*azt azt*oireet rotu*oireet /  dist=poi  link=log type3; 
run;

proc catmod order=data; weight lkm;
        model azt*oireet*rotu = _response_ / pred=freq noiter;
        loglin azt|oireet|rotu @2;
run;

/*************************/
data Dpenalty;
input Trotu $ Urotu $ RKuo $ frekv @@;
cards;
Valk Valk K 19 Valk Valk Ei 132
Valk Musta K 0 Valk Musta Ei 9
Musta Valk K 11  Musta Valk Ei 52
Musta Musta K 6  Musta Musta Ei 97
;
proc genmod order=data; class Trotu  Urotu  RKUO;
   model  frekv = Trotu Urotu RKUO Trotu*Urotu Trotu*RKUO Urotu*RKUO
/ dist=poi link=log type3;
run;

/* Alkoholi, savukkeet ja marihuana  */
/* ********************************  */ 
data drugs;
input alkoholi $ savukkeet $ marihuana $ lkm @@; 
cards;
Kylla Kylla Kylla 911   Kylla Kylla Ei 538   Kylla Ei Kylla 44   Kylla Ei Ei 456
 Ei   Kylla Kylla   3    Ei   Kylla Ei 43     Ei   Ei Kylla  2   Ei    Ei Ei 279
;
/* (XY,XZ,YZ) */
proc genmod order=data; class alkoholi savukkeet marihuana;
   model  lkm = alkoholi savukkeet marihuana 
		alkoholi*savukkeet alkoholi*marihuana savukkeet*marihuana / 
	      dist=poi link=log predicted type3;
proc catmod order=data; weight lkm;
   model alkoholi*savukkeet*marihuana = _response_ / pred=freq noiter; 
loglin alkoholi|savukkeet alkoholi|marihuana savukkeet|marihuana;
run;
/* (XZ,YZ) */
proc genmod; class alkoholi savukkeet marihuana;
   model  lkm = alkoholi savukkeet marihuana alkoholi*marihuana savukkeet*marihuana/ 
	     dist=poi link=log predicted type3;
proc catmod; weight lkm;
   model alkoholi*savukkeet*marihuana = _response_ / pred=freq noiter; 
loglin alkoholi|marihuana savukkeet|marihuana;
run;

/* 10. Harjoitukset */
/* 10.2  */
data drugs;
input alkoholi $ savukkeet $ marihuana $ lkm @@; 
cards;
Kylla Kylla Kylla 911   Kylla Kylla Ei 538   Kylla Ei Kylla 44   Kylla Ei Ei 456
 Ei   Kylla Kylla   3    Ei   Kylla Ei 43     Ei   Ei Kylla  2   Ei    Ei Ei 279
;
/* (AS,SM) */
proc genmod order=data; class alkoholi savukkeet marihuana;
   model  lkm = alkoholi savukkeet marihuana 
		alkoholi*savukkeet savukkeet*marihuana / 
	      dist=poi link=log predicted;
run;

/* 10.5 Aspiriini, vatsahaava ja pohjukaissuolen haava  */

data vhaava;
input Sijainti $ Tauti $ Aspirin $ frekv @@; 
cards;
Vatsa On On 25 Vatsa On Ei 39 Vatsa Ei On 6 Vatsa Ei Ei 62
Pohjs On On  8 Pohjs On Ei 49 Pohjs Ei On 8 Pohjs Ei Ei 53
;
proc genmod order=data; class Sijainti  Tauti  Aspirin;
   model  frekv = Sijainti Tauti Aspirin Sijainti*Tauti Sijainti*Aspirin Tauti*Aspirin
/ dist=poi link=log type3;
run;