Capitolo 7 La verifica delle assunzioni di base

Nel momento in cui eseguiamo test d’ipotesi nell’ambito di un modello lineare, assumiamo implicitamente che i dati rispondano ai seguenti requisiti:

  1. il modello è corretto;
  2. la risposta osservata è una funzione del modello più o meno l’errore sperimentale;
  3. l’ errore sperimentale è indipendente dal modello;
  4. gli errori sono normalmente distribuiti, con media zero e varianze omogenee;
  5. gli errori rilevati in esperimenti ripetuti sono tra loro indipendenti.
  6. non ci sono osservazioni aberranti;

Il rispetto di alcune di queste assunzioni di base è garantito dal disegno sperimentale, con il quale, attraverso la randomizzazione, si cerca di evitare ogni possibile confusione tra effetti deterministici ed errore sperimentale e si fa in modo che le osservazioni siano indipendenti tra di loro. Eventuali vincoli alla randomizzazione introdotti in fase di disegno (es. blocking) vengono gestiti con l’adozione di un modello ANOVA che ne riconosca l’esistenza (vderemo, ad esempio, l’ANOVA a blocchi randomizzati).

Ci sono tuttavia alcune assunzioni di base che possono essere verificate solo a posteriori, cioè dopo aver effettuato la stima dei parametri. Queste assunzioni riguardano la struttura dei residui, che debbono essere normalmente distribuiti, omoscedastici e non debbono contenere valori aberranti. Ovviamente, per ottenere i residui, devo prima procedere al fitting del modello, individuare i valori attesi e quantificare il loro scostamento rispetto alle osservazioni.

Il rispetto delle assunzioni di base è importante, perché ogni deviazione può inficiare la validità dei test d’ipotesi e delle inferenze in genere, modificando il livello di significatività.

7.1 Procedure diagnostiche

Per diagnosticare eventuali patologie dei dati, possiamo utilizzare diverse procedure, così classificate:

  1. metodi grafici
  2. metodi formali, tramite test d’ipotesi

7.2 Analisi grafica dei residui

La gran parte dei pre-requisiti fondamentali di un dataset riguardano la struttura dei residui e, di conseguenza, l’ispezione grafica di questi ultimi, eventualmente accompagnata da semplici strumenti algebrici, possono permetterci di evidenziare la gran parte delle ’patologie’ di cui soffrono i dati sperimentali. Si può affermare che l’ispezione dei residui è uno strumento diagnostico fondamentale, il cui impiego dovrebbe rientrare tra le metodiche di routine per ogni elaborazione statistica dei dati.

Ricordo che i residui sono gli scostamenti tra i valori osservati e quelli attesi sulla base del modello in studio; hanno sempre media 0, mentre la deviazione standard cambia di volta in volta.

7.2.1 Grafico dei residui contro i valori attesi

Il metodo grafico più utilizzato per l’analisi dei residui consiste nel plottare i residui verso i valori attesi. Se non vi sono problemi, i punti nel grafico dovrebbero essere distribuiti in modo assolutamente casuale, come in figura 7.1.

Grafico dei residui verso gli attesi: non è visibile nessuna deviazione rispetto agli assunti di base dei modelli lineari

Figure 7.1: Grafico dei residui verso gli attesi: non è visibile nessuna deviazione rispetto agli assunti di base dei modelli lineari

Le osservazioni aberranti (outliers) sono chiaramente indicate nel grafico dei residui come punti isolati rispetto agli altri (Figura 7.2).

Grafico dei residui verso gli attesi: presenza di un dato aberrante

Figure 7.2: Grafico dei residui verso gli attesi: presenza di un dato aberrante

L’eterogeneità delle varianze è invece indicata da residui che si allargano o si stringono procedendo verso i margini del grafico (Figura 7.3), facendo emergere una sorta di proporzionalità tra media e varianza.

Grafico dei residui verso gli attesi: esempio di eteroscedasticità

Figure 7.3: Grafico dei residui verso gli attesi: esempio di eteroscedasticità

A volte la relazione causa effetto non è lineare o, comunque, il modello devia sistematicamente dai dati osservati. Di conseguenza, i residui mostrano un evidente pattern legato al segno. Ad esempio, nella figura 7.4 i residui sono tendenzialmente negativi per bassi valori attesi e positivi per alti valori.

