1 Analisi della covarianza: esempio introduttivo

2 Da completare con le formule

Abbiamo \(k\) gruppi, ciascuno di \(n_j\) osservazioni. La matrice \(\mathbf{X}\) (nell’impostazione a rango pieno e ipotizzando rette parallele) ha \(n=\sum_{j=1}^{k}n_j\) righe e \(k+1\) colonne, e il vettore della variabile di risposta è al solito una colonna di \(n\) elementi.

\[ \mathbf{y}= \left(\begin{array}{c} y_{1,1}\\ \dots \\ y_{n_1,1}\\ y_{1,2}\\ \dots \\ y_{n_2,2}\\ \dots \\ y_{i,j}\\ \dots \\ y_{1,k}\\ \dots \\ y_{n_k,k}\\ \end{array}\right) \quad \mathbf{X}= \left(\begin{array}{ccccccc} 1&0&&0&&0&z_{1,1}\\ \dots & \dots & \dots & \dots & \dots & \dots & \dots \\ 1&0& \dots & \dots & \dots &0&z_{n_1,1}\\ 0&1& \dots & \dots & \dots &0&z_{1,2}\\ \dots & \dots & \dots & \dots & \dots & \dots & \dots \\ 0&1& \dots & \dots & \dots &0&z_{n_2,2}\\ \dots & \dots & \dots & \dots & \dots & \dots & \dots \\ 0&0& \dots &1& \dots &0&z_{i,j}\\ \dots & \dots & \dots & \dots & \dots & \dots & \dots \\ 0&0& \dots &0& \dots &1&z_{1,k}\\ \dots& \dots& \dots & \dots& \dots & \dots & \dots \\ 0&0& \dots &0& \dots &1&z_{n_k,k}\\ \end{array}\right) \] e il vettore dei parametri \(\boldsymbol{ \beta}\) è: \[ \boldsymbol{ \beta}= \left\{\alpha_1,\alpha_2,\dots,\alpha_k,\beta\right\}^{\mathsf{T}} \]

Si vede che: \[ \mathbf{X}^{\mathsf{T}}\mathbf{X}= \left(\begin{array}{cccc} n_{1}&0&&0\\ & \dots &&0\\ &&n_{k}&0\\ 0&0&0&\sum_{j=1}^{k}\sum_{i=1}^{n_j}z_{ij}^{2}\\ \end{array}\right) \]

\[ \mathbf{X}^{\mathsf{T}}\mathbf{y}= \left(\begin{array}{cc} \sum_{i=1}^{n_1}y_{i1}&1\\ \dots & \dots \\ \sum_{i=1}^{n_k}y_{ik}&k\\ \sum_{j=1}^{k}\sum_{i=1}^{n_j}y_{ij}z_{ij}&k+1\\ \end{array}\right) \]

Per cui le stime di massima verosimiglianza sono: \[ \hat{\alpha}_{j} =M_{j}, j=1,2,\dots,k; \] \[ \hat{\beta}=\mathbf{b}= \frac{\sum_{j=1}^{k}\sum_{i=1}^{n_j}y_{ij}z_{ij}}{\sum_{j=1}^{k}\sum_{i=1}^{n_j}z_{ij}^{2}} \]

stima della pendenza comune:

\(\hat{\beta}=\mathbf{b}\) è una media ponderata dei \(\hat{\beta}_{j}\) dei singoli gruppi.

La devianza residua è:

\[ R(\hat{\theta}) =\sum_{j=1}^{k}\sum_{i=1}^{n_j}(y_{ij}-M_{j}-b^{*}z_{ij})^{2}=\sum_{j=1}^{k}\sum_{i=1}^{n_j}(y_{ij}-M_{j})^{2}-b^{2}\sum_{j=1}^{k}\sum_{i=1}^{n_j}z_{ij}^{2} \] \[ =\sum_{j=1}^{k}\sum_{i=1}^{n_j}(y_{ij}-M_{j})^{2}-[\sum_{j=1}^{k}\sum_{i=1}^{n_j}(y_{ij}-M_{j})z_{ij}]^{2}/\sum_{j=1}^{k}\sum_{i=1}^{n_j}z_{ij}^{2}= \]

=\(DevInt(\mathbf{y})-[Codev.Int(Z,\mathbf{y})]^{2}/DevInt(Z)\)

2.1 Inizializzazione:

vengono caricati i packages necessari per realizzare questo documento

Data Frame Summary

firms

Dimensions: 286 x 7
Duplicates: 0

Variable

Stats / Values

Freqs (% of Valid)

Graph

codice [factor]

1. azienda1 2. azienda10 3. azienda100 4. azienda101 5. azienda102 6. azienda103 7. azienda104 8. azienda105 9. azienda106 10. azienda107 [ 276 others ]

1	(	0.3%	)
1	(	0.3%	)
1	(	0.3%	)
1	(	0.3%	)
1	(	0.3%	)
1	(	0.3%	)
1	(	0.3%	)
1	(	0.3%	)
1	(	0.3%	)
1	(	0.3%	)
276	(	96.5%	)

PROV [factor]

1. AG 2. CL 3. CT 4. EN 5. ME 6. PA 7. RG 8. SR 9. TP

19	(	6.6%	)
5	(	1.8%	)
71	(	24.8%	)
4	(	1.4%	)
30	(	10.5%	)
89	(	31.1%	)
17	(	5.9%	)
26	(	9.1%	)
25	(	8.7%	)

FATTURATO [integer]

Mean (sd) : 9666 (16405.5) min < med < max: 2017 < 4035 < 149161 IQR (CV) : 5604.5 (1.7)

