Calcolo delle Probabilita'e Statistica Matematica 2
A.A. 2010/11
LABORATORIO DI SAS – A. MICHELETTI
LETTURA DATI E STATISTICA DESCRITTIVA
LETTURA DATI
Dati scritti nel corpo del programma:
ESEMPIO1.SAS
data class;
input nome$ alt peso age;
cards;
Alfredo 168 70.0 14
Alice
156 58.0 13
Barbara 165 55.5 13
Carla
162 54.5 14
Enrico
163 70.0 14
Giorgio 170 62.5 13
Carlo
168 60.5 12
Luisa
158 52.0 12
Franco
164 60.0 14
Sandra
155 48.0 13
;
proc print data=class;
run;
SAS System
Oss
1
2
3
4
5
6
7
8
9
10
nome
Alfredo
Alice
Barbara
Carla
Enrico
Giorgio
Carlo
Luisa
Franco
Sandra
13:13 Tuesday, October 12, 2010
alt
peso
age
168
156
165
162
163
170
168
158
164
155
70.0
58.0
55.5
54.5
70.0
62.5
60.5
52.0
60.0
48.0
14
13
13
14
14
13
12
12
14
13
Dati scritti in un file testo separato
IMPDATITESTO.SAS
data class2;
filename dati 'E:\didattica\aa10-11\SAS\file_lezione1\studenti.txt';
infile dati;
input nome$ alt peso age;
cards;
Dati da prelevare solo da alcune colonne di un file di testo
1
data lavoro1;
filename dati 'E:\didattica\aa10-11\SAS\dati\Laur92fs.txt';
infile dati;
input tipolaurea 1-1 regione 2-3 mesi 41-42;
cards;
Dati scritti in un file di formato database (Excel, Access, etc.). Questo file
si produce anche usando la modalita’ interattiva di caricamento dei dati.
IMPORTDATI.SAS
PROC IMPORT OUT= Tmp1.studenti
DATAFILE= "C:\user\statmat1\studenti1.xls"
DBMS=EXCEL2000 REPLACE;
GETNAMES=YES;
RUN;
PROCEDURA MEANS: INDICI STATISTICI DESCRITTIVI
Useremo spesso come dataset di riferimento i file del tipo
telefonate_gen2008
In esso sono riportati i dati relativi alle telefonate ricevute dal centralino del servizio 118 di Milano nel 2008, in questo caso nel mese di
gennaio. Ne useremo altri relativi ad altri mesi dell’anno.
In tutti, le variabili riportate sono:
DT_INLINE= data e ora di arrivo della telefonata. Se non ci sono operatori liberi, la chiamata viene messa in attesa
DT_INCALL=data e ora in cui risponde un operatore. I dati persi si riferiscono a telefonate in cui il chiamante ha riappeso prima che un
operatore rispondesse
DT_ENDCALL=data e ora di fine telefonata, identificata con l’istante in cui il chiamante riappende.
GIORNO (del mese), MESE, Weekday (giorno della settimana),ORA, Minuti, Secondi di inizio chiamata: estratti da DT_INLINE
ImpegnoOper = tempo passato al telefono con un operatore = DT_ENDCALL-DT_INCALL; espresso in secondi
Dur_Telef= durata della telefonata = DT_ENDCALL – DT_INLINE; espresso in secondi
DS_COMUNE= comune in cui e’ avvenuto l’infortunio per cui e’ richiesto l’intervento del 118
ID_CODICE=codice di urgenza dell’infortunio: R=rosso, G=giallo, V=verde, B=bianco, X=codice errato
DS_MOTIVO=motivo dell’infortunio
DS_DET_MOTIVO=motivo specifico dell’infortunio
MEANS.SAS
data Tmp1.telefonate_gen08;
proc means;
var impegnoOper Dur_Telef;
run;
SAS
18:55 Thursday, October 14, 2010
La procedura MEANS
Variabile
Etichetta
N
Media
Dev std
Minimo
Massimo
ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒ
ImpegnoOper
ImpegnoOper
54611
60.2167512
54.4911097
0
1997.00
Dur_Telef
Dur_Telef
65534
68.1801202
59.8549488
0
2000.00
ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒ
MEANS2.SAS
proc means data=mylib.telefonate_gen2008 n max min range mean std median kurtosis skewness
maxdec=2;;
var impegnoOper Dur_Telef;
run;
quit;
SAS
10:22 Friday, October 15, 2010
1
La procedura MEANS
Variabile
Etichetta
N
Media
Dev std
Minimo
Massimo
ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒ
ImpegnoOper
ImpegnoOper
48916
60.7105446
54.8978866
0
1997.00
Dur_Telef
Dur_Telef
58770
68.6420112
60.2828170
0
2000.00
ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒ
SAS
10:22 Friday, October 15, 2010
2
La procedura MEANS
Variabile
Etichetta
N
Massimo
Minimo
Range
Media
Dev std
ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒ
ImpegnoOper ImpegnoOper 48916
1997.00
0.00
1997.00
60.71
54.90
Dur_Telef
Dur_Telef
58770
2000.00
0.00
2000.00
68.64
60.28
ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒ
Variabile
Etichetta
Mediana
Curtosi
Skewness
ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒ
ImpegnoOper ImpegnoOper
51.00
60.26
4.15
Dur_Telef
Dur_Telef
59.00
35.84
2.98
ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒ
MEANS3.SAS
proc means data=mylib.telefonate_gen2008 n max min range mean std median kurtosis skewness
maxdec=2;;
var impegnoOper Dur_Telef;
class ID_codice;
run;
quit;
SAS
10:22 Friday, October 15, 2010
La procedura MEANS
N.
ID_CODICE
oss Variabile
Etichetta
N
Massimo
Minimo
Range
ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒ
B
11638 ImpegnoOper ImpegnoOper 11638
1326.00
1.00
1325.00
Dur_Telef
Dur_Telef
11638
1353.00
5.00
1348.00
G
10462
ImpegnoOper
ImpegnoOper
10462
699.00
1.00
698.00
3
Dur_Telef
Dur_Telef
10462
709.00
5.00
704.00
R
1734
ImpegnoOper
Dur_Telef
ImpegnoOper
Dur_Telef
1734
1734
401.00
450.00
0.00
5.00
401.00
445.00
V
9067
ImpegnoOper
Dur_Telef
ImpegnoOper
Dur_Telef
9067
9067
620.00
646.00
1.00
8.00
619.00
638.00
X
8507
ImpegnoOper ImpegnoOper
8507
1997.00
0.00
1997.00
Dur_Telef
Dur_Telef
8507
2000.00
2.00
1998.00
ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒ
N.
ID_CODICE
oss Variabile
Etichetta
Media
Dev std
Mediana
Curtosi
ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒ
B
11638 ImpegnoOper ImpegnoOper
68.26
67.50
50.00
35.85
Dur_Telef
Dur_Telef
88.73
70.94
72.00
29.38
G
10462
ImpegnoOper
Dur_Telef
ImpegnoOper
Dur_Telef
66.60
86.52
37.78
42.76
59.00
79.00
19.82
12.48
R
1734
ImpegnoOper
Dur_Telef
ImpegnoOper
Dur_Telef
78.87
98.29
39.83
44.29
70.00
90.00
7.96
6.13
V
9067
ImpegnoOper
Dur_Telef
ImpegnoOper
Dur_Telef
81.84
103.07
46.43
51.11
71.00
93.00
13.35
9.69
X
8507
ImpegnoOper ImpegnoOper
27.28
48.91
11.00
327.26
Dur_Telef
Dur_Telef
42.49
53.18
26.00
229.73
ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒ
N.
ID_CODICE
oss Variabile
Etichetta
Skewness
ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒ
B
11638 ImpegnoOper ImpegnoOper
4.08
Dur_Telef
Dur_Telef
3.56
G
10462
ImpegnoOper
Dur_Telef
ImpegnoOper
Dur_Telef
2.69
2.01
R
1734
ImpegnoOper
Dur_Telef
ImpegnoOper
Dur_Telef
1.91
1.57
V
9067 ImpegnoOper ImpegnoOper
2.55
ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒ
SAS
10:22 Friday, October 15, 2010
La procedura MEANS
N.
ID_CODICE
oss Variabile
Etichetta
Skewness
ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒ
V
9067 Dur_Telef
Dur_Telef
2.08
X
8507
ImpegnoOper ImpegnoOper
10.82
Dur_Telef
Dur_Telef
8.59
ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒ
4
PROCEDURA FREQ: CALCOLO FREQUENZE E ISTOGRAMMI (senza grafici)
Studiamo dapprima la distribuzione di una variabile discreta:
FREQ1.sas
/* analizzo la distribuzione delle ore del giorno:
ci sono orari con chiamate piu' frequenti?*/
proc freq data=mylib.telefonate_gen2008 ;
TABLES ora /out=work.contotel; /*salvo l'output nel dataset work.contotel*/
run;
SAS
10:22 Friday, October 15, 2010
23
La procedura FREQ
ORA
Frequenza
Percentuale
ORA
Frequenza
Percentuale
cumulativa
cumulativa
ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒ
0
1819
3.10
1819
3.10
1
1671
2.84
3490
5.94
2
1454
2.47
4944
8.41
3
1294
2.20
6238
10.61
4
1122
1.91
7360
12.52
5
1145
1.95
8505
14.47
6
1292
2.20
9797
16.67
7
1685
2.87
11482
19.54
8
2169
3.69
13651
23.23
9
2949
5.02
16600
28.25
10
3344
5.69
19944
33.94
11
3279
5.58
23223
39.52
12
3113
5.30
26336
44.81
13
3080
5.24
29416
50.05
14
3083
5.25
32499
55.30
15
2980
5.07
35479
60.37
16
3009
5.12
38488
65.49
17
3062
5.21
41550
70.70
18
3266
5.56
44816
76.26
19
3231
5.50
48047
81.75
20
3321
5.65
51368
87.41
21
2951
5.02
54319
92.43
22
2347
3.99
56666
96.42
23
2104
3.58
58770
100.00
Analizziamo ora una variabile continua: prima di applicare la procedura freq dovremo
“discretizzarla” attraverso la procedura format:
FREQ2.sas
/* analizzo la distribuzione dell' impegno operatore. Le durate sono in secondi*/
proc format; /* definisce range di valori:
serve per raggruppare i dati,
associando un'etichetta a ogni gruppo*/
value durata low-30='A meno di mezzo minuto' /*nell'output elenca le etichette in ordine
lessicografico */
30-60='B mezzo-un minuto'
60-120='C uno-due minuti'
120-300='D due-cinque minuti'
300-high='E più di 5 minuti';
run;
proc freq data=mylib.telefonate_gen2008 order=formatted;
TABLES ImpegnoOper;
format ImpegnoOper durata.;
run;
quit;
SAS
La procedura FREQ
ImpegnoOper
Frequenza
Percentuale
ImpegnoOper
Frequenza
Percentuale
cumulativa
cumulativa
ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒ
A meno di mezzo minuto
14233
29.10
14233
29.10
B mezzo-un minuto
14620
29.89
28853
58.98
C uno-due minuti
15186
31.05
44039
90.03
D due-cinque minuti
4616
9.44
48655
99.47
E più di 5 minuti
261
0.53
48916
100.00
Frequenza mancanti = 9854
Vediamo ora come si possono effettuare analisi separate per classi diverse, dove le classi sono
i valori assunti da una variabile discreta. In tal caso si usa prima la procedura sort per
ordinare i dati secondo valori crescenti di tale variabile discreta, e poi si usa lo statement
“by” all’interno della procedura freq.
FREQ3.SAS
/* analizzo la distribuzione delle ore del giorno divise per codici
ci sono orari con codici piu' frequenti?*/
/*prima devo ordinare i dati per valori crescenti della variabile ID_CODICE*/
proc sort data=tmp1.telefonate_gen2008 out=work.ordinati;
by ID_CODICE;
run;
proc freq data=work.ordinati;
TABLES ora /out=work.contocodice; /*salvo l'output nel dataset work.contocodice*/
by ID_CODICE;
run;
quit;
OUTPUT:
SAS
----------------------------------------- ID_CODICE=
17:33 Monday, April 4, 2011
1
------------------------------------------
La procedura FREQ
ORA
Frequenza
Percentuale
ORA
Frequenza
Percentuale
cumulativa
cumulativa
ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒ
0
457
2.63
457
2.63
1
398
2.29
855
4.92
2
356
2.05
1211
6.98
3
334
1.92
1545
8.90
4
274
1.58
1819
10.48
5
294
1.69
2113
12.17
6
339
1.95
2452
14.12
7
403
2.32
2855
16.44
8
512
2.95
3367
19.39
9
754
4.34
4121
23.74
10
958
5.52
5079
29.25
11
1063
6.12
6142
35.38
12
1064
6.13
7206
41.50
13
988
5.69
8194
47.20
14
1024
5.90
9218
53.09
15
924
5.32
10142
58.41
16
969
5.58
11111
64.00
17
978
5.63
12089
69.63
18
1159
6.68
13248
76.30
19
1084
6.24
14332
82.55
20
1105
6.36
15437
88.91
21
830
4.78
16267
93.69
22
570
3.28
16837
96.98
23
525
3.02
17362
100.00
SAS
17:33 Monday, April 4, 2011
2
----------------------------------------- ID_CODICE=B -----------------------------------------La procedura FREQ
ORA
Frequenza
Percentuale
ORA
Frequenza
Percentuale
cumulativa
cumulativa
ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒ
0
450
3.87
450
3.87
1
387
3.33
837
7.19
2
338
2.90
1175
10.10
3
301
2.59
1476
12.68
4
251
2.16
1727
14.84
5
263
2.26
1990
17.10
6
244
2.10
2234
19.20
7
306
2.63
2540
21.83
8
404
3.47
2944
25.30
9
522
4.49
3466
29.78
10
633
5.44
4099
35.22
11
555
4.77
4654
39.99
12
542
4.66
5196
44.65
13
508
4.37
5704
49.01
14
525
4.51
6229
53.52
15
471
4.05
6700
57.57
16
499
4.29
7199
61.86
17
519
4.46
7718
66.32
18
548
4.71
8266
71.03
19
662
5.69
8928
76.71
20
816
7.01
9744
83.73
21
742
6.38
10486
90.10
22
631
5.42
11117
95.52
23
521
4.48
11638
100.00
SAS
17:33 Monday, April 4, 2011
3
----------------------------------------- ID_CODICE=G -----------------------------------------La procedura FREQ
ORA
Frequenza
Percentuale
ORA
Frequenza
Percentuale
cumulativa
cumulativa
ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒ
0
310
2.96
310
2.96
1
299
2.86
609
5.82
2
290
2.77
899
8.59
3
224
2.14
1123
10.73
4
204
1.95
1327
12.68
5
194
1.85
1521
14.54
6
252
2.41
1773
16.95
7
341
3.26
2114
20.21
8
513
4.90
2627
25.11
9
646
6.17
3273
31.28
10
645
6.17
3918
37.45
11
631
6.03
4549
43.48
12
603
5.76
5152
49.24
13
556
5.31
5708
54.56
14
514
4.91
6222
59.47
15
499
4.77
6721
64.24
16
486
4.65
7207
68.89
17
506
4.84
7713
73.72
18
552
5.28
8265
79.00
19
499
4.77
8764
83.77
20
498
4.76
9262
88.53
21
452
4.32
9714
92.85
22
384
3.67
10098
96.52
23
364
3.48
10462
100.00
SAS
17:33 Monday, April 4, 2011
4
----------------------------------------- ID_CODICE=R -----------------------------------------La procedura FREQ
ORA
Frequenza
Percentuale
ORA
Frequenza
Percentuale
cumulativa
cumulativa
ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒ
0
70
4.04
70
4.04
1
64
3.69
134
7.73
2
35
2.02
169
9.75
3
36
2.08
205
11.82
4
36
2.08
241
13.90
5
36
2.08
277
15.97
6
50
2.88
327
18.86
7
53
3.06
380
21.91
8
99
5.71
479
27.62
9
109
6.29
588
33.91
10
118
6.81
706
40.72
11
86
4.96
792
45.67
12
94
5.42
886
51.10
13
80
4.61
966
55.71
14
15
16
17
18
19
20
21
22
23
79
85
73
71
83
103
74
71
71
58
4.56
4.90
4.21
4.09
4.79
5.94
4.27
4.09
4.09
3.34
SAS
1045
1130
1203
1274
1357
1460
1534
1605
1676
1734
60.27
65.17
69.38
73.47
78.26
84.20
88.47
92.56
96.66
100.00
17:33 Monday, April 4, 2011
5
----------------------------------------- ID_CODICE=V -----------------------------------------La procedura FREQ
ORA
Frequenza
Percentuale
ORA
Frequenza
Percentuale
cumulativa
cumulativa
ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒ
0
237
2.61
237
2.61
1
268
2.96
505
5.57
2
201
2.22
706
7.79
3
204
2.25
910
10.04
4
188
2.07
1098
12.11
5
167
1.84
1265
13.95
6
197
2.17
1462
16.12
7
316
3.49
1778
19.61
8
397
4.38
2175
23.99
9
589
6.50
2764
30.48
10
594
6.55
3358
37.04
11
541
5.97
3899
43.00
12
460
5.07
4359
48.08
13
512
5.65
4871
53.72
14
493
5.44
5364
59.16
15
482
5.32
5846
64.48
16
457
5.04
6303
69.52
17
448
4.94
6751
74.46
18
400
4.41
7151
78.87
19
427
4.71
7578
83.58
20
407
4.49
7985
88.07
21
436
4.81
8421
92.88
22
343
3.78
8764
96.66
23
303
3.34
9067
100.00
SAS
17:33 Monday, April 4, 2011
6
----------------------------------------- ID_CODICE=X -----------------------------------------La procedura FREQ
ORA
Frequenza
Percentuale
ORA
Frequenza
Percentuale
cumulativa
cumulativa
ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒ
0
295
3.47
295
3.47
1
255
3.00
550
6.47
2
234
2.75
784
9.22
3
195
2.29
979
11.51
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
169
191
210
266
244
329
396
403
350
436
448
519
525
540
524
456
421
420
348
333
1.99
2.25
2.47
3.13
2.87
3.87
4.65
4.74
4.11
5.13
5.27
6.10
6.17
6.35
6.16
5.36
4.95
4.94
4.09
3.91
1148
1339
1549
1815
2059
2388
2784
3187
3537
3973
4421
4940
5465
6005
6529
6985
7406
7826
8174
8507
13.49
15.74
18.21
21.34
24.20
28.07
32.73
37.46
41.58
46.70
51.97
58.07
64.24
70.59
76.75
82.11
87.06
91.99
96.09
100.00
ESERCIZIO 1
Utilizzando le procedure di statistica descrittiva viste, dare una prima risposta alle seguenti domande:
• Il processo di arrivo di telefonate al centralino del 118 (DT_INLINE) si può considerare un processo di Poisson
a intensità costante?
• Ciascuno dei sottoprocessi degli arrivi di telefonate divisi per codice di emergenza, si può considerare un
processo di Poisson a intensità costante?
• Ciascuno dei sottoprocessi degli arrivi di telefonate divisi per motivi (DS_MOTIVO), si può considerare un
processo di Poisson a intensità costante?
MANIPOLAZIONE DI DATASETS
Vediamo come si possono combinare più datasets fra loro
CREO_DATASETS.sas
data dati1;
input anno var1 var2;
cards;
2000 1 11
2001 2 22
2003 3 33
2004 4 44
;
data dati2;
input anno var1 var2;
cards;
2000 5 55
2001 6 66
2002 7 77
2003 8 88
;
/*unisco i due dataset mettendoli uno sopra l'altro */
/*ATTENZIONE: il numero di variabili e la loro collocazione nei due datasets deve essere
uguale*/
data combined;
set dati1 dati2;
run;
/* affianco i due datasets */
/*prima creo due nuove variabili con nomi diversi dalle precedenti,
salvandole in un nuovo dataset*/
data dati22;
set dati2;
var3=var1;
var4=var2;
run;
/*ora affianco dati1 e dati22 */
data combined2;
merge dati22 dati1;
run;
/*ora li affianco ma considerando i valori diversi assunti da "anno"*/
data combined3;
merge dati22 dati1;
by anno;
run; /*N.B.: poichè non ho eliminato var1 e var2 da dati22, mantiene il loro
contenuto dove dati1 e' vuoto, ossia in corrispondenza dell'anno 2002.
Rinominiamo var1 e var2 in dati2 e rieseguiamo: */
data dati222;
set dati2 (rename=(var1=var3 var2=var4));
run;
data combined33;
merge dati1 dati222;
by anno;
run;
/*se voglio aggiornare alcuni valori di un dataset con quelli contenuti in un altro: */
data dati3;
input anno var1 var2;
cards;
2000 8 88
2004 9 99
;
/*aggiorno dati1 con dati3 */
data updated;
update dati1 dati3;
by anno;
run;
/*creo un nuovo dataset contenente solo alcuni dati e con una nuova variabile*/
data new;
set dati1;
where var1 > 1;
if anno < 2002 then var3=1;
else var3=0;
run;
Istogrammi e torte
GRAPH1.SAS
/*disegno un istogramma e una torta del numero di telefonate per ora */
proc gchart data=tmp1.telefonate_gen2008;
vbar ORA / type=percent midpoints=(0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
21 22 23);
pie ORA / fill=solid type=percent other=5 midpoints=(0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
15 16 17 18 19 20 21 22 23);
pie ORA /fill=solid type=percent other=5 slice=arrow midpoints=(0 1 2 3 4 5 6 7 8 9 10
11 12 13 14 15 16 17 18 19 20 21 22 23);
run;
quit;
PERCENTUALE
6
5
4
3
2
1
0
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
ORA
PERCENTUALE di ORA
15
5.07%
14
5.25%
16
5.12%
17
5.21%
13
12
5.24%5.30%
PERCENTUALE di ORA
11
5.58%
10
5.69%
9
5.02%
18
5.56%
19
5.50%
20
5.65%
21
5.02%
ALTRO
30.80%
14
15
5.25%
5.07%
16
5.12%
17
5.21%
18
5.56%
19
5.50%
20
5.65%
21
5.02%
13
12
5.30%
5.24%
11
5.58%
10
5.69%
9
5.02%
ALTRO
30.80%
Gli stessi grafici si possono produrre anche usando l’ opzione DISCRETE, che consente di disegnare una barra o una fetta per
ogni valore della variabile considerata:
GRAPH1_bis.sas
proc gchart data=tmp1.telefonate_gen2008;
vbar ORA / discrete type=percent;
pie ORA / discrete fill=solid type=percent other=5;
pie ORA /discrete fill=solid type=percent other=5 slice=arrow;
run;
quit;
Vediamo ora come usare titoli e come dividere i grafici in sottogruppi:
GRAPH2.SAS
data tmp1.telefonate_gen2008_2; /*creo nuovo dataset*/
set tmp1.telefonate_gen2008; /*copio il vecchio al suo interno ma filtrando i dati*/
where (ID_CODICE eq 'V' or ID_CODICE eq 'G' or ID_CODICE eq 'R') and DS_Comune ne
'#VEDINOTE'; /*copio solo i records con codice significativo e destinazione nota*/
if DS_Comune = 'MILANO' then zona = 'MILANO
'; /*creo una nuova variabile: zona*/
else zona = 'PROVINCIA';
goptions device=win;
title1 'LOCALIZZAZIONE INFORTUNI GENNAIO';
pattern v=s c=pink;
proc gchart data=tmp1.telefonate_gen2008_2;
vbar zona / type=percent; /*disegno istogramma di zona*/
run;
title2 'istogrammi per codice';
proc gchart data=tmp1.telefonate_gen2008_2;
hbar zona / type=percent group=ID_codice autoref g100; /*disegno un istogramma diverso
per ogni codice*/
run;
LOCALIZZAZIONE INFORTUNI GENNAIO
PERCENTUALE
60
50
40
30
20
10
0
MILANO
PROVINCIA
zona
LOCALIZZAZIONE INFORTUNI GENNAIO
istogrammi per codice
ID_CODICE
G
R
V
FREQ.
CUM.
FREQ.
PCT.
CUM.
PCT.
MILANO
6057
6057
57.90
57.90
PROVINCIA
4405
10462
42.10
100.00
MILANO
902
902
52.02
52.02
PROVINCIA
832
1734
47.98
100.00
MILANO
5427
5427
59.85
59.85
PROVINCIA
3640
9067
40.15
100.00
zona
0
10
20
30
40
50
60
PERCENTUALE
Generiamo ora grafici separate per i 3 tipi di codice
GRAPH3.SAS
/* per dividere i grafici in piu' finestre si usa "by" ma bisogna prima ordinare i dati
rispetto alla variabile usata nello statement by */
proc sort data=tmp1.telefonate_gen2008_2 out=tmp1.telefonate_gen2008_sorted;
by ID_codice;
run;
goptions device=win;
title1 color=blue "Infortuni gennaio";
title2; /*cancella tutti i titoli dal secondo in poi. Se non lo faccio restano attivi per tutta
la sessione*/
pattern c=blue;
proc gchart data=tmp1.telefonate_gen2008_sorted;
vbar zona / type=percent ;
by ID_codice; /*disegno un istogramma diverso per ogni codice*/
run;
/*ora affianco gli esiti di codici diversi nello stesso */
proc gchart data=tmp1.telefonate_gen2008_2;
hbar zona / type=percent subgroup=ID_codice;
run;
quit;
Infortuni gennaio
Infortuni gennaio
Infortuni gennaio
ID_CODICE=V
ID_CODICE=R
ID_CODICE=G
PERCENTUALE
PERCENTUALE
PERCENTUALE
60
60
60
50
50
50
40
40
40
30
30
30
20
20
20
10
10
10
0
0
MILANO
0
MILANO
PROVINCIA
PROVINCIA
R
zona
MILANO
zona
PROVINCIA
V
ID_CODICE
zona
ID_CODICE
Infortuni gennaio
zona
MILANO
PROVINCIA
0
10
20
30
40
50
FREQ.
CUM.
FREQ.
PCT.
CUM.
PCT.
12386
12386
58.25
58.25
8877
21263
41.75
100.00
60
PERCENTUALE
ID_CODICE
G
R
V
GRAPH4.SAS
/*disegno un istogramma e una torta dell'impegno operatore */
title1 color=black "Impegno operatore";
pattern c=cyan;
proc gchart data=tmp1.telefonate_gen2008;
vbar ImpegnoOper / type=percent MIDPOINTS= (30 60 90 120 150 180 210 240 270 300 330);
pie ImpegnoOper/ coutline=black woutline=2 percent=inside slice=outside other=5
MIDPOINTS= (60 120 180 240 300);
run;
quit;
Impegno operatore
Impegno operatore
FREQUENZA di ImpegnoOper
PERCENTUALE
50
40
60
39280
30
80.30%
20
5.30%
10
ALTRO
2592
14.40%
0
30
60
90
120
150
180
210
240
270
300
330
120
7044
ImpegnoOper
ESERCIZIO 2
Si considerino i dati sulle telefonate al servizio 118 di Milano. Attraverso procedure numeriche e grafiche di analisi
descrittiva si dia una prima risposta alle seguenti domande:
1. Nelle telefonate del mese di gennaio ci sono differenze significative fra i vari codici di emergenza nella durata
dei “tempi morti”, ossia la differenza: durata telefonata – impegno operatore?
2. Si considerino le sole telefonate relative ad INCIDENTI STRADALI.
a. La distribuzione della tipologia specifica di incidente (auto-auto, auto-moto, etc.) .) di gennaio e di
maggio è diversa?
b. La distribuzione della tipologia specifica di incidente è diversa fra i vari giorni della settimana?
3. Si considerino le sole telefonate relative ad eventi di tipo MEDICO ACUTO.
a. La distribuzione della tipologia specifica di infortunio (cardiocircolatoria, respiratoria, etc.) di gennaio e
di maggio è diversa?
b. La distribuzione della tipologia specifica di infortunio è diversa fra i vari giorni della settimana?
PROCEDURA UNIVARIATE:
INDICI DESCRITTIVI E TEST SULLA MEDIA
UNIVAR.SAS
proc univariate data=tmp1.telefonate_gen2008_2;
var Tmorto;
histogram Tmorto/normal;
qqplot Tmorto;
run;
La procedura UNIVARIATE
Variabile: Tmorto
Momenti
N
Media
Deviazione std
Skewness
SS non corretta
Coeff variaz
48916
19.3736201
20.1530228
1.65336041
38226542
104.02301
Somma dei pesi
Somma delle osservazioni
Varianza
Curtosi
SS corretta
Errore std media
48916
947680
406.144326
2.74586304
19866549.7
0.09112021
Misure statistiche di base
Posizione
Media
Mediana
Moda
Variabilità
19.37362
11.00000
3.00000
Deviazione std
Varianza
Range
Range interquartile
20.15302
406.14433
172.00000
23.00000
Test di posizione: Mu0=0
Test
-Statistica-
-----P-value------
T di Student
Segno
Rango con segno
t
M
S
Pr > |t|
Pr >= |M|
Pr >= |S|
212.6161
24457
5.9816E8
Quantili (Definizione 5)
Quantile
100% Max
99%
95%
90%
75% Q3
50% Mediana
25% Q1
10%
5%
1%
0% Min
Stima
172
84
64
50
27
11
4
3
3
2
0
<.0001
<.0001
<.0001
17:13 Sunday, November 14, 2010
14
La procedura UNIVARIATE
Variabile: Tmorto
Osservazioni estreme
---Inferiori---
---Superiori---
Valore
Oss
Valore
Oss
0
0
1
1
1
30367
10286
56710
54193
52878
152
154
154
163
172
13618
1294
33509
1298
1295
Valori mancanti
Valore
mancante
Conteggio
.
9854
---Percentuale di--Tutte
Oss
le oss
mancanti
16.77
100.00
17:13 Sunday, November 14, 2010
La procedura UNIVARIATE
Distribuzione normale adattata per Tmorto
Parametri per la distribuzione Normale
Parametro
Simbolo
Media
Dev std
Mu
Sigma
Stima
19.37362
20.15302
Test della bontà di adattamento per la distribuzione Normale
Test
-------Statistica-------
---------P-value---------
Kolmogorov-Smirnov
Cramer-von Mises
Anderson-Darling
D
W-Qu
A-Qu
Pr > D
Pr > W-Qu
Pr > A-Qu
0.19311
589.08988
3309.25883
Quantili per la distribuzione Normale
Percentuale
1.0
5.0
10.0
25.0
50.0
75.0
-------Quantile-----Osservato
Stimato
2.00000
3.00000
3.00000
4.00000
11.00000
27.00000
-27.50932
-13.77515
-6.45352
5.78061
19.37362
32.96663
<0.010
<0.005
<0.005
15
90.0
95.0
99.0
50.00000
64.00000
84.00000
45.20076
52.52239
66.25656
distribuzione medico acuto gen-mag
30
Percentuale
25
20
15
10
5
0
0
9
18
27
36
45
54
63
72
81
90
99
108
117
126
135
144
153
162
171
Tmorto
175
150
Tmorto
125
100
75
50
25
0
-6
-4
-2
0
2
4
6
Quantili normali
La distribuzione di Tmorto non e’ gaussiana. Verifichiamo la robustezza del test sulla media provando a testare l’ipotesi che la
media sia uguale a 19, ossia circa il valore della media aritmetica dei dati. Inoltre sovrapponiamo all’istogramma un grafico di una
gaussiana con media e varianza uguali a quelle campionarie:
UNIVAR2.SAS
proc univariate data=tmp1.telefonate_gen2008_2 Mu0=19;
var Tmorto;
histogram Tmorto/normal (mu=est sigma=est);
run;
La procedura UNIVARIATE
Variabile: Tmorto
Momenti
N
Media
Deviazione std
Skewness
SS non corretta
Coeff variaz
48916
19.3736201
20.1530228
1.65336041
38226542
104.02301
Somma dei pesi
Somma delle osservazioni
Varianza
Curtosi
SS corretta
Errore std media
48916
947680
406.144326
2.74586304
19866549.7
0.09112021
Misure statistiche di base
Posizione
Media
Mediana
Moda
Variabilità
19.37362
11.00000
3.00000
Deviazione std
Varianza
Range
Range interquartile
20.15302
406.14433
172.00000
23.00000
Test di posizione: Mu0=19
Test
-Statistica-
-----P-value------
T di Student
Segno
Rango con segno
t
M
S
Pr > |t|
Pr >= |M|
Pr >= |S|
4.100299
-6984
-9.714E7
<.0001
<.0001
<.0001
Quantili (Definizione 5)
Quantile
100% Max
99%
95%
90%
75% Q3
50% Mediana
25% Q1
10%
5%
1%
0% Min
Stima
172
84
64
50
27
11
4
3
3
2
0
17:13 Sunday, November 14, 2010
La procedura UNIVARIATE
Variabile: Tmorto
Osservazioni estreme
---Inferiori---
---Superiori---
Valore
Oss
Valore
Oss
0
30367
152
13618
20
0
1
1
1
10286
56710
54193
52878
154
154
163
172
1294
33509
1298
1295
Valori mancanti
Valore
mancante
Conteggio
.
9854
---Percentuale di--Tutte
Oss
le oss
mancanti
16.77
100.00
17:13 Sunday, November 14, 2010
La procedura UNIVARIATE
Distribuzione normale adattata per Tmorto
Parametri per la distribuzione Normale
Parametro
Simbolo
Media
Dev std
Mu
Sigma
Stima
19.37362
20.15302
Test della bontà di adattamento per la distribuzione Normale
Test
-------Statistica-------
---------P-value---------
Kolmogorov-Smirnov
Cramer-von Mises
Anderson-Darling
D
W-Qu
A-Qu
Pr > D
Pr > W-Qu
Pr > A-Qu
0.19311
589.08988
3309.25883
Quantili per la distribuzione Normale
Percentuale
1.0
5.0
10.0
25.0
50.0
75.0
90.0
95.0
99.0
-------Quantile-----Osservato
Stimato
2.00000
3.00000
3.00000
4.00000
11.00000
27.00000
50.00000
64.00000
84.00000
-27.50932
-13.77515
-6.45352
5.78061
19.37362
32.96663
45.20076
52.52239
66.25656
<0.010
<0.005
<0.005
21
30
Percentuale
25
20
15
10
5
0
0
9
18
27
36
45
54
63
72
81
90
99
108
117
126
135
144
153
162
171
Tmorto
Proviamo ora a utilizzare lo strumento SAS ODS Graphics (Novità della versione 9.2 di SAS). Esso consente di ottenere output
grafici e di testo in molti formati: rtf, html, pdf, ….
Univar_ODS.sas
ODS graphics on;
ods listing close;
ods rtf path="C:\Users\Alex\Desktop\SAS\ODSout";
proc univariate data=tmp1.telefonate_gen2008_2;
var Tmorto;
histogram Tmorto/normal (mu=est sigma=est);
inset normal(mean std) skewness kurtosis; /*crea nel grafico
un riquadro con le stat. specificate*/
run;
ods rtf close;
ods listing;
ods graphics off;
L’output (sul file sasrtf.rtf) e’ il seguente:
Momenti
N
48916 Somma dei pesi
48916
Media
19.3736201 Somma delle osservazioni
947680
Deviazione std
20.1530228 Varianza
406.144326
Skewness
1.65336041 Curtosi
2.74586304
SS non corretta
38226542 SS corretta
19866549.7
Coeff variaz
104.02301 Errore std media
0.09112021
Misure statistiche di base
Posizione
Media
Variabilità
19.37362 Deviazione std
20.15302
Mediana 11.00000 Varianza
Moda
406.14433
3.00000 Range
172.00000
Range interquartile
23.00000
Test di posizione: Mu0=0
Test
Statistica
T di Student
t
Segno
M
P-value
212.6161 Pr > |t|
24457 Pr >= |M| <.0001
5.9816E8 Pr >= |S|
Rango con segno S
Quantili (Definizione 5)
Quantile
Stima
172
100% Max
99%
84
95%
64
90%
50
75% Q3
27
50% Mediana
11
25% Q1
4
10%
3
5%
3
1%
2
0% Min
0
Osservazioni estreme
Inferiori
Valore
<.0001
Superiori
Oss Valore
Oss
0 30367
152 13618
0 10286
154
1 56710
154 33509
1 54193
163
1298
1 52878
172
1295
1294
<.0001
Valori mancanti
Percentuale di
Valore
Tutte
Oss
mancante Conteggio le oss mancanti
.
9854
16.77
100.00
26
Parametri per la distribuzione
Normale
Parametro Simbolo
Stima
Media
Mu
19.37362
Dev std
Sigma
20.15302
Test della bontà di adattamento per la distribuzione Normale
Test
Statistica
P-value
0.19311 Pr > D
Kolmogorov-Smirnov D
Cramer-von Mises
W-Qu
Anderson-Darling
A-Qu
<0.010
589.08988 Pr > W-Qu <0.005
3309.25883 Pr > A-Qu
<0.005
Quantili per la distribuzione
Normale
Quantile
Percentuale Osservato
Stimato
1.0
2.00000 -27.50932
5.0
3.00000 -13.77515
10.0
3.00000
-6.45352
25.0
4.00000
5.78061
50.0
11.00000
19.37362
75.0
27.00000
32.96663
90.0
50.00000
45.20076
95.0
64.00000
52.52239
99.0
84.00000
66.25656
Confrontiamo ora la variabile ImpegnoOper dei mesi di maggio e agosto per vedere se ci sono differenze nella
distribuzione.
Concentriamoci sulla parte grafica, utilizzando ODS:
UNIVAR_ODS2.sas
/*creo un dataset contenente i dati di maggio e agosto */
data tmp1.combined_genmag;
set tmp1.telefonate_mag2008 tmp1.telefonate_ago2008;
run;
/*affianchiamo due grafici per mese col comando class
per
ODS
ods
ods
vedere se ci sono differenze sull'impegno operatore sui due mesi*/
27
graphics on;
listing close;
rtf path="C:\Users\Alex\Desktop\SAS\ODSout"
file='univar_ods2.rtf'; /*così l'output è su questo file*/
title 'Impegno operatore maggio-agosto';
proc univariate data=tmp1.combined_genmag;
var ImpegnoOper;
class mese;
histogram ImpegnoOper/normal (mu=est sigma=est) kernel(color=red);
inset normal(mean std) skewness kurtosis /position=ne; /*crea nel grafico
un riquadro con le stat. specificate e lo mette a nord-est della figura*/
run;
ods rtf close;
ods listing;
ods graphics off;
Oltre alle analisi statistiche, otteniamo il seguente output:
La linea tratteggiata rossa è una Kernel density estimation della p.d.f. della variabile aleatoria: è un metodo
alternativo all’istogramma per rappresentare una stima della p.d.f.; la si effettua per distribuzioni CONTINUE.
Essa consente di visualizzare meglio le differenze con altre distribuzioni continue come la normale.
Si noti che stranamente ci sono dei dati negativi (errori nel database!).
Proviamo ora a sovrapporre grafici della distribuzione gamma, e fissiamo i midpoints degli istogrammi:
UNIVAR_ODS3.sas
/*creo un dataset ripulito dai dati "strani" */
data tmp1.combined_genmag_2;
set tmp1.combined_genmag;
where ImpegnoOper>0;
run;
/*fisso i midpoints e sovrappongo una gamma*/
ODS graphics on;
ods listing close;
ods rtf path="C:\Users\Alex\Desktop\SAS\ODSout"
file='univar_ods3.rtf'; /*così l'output è su questo file*/
title 'Impegno operatore maggio-agosto';
proc univariate data=tmp1.combined_genmag_2;
var ImpegnoOper;
class mese;
histogram ImpegnoOper/gamma (alpha=est sigma=est) kernel(color=red)
midpoints = 30 60 90 120 150 180 210 240;
histogram ImpegnoOper/exponential (theta=0 sigma=est) kernel(color=red)
midpoints = 30 60 90 120 150 180 210 240;
run;
ods rtf close;
ods listing;
ods graphics off;
Dei risultati riportiamo qui solo i grafici e le statistiche relative ai test non
parametrici:
28
29
Test della bontà di adattamento per la distribuzione Gamma
Test
Statistica
Kolmogorov-Smirnov D
Cramer-von Mises
P-value
0.0499011 Pr > D
<0.001
W-Qu 29.1419881 Pr > W-Qu <0.001
30
Test della bontà di adattamento per la distribuzione
Esponenziale
Test
Statistica
Kolmogorov-Smirnov D
ESERCIZIO 3
P-value
0.098396 Pr > D
<0.001
Cramer-von Mises
W-Qu 138.655520 Pr > W-Qu <0.001
Anderson-Darling
A-Qu
717.430049 Pr > A-Qu
<0.001
Si consideri il SAS Data set tel_gen08_tattesa , contenente i dati relativi alle telefonate arrivate al 31
118 di Milano nel gennaio 2008, a cui sono state aggiunte la variabili
tattesa, contenente il tempo di attesa fra ciascuna telefonata e la precedente
tattesa_cum, contenente il tempo di attesa cumulato (ogni riga contiene la somma di tutte le righe
precedenti).
Nel dataset i dati sono ordinati per DT_INLINE crescente.
Al fine di verificare se il processo di arrivo delle telefonate sia un processo di Poisson, verificare se
1. I tempi di attesa fra una telefonata e la successiva abbiano distribuzione esponenziale
2. I tempi di attesa fra una telefonata e la successiva NEI SOLI WEEKEND abbiano distribuzione
esponenziale
FACOLTATIVI (aiutarsi con Excel per la creazione dei datasets):
3. I tempi di attesa fra un codice rosso e il successivo abbiano distribuzione esponenziale
4. I tempi di attesa fra un codice giallo e il successivo abbiano distribuzione esponenziale
PROCEDURA TTest:
Confronto di medie
32
Test T tra campioni indipendenti
Durante un esperimento un gruppo di animali da laboratorio viene
diviso in due sottogruppi, di cui uno viene trattato con un
anestetico e l’altro non viene trattato (1=trattato, 0=non trattato).
Si misura in ogni animale il valore del flusso sanguigno nei reni e
ci si chiede se l’anestetico influenzi o meno tale valore.
animale
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
flusso_sangue
2,35
2,55
1,95
2,79
3,21
2,97
3,44
2,58
2,66
2,31
3,43
2,37
1,82
2,98
2,53
2
1,71
2,22
2,71
1,83
2,14
3,72
2,1
2,58
1,32
3,7
1,59
2,07
2,15
2,05
TTEST_INDIP.SAS
trattato
0
0
0
1
0
1
1
0
1
1
1
1
0
0
0
0
0
0
0
1
1
1
1
0
1
0
0
1
1
1
ODS graphics on; /*attivo la grafica*/
ods listing close;
ods rtf path="C:\Users\Alex\Desktop\SAS\ODSout"
file='ttest_indip.rtf'; /*così l'output è su questo file*/
proc ttest data=tmp1.anestetico alpha=0.01 /*livello di sign.
per gli intervalli di fiducia*/
ci=equal /*usa due code uguali per costruire int. di fiducia
per la varianza; oppure: ci=UMPU o ci=none */
cochran /*calcola il T Test anche per varianze non uguali*/
H0=0;
class trattato;
var flusso_sangue;
run;
ods rtf close;
ods listing;
Otteniamo il seguente output:
trattato
N Media Dev std Err std Minimo Massimo
0
15
2.4320
0.5815
0.1501
1.5900
3.7000
1
15
2.4900
0.6687
0.1727
1.3200
3.7200
-0.0580
0.6266
0.2288
Diff (1-2)
Media
Media CL al
99%
Dev std
0
2.4320
1.9851 2.8789
0.5815 0.3888 1.0778
1
2.4900
1.9760 3.0040
0.6687 0.4471 1.2395
-0.0580 -0.6902 0.5742
0.6266 0.4643 0.9392
trattato
Metodo
Diff (1-2) Aggregati
Dev std CL
al 99%
Diff (1-2) Satterthwaite -0.0580 -0.6911 0.5751
Metodo
Varianze
DF Valore t Pr > |t|
Aggregati
Uguale
28
-0.25
0.8017
Satterthwaite Diverso
27.47
-0.25
0.8018
Uguaglianza di varianze
Metodo
Folded F
DF num DF den
14
14
Valore
F Pr > F
1.32
0.6081
33
34
35
Per verificare la normalita’ dei dati, applico anche una procedura univariate, facendo calcolare i test
non parametrici di fitting con la normale:
Test_normality_indip.sas
proc univariate data=tmp1.anestetico; /*verifico con test non parametrici
se la distribuzione dei dati è gaussiana*/
var flusso_sangue;
histogram flusso_sangue/normal (mu=est sigma=est);
class trattato;
run;
quit;
Riporto solo i risultati dei test non parametrici:
Trattato=0
Test della bontà di adattamento per la distribuzione Normale
Test
-------Statistica-------
---------P-value---------
Kolmogorov-Smirnov
Cramer-von Mises
Anderson-Darling
D
W-Qu
A-Qu
Pr > D
Pr > W-Qu
Pr > A-Qu
0.13287678
0.03683466
0.24120398
>0.150
>0.250
>0.250
Trattato=1
Test della bontà di adattamento per la distribuzione Normale
Test
-------Statistica-------
---------P-value---------
Kolmogorov-Smirnov
Cramer-von Mises
Anderson-Darling
D
W-Qu
A-Qu
Pr > D
Pr > W-Qu
Pr > A-Qu
0.17121131
0.07951635
0.45065022
>0.150
0.202
0.242
In entrambi i casi possiamo accettare l’ipotesi di normalita’, sia dei trattati che dei non trattati, e
possiamo quindi applicare il t-test anche se il campione non e’ molto numeroso.
36
Test T tra campioni accoppiati
Supponiamo ora che sugli stessi animali venga misurato il flusso del sangue prima e dopo aver loro
somministrato l’anestetico. Si vuole verificare se il flusso subisca un cambiamento significativo.
animale
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
pre_anestetico dopo_anestetico
2,35
2
2,55
1,71
1,95
2,22
2,79
2,71
3,21
1,83
2,97
2,14
3,44
3,72
2,58
2,1
2,66
2,58
2,31
1,32
3,43
3,7
2,37
1,59
1,82
2,07
2,98
2,15
2,53
2,05
In tal caso i due campioni non sono indipendenti e si dovra’ eseguire un T test accoppiato
TTEST_ACCOPP.SAS
ODS graphics on; /*attivo la grafica*/
ods listing close;
ods rtf path="C:\Users\Alex\Desktop\SAS\ODSout"
file='ttest_accopp.rtf'; /*così l'output è su questo file*/
proc ttest data=tmp1.anestetico_accoppiati alpha=0.05;
paired pre_anestetico*dopo_anestetico;
run;
ods rtf close;
ods listing;
N Media Dev std Err std Minimo Massimo
15 0.4033
0.5358
0.1383
-0.2800
1.3800
Media CL al
Media
95%
Dev std
0.4033 0.1066 0.7000
Dev std CL
al 95%
0.5358 0.3923 0.8450
DF Valore t Pr > |t|
14
2.92
0.0113
Grafici per studiare le discrepanze fra i due campioni:
37
38
Grafici per studiare la normalità dei dati:
39
Questi grafici rivelano già qualitativamente un buon adattamento con la normale. Anche qui, per
sicurezza, dato che il campione è piccolo e non si può applicare il TLC, si potrebbero effettuare dei test
non parametrici di adattamento della variabile differenza con la normale:
test_normality_accopp.sas
data work.testdiff;
set tmp1.anestetico_accoppiati;
differenza=pre_anestetico-dopo_anestetico;
run;
proc univariate data=work.testdiff; /*verifico con test non parametrici
se la distribuzione dei dati è gaussiana*/
var differenza;
histogram differenza/normal (mu=est sigma=est);
run;
quit;
Risultati dei test di adattamento:
Test della bontà di adattamento per la distribuzione Normale
Test
-------Statistica-------
---------P-value---------
Kolmogorov-Smirnov
Cramer-von Mises
Anderson-Darling
D
W-Qu
A-Qu
Pr > D
Pr > W-Qu
Pr > A-Qu
0.15897741
0.06927057
0.48401877
>0.150
>0.250
0.203
40
Di nuovo i p-value sono alti, per cui possiamo accettare l’ipotesi di normalità dei dati.
ESERCIZIO 4
Si considerino i dati sulle telefonate al servizio 118 di Milano.
Attraverso opportuni test di ipotesi si dia una risposta alle
seguenti domande:
1. Nelle telefonate del mese di gennaio ci sono differenze
significative fra i tempi medi di “impegno operatore” fra i
codici rossi e i codici verdi?
2. Si considerino le sole telefonate relative ad INCIDENTI
STRADALI:
a. Il numero medio giornaliero di incidenti di tipo auto-moto
o caduta da moto e’ significativamente diverso nei mesi di
gennaio e di maggio?
b. Il numero medio giornaliero di incidenti di gennaio di tipo
auto-moto o caduta da moto e’ significativamente diverso
nei giorni lavorativi e nei weekend?
3. Si consideri il SAS Data set
tel_gen08_tattesa, dell’esercizio
precedente.
a) I tempi di attesa fra una telefonata e l’altra hanno
medie significativamente diverse nei giorni lavorativi e
nei weekend?
b) I tempi attesa fra una telefonata e l’altra nei giorni
lavorativi hanno medie significativamente diverse nelle
ore diurne (fra le 7 e le 18) e notturne (fra le 19 e le
6)?
TEST CHI QUADRO
Procedura FREQ:tabelle di contingenza e confronti di proporzioni
Utilizzo i dati delle telefonate al 118 nel mese di Gennaio. Voglio verificare se vi sia una differenza
significativa fra la distribuzione dei motivi di chiamata fra giorni lavorativi e week end.
Programma chisq.sas
data tmp1.telefonate_gen2008_3;
set tmp1.telefonate_gen2008;
where DS_MOTIVO and DS_MOTIVO ne 'INC FERROVIA'; /*seleziona solo le righe in cui
DS_MOTIVO non e' vuoto.
Gli zeri e i vuoti sono FALSE per where, gli altri valori sono TRUE.
Escludo inoltre gli incidenti ferrovia (e' uno solo)*/
if weekday='sab' or weekday='dom' then daylabel='festivo';
else daylabel='feriale'; /* creo una nuova variabile che mi dice se il giorno
e' lavorativo o un week end */
run;
ODS graphics on; /*attivo la grafica*/
ods listing close;
ods rtf path="C:\Users\Alex\Desktop\SAS\ODSout"
file='chisq.rtf'; /*così l'output è su questo file*/
proc freq data=tmp1.telefonate_gen2008_3;
tables DS_MOTIVO*daylabel / chisq expected deviation; /*effettuo
un test chi quadro, plottando anche le frequenze attese
e le differenze f_osservate - f_attese*/
run;
ods rtf close;
ods listing;
41
42
Tabella di DS_MOTIVO per daylabel
DS_MOTIVO(DS_MOTIVO)
Frequenza
Atteso
Deviazione
Percentuale
Pct riga
Pct col
ALTRO
ANIMALI
CADUTA
EVENTO DI MASSA
EVENTO VIOLENTO
INC INFORTUNIO
INC STRADALE
daylabel
feriale
festivo
Totale
42
40.347
1.6534
0.17
77.78
0.23
12
13.653
-1.653
0.05
22.22
0.19
54
12
9.7131
2.2869
0.05
92.31
0.07
1
3.2869
-2.287
0.00
7.69
0.02
1867
1802.1
64.851
7.63
77.40
10.21
545
609.85
-64.85
2.23
22.60
8.81
2412
31
32.875
-1.875
0.13
70.45
0.17
13
11.125
1.875
0.05
29.55
0.21
44
298
323.52
-25.52
1.22
68.82
1.63
135
109.48
25.52
0.55
31.18
2.18
362
377.32
-15.32
1.48
71.68
1.98
143
127.68
15.316
0.58
28.32
2.31
1488
1464.4
23.567
6.08
75.92
8.14
472
495.57
-23.57
1.93
24.08
7.63
0.22
13
0.05
9.86
0.18
433
1.77
505
2.06
1960
8.01
43
Tabella di DS_MOTIVO per daylabel
DS_MOTIVO(DS_MOTIVO)
Frequenza
Atteso
Deviazione
Percentuale
Pct riga
Pct col
INTOSSICAZIONE
MEDICO ACUTO
NON NOTO
PREVENZIONE
SOCCORSO PERSONA
Totale
daylabel
feriale
festivo
Totale
380
394.5
-14.5
1.55
71.97
2.08
148
133.5
14.5
0.60
28.03
2.39
528
13315
13344
-29.27
54.41
74.55
72.83
4545
4515.7
29.274
18.57
25.45
73.46
17860
264
275.7
-11.7
1.08
71.54
1.44
105
93.298
11.702
0.43
28.46
1.70
369
21
21.668
-0.668
0.09
72.41
0.11
8
7.3324
0.6676
0.03
27.59
0.13
203
196.5
6.497
0.83
77.19
1.11
60
66.497
-6.497
0.25
22.81
0.97
18283
74.72
6187
25.28
2.16
72.99
1.51
29
0.12
263
1.07
24470
100.00
Statistiche per la tabella di DS_MOTIVO per daylabel
tatistica
DF Valore
Prob
Chi-quadrato
11 29.2284 0.0021
Chi-quadrato rapp verosim
11 29.5326 0.0019
Chi-quadrato MH
1
4.5698 0.0325
Coefficiente Phi
0.0346
Coefficiente di contingenza
0.0345
V di Cramer
0.0346
Dimensione campione = 24470
Posso quindi dire che c’e’ una differenza significativa, se uso un alpha=0.002 o maggiore. Dagli istogrammi si evidenzia che la differenza principale sta nella classe “Medico Acuto”: le chiamate per motivi
medici sono piu’ frequenti nei giorni feriali che nei festivi.
Applichiamo ora il test chi quadro a una variabile continua, che andra’ pertanto discretizzata con un “format”:
verifichiamo se la distribuzione della variabile “impegno operatore” sia significativamente diversa fra i giorni feriali o
festivi.
Programma chisq_cont.sas
data work.telefonate_gen2008_4;
set tmp1.telefonate_gen2008;
if weekday='sab' or weekday='dom' then daylabel='festivo';
else daylabel='feriale'; /* creo una nuova variabile che mi dice se il giorno
e' lavorativo o un week end */
run;
proc format; /* definisce range di valori:
serve per raggruppare o discretizzare i dati,
associando un'etichetta a ogni gruppo*/
value durata low-30='A meno di mezzo minuto' /*nell'output elenca le etichette in ordine
lessicografico */
30-60='B mezzo-un minuto'
60-120='C uno-due minuti'
120-300='D due-cinque minuti'
300-high='E più di 5 minuti';
run;
ODS graphics on; /*attivo la grafica*/
ods listing close;
ods rtf path="C:\Users\Alex\Desktop\SAS\ODSout"
file='chisq_cont.rtf'; /*così l'output è su questo file*/
proc freq data=work.telefonate_gen2008_4;
tables daylabel*ImpegnoOper / chisq expected deviation cellchi2
plots=none; /*effettuo
un test chi quadro, plottando anche le frequenze attese,
le differenze f_osservate - f_attese, il contributo alla statistica
chi2 dato da ciascuna cella.*/
format ImpegnoOper durata.;
run;
ods rtf close;
ods listing;
OUTPUT:
Tabella di daylabel per ImpegnoOper
daylabel
ImpegnoOper(ImpegnoOper)
Frequenza
Atteso
Deviazione
A
ChiQuadrato cella
meno
Percentuale
di B mezzoPct riga
mezzo
un
Pct col
minuto
minuto
feriale
festivo
Totale
C unodue
minuti
D due- E più
cinque
di 5
minuti minuti Totale
10189
10353
-163.9
2.596
20.83
28.64
71.59
10614
10634
-20.44
0.0393
21.70
29.83
72.60
11177
11046
130.86
1.5502
22.85
31.41
73.60
3402
3357.6
44.369
0.5863
6.95
9.56
73.70
199
189.85
9.1513
0.4411
0.41
0.56
76.25
35581
4044
3880.1
163.94
6.9267
8.27
30.33
28.41
4006
3985.6
20.439
0.1048
8.19
30.04
27.40
4009
4139.9
-130.9
4.1363
8.20
30.06
26.40
1214
1258.4
-44.37
1.5644
2.48
9.10
26.30
62
71.151
-9.151
1.177
0.13
0.46
23.75
13335
14233
29.10
14620
29.89
15186
31.05
4616
9.44
261
0.53
48916
100.00
Frequenza mancanti = 9854
Statistiche per la tabella di daylabel per ImpegnoOper
Statistica
DF Valore
Prob
Chi-quadrato
4 19.1222 0.0007
Chi-quadrato rapp verosim
4 19.1417 0.0007
Chi-quadrato MH
1 15.2516 <.0001
Coefficiente Phi
0.0198
Coefficiente di contingenza
0.0198
V di Cramer
0.0198
72.74
27.26
Dimensione effettiva campione = 48916
Frequenza mancanti = 9854
WARNING: 17% dei dati è mancante.
Dunque ci sono differenze significative fra giorni feriali e festivi e le differenze maggiori sono dovute a un aumento
rispetto ai valori attesi delle telefonate festive che durano meno di mezzo minuto, con conseguente diminuzione delle
telefonate piu’ lunghe.
50
SAS
La procedura FREQ
ID_CODICE=' '
Confronto fra una proporzione sperimentale e una teorica CASO 1: utilizzo dell’approssimazione asintotica della binomiale con la normale
Mi chiedo se la proporzione di telefonate che durano fino a due minuti sia pari al 80% del totale, indistinto rispetto al
giorno della settimana, ossia verifico
H0: prop. telefonate con durata fino a 2 min = 0.8
prop. telefonate con durata superiore a 2 min = 0.2
contro tutte le alternative.
Programma conf_prop_bin1.sas
data work.telefonate_gen2008_5;
set tmp1.telefonate_gen2008;
if ImpegnoOper le 120 then telefonata='corta';
else telefonata='lunga'; /* creo una nuova variabile che mi indica la durata
la durata della telefonata */
run;
ODS graphics on; /*attivo la grafica*/
ods listing close;
ods rtf path="C:\Users\Alex\Desktop\SAS\ODSout"
file='bin1.rtf'; /*così l'output è su questo file*/
proc freq data=work.telefonate_gen2008_5;
tables telefonata / binomial (p=0.8) plots(only)=freqplot(type=dot scale=percent);
run;
ods rtf close;
ods listing;
OUTPUT:
Frequenza Percentuale
telefonata Frequenza Percentuale cumulativa cumulativa
corta
53893
91.70
53893
91.70
lunga
4877
8.30
58770
100.00
51
SAS
La procedura FREQ
ID_CODICE=' '
Proporzione binomiale per telefonata =
corta
Proporzione
0.9170
ASE
0.0011
Limite conf inferiore al 95%
0.9148
Limite conf superiore al 95%
0.9192
Limiti conf esatti
Limite conf inferiore al 95%
0.9148
Limite conf superiore al 95%
0.9192
52
SAS
La procedura FREQ
ID_CODICE=' '
Test di H0: Proporzione =
0.8
ASE con H0
Z
0.0016
70.9188
Pr unilaterale > Z
<.0001
Pr bilaterale > |Z|
<.0001
Dimensione campione = 58770
Dunque la proporzione di telefonate di durata inferiore ai due minuti è significativamente maggiore di 0.8. Si noti però
che a causa della grande dimensione del campione gli intervalli di fiducia asintotico (con approssimazione della
binomiale con la normale) ed esatto sono indistinguibili.
Eseguiamo lo stesso test sulle telefonate divise per codice:
Programma conf_prop_bin2.sas
data work.telefonate_gen2008_5;
set tmp1.telefonate_gen2008;
if ImpegnoOper le 120 then telefonata='corta';
else telefonata='lunga'; /* creo una nuova variabile che mi indica la durata
la durata della telefonata */
run;
proc sort data=work.telefonate_gen2008_5;
by ID_codice;
run;
ODS graphics on; /*attivo la grafica*/
ods listing close;
ods rtf path="C:\Users\Alex\Desktop\SAS\ODSout"
file='bin2.rtf'; /*così l'output è su questo file*/
proc freq data=work.telefonate_gen2008_5;
tables telefonata / binomial (p=0.8) plots=none;
by id_codice;
run;
ods rtf close;
ods listing;
OUTPUT:
53
SAS
La procedura FREQ
ID_CODICE=' '
Frequenza Percentuale
telefonata Frequenza Percentuale cumulativa cumulativa
corta
16823
96.90
16823
96.90
lunga
539
3.10
17362
100.00
Proporzione binomiale per telefonata =
corta
Proporzione
0.9690
ASE
0.0013
Limite conf inferiore al 95%
0.9664
Limite conf superiore al 95%
0.9715
Limiti conf esatti
Limite conf inferiore al 95%
0.9663
Limite conf superiore al 95%
0.9715
Test di H0: Proporzione =
0.8
ASE con H0
Z
0.0030
55.6559
Pr unilaterale > Z
<.0001
Pr bilaterale > |Z|
<.0001
Dimensione campione = 17362
05:21 martedì, giugno 07, 2011 54
SAS
La procedura FREQ
ID_CODICE=B
Frequenza Percentuale
telefonata Frequenza Percentuale cumulativa cumulativa
corta
10033
86.21
10033
86.21
lunga
1605
13.79
11638
100.00
Proporzione binomiale per telefonata =
corta
Proporzione
0.8621
ASE
0.0032
Limite conf inferiore al 95%
0.8558
Limite conf superiore al 95%
0.8684
Limiti conf esatti
Limite conf inferiore al 95%
0.8557
Limite conf superiore al 95%
0.8683
Test di H0: Proporzione =
0.8
ASE con H0
Z
0.0037
16.7455
Pr unilaterale > Z
<.0001
Pr bilaterale > |Z|
<.0001
Dimensione campione = 11638
05:21 martedì, giugno 07, 2011 55
SAS
La procedura FREQ
ID_CODICE=G
Frequenza Percentuale
telefonata Frequenza Percentuale cumulativa cumulativa
corta
9646
92.20
9646
92.20
lunga
816
7.80
10462
100.00
Proporzione binomiale per telefonata =
corta
Proporzione
0.9220
ASE
0.0026
Limite conf inferiore al 95%
0.9169
Limite conf superiore al 95%
0.9271
Limiti conf esatti
Limite conf inferiore al 95%
0.9167
Limite conf superiore al 95%
0.9271
Test di H0: Proporzione =
0.8
ASE con H0
Z
0.0039
31.1975
Pr unilaterale > Z
<.0001
Pr bilaterale > |Z|
<.0001
Dimensione campione = 10462
05:21 martedì, giugno 07, 2011 56
SAS
La procedura FREQ
ID_CODICE=R
Frequenza Percentuale
telefonata Frequenza Percentuale cumulativa cumulativa
corta
1505
86.79
1505
86.79
lunga
229
13.21
1734
100.00
Proporzione binomiale per telefonata =
corta
Proporzione
0.8679
ASE
0.0081
Limite conf inferiore al 95%
0.8520
Limite conf superiore al 95%
0.8839
Limiti conf esatti
Limite conf inferiore al 95%
0.8511
Limite conf superiore al 95%
0.8835
Test di H0: Proporzione =
0.8
ASE con H0
0.0096
Z
7.0723
Pr unilaterale > Z
<.0001
Pr bilaterale > |Z|
<.0001
Dimensione campione = 1734
05:21 martedì, giugno 07, 2011 57
SAS
La procedura FREQ
ID_CODICE=V
Frequenza Percentuale
telefonata Frequenza Percentuale cumulativa cumulativa
corta
7724
85.19
7724
85.19
lunga
1343
14.81
9067
100.00
Proporzione binomiale per telefonata =
corta
Proporzione
0.8519
ASE
0.0037
Limite conf inferiore al 95%
0.8446
Limite conf superiore al 95%
0.8592
Limiti conf esatti
Limite conf inferiore al 95%
0.8444
Limite conf superiore al 95%
0.8591
Test di H0: Proporzione =
0.8
ASE con H0
Z
0.0042
12.3502
Pr unilaterale > Z
<.0001
Pr bilaterale > |Z|
<.0001
Dimensione campione = 9067
Frequenza Percentuale
Frequenza Percentuale cumulativa cumulativa
corta
8162
95.94
8162
95.94
lunga
345
4.06
8507
100.00
Proporzione binomiale per telefonata =
corta
Proporzione
0.9594
ASE
0.0021
Limite conf inferiore al 95%
0.9553
Limite conf superiore al 95%
0.9636
Limiti conf esatti
Limite conf inferiore al 95%
0.9550
Limite conf superiore al 95%
0.9635
Test di H0: Proporzione =
0.8
ASE con H0
Z
0.0043
36.7654
Pr unilaterale > Z
<.0001
Pr bilaterale > |Z|
<.0001
Dimensione campione = 8507
_______________________________________________________________________________
Le percentuali sono ancora tutte significativamente maggiori di 0.8, anche se per alcuni codici esse scendono. Si noti
inoltre che a causa delle piu’ piccole dimensioni dei singoli campioni, gli intervalli di fiducia con l’approssimazione
normale e la distribuzione binomiale esatta cambiano, anche se di poco.
CASO 2: utilizzo del test chi quadro per il confronto fra una proporzione teorica e una sperimentale
Supponiamo ad esempio di voler verificare se l’arrivo di telefonate al centralino del 118 nel mese di gennaio segue un processo di Poisson a
parametro costante.
Programma Confr_distr_teorica1.sas
Estraggo la variabile di interesse e la discretizzo
/*conto il numero di telefonate al giorno in gennaio */
proc freq data=tmp1.telefonate_gen2008;
tables giorno/out=work.freq_giorno_gen;
run;
Calcolo alcune statistiche descrittive che mi sono utili per l’analisi, in particolare per determinare il parametro della
distribuzione di Poisson di riferimento
/*calcolo la media dei conteggi*/
proc means data=work.freq_giorno_gen mean min max var;
var count/;
run;
/* verifico ora se il numero di telefonate al giorno ha
distribuzione di Poisson a parametro uguale alla media calcolata
(circa 1896)*/
/*calcolo le proporzioni sperimentali, discretizzando il range di
count*/
data work.prop;
set work.freq_giorno_gen;
/*discretizzo la variabile count */
if count<=1870 then partition='I1';
if count>1870 and count<=1900 then partition='I2';
if count>1900 and count<=1930 then partition='I3';
if count>1930 then partition='I4';
run;
Ora calcolo le proporzioni teoriche che dovrei osservare se fosse vera
H0: i dati hanno distribuz. di Poisson a parametro 1896
per la discretizzazione scelta. Leggero’ i risultati nel dataset work.prop_teorica.
ATTENZIONE: devo esprimere le proporzioni come percentuali intere, e la loro somma deve essere 100. Arrotondare
in modo opportuno! (ceil arrotonda verso –inf e floor verso +inf ).
data work.prop_teorica;
/* Calcolo le proporzioni teoriche in un dataset diverso */
/* sintassi: prob=CDF('POISSON',n,l) con l=media della Poisson, n=punto in
cui calcolo la cdf.
Ossia: CDF('POISSON',n,l) calcola la probabilità che una v.a. di Poisson a
parametro l sia <=n */
array p{4} p1 p2 p3 p4;
p{1}=floor(CDF('POISSON',1870,1896)*100);
p{2}=floor(CDF('POISSON',1900,1896)*100)-floor(CDF('POISSON',1870,1896)*100);
p{3}=floor(CDF('POISSON',1930,1896)*100)-floor(CDF('POISSON',1900,1896)*100);
p{4}=100-floor(CDF('POISSON',1930,1896)*100);
run;
Dopo aver letto i risultati dal dataset work.prop_teorica posso applicare la procedura freq ai dati sperimentali con le
percentuali opportune
ODS graphics on; /*attivo la grafica*/
ods listing close;
ods rtf path="C:\Users\Alex\Desktop\SAS\ODSout"
file='prop_teo.rtf'; /*così l'output è su questo file*/
proc freq data=work.prop;
tables partition/ chisq testp=(27,27,24,22)
plots(only)=deviationplot(type=dot);
run;
ods rtf close;
ods listing;
OUTPUT:
Test sulla Frequenza Percentuale
partition Frequenza Percentuale percentuale cumulativa cumulativa
I1
16
51.61
27.00
16
51.61
I2
4
12.90
27.00
20
64.52
I3
3
9.68
24.00
23
74.19
I4
8
25.81
22.00
31
100.00
Test del chi-quadrato
per proporzioni specificate
Chi-quadrato
DF
Pr > ChiQuadr
12.0909
3
0.0071
Dimensione campione = 31
Dunque la distribuzione dei dati e’ significativamente diversa da una Poisson(1896).
CASO 3: utilizzo del test chi quadro per tabelle di contingenza con conteggi già forniti
Vediamo ora come si può effettuare un test chi quadro quando i dati contengono già i conteggi delle celle della tabella di contingenza. Si
consideri il seguente esempio, i cui dati si riferiscono a un sondaggio effettuato in Italia in un campione di bambini delle scuole elementari (file
animali_bambini.xls)
Censimento a scuola: scuole elementari
Quale animale/i possiedi?
Dimensione Comune: Metropolitano
Animali posseduti
Gatto e
Cane e
Cane e
Cane e
Nessuno Altro Uccellino Pesciolino Gatto pesciolino Cane uccellino pesciolino gatto
Num.
Num.
Num.
Num.
Num.
Num.
Num.
Num.
Num.
Num.
Ripartizione
Nord Ovest
Nord Est
Centro
Sud
Isole
Totale
1164
147
376
683
604
2974
516
90
201
320
253
1380
77
17
41
149
145
429
328
37
93
192
56
706
236
36
81
59
62
474
37
12
8
2
5
64
251
27
105
215
255
853
13
4
8
30
51
106
44
4
17
25
16
106
39
6
16
15
37
113
Dimensione Comune: Non metropolitano
Animali posseduti
Gatto e
Cane e
Cane e
Cane e
Nessuno Altro Uccellino Pesciolino Gatto pesciolino Cane uccellino pesciolino gatto
Num.
Num.
Num.
Num.
Num.
Num.
Num.
Num.
Num.
Num.
Ripartizione
Nord Ovest
Nord Est
Centro
Sud
Isole
Totale
3399 3779
2171 2523
1621 2490
6140 4425
2132 1992
15463 15209
430
323
225
942
522
2442
1148
537
415
1110
388
3598
925
675
576
942
423
3541
219
140
124
116
55
654
1329
865
817
2729
1314
7054
171
148
99
294
167
879
267
179
120
270
102
938
472
329
372
639
259
2071
Vogliamo verificare se la distribuzione degli animali posseduti sia diversa fra bambini che abitano in centri metropolitani o non metropolitani.
Considereremo pertanto solo l’ultima riga di entrambe le tabelle precedenti, e rielaboriamo i dati in modo da poter essere letti da SAS
adeguatamente, dato che i dati forniscono gia’ i conteggi da inserire nelle celle della tabella di contingenza che analizzeremo.
Riscriviamo i dati come segue (foglio 2 dello stesso file excel di prima):
Animali
Nessuno
Altro
Uccellino
dim_comune
metropolitano
metropolitano
metropolitano
frequenze
2974
1380
429
Pesciolino
Gatto
Gatto e
pesciolino
Cane
metropolitano
metropolitano
706
474
metropolitano
metropolitano
64
853
Cane e
uccellino
Cane e
pesciolino
Cane e gatto
Nessuno
Altro
Uccellino
Pesciolino
Gatto
Gatto e
pesciolino
Cane
Cane e
uccellino
Cane e
pesciolino
Cane e gatto
metropolitano
106
metropolitano
106
metropolitano
non
metropolitano
non
metropolitano
non
metropolitano
non
metropolitano
non
metropolitano
non
metropolitano
non
metropolitano
non
metropolitano
non
metropolitano
non
metropolitano
113
15463
15209
2442
3598
3541
654
7054
879
938
2071
Questi sono i dati che leggeremo con SAS, cosi memorizzeremo 3 variabili:
animali, contenente il tipo di animali posseduti
dim_comune, che indica se il comune e’ metropolitano o no
frequenze, che indica i conteggi da inserire nella tabella di contingenza.
Eseguiamo a questo punto il programma SAS seguente:
Chi_quadro_bambini_animali.sas
/* leggo i dati */
PROC IMPORT OUT= WORK.bambini
DATAFILE= "C:\Users\Alex\Desktop\SAS\dati\animali_elementari.xls"
DBMS=EXCEL REPLACE;
SHEET="foglio2$";
GETNAMES=YES;
RUN;
proc freq data=work.bambini;
weight frequenze; /*mi dice che in frequenze e' contenuto il conteggio osservato in ogni
cella */
tables animali*dim_comune / chisq expected cellchi2;
run;
OUTPUT:
La procedura FREQ
Tabella di Animali per dim_comune
Animali(Animali)
dim_comune(dim_comune)
Frequenza
‚
Atteso
‚
ChiQuadrato cella‚
Percentuale
‚
Pct riga
‚
Pct col
‚metropol‚non metr‚ Totale
‚itano
‚opolitan‚
‚
‚o
‚
ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒˆ
Altro
‚
1380 ‚ 15209 ‚ 16589
‚
2024 ‚ 14565 ‚
‚ 204.9 ‚ 28.472 ‚
‚
2.34 ‚ 25.75 ‚ 28.09
‚
8.32 ‚ 91.68 ‚
‚ 19.15 ‚ 29.33 ‚
ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒˆ
Cane
‚
853 ‚
7054 ‚
7907
‚ 964.71 ‚ 6942.3 ‚
‚ 12.935 ‚ 1.7975 ‚
‚
1.44 ‚ 11.94 ‚ 13.39
‚ 10.79 ‚ 89.21 ‚
‚ 11.84 ‚ 13.60 ‚
ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒˆ
Cane e gatto
‚
113 ‚
2071 ‚
2184
‚ 266.46 ‚ 1917.5 ‚
‚ 88.384 ‚ 12.282 ‚
‚
0.19 ‚
3.51 ‚
3.70
‚
5.17 ‚ 94.83 ‚
‚
1.57 ‚
3.99 ‚
ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒˆ
Cane e pesciolin ‚
106 ‚
938 ‚
1044
o
‚ 127.38 ‚ 916.62 ‚
‚ 3.5871 ‚ 0.4985 ‚
‚
0.18 ‚
1.59 ‚
1.77
‚ 10.15 ‚ 89.85 ‚
‚
1.47 ‚
1.81 ‚
ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒˆ
Cane e uccellino ‚
106 ‚
879 ‚
985
‚ 120.18 ‚ 864.82 ‚
‚ 1.6724 ‚ 0.2324 ‚
‚
0.18 ‚
1.49 ‚
1.67
‚ 10.76 ‚ 89.24 ‚
‚
1.47 ‚
1.70 ‚
ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒˆ
Totale
7205
51849
59054
12.20
87.80
100.00
(Continua)
La procedura FREQ
Tabella di Animali per dim_comune
Animali(Animali)
dim_comune(dim_comune)
Frequenza
‚
Atteso
‚
ChiQuadrato cella‚
Percentuale
‚
Pct riga
‚
Pct col
‚metropol‚non metr‚ Totale
‚itano
‚opolitan‚
‚
‚o
‚
ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒˆ
Gatto
‚
474 ‚
3541 ‚
4015
‚ 489.86 ‚ 3525.1 ‚
‚ 0.5134 ‚ 0.0713 ‚
‚
0.80 ‚
6.00 ‚
6.80
‚ 11.81 ‚ 88.19 ‚
‚
6.58 ‚
6.83 ‚
ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒˆ
Gatto e pescioli ‚
64 ‚
654 ‚
718
no
‚ 87.601 ‚ 630.4 ‚
‚ 6.3585 ‚ 0.8836 ‚
‚
0.11 ‚
1.11 ‚
1.22
‚
8.91 ‚ 91.09 ‚
‚
0.89 ‚
1.26 ‚
ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒˆ
Nessuno
‚
2974 ‚ 15463 ‚ 18437
‚ 2249.4 ‚ 16188 ‚
‚ 233.38 ‚ 32.431 ‚
‚
5.04 ‚ 26.18 ‚ 31.22
‚ 16.13 ‚ 83.87 ‚
‚ 41.28 ‚ 29.82 ‚
ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒˆ
Pesciolino
‚
706 ‚
3598 ‚
4304
‚ 525.12 ‚ 3778.9 ‚
‚ 62.307 ‚ 8.6582 ‚
‚
1.20 ‚
6.09 ‚
7.29
‚ 16.40 ‚ 83.60 ‚
‚
9.80 ‚
6.94 ‚
ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒˆ
Uccellino
‚
429 ‚
2442 ‚
2871
‚ 350.28 ‚ 2520.7 ‚
‚ 17.69 ‚ 2.4582 ‚
‚
0.73 ‚
4.14 ‚
4.86
‚ 14.94 ‚ 85.06 ‚
‚
5.95 ‚
4.71 ‚
ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒˆƒƒƒƒƒƒƒƒˆ
Totale
7205
51849
59054
12.20
87.80
100.00
3
La procedura FREQ
Statistiche per la tabella di Animali per dim_comune
Statistica
DF
Valore
Prob
ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒ
Chi-quadrato
9
719.5111
<.0001
Chi-quadrato rapp verosim
9
742.6204
<.0001
Chi-quadrato MH
1
591.8498
<.0001
Coefficiente Phi
0.1104
Coefficiente di contingenza
0.1097
V di Cramer
0.1104
Dimensione campione = 59054
I risultati indicano che vi è una distribuzione significativamente diversa fra comuni metropolitani e non metropolitani, e le maggiori differenze
si rilevano nell’alto numero di bambini di comuni metropolitani che non hanno animali o ne hanno uno piccolo (pesciolino o uccellino) o di
bambini che posseggono animali diversi da quelli elencati nel sondaggio.
ESERCIZIO 5
Si consideri il SAS Data set
tel_gen08_tattesa dell’Esercizio 3.
a) I tempi di attesa fra una telefonata e l’altra hanno distribuzioni
significativamente diverse nei giorni lavorativi e nei weekend?
b) I tempi attesa fra una telefonata e l’altra nei giorni lavorativi
hanno distribuzioni significativamente diverse nelle ore diurne
(fra le 7 e le 18) e notturne (fra le 19 e le 6)?
ESERCIZIO 6
Si consideri il problema affrontato nell’esempio del programma
Confr_distr_teorica1.sas
Se si cambia la discretizzazione della variabile count (conteggio del numero di
telefonate al giorno) cambiano i risultati del test chi quadro?
Esercizio 7
Si consideri il file di dati excel Colazione_bambini.xls, che si riferisce a un
sondaggio effettuato presso alcune scuole elementari italiane sulle abitudini
alimentari dei bambini.
Rielaborando opportunamente le tabelle presenti nei 3 fogli excel del file, si
risponda alla seguente domanda:
ci sono differenze significative di alimentazione fra bambini che vivono in
centri metropolitani o in centri non metropolitani?
LEZIONE 9:
REGRESSIONE e ANALISI DELLA VARIANZA
Regressione: Procedura Reg
Utilizziamo il SAS data set Milano_covar.sas7bdat, contenente il numero di interventi giornalieri effettuati dal servizio
118 di Milano dal 2005 al 2007 sul solo comune di Milano, sia totale che diviso per tipologia o sotto-tipologia di
infortunio. Nel data set sono anche riportati i valori di alcune variabili climatiche (temperatura, pressione, etc.) e di
alcuni inquinanti, nonche’ l’incidenza della sindrome influenzale, misurata come n. di casi segnalati dai medici di base
ogni 1000 assistiti.
Vogliamo dapprima vedere se gli interventi di tipo medico acuto abbiano una dipendenza lineare (usiamo una retta)
dalla temperatura.
REGR1.SAS
ODS graphics on; /*attivo la grafica*/
ods listing close;
ods rtf path="C:\Users\Alex\Desktop\SAS\ODSout"
file='reg1.rtf'; /*così l'output è su questo file*/
proc reg data=tmp1.Milano_covar;
model totale_medico_acuto=temperatureavgf;
run;
ods rtf close;
ods listing;
OUTPUT:
Numero osservazioni lette
1094
Numero osservazioni usate
1046
Numero osservazioni con valori mancanti
48
Analisi della varianza
Somma dei
Media Valore
quadrati quadratica
F Pr > F
Origine
DF
Modello
1
108132
108132
1044
949423
909.40892
Totale corretto 1045
1057555
Errore
Radice dell'MSE
Media dip.
Coeff Var
30.15641 R-quadro
118.90 <.0001
0.1022
255.71033 R-quadro corr 0.1014
11.79319
Stime dei parametri
Variabile
Etichetta
Intercept
Intercept
TemperatureAvgF TemperatureAvgF
Stima dei
Errore
DF parametri standard Valore t Pr > |t|
1
275.18809
2.01497
136.57
<.0001
1
-1.28939
0.11825
-10.90
<.0001
7
7
Vogliamo ora verificare se vi sia una dipendenza del numero di interventi di tipo medico acuto da
alcune variabili climatiche, alcuni inquinanti, l’incidenza dell’influenza e il giorno_n (ossia se ci
sia evidenza di un trend crescente o decrescente sui 3 anni). Vogliamo inoltre selezionare solo le
variabili piu’ significative fra quelle elencate, attraverso 4 diversi tipi di metodi di selezione delle
variabili: forward, backward, stepwise, maxR.
REGR2.SAS
ODS graphics on; /*attivo la grafica*/
ods listing close;
ods rtf path="C:\ODSout"
file='reg2.rtf'; /*si vedano i risultati su questo file*/
proc reg data=tmp1.Milano_covar;
model totale_medico_acuto=temperatureavgf
humidityavg
pressure_mean
PrecipitationSumIn
giorno_n
INFLUENZA_incidenza_giornaliera
mediadimedia_di_ozono
CO NO NO2 PM10
/ selection=forward slentry=0.1;
model totale_medico_acuto=temperatureavgf
humidityavg
pressure_mean
PrecipitationSumIn
giorno_n
7
INFLUENZA_incidenza_giornaliera
mediadimedia_di_ozono
CO NO NO2 PM10
/ selection=backward slstay=0.1;
model totale_medico_acuto=temperatureavgf
humidityavg
pressure_mean
PrecipitationSumIn
giorno_n
INFLUENZA_incidenza_giornaliera
mediadimedia_di_ozono
CO NO NO2 PM10
/ selection=stepwise slentry=0.05
slstay=0.1;
model totale_medico_acuto=temperatureavgf
humidityavg
pressure_mean
PrecipitationSumIn
giorno_n
INFLUENZA_incidenza_giornaliera
mediadimedia_di_ozono
CO NO NO2 PM10
/ selection=maxr;
run;
ods rtf close;
ods listing;
Per l’OUTPUT si veda il file reg2.rtf
7
ANALISI DELLA VARIANZA: Procedura ANOVA
In due laboratori vengono eseguite delle misurazioni su un farmaco in 12 giorni diversi, con i
seguenti risultati (FILE laboratori1.xls):
misura
99,57
98,33
100,13
100,29
98,85
101,19
101,19
99,96
100,33
100,17
99,81
100,72
99,91
102,68
100,36
100,61
101,57
100,56
100,40
99,67
100,79
99,16
101,21
102,12
Lab
1
1
1
1
1
1
1
1
1
1
1
1
2
2
2
2
2
2
2
2
2
2
2
2
giorno
1
2
3
4
5
6
7
8
9
10
11
12
1
2
3
4
5
6
7
8
9
10
11
12
Verifichiamo se la misura e’ influenzata o meno dal laboratorio in cui e’ stata eseguita
ANOVA1.SAS
ODS graphics on; /*attivo la grafica*/
ods listing close;
ods rtf path="C:\Users\Alex\Desktop\SAS\ODSout"
file='anova1.rtf'; /*così l'output è su questo file*/
proc anova data=tmp1.Laboratori1;
class Lab; /*variabile categorica*/
model misura=Lab;
run;
ods rtf close;
ods listing;
OUTPUT:
7
Informazioni sui
livelli di
classificazione
Classe Livelli Valori
Lab
2 12
Numero osservazioni lette
1529
Numero osservazioni usate
Origine
DF
Modello
1
24
Somma dei
Media Valore
quadrati quadratica
F Pr > F
3.03305490
3.03305490
Errore
22 19.10676444
0.86848929
Totale corretto
23 22.13981934
3.49
0.0750
Radice
R-quadro Var coeff dell'MSE Media di misura
0.136995
0.928214
0.931928
100.4001
Media Valore
Origine DF SS Anova quadratica
F Pr > F
Lab
1 3.03305490
3.03305490
3.49
0.0750
7
Se cerchiamo di verificare se la misura e’ influenzata o meno dal giorno in cui e’ stata eseguita:
ANOVA2.SAS
ODS graphics on; /*attivo la grafica*/
ods listing close;
ods rtf path="C:\Users\Alex\Desktop\SAS\ODSout"
file='anova2.rtf'; /*così l'output è su questo file*/
proc anova data=tmp1.Laboratori1;
class giorno; /*variabile categorica*/
model misura=giorno;
run;
ods rtf close;
ods listing;
7
OUTPUT:
Informazioni sui livelli di
classificazione
Classe Livelli Valori
12 1 2 3 4 5 6 7 8 9 10 11 12
giorno
Numero osservazioni lette
1529
Numero osservazioni usate
24
Somma dei
Media Valore
DF
quadrati quadratica
F Pr > F
Origine
Modello
11
5.71709639
0.51973604
Errore
12 16.42272294
1.36856025
Totale corretto
23 22.13981934
0.38
0.9402
Radice
R-quadro Var coeff dell'MSE Media di misura
0.258227
1.165193
1.169855
100.4001
Media Valore
Origine DF SS Anova quadratica
F Pr > F
giorno
11 5.71709639
0.51973604
0.38
0.9402
7
Gli esperimenti vengono ora eseguiti in soli 3 giorni diversi, sempre in due laboratori:
misura
99,57
98,33
100,13
100,29
98,85
101,19
101,19
99,96
100,33
100,17
99,81
100,72
99,91
102,68
100,36
100,61
101,57
Lab
1
1
1
1
1
1
1
1
1
1
1
1
2
2
2
2
2
giorno
1
1
1
1
2
2
2
2
3
3
3
3
1
1
1
1
2
7
100,56
100,40
99,67
100,79
99,16
101,21
102,12
2
2
2
2
2
2
2
2
2
2
3
3
3
3
Vediamo se c’e’ una dipendenza sia dal laboratorio che dal giorno e dalla loro “interazione”
ANOVA3.SAS
ODS graphics on; /*attivo la grafica*/
ods listing close;
ods rtf path="C:\Users\Alex\Desktop\SAS\ODSout"
file='anova3.rtf'; /*così l'output è su questo file*/
proc anova data=tmp1.Laboratori2;
class lab giorno; /*variabile categorica*/
model misura=lab giorno lab*giorno;
run;
ods rtf close;
ods listing;
OUTPUT:
Informazioni sui
livelli di
classificazione
Classe Livelli Valori
Lab
2 12
giorno
3 123
Numero osservazioni lette
1529
Numero osservazioni usate
Origine
DF
Modello
5
24
Somma dei
Media Valore
quadrati quadratica
F Pr > F
4.60262472
0.92052494
Errore
18 17.53719462
0.97428859
Totale corretto
23 22.13981934
0.94
0.4762
7
Radice
R-quadro Var coeff dell'MSE Media di misura
0.207889
Origine
0.983127
0.987061
100.4001
Media Valore
DF SS Anova quadratica
F Pr > F
Lab
1 3.03305490
3.03305490
3.11
0.0946
giorno
2 0.37837044
0.18918522
0.19
0.8252
Lab*giorno
2 1.19119938
0.59559969
0.61
0.5535
ESERCIZIO 8
Si considerino i dati sugli interventi del servizio 118 di Milano
riportati nel sas data set Milano_covar.sas7bdat. Esso contiene
sia il conteggio del numero di interventi giornalieri effettuati
dal 2005 al 2007, diviso per tipologie, sia i dati relativi ad
alcune variabili esterne (climatiche, sindrome influenzale,
inquinanti), sia alcune variabili dicotomiche, che indicano se il
giorno considerato era festivo, una festa scolastica, se ha
nevicato, etc.
Attraverso un modello di regressione o di analisi della varianza
si dia una risposta alle seguenti domande:
1. C’e’ una dipendenza fra il numero giornaliero di incidenti
stradali e le variabili climatiche? Quali variabili
climatiche sono piu’ significative nello spiegare la
variabilità degli incidenti stradali?
2. C’e’ una dipendenza degli incidenti stradali dai giorni
festivi, dalle festività nazionali e dalle vacanze
scolastiche?
8
Matrici e vettori in SAS: SAS/IML
La creazione e manipolazione di matrici e vettori in SAS si effettua attraverso il SAS/IML
(Interactive Matrix Language). Il SAS/IML può essere utilizzato, come altri software matematicoscientifici come Matlab, anche per effettuare simulazioni e Calcolo scientifico in genere.
Vediamo un primo esempio di semplice programma:
Esempio1_iml.sas
proc iml; /* attivo iml=interactive matrix language: manipola
matrici e vettori come Matlab */
A={ 1 2,
3 4}; /*creazione di una matrice A 2x2 */
print A; /*stampa A sulla finestra di output*/
B={1 -2 6,
22 9 -1}; /*creazione di una matrice B 2x3 */
print B;
C=A*B;
/*prodotto di matrici*/
D=inv(A)*B; /*inv(A) calcola l'inversa */
print C; /*stampo i risultati sulla finestra di output */
print D;
show names; /*mostra dimensioni e tipologia di tutte le variabili
in memoria*/
quit;
Eseguendo il programma, sulla finestra di output compaiono i
risultati:
A
1
3
2
4
B
1
22
-2
9
6
-1
C
8
45
91
16
30
4
14
D
20
-9.5
13
-7.5
-13
9.5
SYMBOL
ROWS
COLS TYPE
SIZE
----------- ------ ---- -----A
2
2 num
8
B
2
3 num
8
C
2
3 num
8
D
2
3 num
8
Number of symbols = 5 (includes those without values)
IML resta attivo fino a che non do il comando quit; SAS/IML puo’ essere usato anche in
modalita’ interattiva. Se chiudo la finestra di output i risultati compaiono sulla finestra “message
LOG”.
Vediamo come si definisce uno scalare, di varie tipologie:
Esempio2_iml.sas
proc iml;
x=12; /* intero */
y=12.34; /* reale */
z=.; /*dato perso o variabile vuota */
a='Hello'; /*variabile carattere */
b="Buon giorno"; /*variabile carattere */
print x y z a b;
Output:
X
Y
12
12.34
Z A
B
. Hello Buon giorno
Vediamo come si creano vettori e matrici: gli elementi si racchiudono fra graffe, le righe sono
separate da virgole
xx={1 2 3 4 5 6}; /*vettore riga*/
print xx;
8
Output:
XX
1
2
3
4
5
6
xx={1, 2, 3, 4, 5, 6}; /*vettore colonna*/
print xx;
Output:
XX
1
2
3
4
5
6
z={1 2, 3 4, 5 6}; /*matrice 3x2 */
print z;
Output:
Z
1
3
5
2
4
6
w=3#z; /*moltiplica lo scalare 3 per la matrice z*/
print w;
Output:
W
3
9
15
6
12
18
ww=z##2; /*eleva al quadrato elemento per elemento*/
print ww;
8
Output:
ww
1
9
25
4
16
36
Matrici con variabili carattere:
aa={abc defg}; /*matrice di var carattere. Trasforma tutto in
upper case*/
print aa;
Output:
AA
ABC
DEFG
bb={'ABC' 'defg'}; /*matrice di var carattere. Usando ' o "
mantiene il case immesso*/
cc={"abc" "DEFG"};
print bb cc;
Output:
BB
CC
ABC defg abc DEFG
store xx; /*memorizza nel dataset work.imlstore i contenuti del vettore xx */
show storage; /*mostra contenuto dello storage */
Output:
Contents of storage library = WORK.IMLSTOR
Matrices:
XX
free xx; /*cancello xx dalla memoria*/
show names;
Output:
Modules:
SYMBOL
ROWS
COLS TYPE
SIZE
------ ------ ------ ---- -----a
2
2 num
8
aa
1
2 char
4
8
b
2
bb
1
c
2
cc
1
d
2
r
2
w
3
ww
3
x
1
y
1
z
3
Number of symbols
2
2
2
2
2
2
2
2
1
1
2
=
num
8
char
4
num
8
char
4
num
8
num
8
num
8
num
8
num
8
num
8
num
8
14 (includes those without values)
load xx; /* lo ricarico */
show names;
Output:
SYMBOL
ROWS
COLS TYPE
SIZE
------ ------ ------ ---- -----a
2
2 num
8
aa
1
2 char
4
b
2
2 num
8
bb
1
2 char
4
c
2
2 num
8
cc
1
2 char
4
d
2
2 num
8
r
2
2 num
8
w
3
2 num
8
ww
3
2 num
8
x
1
1 num
8
xx
6
1 num
8
y
1
1 num
8
z
3
2 num
8
Number of symbols = 14 (includes those without values)
8
Functions as Expressions
Matrices can be created as a result of a function call. Scalar functions such as LOG
or SQRT operate on each element of a matrix, while matrix functions such as INV or
RANK operate on the entire matrix.
For example, the statement
a=sqrt(b);
assigns the square root of each element of B to the corresponding element of A.
The statement
y=inv(x);
calls the INV function to compute the inverse matrix of X and assign the results to
Y.
The statement
r=rank(x);
creates a matrix R with elements that are the ranks of the corresponding elements of
X. (rank = ordine lessicografico crescente)
Control Statements
SAS/IML software has a set of statements for controlling program execution. Control
statements direct the flow of execution of statements in IML. With them, you can
define DO-groups and modules (also known as subroutines) and route execution of
your program. Some control statements are described as follows.
Statements Action
DO, END group statements
iterative DO, END define an iteration loop
GOTO, LINK transfer control (attenzione: evitarli se possibile per non causare loop)
IF-THEN/ELSE routes execution conditionally
PAUSE instructs a module to pause during execution
QUIT ends a SAS/IML session
RESUME instructs a module to resume execution
RETURN returns from a LINK statement or a CALL module
RUN executes a module
START, FINISH define a module
STOP, ABORT stop execution of an IML program
Alcuni di questi comandi saranno utilizzati nei programmi successivi.
Per ulteriori dettagli su SAS/IML si veda la guida per l’utente iml_ug_7306.pdf (disponibile su ariel
o sul sito del SAS Institute, alla voce “documentation”).
8
Permutation Tests
Riprendiamo la domanda a cui abbiamo cercato di rispondere nell’Esercizio 4, riferita alle telefonate al
servizio 118 di Milano:
Nelle telefonate del mese di gennaio ci sono differenze significative fra i
tempi medi di “impegno operatore” fra i codici rossi e i codici verdi?
Nell’esercizio 4 avevamo applicato un t-test per il confronto di medie fra campioni indipendenti. Proviamo
ora a risolvere lo stesso problema, ossia effettuare un confronto di medie, utilizzando un permutation test.
Assumeremo che le distribuzioni dei due campioni “tempi di impegno operatore con codici rossi” e “tempi di
impegno operatore con codici verdi” differiscano al piu’ per la media (si dice che stiamo effettuando un
location test) e verifichiamo l’ipotesi
H0: media codici rossi = media codici verdi
contro tutte le alternative.
Permutation1.sas
/*DOMANDA da Esercizio 4: Nelle telefonate del mese di gennaio ci sono
differenze
significative fra i tempi medi di “impegno operatore” fra i
codici rossi e i codici verdi?*/
/* Filtro i dati in due diversi dataset, contenenti i due campioni di
interesse*/
data work.filtro1;
set tmp1.telefonate_gen2008;
where ID_CODICE='R';
run;
data work.filtro2;
set tmp1.telefonate_gen2008;
where ID_CODICE='V';
run;
proc iml;
nsamples=500; /*numero di permutazioni da effettuare*/
*Estraggo i vettori di interesse, che userò per confrontare le medie;
*primo campione: lo salvo nel vettore "rossi";
edit work.filtro1;
read all var {ImpegnoOper} into rossi;
nx=nrow(rossi);
*secondo campione: lo salvo nel vettore "verdi";
edit work.filtro2;
read all var {ImpegnoOper} into verdi;
*function che effettua la permutazione dei dati;
start RandPerm(N);
u = 1:N;
/** allocate u **/
call randgen(u, "Uniform"); /** fill with random numbers U(0,1) **/
return(rank(u));
/** return rank (order) of elements **/
finish RandPerm;
8
*function che permuta le righe della matrice o vettore M;
start PermuteRows(M);
n = nrow(M);
/** how many rows? **/
r = RandPerm(n);
/** get a random permutation **/
return ( M[r] );
/** return matrix with permuted rows **/
finish PermuteRows;
n1=nrow(rossi);
n2=nrow(verdi);
media_rossi=(1/n1)*rossi[1:n1][+];
media_verdi=(1/n2)*verdi[1:n2][+]; /*ho calcolato le medie sui dati inziali */
* print media_rossi;
* print media_verdi;
differenza_originale=media_rossi - media_verdi;
print differenza_originale;
tutti=rossi // verdi; /*sovrappongo i due vettori in un unico vettore chiamato
"tutti" */
nrighe=nrow(tutti);
differenza= J(nsamples,1,0); /*crea una matrice con nsamples righe, 1 colonna,
con tutti gli elementi =0*/
* creo gli nsamples campioni permutati;
do nperm=1 to nsamples;
M= PermuteRows(tutti); *permuto il vettore di tutti i dati insieme;
red=M[1:n1]; * i primi n1 sono rossi: li salvo in red;
green=M[n1+1:nrighe]; * gli altri sono verdi: li salvo in green;
media_red=(1/n1)*red[1:n1][+];
media_green=(1/n2)*green[1:n2][+];
differenza[nperm]=media_red - media_green; * ho calcolato medie e loro
differenza sul campione permutato;
end;
*print differenza;
*ricerca del p-value: conto numero di componenti di differenza maggiori in
modulo di differenza_originale;
check=abs(differenza)>=abs(differenza_originale); /*vettore che contiene 1 e 0 a
seconda che la disuguaglianza sia soddisfatta o no */
*print check;
pvalue1=check[1:nsamples][+];
pvalue=pvalue1/nsamples;
print pvalue;
quit;
Ovviamente a causa dell’aleatorietà delle permutazioni scelte troveremo risultati diversi in
esecuzioni diverse, anche se differiscono di poco. Riportiamo i risultati di alcune esecuzioni
differenza_originale
-2.975944
Esecuzione 1
Pvalue 0.022
Esecuzione 2
Pvalue 0.02
Esecuzione 3
Pvalue 0.006
Esecuzione 4
Pvalue 0.01
...
8
Se invece aumento a nperm=1000, ottengo minore variabilità nel calcolo dei p-value
Esecuzione 1
pvalue 0.012
Esecuzione 2
Pvalue 0.013
Esecuzione 3
Pvalue 0.011
Esecuzione 4
Pvalue 0.008
Quindi se i costi computazionali non sono troppo alti, conviene aumentare ad almeno 1000 il
numero di permutazioni da effettuare.
8
BOOTSTRAP
Utilizziamo il bootstrap per calcolare un intervallo di fiducia per il numero medio giornaliero di
telefonate al servizio 118 nel mese di gennaio.
Bootstrap.sas
Dapprima effettuo il conteggio del numero di telefonate al giorno e ne faccio un istogramma:
/*conto il numero di telefonate al giorno in gennaio e ne faccio un istogramma
*/
proc freq data=tmp1.telefonate_gen2008;
tables giorno/out=work.freq_giorno_gen;
run;
proc gchart data=work.freq_giorno_gen;;
vbar count/ midpoints=(1700 1800 1900 2000 2100 2200);
run;
Output:
FREQUENZA
12
11
10
9
8
7
6
5
4
3
2
1
0
1700
1800
1900
2000
2100
2200
Conteggio di frequenza
La distribuzione evidentemente non e’ normale e il campione non e’ molto grande, quindi non so se
posso applicare il Teorema Limite Centrale.
Applico allora il bootstrap per il calcolo di un intervallo di fiducia per la media dei dati. Effettuero’
poi un confronto con l’intervallo ottenuto usando il TLC e l’approssimazione con la normale.
/*applico il bootstrap*/
proc iml;
nsamples=500; /*numero di repliche bootstrap da generare*/
alpha=0.05; /*significativita' scelta: costruisco intervallo a livello (12alpha).
NB: scegliere alpha e nsamples in modo che il loro prodotto sia un numero
intero*/
*Estraggo il vettori di interesse: Tmorto, che userò per confrontare le medie;
* campione: lo salvo nel vettore "tempi_morti";
edit work.freq_giorno_gen;
read all var {count} into campione;
close work.freq_giorno_gen;
nx=nrow(campione);
9
*function che genera un numero casuale intero uniforme in [1, N];
start RandIntero(N);
call randgen(u, "Uniform"); /** genero un numero casuale U(0,1) **/
u=(u*N)+1;
return(floor(u));
/** ritorna il valore intero di u **/
finish RandIntero;
*function che genera una replica bootstrap di x;
start boots(x);
n = nrow(x);
/** how many rows? **/
y=j(n,1,0); /*inizializzo y */
do i=1 to n;
r = RandIntero(n);
/** seleziono una componente del vettore di
dati di partenza **/
y[i]=x[r];
end;
return (y);
/** return the bootstrap sample **/
finish boots;
/*esempio di esecuzione di boots:
x={1, 2, 3, 4, 5};
print x;
y=boots(x);
print y;
*/
* calcolo lo stimatore sul campione di partenza;
mean_ini=(1/nx)*campione[1:nx][+];
print mean_ini;
media=j(nsamples,1,0); *inizializzo il vettore media: contiene le repliche
bootstrap dello stimatore;
do i=1 to nsamples;
z=boots(campione);
media[i]=(1/nx)*z[1:nx][+]; *calcolo la media di z= i-esima replica
bootstrap dello stimatore;
end;
*calcolo l'intervallo di fiducia: salvo le repliche bootstrap nel nuovo dataset
work.repliche;
create work.repliche from media [colname='media']; *create genera un data set a
partire da una matrice creata con IML.
colname chiama la colonna media invece di col1;
append from media; *riempio le variabili usando il vettore media;
close work.repliche; *chiudo il dataset;
sort work.repliche by media; *ordino il dataset;
*rileggo dal dataset il vettore delle repliche bootstrap ordinato;
edit work.repliche;
read all var {media} into media_ord;
close work.repliche;
k=alpha*nsamples; *calcolo l'indice degli estremi dell'intervallo;
Theta_lo=media_ord[k]; *estraggo gli estremi;
Theta_up=media_ord[nsamples-k];
print Theta_lo Theta_up;
*calcolo per confronto un intervallo usando il TLC;
9
*calcolo la varianza camp. e la deviazione standard dei dati;
mm=j(nx,1,mean_ini);
vv=(campione-mm)##2;
var_ini=(1/(nx-1))*vv[1:nx][+];
se_ini=sqrt(var_ini);
*estremi dell'intervallo gaussiano;
MeanN_lo=mean_ini-1.64*(se_ini/sqrt(nx));
MeanN_up=mean_ini+1.64*(se_ini/sqrt(nx));
print MeanN_lo MeanN_up;
quit;
OUTPUT:
mean_ini: media stimata sul campione di partenza
1895.8065
Estremi intervallo bootstrap: Theta_lo
Theta_up
1843.9355 1955.4194
Estremi intervallo gaussiano: MeanN_lo
MeanN_up
1840.9031 1950.7098
Esecuzioni diverse danno esiti diversi per quanto riguarda il solo intervallo bootstrap. In generale
pero’ si ottengono intervalli bootstrap un po’ piu’ larghi di quelli gaussiani e spostati piu’ a sinistra.
Cio’ potrebbe evidenziare la presenza di un bias negli intervalli gaussiani, asintoticamente tendente
a 0, ma non corretto qui a causa della piccola dimensione del campione.
9