Grafico dei residui verso gli attesi: esempio di lack of fit (nonlinearità)

Figure 7.4: Grafico dei residui verso gli attesi: esempio di lack of fit (nonlinearità)

7.2.2 QQ-plot

Il grafico dei residui verso i valori attesi non è in grado di evidenziare problemi di non-normalità. A questo fine, risulta molto utile un QQ-plot (quantile-quantile plot), nel quale i residui standardizzati vengono plottati contro i rispettivi percentili di una distribuzione normale standardizzata. Se la distribuzione dei dati è effettivamente normale, le due grandezze (i percentili nel campione e i percentili di una distribuzione normale standardizzata) dovrebbero essere largamente coincidenti; ad esempio, se la distribuzione è normale i residui negativi dovrebbero essere più o meno ugualmente numerosi di quelli positivi e non ci dovrebbero essere residui troppo alti in valore asoluto, soprattutto se il dataset è piccolo (infatti i valori estremi sono rari e dovrebbero comparire di rado). Se questo è vero, il QQ-plot dovrebbe essere costituito da punti che giacciono lungo la bisettrice del primo e del terzo quadrante (7.5)

QQ-plot per un dataset normalmente distribuito

Figure 7.5: QQ-plot per un dataset normalmente distribuito

Le deviazioni più diffuse dalla normalità sono relative alla simmetria (skewness) e alla curtosi della popolazione. In particolare, può capitare che i residui abbiano asimmetria positiva (right-skewed: la media è maggiore della mediana), così che quelli negativi sono più numerosi di quelli positivi, ma questi ultimi sono mediamente di più elevato valore assoluto. Ad esempio, una distribuzione log-normale centrata è right-skewed ed, estraendo da questa una serie di valori, il QQ-plot si presenta come in figura 7.6.

QQ-plot per un dataset con asimmetria positiva (right-skewed)

Figure 7.6: QQ-plot per un dataset con asimmetria positiva (right-skewed)

Al contrario, in una distribuzione left-skewed (asimmetria negativa), la media è minore della mediana e, di conseguenza, i residui positivi sono più numerosi, ma di valore assoluto più basso che non nella distribuzione normale corrispondente. Ad esempio, una distribuzione beta, in certe condizioni ( traslazione, alto \(a\) e basso \(b\)), è left-skewed e i valori campionati presentano un QQ-plot con un andamento tipico (7.7):

QQ-plot per un dataset con asimmetria negativa (left-skewed)

Figure 7.7: QQ-plot per un dataset con asimmetria negativa (left-skewed)

Per quanto riguarda la curtosi, è necessario osservare le code della distribuzione: se queste sono più alte di una distribuzione normale parliamo di distribuzione platicurtica, mentre se sono più basse parliamo di distribuzione leptocurtica. Ad esempio, una distribuzione di t di Student con pochi gradi di libertà è platicurtica ed il QQ-plot mostra l’andamento indicato in figura 7.8.

QQ-plot per un dataset con distribuzione platicurtica

Figure 7.8: QQ-plot per un dataset con distribuzione platicurtica

Al contrario, una distribuzione uniforme in un intervallo ristretto è tipicamente leptocurtica e presenta un QQ-plot come quello riportato in Figura 7.9.

QQ-plot per un dataset con distribuzione leptocurtica

Figure 7.9: QQ-plot per un dataset con distribuzione leptocurtica

7.3 Strumenti diagnostici formali

La valutazioni precedentemente esposte sono di tipo grafico e sono considerate sufficientemente robuste per la maggior parte delle situazioni. Tuttavia, esistono anche test statistici che consentono di saggiare l’ipotesi nulla di ’assenza di deviazioni’. Con questi test, basta guardare il ’p-value’: se questo è inferiore a 0.05 l’ipotesi nulla deve essere rifiutata e può essere necessario intraprendere azioni correttive.

Per l’omogeneità delle varianze veniva utilizzato il test di Bartlett, che, tuttavia, è oggi caduto in disuso, data la sua sensibilità alla non-normalità dei residui, quasi sempre presente, insieme all’eteroscedasticità. Oggi si preferisce utilizzare il test di Levene, che consiste nel sottoporre ad ANOVA i residui in valore assoluto, al posto dei dati osservati. Per ogni trattamento, i residui hanno media zero, ma se vengono presi in valore assoluto, hanno medie più alte quando la loro varianza è alta. Per esempio, possiamo prendere due campioni centrati, con media zero a varianza rispettivamente pari a 1 e 4:

A <- c(-1, 0, 1); B <- c(-4, 0, 4)
mean(A); mean(B)
## [1] 0
## [1] 0
var(A); var(B)
## [1] 1
## [1] 16

Se prendiamo i valori assoluti, la media del primo campione è 2/3, mentra la media del secondo campione è 8/3. Se questa differenza è significativa, essa produce il rifiuto dell’ipotesi nulla nel test F di ANOVA e conferma l’eterogeneità delle varianze. Il test di Levene, in R, si può eseguire con la funzione ‘leveneTest()’ nel package ‘car’.

res <- c(A, B)
tratt <- rep(c("A", "B"), each = 3)
model <- lm(res ~ factor(tratt))
anova(model)
## Analysis of Variance Table
## 
## Response: res
##               Df Sum Sq Mean Sq F value Pr(>F)
## factor(tratt)  1      0     0.0       0      1
## Residuals      4     34     8.5
car::leveneTest(res ~ factor(tratt), center=mean)
## Levene's Test for Homogeneity of Variance (center = mean)
##       Df F value Pr(>F)
## group  1  2.1176 0.2193
##        4

Il test di Levene può anche essere effettuato considerando gli scarti rispetto alla mediana (invece che alla media), in modo da ottenere un test più robusto nei confronti degli outliers.

Per quanto riguarda le deviazioni dalla normalità, può essere utilizzato il test di Shapiro-Wilks. Per esempio, nel caso di un dataset ottenuto da una distribuzione uniforme (quindi non-normale), il test di Shapiro porta ai risultati sotto indicati.

shapiro.test(runif(100, min = -2, max = 2))
## 
##  Shapiro-Wilk normality test
## 
## data:  runif(100, min = -2, max = 2)
## W = 0.94547, p-value = 0.0004227

7.4 Risultati contraddittori

La valutazione degli assunti di base è un passo fondamentale nell’analisi dei dati sperimentali è non deve mai essere dimenticata. Tuttavia, l’esperienza insegna che, nella pratica, è molto facile incontrare situazioni dubbie, nelle quali i risultati ottenuti con le diverse tecniche diagnostiche appaiono contraddittori e difficili da interpretare. Come comportarsi in questi casi? A mio parere bisogna sempre ricordare che la ’verità vera’ ci sfugge e, di conseguenza, tutte le valutazioni debbono essere sempre condotte con il massimo ’buon senso’.

Un aspetto importante da considerare è la tipologia del dato: conteggi e proporzioni difficilmente sono normalmente distribuiti ed omoscedastici e, con questi dati, la prudenza non è mai troppa, quando si tratta di impiegare modelli lineari. Allo stesso modo, è necessaria una grande prudenza quando si analizzano variabili quantitative dove la differenza tra le diverse tesi è molto grande, orientativamente più di un ordine di grandezza. Con questi dati, l’assunzione di omogeneità delle varianze è quasi sempre violata.

7.5 ‘Terapia’

Se le procedure diagnostiche hanno evidenziato deviazioni rispetto agli assunti di base, è necessario valutare se e come intraprendere azioni correttive. Ovviamente, la ‘terapia’ cambia al cambiare della ‘patologia’.

7.5.1 Correzione/Rimozione degli outliers

In presenza di outliers, la ‘terapia’ più opportuna è, banalmente, quella di rimuoverli, ottenendo così un dataset ‘sbilanciato’ (diverso numero di repliche per trattamento). Oggigiorno, trattare un dataset sbilanciato non costituisce un problema, ovviamente se si utilizzano le metodiche di analisi opportune. Qualche anno fa, al contrario, si cercava di evitare lo sbilanciamento a tutti i costi, utilizzando tecniche di imputing per l’immissione di valori ’ragionevoli’ in sostituzione di quelli mancanti/aberranti. Con le dovute eccezioni, le tecniche di imputing sembrano oggi abbastanza obsolete.