282 distinct values

DIPENDEN [integer]

Mean (sd) : 29.2 (59.7) min < med < max: 1 < 13 < 669 IQR (CV) : 25 (2)

74 distinct values

UTILE [integer]

Mean (sd) : -154.7 (4084.7) min < med < max: -67610 < 16 < 8218 IQR (CV) : 105.8 (-26.4)

210 distinct values

COSTOLAV [integer]

Mean (sd) : 1274.9 (3079.1) min < med < max: 5 < 455.5 < 28915 IQR (CV) : 914 (2.4)

261 distinct values

ISTAT_91 [factor]

1. agricolt 2. alberg_p 3. commerc 4. costruzi 5. cred_ass 6. energia 7. manifatt 8. serv_p_a 9. tras_tel

3	(	1.1%	)
11	(	3.9%	)
125	(	43.7%	)
41	(	14.3%	)
3	(	1.1%	)
6	(	2.1%	)
70	(	24.5%	)
7	(	2.5%	)
20	(	7.0%	)

Generated by summarytools 0.9.6 (R version 4.0.2)
2020-09-16

2.2 Il problema generale

'data.frame':   286 obs. of  7 variables:
 $ codice   : Factor w/ 286 levels "azienda1","azienda10",..: 1 112 210 221 232 243 254 265 276 2 ...
 $ PROV     : Factor w/ 9 levels "AG","CL","CT",..: 3 8 3 3 6 3 6 6 3 6 ...
 $ FATTURATO: int  149161 109433 90504 90416 79919 72705 59578 58318 50586 49338 ...
 $ DIPENDEN : int  33 472 669 348 105 207 39 119 48 60 ...
 $ UTILE    : int  0 69 -1209 68 470 59 56 -67610 160 1833 ...
 $ COSTOLAV : int  1368 24112 24166 16871 7326 6321 1609 7461 1930 845 ...
 $ ISTAT_91 : Factor w/ 9 levels "agricolt","alberg_p",..: 4 7 4 4 9 3 7 7 3 3 ...

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    5.0   198.5   455.5  1274.9  1112.5 28915.0

Plot dev stand. vs. medie all’interno dei gruppi per la sola variabile COSTOLAV ¹

Dal momento che sembra esserci un legame diretto fra medie e deviazione standard all’interno dei gruppi, vediamo, almeno dal punto di vista esplorativo, come appare la relazione fra i logaritmi delle variabili ²

Da questo momento in poi consideriamo i logaritmi in base 10 delle variabili COSTOLAV (variabile di risposta) e DIPENDEN (variabile concomitante)

2.3 analisi della varianza semplice ad una via

Vediamo intanto il modello q1, con la variabile di risposta y e la sola variabile qualitativa as.factor(z) (senza variabile concomitante)


Call:
lm(formula = y ~ as.factor(z))

Residuals:
     Min       1Q   Median       3Q      Max 
-2.17957 -0.31957  0.02198  0.32926  1.53228 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)            2.5268     0.3192   7.916 5.93e-14 ***
as.factor(z)alberg_p   0.4002     0.3601   1.111    0.267    
as.factor(z)commerc   -0.1544     0.3230  -0.478    0.633    
as.factor(z)costruzi   0.3517     0.3307   1.064    0.288    
as.factor(z)cred_ass  -0.1912     0.4514  -0.424    0.672    
as.factor(z)energia    0.3535     0.3909   0.904    0.367    
as.factor(z)manifatt   0.4002     0.3260   1.228    0.221    
as.factor(z)serv_p_a   0.2247     0.3815   0.589    0.556    
as.factor(z)tras_tel   0.4020     0.3423   1.174    0.241    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.5529 on 277 degrees of freedom
Multiple R-squared:  0.1929,    Adjusted R-squared:  0.1696 
F-statistic: 8.274 on 8 and 277 DF,  p-value: 4.697e-10

2.4 analisi della covarianza semplificata ad una via

Vediamo intanto il modello m0, con la variabile di risposta y in funzione della variabile qualitativa as.factor(z) e della variabile quantitativa x (variabile concomitante)

Abbiamo sostanzialmente adattato un modello con * k rette parallele *

modello m con variabile concomitante x senza fattore di classificazione * retta unica per tutti i gruppi *

2.5 Modello moltiplicativo

Vediamo intanto il modello m1, con la variabile di risposta y in funzione del prodotto della variabile qualitativa as.factor(z) e della variabile quantitativa x * (k rette distinte nei k gruppi) *

inserire formule *

2.6 confronto fra il modello m0 e il modello m1

Praticamente stiamo effettuando un test sull’interazione fra \(x\) e \(z\) ossia stiamo saggiando l’ipotesi di parallelismo fra le k rette

Analysis of Variance Table

Model 1: y ~ x + as.factor(z)
Model 2: y ~ x * as.factor(z)
  Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
1    276 29.260                              
2    268 27.663  8    1.5967 1.9336 0.05526 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

2.7 Confronto fra il modello m e il modello m0

stiamo saggiando l’ipotesi che non vi sia un effetto del gruppo z ossia stiamo saggiando l’ipotesi che le k rette coincidano con un unica retta generale

Analysis of Variance Table

Model 1: y ~ x
Model 2: y ~ x + as.factor(z)
  Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
1    284 31.297                              
2    276 29.260  8    2.0366 2.4013 0.01613 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

in realtà andrebbe fatta sui residui delle regressioni interne ai gruppi.↩︎
Converrà in questo caso utilizzare logaritmi in base 10 per rendere più leggibili le trasformazioni logaritmiche.↩︎

Analisi della covarianza

Marcello Chiodi

2020