17/11/2013
Analisi della regressione
L’analisi della regressione è una tecnica statistica per modellare e
investigare le relazioni tra due (o più) variabili.
Nella tavola è riportata la % di purezza di ossigeno,
rilasciata in un processo di distillazione chimica, e il
livello di idrocarbonio, presente nel condensatore principale di unità di distillazione.
Dati
salvati
in
un file
Osservazioni Liv.Idrocarbonio Purezza
1
0,99
90,01
2
1,02
89,05
3
1,15
91,43
4
1,29
93,74
5
1,46
96,73
6
1,36
94,45
7
0,87
87,59
8
1,23
91,77
9
1,55
99,42
10
1,4
93,65
11
1,19
93,54
12
1,15
92,52
13
0,98
90,56
14
1,01
89,54
15
1,11
89,85
16
1,2
90,39
17
1,26
93,25
18
1,32
93,41
19
1,43
94,98
20
0,95
87,33
2
1
17/11/2013
Con questa concatenazione di functions è possibile far leggere i dati in
un file memorizzati sotto forma di tabella e assegnarli ad una variabile.
> dati<-matrix(scan("C:/Programmi/R/R-3.0.2/regressionedati.R",
+ n=20*2),20,2,byrow=TRUE)
Read 40 items
Un primo modo di studiare che tipo
di relazione sussiste tra le popolazioni che hanno generato i dati è
il diagramma di dispersione (scatter
plot).
> datix<-dati[,1]
> datiy<-dati[,2]
> plot(datix,datiy,xlab='purezza
+ ossigeno',ylab='% di idrocarbonio')
3
CONSIDERAZIONI SULLO SCATTER DIAGRAM
Le correlazioni, possibilmente positiva o possibilmente negativa, sono
In forse quando i punti rappresentativi delle coppie di dati, pur disponendosi attorno ad una delle due diagonali del diagramma, presentano
una dispersione piuttosto accentuata tale da far presumere l’esistenza
di altre cause che intervengono a determinare l'effetto studiato.
2
17/11/2013
C’è un modo per misurare il grado di correlazione
tra due variabili aleatorie?
D efin iz io n e
Se µ X = E [X
]
e µ Y = E [Y ] , si d efin isce co va ria n za d i
X e Y , la q u a n tità
C o v (X ,Y ) = E [( X − µ X )(Y − µ Y )].
cov( X , Y ) = E[ XY ] − µ X µY
Teorema : Se X e Y sono indipendenti, Cov( X , Y ) = 0.
Il viceversa non vale.
Controesempio
X
p ( x)
−1
1
3
0
1
3
1
2
1 ⇒Y = X
⇒ Cov ( X , Y ) = 0
3
5
Teorema : Var ( X ± Y ) = Var ( X ) + Var (Y ) ± 2 cov( X , Y )
Definizione
La correlazione tra le variabili aleatorie X e Y è la quantità:
ρ=
σ
cov( X , Y )
= XY
Var ( X )Var (Y ) σ X σ Y
Se la covarianza tra due variabili aleatorie è positiva, negativa o nulla,
anche la correlazione sarà positiva, negativa o nulla.
Teorema : La correlazione tra le variabili aleatorie X e Y gode
della seguente proprietà: -1 ≤ ρ ≤ 1
Teorema : Se ρ = ±1 ⇒ P(Y = aX ± b) = 1
6
3
17/11/2013
Il coefficiente di correlazione non è una misura generale della relazione
tra due variabili, ma esprime solo il grado di linearità della correlazione
in un grafico a dispersione.
Gli outliers possono modificare significativamente il valore
del coefficiente di correlazione.
4
17/11/2013
In R la function per
calcolare il coefficiente
di correlazione è
cor()
Bisogna specificare in
input il tipo di coefficiente
da calcolare: in questo
caso quello di Pearson.
> cor(datix, datiy, method="pearson")
[1] 0.9367154
5
17/11/2013
Se si immagina che la relazione tra X e Y sia lineare, allora bisogna
ricercare i coefficienti della trasformazione
Y = mX + b + ε
• m, b coefficienti di regressione
• X (livello di idrocarbonio) var. aleat. indipendente
• Y (purezza dell'ossigeno) var. aleat. dipendente
• ε (errore casuale), E[ε ] = 0, Var[ε ] = σ 2 .
Come si calcolano i coefficienti?
IL METODO DEI MINIMI QUADRATI
Minimizzare la distanza tra i punti delle osservazioni e la retta stessa.
11
residui
Si cerca il minimo della funzione 

n
2
 rispetto m e b
L ( a, b ) = ∑ ( yi − mxi − b )

i =1

6
17/11/2013
La funzione di R per stimare un modello di regressione lineare è lm
(linear models). L’output di tale funzione è molto complesso: in generale è opportuno salvarlo in una variabile di classe per poi estrarne le
quantità necessarie all’analisi dei dati
>result<-lm( datiy ~ datix)
>summary(result)
Residuals:
Min
1Q
Median
3Q Max
-1.83029 -0.73334 0.04497 0.69969 1.96809
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 74.283
1.593
46.62
< 2e-16 ***
datix
14.947
1.317
11.35 1. 23e-09 ***
--Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
13
Pertanto il modello lineare stimato è
y=
y = 74.2 + 14.9 x
Possiamo sovrapporre la retta al grafico
con il comando:
> abline(lm(datiy ~ datix))
14
7
17/11/2013
>summary(result)
Residuals:
Min
1Q
Median
3Q Max
-1.83029 -0.73334 0.04497 0.69969 1.96809
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 74.283
1.593
46.62
< 2e-16 ***
datix
14.947
1.317
11.35 1. 23e-09 ***
--Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
L’errore standard misura il grado di precisione della stima dei parametri.
15
>summary(result)
Residuals:
Min
1Q
Median
3Q Max
-1.83029 -0.73334 0.04497 0.69969 1.96809
Coefficients:
Estimate Std. Error t value Pr(>|t|)
< 2e-16 ***
(Intercept) 74.283
1.593
46.62
datix
14.947
1.317
11.35 1.23e-09 ***
--Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
E’ l’output di un test per verificare se i valori dei coefficienti sono nulli.
La legenda associa alle stelle un valore numerico: quando questi valori numerici sono inferiori a 0.05, allora si può ritenere non vera l’ipotesi che i coefficienti sono nulli.
16
8
17/11/2013
>summary(result)
…Multiple R-squared: 0.8774, Adjusted R-squared: 0.8706
R^2 è noto nella letteratura come coefficiente di determinazione.
> rho<-cor(datix, datiy, method="pearson")
> rho^2
[1] 0.8774357
Spiega quanta parte della variabilità della Y è imputabile alla variabilità
della X, in questo caso l’87%.
Il valore numerico 0.8706 è ottenuto usando una formula corretta e
più robusta.
I valori vengono riportati entrambi perché grosse discrepanze segnalano
la presenza di valori anomali.
17
> summary(result)
F-statistic: 128.9 on 1 and 18 DF, p-value: 1.227e-09
Questo è un test per stabilire se è valido il modello con predittore
Y = mX + b + ε oppure quello senza predittore Y = b + ε
L’esito del test si legge analogamente a prima.
Si osservi che questo test equivale al test per riconoscere se m = 0.
18
9
17/11/2013
ANALISI DEI RESIDUI
> summary(result)
Residuals:
Min
1Q
Median
3Q
Max
-1.83029 -0.73334 0.04497 0.69969 1.96809
Residual standard error: 1.087 on 18 degrees of freedom
Si definiscono residui le distanze tra i valori della variabile dipendente
ottenuti mediante la retta di regressione e quelli osservati.
ei = yi − yˆi
Anche in questo caso l’errore standard misura la variabilità dell’errore
stocastico ε inserito nel modello
Perchè il modello sia valido è necessario provare che ε ∼ N(0,σ 2 )19
Per ora usiamo due strumenti grafici. Il grafico dei residui e il qqnorm.
> str(result)
List of 12
$ coefficients : Named num [1:2] 74.3 14.9
..- attr(*, "names")= chr [1:2] "(Intercept)" "datix"
$ residuals : Named num [1:20] 0.9287 -0.4797 -0.0429 0.1744 0.6234 ...
..- attr(*, "names")= chr [1:20] "1" "2" "3" "4" ...
$ effects : Named num [1:20] -412.154 12.334 -0.207 -0.084 0.25 ...
..- attr(*, "names")= chr [1:20] "(Intercept)" "datix" "" "" ...
$ rank
: int 2
$ fitted.values: Named num [1:20] 89.1 89.5 91.5 93.6 96.1 ...
..- attr(*, "names")= chr [1:20] "1" "2" "3" "4" ...
$ assign
: int [1:2] 0 1
$ qr
:List of 5
..$ qr : num [1:20, 1:2] -4.472 0.224 0.224 0.224 0.224 ...
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$ : chr [1:20] "1" "2" "3" "4" ...
.. .. ..$ : chr [1:2] "(Intercept)" "datix"
.. ..- attr(*, "assign")= int [1:2] 0 1
..$ qraux: num [1:2] 1.22 1.17…
10
17/11/2013
> qqnorm(result$residuals)
>
I punti si distribuiscono lungo
una retta.
> plot(1:20,result$residuals,type='p',
+ main='Grafico dei residui',
+ xlab='residui')
> abline(h=0, lwd=2)
La dispersione non presenta
particolari patterns.
ESEMPI DI GRAFICI DI RESIDUI
Omoschedasticità=
varianza
costante per
l’errore
Eteroschedasticità=
varianza
diversa per
l’errore
11
17/11/2013
La relazione tra due campioni casuali non è sempre lineare. A volte è
possibile trasformare opportunamente le scale di misurazione perché
il modello di regressione lineare sia sufficiente.
Il seguente dataset rappresenta il peso corporeo espresso in kg e il
peso del cervello espresso in gr. per 62 diverse specie di mammiferi.
> cervello<-matrix(scan("C:/Programmi/R/R-3.0.2/cervello.txt",
+ skip=1,n=62*2),62,2,byrow=TRUE)
Read 124 items
Effettuando un diagramma
di dispersione dei dati, si
osserva un andamento molto
lontano da quello lineare.
Provando a rappresentare i dati in scala logaritmica, si ottiene…
> plot(log(cervello[,1]),log(cervello[,2]),type='p',main='Scatter-plot
+ Scala Log',xlab='Peso corporeo',ylab='Peso cervello')
> abline(lm(log(cervello[,2]) ~ log(cervello[,1])) )
>
12
17/11/2013
Torniamo all’esercizio sul fitting dei dati:
Esempio: si misura il tempo di vita di 20 lampadine, che è stato riportato
di seguito. Stimare la funzione guasto, la funzione di affidabilità, la densità di guasto, il tempo medio di vita. Determinare quale modello teorico
potrebbe descrivere il tempo di vita.
10,45
6,10
2,30 19,71 49,81 46,89 69,08
1,01 0,71 3,89 5,38 0,37
0,32 11,36 43,20 3,24
7,29 9,12 17,52 9,60
Usando i probability plotting papers, avevamo stabilito che il modello
di Weibull si adattava bene ai dati. Usando la function fitdistr()
> stima<-fitdistr(dati,"weibull")
> str(stima)
Usando la regressione
sui dati del probability
plotting paper?
List of 5
$ estimate: Named num [1:2] 0.764 13.497
..- attr(*, "names")= chr [1:2] "shape" "scale"
$ sd : Named num [1:2] 0.133 4.176
..- attr(*, "names")= chr [1:2] "shape" "scale"
…
E’ necessario costruire il modello di regressione tra
> source('datiprimoesercizio.r')
> x<-log(sort(dati))
> index<-seq(1,20,1)
> cdf<-(index-0.3)/20.4
> z<-log(log(1/(1-cdf)))
Con fitdistr()
$ estimate: 0.764 13.497
"shape" "scale”
Per ottenere i coefficienti è sufficiente che:
> result<-lm(z ~ x)
> result
Coefficients:
(Intercept)
x
-1.8988
0.7276
Ovviamente questi non sono i parametri
di scala e di forma, perché nella retta di
regressione associata al probability plotting paper è
z = β x − β log α ⇒
m = β e b = − β log α
ossia β = 0.7276 e α =13.59464
13
17/11/2013
Un caso in cui l’uso Del Modello di regressione
è necessario
L’esercizio che segue si riferisce a un campione casuale di 90 tempi di
guasto.
I tempi sono ripartiti in modalità e viene fornita la tabella seguente.
Rank
Tempo guasto
# guasti
Guasti cum.
Median rank %
1
600
23
23
25.1
2
1200
9
32
35.1
3
1800
2
34
37.3
4
2400
3
37
40.6
5
3000
11
48
52.8
6
3600
6
54
59.4
7
4200
3
57
62.7
8
4800
6
63
69.4
9
5400
4
67
73.8
10
6600
3
70
77.1
11
6900
4
74
81.5
12
7200
8
82
90.4
13
7800
3
85
93.7
14
8400
0
85
93.7
15
9000
1
86
94.8
16
9600
2
88
97.0
17
10200
2
90
99.2
14
17/11/2013
Un altro modo per
leggere le tabelle di
dati…
In R, aprire un file dove
riportare i dati e salvarlo con estensione .csv.
Questi file, una volta
letti con il comando
read.csv() producono
delle matrici con
la stessa struttura di
quella usata nel file
> guasti<-read.csv("C:/Programmi/R/guasti.csv",header=F,sep=",")
Per costruire il probability plotting paper relativo al modello di Weibull
> guasti
V1 V2
1 600 23
2 1200 9
3 1800 2
4 2400 3
5 3000 11
6 3600 6
7 4200 3
8 4800 6
9 5400 4
10 6600 3
11 6900 4
12 7200 8
13 7800 3
14 8400 0
15 9000 1
16 9600 2
17 10200 2
> x=log(guasti[,1])
> cdf=(cumsum(guasti[,2])-0.3)/(90+0.4)
> y=log(log(1/(1-cdf)))
Grafico
Non
Corretto
15
17/11/2013
Anche ricostruendo il campione casuale a partire dalla tabella…
> dati<-c(rep(guasti[1,1],guasti[1,2]))
> dati
[1] 600 600 600 600 600 600 600 600 600 600 600 600 600 600 600 600 600 600 600
[20] 600 600 600 600
> for (i in 2:17) {
+ dati<-c(dati, rep(guasti[i,1],guasti[i,2]))
+}
> dati
[1] 600 600 600 600 600 600 600 600 600 600 600 600
[13] 600 600 600 600 600 600 600 600 600 600 600 1200
[25] 1200 1200 1200 1200 1200 1200 1200 1200 1800 1800 2400 2400
[37] 2400 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000
[49] 3600 3600 3600 3600 3600 3600 4200 4200 4200 4800 4800 4800
[61] 4800 4800 4800 5400 5400 5400 5400 6600 6600 6600 6900 6900
[73] 6900 6900 7200 7200 7200 7200 7200 7200 7200 7200 7800 7800
[85] 7800 9000 9600 9600 10200 10200
>
Il grafico che ne risulta non consente di ritenere adeguato il modello
di Weibull…
…perché vi sono molti dati che si ripetono
16
17/11/2013
Nel grafico originale invece riconosciamo tre andamenti:
Rank
Tempo guasto
1
600
# guasti
23
Guasti cum.
23
Median rank %
25.1
2
1200
9
32
35.1
3
1800
2
34
37.3
4
2400
3
37
40.6
5
3000
11
48
52.8
6
3600
6
54
59.4
7
4200
3
57
62.7
8
4800
6
63
69.4
9
5400
4
67
73.8
10
6600
3
70
77.1
11
6900
4
74
81.5
12
7200
8
82
90.4
13
7800
3
85
93.7
14
8400
0
85
93.7
15
9000
1
86
94.8
16
9600
2
88
97.0
17
10200
2
90
99.2
17
17/11/2013
Num.prova
Tempo guasto
# guasti
Guasti cum.
Median rank %
1
600
23
23
60.7
2
1200
9
32
84.8
3
1800
2
34
90.1
4
2400
3
37
98.1
5
3000
11
48
28.6
6
3600
6
54
44.7
7
4200
3
57
52.7
8
4800
6
63
68.7
9
5400
4
67
79.4
10
6600
3
70
87.4
11
6900
4
74
98.1
12
7200
8
82
47.0
13
7800
3
85
65.3
14
8400
0
85
65.3
15
9000
1
86
71.3
16
9600
2
88
83.6
17
10200
2
90
95.8
Num.prova
Tempo guasto
# guasti
Guasti cum.
Median rank %
1
600
23
23
60.7
2
1200
9
32
84.8
3
1800
2
34
90.1
4
2400
3
37
98.1
5
3000
11
11
28.6
6
3600
6
17
44.7
7
4200
3
20
52.7
8
4800
6
26
68.7
9
5400
4
30
79.4
10
6600
3
33
87.4
11
6900
4
37
98.1
12
7200
8
8
47.0
13
7800
3
11
65.3
14
8400
0
11
65.3
15
9000
1
12
71.3
16
9600
2
14
83.6
17
10200
2
16
95.8
18
17/11/2013
> dati1<-dati[1:37]
> dati2<-dati[38:74]
> dati3<-dati[74:90]
> library(MASS)
> stima1<-fitdistr(dati1,'weibull')
Error in fitdistr(dati1, "weibull") : optimization failed
In addition: Warning messages:
1: In densfun(x, parm[1], parm[2], ...) : NaNs produced
2: In densfun(x, parm[1], parm[2], ...) : NaNs produced
Il motivo dell’errore è che nella formula che stima uno dei due parametri della Weibull, compaiono le distanze tra i dati al denominatore.
Si può ovviare all’inconveniente usando la regressione lineare.
Una procedura più corretta consiste nell’uso dei cosiddetti modelli lineari generalizzati. Tuttavia ci “accontentiamo” di costruire un modello di regressione per ciascuno dei 3 sottogruppi trovati.
Per il primo gruppo:
> guasti1<-guasti[1:4,1]
> x<-log(sort(guasti1))
> fas<-guasti[1:4,2]
> cdf<-(cumsum(fas)-0.3)/37.4
> z<-log(log(1/(1-cdf)))
> result<-lm(z ~ x)
> result
Coefficients:
(Intercept)
x
-6.352
0.980
Con la trasformazione vista prima:
β =0.980 e α = 653.0363
19
17/11/2013
Per il secondo gruppo:
> guasti2<-guasti[5:11,1]
> x<-log(sort(guasti2))
> fas<-guasti[5:11,2]
> cdf<-(cumsum(fas)-0.3)/37.4
> z<-log(log(1/(1-cdf)))
> result<-lm(z ~ x)
> result
Coefficients:
(Intercept)
x
-22.249
2.643
Con la trasformazione vista prima:
β =2.643 e α = 4528.226
>
Per il terzo gruppo:
> guasti3<-guasti[12:17,1]
> x<-log(sort(guasti3))
> fas<-guasti[12:17,2]
> cdf<-(cumsum(fas)-0.3)/37.4
> z<-log(log(1/(1-cdf)))
> result<-lm(z ~ x)
> result
Coefficients:
(Intercept)
x
-20.791
2.183
Con la trasformazione vista prima:
β =2.183 e α = 13436.44
>
20
17/11/2013
L’esempio mostra come per questo componente le informazioni circa il
tasso di guasto siano ripartite lungo l’arco temporale. La affidabilità totale è una somma di affidabilità, ciascuna si riferisce a un periodo della
vita, secondo dei pesi che restituiscono il peso di quel tasso di guasto
lungo l’arco temporale
R (t ) = a1 R1 (t ) + a2 R2 (t ) + a3 R3 (t )
con a1 + a2 + a3 = 1
f (t ) = a1 f1 (t ) + a2 f 2 (t ) + a3 f3 (t )
Iperesponenziale
MISTURE
F (t ) = a1 F1 (t ) + a2 F2 (t ) + a3 F3 (t )
Nell’esempio è:
β1
β2
 t 
 t 
 t 
37
37
16
R(t ) = exp  −  + exp  −  + exp  − 
90
90
90
 α1 
 α2 
 α3 
β3
21