Vorrei porre l’attenzione sul fatto che i dati aberranti sono dati molto ‘influenziali’, nel senso che possono influenzare in modo molto marcato il risultato dell’analisi. Pertanto, se è sbagliato correggerli arbitrariamente, senza aver prima accertato che siano effettivamente frutto di errore, è altrettanto sbagliato lasciarli nel dataset. Ovviamente, la correzione non può che riguardare una larga minoranza dei dati sperimentali raccolti (uno o due dati), altrimenti si dovrà necessariamente pensare di ripetere l’esperimento.

7.5.2 Correzione del modello

Abbiamo visto che il modello impiegato potrebbe non essere adatto a descrivere il dataset (mancanza di adattamento). Gli effetti potrebbero non essere addittivi, ma moltiplicativi, oppure potrebbero essere non-lineari. Potrebbero essere presenti asintoti che il nostro modello non è in grado di descrivere, oppure la crescita/decrescita osservata potrebbero essere più lente/veloci di quanto la nostra funzione sia in grado di modellizzare. Per tutti questi casi, ovviamente, la strategia più consigliabile è quella di utilizzare un diverso modello, capace di descrivere meglio le osservazioni sperimentali.

7.5.3 Trasformazione della variabile indipendente

A volte la variabile dipendente non è qualitativa, bensì quantitativa. Vedremo meglio questo aspetto nei prossimi capitoli, quando parleremo di regressione. Tuttavia, anticipiamo che, quando la variabile indipendente è quantitativa, essa può essere sottoposta ad opportuna trasformazione. Ad esempio, se i dati mostrano un andamento esponenziale, la trasformazione della variabile indipendente in logaritmo può portare al sensibile miglioramento del fitting con un’equazione lineare.

7.5.4 Impiego di metodiche statistiche avanzate

In presenza di deviazioni sistematiche rispetto agli assunti di base, la cosa più logica sarebbe quella di chiedersi perché il dataset è non-normale e/o eteroscedastico ed incorporare queste informazioni nel modello. Ad esempio, un conteggio potrebbe seguire la distribuzione di Poisson, oppure una serie di proporzioni potrebbero seguire la distribuzione binomiale. In questi casi sarebbe opportuno utilizzare modelli lineari generalizzati (GLiM), basati non sulla distribuzione normale, ma su altre distribuzioni, più adatte per il fenomeno in studio. Allo stesso modo, l’eterogeneità delle varianze può essere incorporata nel modello, utilizzando tecniche dei minimi quadrati generalizzati (GLS). In altri casi, quando non si riescono a rispettare le assunzioni di base dei modelli, si può ricorrere a metodiche statistiche che ne fanno di menoe che, pertanto, sono dette metodiche ‘non-parametriche’.

In questo libro, non tratteremo né GLiM, né i GLs, né le metodiche non parametriche. Per quello che riguarda GLiM e GLS, si tratta di metodiche che richiedono un corso di livello più avanzato; per quanto riguarda le metodiche non-parametriche, di esse non parleremo, per una questione di preferenze personali: a mio modo di vedere, utilizzare metodiche non-parametriche è come rinuciare in partenza a comprendere le basi biologiche del fenomeno e le intrinseche relazioni causa-effetto che sussistono nella realtà.

7.5.5 Trasformazioni stabilizzanti

Una strategia empirica, ma molto seguita in pratica e caratterizzata da un’efficacia non disprezzabile, è quella di ricorrere alle trasformazioni correttive. Con questa tecnica, si va a cercare una metrica sulla quale le assunzioni di base dell’ANOVA siano rispettate e si esegue l’elaborazione dei dati trasformati invece che di quelli non trasformati.

Per i conteggi e per l’eterogeneità delle varianze di variabili continue, si consiglia la trasformazione in radice quadrata o in logaritmo, scegliendo in base a quella che consente la miglior distribuzione dei residui. Per le proporzioni, taluni consigliano la trasformazione nell’arcoseno della radice del valore, che è implementata nel package ‘aomisc’, nella funzione ‘angularTransform()’. Questa funzione rivceve come input un valore percentuale e resituisce l’arcoseno della radice quadrata della proporzione corrispondente.

aomisc::angularTransform(c(26, 47, 25, 28, 24))
## [1] 30.65730 43.28009 30.00000 31.94806 29.33387

Per evitare di schegliere la trasformazione ‘al buio’, si può impiegare la procedura suggerita da Box e Cox (1964), basata su una ‘famiglia di trasformazioni’, definita come segue:

\[ W = \left\{ \begin{array}{ll} \frac{Y^\lambda}{\lambda - 1} & \quad \textrm{if} \,\,\, \lambda \neq 0 \\ \log(Y) & \quad \textrm{if} \,\,\, \lambda = 0 \end{array} \right. \]

dove W è la variabile trasformata, Y è la variabile originale e \(\lambda\) è il parametro che definisce la trasformazione. In particolare, osserviamo che se \(\lambda\) è pari ad 1 i dati non sono, di fatto, trasformati, se è pari a 0.5 abbiamo una trasformazione equivalente alla radice quadrata, se è pari a 0 abbiamo la trasformazione logaritmica, se è pari a -1 abbiamo una trasformazione equivalente al reciproco.

La scelta del valore \(\lambda\) può essere effettuata in modo empirico, confrontando la verosimiglianza di modelli basati su trasformazioni diverse. In R, questa procedura è automatizzata nella funzione ‘boxcox()’, disponibile nel package MASS è verrà illustrata nell’esempio succesivo.

7.6 Esempio

Proviamo ad analizzare il dataset ‘insects’, disponibile nel package ‘aomisc’. Si tratta di un dataset nel quale quindici piante sono state trattate con tre insetticidi diversi in modo completamente randomizzato, scegliendo cinque piante a caso per insetticida. Alcune settimane dopo il trattamento è stato rilevato il numero di insetti presenti su ciascuna pianta. Lasciando da parte le statistiche descrittive, eseguiamo subito l’ANOVA per questo dataset.

library(aomisc)
data(insects)
head(insects)
##   Insecticide Rep Count
## 1           A   1   448
## 2           A   2   906
## 3           A   3   484
## 4           A   4   477
## 5           A   5   634
## 6           B   1   211
mod <- lm(Count ~ Insecticide, data = insects)

Dopo aver effettuato l’ANOVA, i grafici dei residui possono essere ottenuti utilizzando la funzione ‘plot()’ applicata al risultato del fitting lineare. L’argomento ‘which’ specifica il tipo di grafico: se utilizziamo: ‘which = 1’ otteniamo il grafico dei residui verso gli attesi, se invece utilizziamo ‘which = 2’ otteniamo il QQ-plot. I due comandi sono:

plot(mod, which = 1)
plot(mod, which = 2)

e forniscono l’output riportato in figura 7.10.

Analisi grafica dei residui per il dataset 'insects.csv'

Figure 7.10: Analisi grafica dei residui per il dataset ‘insects.csv’

Questo dataset mostra una chiara eteroscedasticità (vedi il grafico di sinistra) e qualche indizio di asimmetria positiva.

car::leveneTest(mod)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  1.0263 0.3878
##       12
shapiro.test(residuals(mod) )
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mod)
## W = 0.87689, p-value = 0.04265

I test di Levene e Shapiro-Wilks confermano che la mancanza di normalità è significativa e, pertanto, scegliamo di impiegare una trasformazione correttiva. Utilizziamo quindi la funzione ‘boxcox()’ per individuare la trasformazione più adatta a correggere la patologia riscontrata. Il comando è:

library(MASS)
boxcox(mod)

e fornisce l’output in figura 7.11.

Scelta del valore ottimale di lambda per la trasformazione di Box e Cox

Figure 7.11: Scelta del valore ottimale di lambda per la trasformazione di Box e Cox

Vediamo che la verosimiglianza del modello raggiunge il massimo valore quando \(\lambda\) è pari a 0.14. I limiti dell’intervallo di confidenza vanno da poco sotto lo 0 a 0.5 circa. Rimanendo nell’ambito dell’intervallo di confidenza, scegliamo il valore \(\lambda = 0\), che corrisponde alla trasformazione logaritmica. Questa scelta è motivata dal fatto che si tratta di una trasformazione molto nota e facilmente comprensibile.

Pertanto, trasformiamo i dati nel logaritmo e ripetiamo l’ANOVA.

mod <- lm(log(Count) ~ Insecticide, data = insects)
par(mfrow = c(1,2))
plot(mod, which = 1)
plot(mod, which = 2)
Analisi grafica dei residui per il dataset 'insects.csv', previa trasformazione logaritmica

Figure 7.12: Analisi grafica dei residui per il dataset ‘insects.csv’, previa trasformazione logaritmica

Vediamo che i dati trasformati non mostrano più alcun sintomo di eteroscedasticità e, di conseguenza, l’ANOVA su questa metrica è totalmente affidabile. Ovviamente, avendo lavorato con il logaritmo dei dati, commentare i risultati potrebbe essere più complicato, in quanto dovremmo calcolare le medie marginali attese su scala logaritmica, come indicato nel codice sottostante.

library(emmeans)
medie <- emmeans(mod, ~Insecticide)
medie
##  Insecticide emmean    SE df lower.CL upper.CL
##  A             6.34 0.178 12     5.96     6.73
##  B             5.81 0.178 12     5.43     6.20
##  C             3.95 0.178 12     3.56     4.34
## 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95

È evidente che presentare le medie su scala logaritmica potrebbe non essere di immediata o facile lettura. Per questo, potremmo pensare di retrotrasformare le medie, utilizzando la trasformazione inversa di quella logaritmica, cioè l’antilogaritmo. Ad esempio, per la prima media:

\[e^{6.34} = 566.7963\]

In questo modo la nostra unità di misura ridiventa quella originale, anche se il valore ottenuto non coincide con la media dei dati originali; in effetti la trasformazione è non lineare e la media dei logaritmi non può coincidere con il logaritmo della media. Possiamo osservare che la media del trattamento A, sulla scala originale, è:

mean(insects[insects$Insecticide == "A","Count"])
## [1] 589.8

e risulta più alta della media retrotrasformata. In realtà, se è vero che i logaritmi sono normalmente distribuiti, la media dei logaritmi (6.34) dovrebbe essere uguale alla mediana (ricordiamo che in una distribuzione normale media e mediana coincidono). La mediana è il valore centrale; dato che la retro-trasformazione è monotona, il valore centrale resta centrale, anche se io retrotrasformo. Quindi la media retrotrasformata è uno stimatore della mediana della popolazione originale, non della sua media. Questo non è uno svantaggio: infatti il QQ-plot suggerisce un’asimmetria positiva (confronta la Figura 7.11 con la Figura 7.6 cosa che è confermata dal fatto che la mediana è minore della media. Se la distribuzione dei dati è asimmetrica, la mediana è un indicatore di tendenza centrale migliore della media, perché meno sensibile ai valori estremi, che sono più frequenti in caso di asimmetria.

Il problema è che, se vogliamo utilizzare la media retrotrasformata, dobbiamo trovare un valore di errore standard per questo stima, con il quale esprimere la sua incertezza. In realtà, anche l’errore standard può essere retrotrasformato, con una tecnica detta metodo ‘delta’, che costituisce un estensione della legge di propagazione degli errori per le trasformazioni non-lineari. È inutile andare nel dettaglio; diciamo solo che la funzione ‘emmeans()’ rende semplicissima l’implementazione del metodo delta, con il comando seguente:

retroMedie <- emmeans(mod, ~Insecticide, transform = "response")
retroMedie
##  Insecticide response     SE df lower.CL upper.CL
##  A              568.6 101.01 12    348.5      789
##  B              335.1  59.54 12    205.4      465
##  C               51.9   9.22 12     31.8       72
## 
## Confidence level used: 0.95

Con questo abbiamo tutto quello che ci serve: stime ed errori standard, che, ovviamente, sono diversi per le diverse medie retrotrasformate, coerentemente con la mancanza di omoscedasticità.


7.7 Per approfondire un po’ …

  1. Ahrens, W. H., D. J. Cox, and G. Budwar. 1990, Use of the arcsin and square root trasformation for subjectively determined percentage data. Weed science 452-458.
  2. Anscombe, F. J. and J. W. Tukey. 1963, The examination and analysis of residuals. Technometrics 5: 141-160.
  3. Box, G. E. P. and D. R. Cox. 1964, An analysis of transformations. Journal of the Royal Statistical Society, B-26, 211-243, discussion 244-252.
  4. D’Elia, A. 2001, Un metodo grafico per la trasformazione di Box-Cox: aspetti esplorativi ed inferenziali. STATISTICA LXI: 630-648.
  5. Saskia, R. M. 1992, The Box-Cox transformation technique: a review. Statistician 41: 169-178.
  6. Weisberg, S., 2005. Applied linear regression, 3rd ed. John Wiley & Sons Inc. (per il metodo ‘delta’)