Test di Fisher

Quando si effettuano dei conteggi nei quali le osservazioni sono organizzate in una tabella di due righe per due colonne, si applicano queste regole:

→ si impiega il test chi-quadrato se le osservazioni sono indipendenti e sono numerose;

→ si impiega il test di Fisher se le osservazioni sono indipendenti e sono poche;

→ si impiega il test di McNemar nel caso di osservazioni non indipendenti (cioè nel caso di dati appaiati).

Quindi bisogna intendersi su cosa significano "numerose" e "poche". Armitage [1] e Snedecor [2] indicano la necessità di impiegare il test di Fisher (detto anche test esatto di Fisher) in alternativa al test chi-quadrato quando:

→ il numero totale di osservazioni è inferiore a 20;

→ il numero totale di osservazioni è compreso tra 20 e 40 e vi sono frequenze attese inferiori a 5.

Le indicazioni che potete trovare su altri test di statistica sono meno vincolanti sui limiti da applicare in termini di numerosità delle osservazioni, mentre concordano tutte nel sottolineare la necessità di applicare il test di Fisher quando in una o più celle vi è un numero ridotto di frequenze attese, le frequenze che dovremmo osservare in teoria e che sono calcolate per ciascuna cella come prodotto delle somme marginali (di riga e di colonna) corrispondenti, diviso per il numero totale dei casi e così ad esempio nella tabella riportata qui sotto:

→ per la cella in alto a sinistra contenente il valore osservato 4 la frequenza attesa sarà (4+16)·(4+1)/(4+16+1+21) = 20·5/42 = 2.38 

→ per la cella in basso a destra con il valore osservato 21 la frequenza attesa sarà (21+1)·(21+16)/(4+16+1+21) = 22·37/42 = 19.38 

Il test di Fisher lo applichiamo all'analisi della relazione che potrebbe intercorrere (e che vogliamo verificare) tra il tipo di allattamento, naturale o artificiale, e la presenza in età successiva nel bambino di malocclusione dentaria. La suzione da una tettarella in gomma non è del tutto fisiologica e potrebbe determinare un qualche effetto sulla morfogenesi delle arcate dentarie che è in atto nel lattante. I dati sono tratti da [3] e sono riportati in questa tabella che contiene quindi le frequenze osservate:

Questo script prevede di inserire i dati manualmente. Copiatelo e incollatelo nella Console di R e premete ↵ Invio:

# TEST DI FISHER - 2 righe · 2 colonne

#

cells <- c(4,16,1,21) # genera l'array cells con i valori numerici contenuti nelle celle

rnames <- c("Allattamento_naturale", "Allattamento_artificiale") # genera l'array rnames con i nomi delle righe

cnames <- c("Denti_normali", "Malocclusione") # genera l'array cnames con i nomi delle colonne

mydata <- matrix(cells, nrow=2, ncol=2, byrow=TRUE, dimnames=list(rnames, cnames)) # genera la matrice dei dati

mydata # mostra i dati

fisher.test(mydata) # esegue il test di Fisher

#

Nelle prime tre righe sono generati con la funzione c():

→ il vettore che contiene i quattro dati (che come si vede devono essere inseriti in sequenza leggendoli da sinistra a destra e dall'alto in basso) che sono salvati nell'oggetto cells;

→ il vettore che contiene i nomi delle righe, salvato nell'oggetto rnames;

→ il vettore che contiene i nomi delle colonne, salvato nell'oggetto cnames.

I tre vettori sono combinati a formare la matrice dei dati mediante la funzione matrix() che impiega gli argomenti che indicano:

→ l'oggetto/vettore contenente i dati (cells);

→ il numero di righe (nrow=2) e il numero di colonne (ncol=2) della matrice che sarà pertanto di due righe per due colonne;

→ la modalità di riempimento della matrice, che deve essere riempita per righe (byrow=TRUE) quindi da sinistra a destra e dall'alto in basso;

→ i nomi da assegnare alle righe e alle colonne (dimnames=list(rnames, cnames)).

La matrice dei dati così costruita viene salvata nell'oggetto mydata, che viene infine mostrato per un controllo finale.

Come si vede la parte principale del lavoro consiste nella corretta immissione e strutturazione dei dati, mentre il calcolo vero e proprio del test viene effettuato in modo molto semplice, con la funzione fisher.test() [4].

Questi sono i risultati:

> mydata # mostra i dati 

                         Denti_normali Malocclusione

Allattamento_naturale                4            16

Allattamento_artificiale             1            21

> fisher.test(mydata) # esegue il test di Fisher

        Fisher's Exact Test for Count Data

data:  mydata

p-value = 0.1745

alternative hypothesis: true odds ratio is not equal to 1

95 percent confidence interval:

   0.4440411 270.1636947

sample estimates:

odds ratio 

  5.059936 

Considerando significativi i valori inferiori al tradizionale valore soglia p=0.05 il test di Fisher fornisce nel nostro caso come risultato un valore di p=0.1745 che non è significativo: pertanto la conclusione è che non abbiamo evidenza statistica del fatto che i due tipi di allattamento comportino differenze nella prevalenza di malocclusione.

La funzione fisher.test() fornisce una analisi aggiuntiva della tabella calcolando anche l'odds ratio. In questo caso un risultato significativamente diverso da 1 deporrebbe per una associazione tra tipo di allattamento e tipo di occlusione dentaria. E un risultato diverso da 1 si avrebbe se i limiti di confidenza dell'odds ratio trovato, che è uguale a 5.059936, non includessero il valore 1. Dato che lo includono (vanno infatti da 0.4440411 a 270.1636947) il risultato non è significativo, conferma la mancanza di associazione tra i fattori presenti nella tabella e conferma anche il valore di p ottenuto con il test di Fisher. Per il tema odds ratio si rimanda ai test di statistica: si trova tra gli altri in Marubini [5], Campbell [6] e Ingelfinger [7].

Abbiamo qui una tipica applicazione del buonsenso in statistica: con un numero totale di osservazioni superiore a 40 (sono in totale 42 osservazioni) avremmo dovuto a rigore applicare il test chi-quadrato, ma è stato applicato il test di Fisher in quanto è stato dato maggior peso al fatto che vi sono celle con un numero ridotto di osservazioni. A conferma di ciò potete eseguire il test chi-quadrato digitando

chisq.test(mydata)

con questo risultato:

> chisq.test(mydata)

        Pearson's Chi-squared test with Yates' continuity correction

data:  mydata

X-squared = 1.1398, df = 1, p-value = 0.2857

Messaggio di avvertimento:

In chisq.test(mydata) : Chi-squared approximation may be incorrect

Il test chi-quadrato riporta una differenza non significativa (p=0.2857) ma un messaggio avverte che l'impiego del test chi-quadrato in questo caso non è appropriato. La ragione? Se digitate

chisq.test(mydata)$expected

ottenete la tabella delle frequenze attese

> chisq.test(mydata)$expected

                         Denti_normali Malocclusione

Allattamento_naturale         2.380952      17.61905

Allattamento_artificiale      2.619048      19.38095

Messaggio di avvertimento:

In chisq.test(mydata) : Chi-squared approximation may be incorrect

nella quale notate sulla sinistra due valori inferiori a 5. Il messaggio "Chi-squared approximation may be incorrect" viene fornito dalla funzione chisq.test() in quanto anche in R la regola "frequenze attese inferiori a 5" rappresenta il razionale per l'impiego del test di Fisher e prevale sul numero delle osservazioni.

In alternativa potete procedere come segue.

Copiate le tre righe riportate qui sotto, incollatele in un editor di file di testo aggiungendo un ↵ Invio al termine dell'ultima riga, e salvate il tutto in C:\Rdati\ sotto forma di un file di testo denominato chi_Fisher.csv (di default i file di testo sono salvati con l'estensione .txt ma noi qui salviamo il file con l'estensione .csv) [8].

Allattamento;Denti_normali;Malocclusione

Allattamento_naturale;4;16

Allattamento_artificiale;1;21

Impiegate quest'altro script che prevede di eseguire il test di Fisher sui dati letti dal file chi_Fisher.csv precedentemente salvato. Copiatelo e incollatelo nella Console di R e premete ↵ Invio:

# TEST DI FISHER - 2 righe · 2 colonne 

#

mydata <- read.table("C:/Rdati/chi_Fisher.csv", header=TRUE, sep=";", row.names="Allattamento") # importa i dati

mydata # mostra i dati

fisher.test(mydata) # esegue il test di Fisher

#

Come vedete lo script è più compatto del precedente, le uniche cose da notare sono gli argomenti della funzione read.table():

→ "C:/Rdati/chi_Fisher.csv" che specifica nome e posizione del file dal quale importare i dati;

→ header=TRUE che indica che nella prima riga del file sono contenuti i nomi delle variabili;

→ sep=";" che specifica il separatore di campo impiegato nel file;

→ row.names="Allattamento" che indica che i nomi delle righe sono contenuti nel campo “Allattamento”.

Salvate entrambi gli script. Se preferite intervenire sullo script per adattarlo di volta in volta ai nuovi dati, potete impiegare il primo dei due. Se volete evitare di mettere mano ogni volta allo script e preferite intervenire sul file di dati, potete impiegare il secondo.

----------

[1] Armitage P. Statistica medica. Giangiacomo Feltrinelli Editore, MIlano, 1979, p.140.

[2] Snedecor GW, Cochran WG. Statistical Methods. The Iowa State University Press, 1980, ISBN 0-8138-1560-6, p. 127.

[3] Armitage P. Statistica medica. Giangiacomo Feltrinelli Editore, MIlano, 1979, pp. 138-140.

[4] Digitate help(fisher.test) nella Console di R per la documentazione della funzione fisher.test().

[5] Bossi A, Cortinovis I, Duca PG, Marubini E. Introduzione alla statistica medica. La Nuova Italia Scientifica, Roma, 1994, ISBN 88-430-0284-8, pp. 282-285.

[6] Campbell MJ, Machin D. Medical Statistics. A Commonsense Approach. John Wiley & Sons, New York, 1993, ISBN 0-471-93764-9, pp. 121-123 e 154-156.

[7] Ingelfinger JA, Mosteller F, Thibodeau LA, Ware JH. Biostatistica in medicina. Macmillan Publishing Co., Inc.; New York, 1983, ISBN 88-7078-065.1, pp. 30-33.

[8] Per eventuali approfondimenti sul formato .csv (comma separated value) che è il formato dati raccomandato per R vedere la sezione Gestione dei dati nella pagina Indice.

Come riportare i numeri in formato fisso

Se nel confronto tra due numeri sono interessato soprattutto al loro ordine di grandezza la notazione scientifica [1] va benissimo. Ma per un confronto puntuale preferisco il formato fisso.

Così, giusto per fare un esempio, quando devo confrontare due valori di probabilità p anziché p=8.327701e-03 e p=6.412771e-02 preferisco leggere p=0.008327701 e  p=0.06412771.

Copiate lo script che segue, incollatelo nella Console di R e premete ↵ Invio:

# RIPORTARE I NUMERI IN FORMATO FISSO o nella notazione scientifica

#

col_1 <- c(0.00003230609, 0.01850436, 0.008327701) # prima colonna

col_2 <- c(0.6991678, 0.00005211903, 0.02457609) # seconda colonna

mydataframe <- data.frame(col_1, col_2) # combina le due colonne in una tabella

#

mydataframe # mostra la tabella

#

options(scipen=999) # esprime i numeri in formato fisso

mydataframe # mostra la tabella

#

options(scipen=0) # ripristina la notazione scientifica

mydataframe # mostra la tabella

#

Come si vede a fronte dell'inserimento in una tabella di numeri espressi in formato fisso (0.00003230609, 0.01850436, 0.008327701, 0.6991678, 0.00005211903, 0.02457609) di default R riporta i numeri in notazione scientifica

> mydataframe # mostra la tabella

         col_1        col_2

1 3.230609e-05 6.991678e-01

2 1.850436e-02 5.211903e-05

3 8.327701e-03 2.457609e-02

ma consente di passare all'espressione dei numeri in formato fisso impiegando all'interno della funzione options() [2] l'argomento scipen=999 

> options(scipen=999) # esprime i numeri in formato fisso

> mydataframe # mostra la tabella

          col_1         col_2

1 0.00003230609 0.69916780000

2 0.01850436000 0.00005211903

3 0.00832770100 0.02457609000

e di ripristinare l'espressione dei numeri in notazione scientifica impiegando all'interno della funzione options() l'argomento scipen=0 

> options(scipen=0) # ripristina la notazione scientifica

> mydataframe # mostra la tabella

         col_1        col_2

1 3.230609e-05 6.991678e-01

2 1.850436e-02 5.211903e-05

3 8.327701e-03 2.457609e-02

Nota bene: in Windows nella cartella C:\Programmi\R\R-n.n.n\etc (dove n.n.n sono i numeri della versione di R installata) si trova il file Rprofile.site [3] che viene eseguito automaticamente quando si lancia R. Nel file Rprofile.site possono essere salvate le options() che si desidera siano impiegate da R di default. Il file può essere elaborato con un comune editor di testo.

Se volete che i numeri siano di default rappresentati in formato fisso, inserite nel file Rprofile.site la riga

options(scipen=999) 

Ovviamente potete, in alternativa, utilizzare all'interno dei vostri script options(scipen=999) ogniqualvolta preferite il formato fisso (soluzione consigliata), ripristinando poi la notazione scientifica con options(scipen=0).

Infine due indicazioni per completare il tema:

→ il razionale per decidere quale deve essere il numero di cifre significative da impiegare nella rappresentazione di un numero è riportato nel post Come stabilire il giusto numero di cifre significative; 

→ per esprimere in pratica i risultati di una elaborazione statistica, che sono generalmente ottenuti con un numero di cifre in eccesso rispetto a quelle significative, vedere il post Come arrotondare i numeri.

----------

[1] La notazione scientifica è una notazione esponenziale, cioè una notazione che consente di scrivere un numero N assegnato come prodotto di un opportuno numero a per la potenza di un altro numero b^k essendo b la base della notazione esponenziale. In particolare la notazione scientifica è una notazione esponenziale in base 10 (b = 10). Tuttavia "notazione esponenziale" non è sinonimo di "notazione scientifica" in quanto esistono notazioni esponenziali che impiegano una base diversa da quella impiegata dalla notazione scientifica: ne sono un esempio la notazione esponenziale in base 2 (b = 2) impiegata in informatica, e la notazione esponenziale in base e (b = e) - essendo il numero di Nepero e = 2.718 281 828 459... - impiegata in campo matematico e in campo scientifico.

[Precisazione aggiuntiva: in R il separatore delle cifre decimali è il punto (.) e come già riportato altrove questa convenzione per ragioni di omogeneità viene adottata non solo negli script ma anche nei file di dati e in tutto il testo di questo sito, quindi anche nell'espressione del numero e di Nepero. Per le norme che regolano la punteggiatura all'interno dei numeri rimando al post Come separare i decimali e come raggruppare le cifre].

[2] Digitate help(options) nella Console di R per scoprire gli altri argomenti della funzione options(). 

[3] Se avete installato la versione a 32 bit di R in Windows trovate il file Rprofile.site nella cartella C:\Programmi(x86)\R\R-n.n.n\etc (dove n.n.n sono sempre i numeri della versione di R installata).

Regressione polinomiale di terzo grado

Supponiamo di avere, per ciascuno di una serie di valori prefissati x (riportati in ascisse), il risultato di altrettante misure y (riportate in ordinate) e con la funzione plot() visualizziamo i dati sotto forma di un grafico xy.

Copiate il codice che segue nella Console di R e premete ↵ Invio:

#

x <- c(20, 40, 60, 80, 100, 120, 140, 160, 180, 200, 220, 240, 260, 280, 300) # valori in ascisse

y <- c(0.695, 0.699, 0.707, 0.722, 0.740, 0.765, 0.790, 0.815, 0.840, 0.865, 0.885, 0.902, 0.920, 0.932, 0.940) # valori in ordinate

plot(y~x, col="black", pch=1, main="Grafico xy", xlab="Valori prefissati (x)", ylab="Risultati delle misure effettuate (y)") # traccia il grafico xy

#

Dalla ispezione del grafico risulta evidente una relazione non lineare.

Questo ci fa escludere la possibilità di descrivere la relazione di funzione che intercorre tra le due variabili mediante una retta. E pone il tema di quale sia la curva che potrebbe meglio descriverla.

La prima cosa da valutare è se esiste un modello deterministico di un evento fisico, chimico-fisico, biologico, econometrico o quant'altro, che si assume descriva adeguatamente la relazione di funzione che intercorre tra le due variabili: in tal caso non resta che applicare l'equazione fornita dal modello. In caso contrario è necessario ricorrere ad un modello empirico: ed è questo che consideriamo essere il caso nostro.

Qualora si debba ricorrere ad un modello empirico, la regola fondamentale è impiegare l'equazione più semplice compatibile con una adeguata descrizione dei dati. Nel nostro caso visto che è evidente dall'ispezione del grafico che l'equazione di una retta 

y = a + b· ​⁠⁠x

calcolata mediante la classica regressione lineare con il metodo dei minimi quadrati sarebbe inadeguata a descrivere la relazione di funzione che lega la y alla x, possiamo pensare di ricorrere ad una equazione appena più complessa, come quella fornita da una regressione polinomiale di secondo grado nella forma

y = a + b· x + c· x²

e calcolata anch'essa con il metodo dei minimi quadrati [1].

Copiate lo script nella Console di R e premete ↵ Invio: 

# REGRESSIONE POLINOMIALE di secondo grado

#

x <- c(20, 40, 60, 80, 100, 120, 140, 160, 180, 200, 220, 240, 260, 280, 300) # valori in ascisse

y <- c(0.695, 0.699, 0.707, 0.722, 0.740, 0.765, 0.790, 0.815, 0.840, 0.865, 0.885, 0.902, 0.920, 0.932, 0.940) # valori in ordinate

#

options(scipen=999) # esprime i numeri in formato fisso

#

# calcola la regressione 

#

polisec <- lm(y~poly(x, degree=2, raw=TRUE)) # calcola la regressione

summary(polisec) # riporta le statistiche della regressione

#

# traccia il grafico

#

newx <- seq(min(x), max(x), (max(x)-min(x))/1000) # valori x in corrispondenza dei quali calcolare il valore del polinomio

newy <- predict(polisec, list(x=newx)) # calcola il valore corrispondente del polinomio

plot(y~x, col="black", pch=1, main="Regressione polinomiale di secondo grado", xlab="Valori prefissati (x)", ylab="Risultati delle misure effettuate (y)") # predispone il grafico e riporta i dati

lines(newx, newy, col="black", lty=2, lwd=1) # traccia il polinomio

#

options(scipen=0) # ripristina la notazione scientifica

#

La prima cosa che notiamo è che il polinomio di secondo grado descrive in pratica una retta.

Inoltre se andate alla quinta riga di codice tra le statistiche generate con la funzione summary(polisec) trovate:

Coefficients:

                                     Estimate   Std. Error t value             Pr(>|t|)    

(Intercept)                      0.6573076923 0.0100102431  65.664 < 0.0000000000000002 ***

poly(x, degree = 2, raw = TRUE)1 0.0009554000 0.0001439488   6.637             0.000024 ***

poly(x, degree = 2, raw = TRUE)2 0.0000001299 0.0000004374   0.297                0.772    

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Mentre l'intercetta a con un p=0.0000000000000002 e il coefficiente di primo grado b del polinomio con un p=0.000024 sono significativi, il coefficiente di secondo grado c con un p=0.772 non lo è affatto. Pertanto l'equazione del polinomio di secondo grado

y = a + b· x + c· x²  

privata del coefficiente c· x² si riduce all'equazione di una retta

y = a + b · x   

cosa che conferma il grafico ottenuto.

Dati che descrivono una curva con un raggio di curvatura omogeneo, senza massimi, senza minimi e senza punti di flesso - cioè senza punti di passaggio tra due curvature che vanno in senso opposto - sono compatibili con il ramo di una parabola descritta appunto da un polinomio di secondo grado. Ma in questo caso i dati presentano un punto di flesso tra una curvatura iniziale rivolta verso l'alto e una curvatura finale rivolta verso il basso. Il polinomio di secondo grado media tra le due, finendo con il fornire una curvatura talmente ridotta da potere essere assimilata ad una retta.

Dobbiamo pertanto vedere cosa accade applicando il modello immediatamente successivo, rappresentato da una regressione polinomiale di terzo grado nella forma

 y = a + b· x + c· x² + d· x³     

Copiate lo script nella Console di R e premete ↵ Invio: 

# REGRESSIONE POLINOMIALE di terzo grado

#

x <- c(20, 40, 60, 80, 100, 120, 140, 160, 180, 200, 220, 240, 260, 280, 300) # valori in ascisse

y <- c(0.695, 0.699, 0.707, 0.722, 0.740, 0.765, 0.790, 0.815, 0.840, 0.865, 0.885, 0.902, 0.920, 0.932, 0.940) # valori in ordinate

#

options(scipen=999) # esprime i numeri in formato fisso

#

# calcola la regressione 

#

politer <- lm(y~poly(x, degree=3, raw=TRUE)) # calcola la regressione

summary(politer) # riporta le statistiche della regressione

#

# traccia il grafico

#

newx <- seq(min(x), max(x), (max(x)-min(x))/1000) # valori x in corrispondenza dei quali calcolare il valore del polinomio

newy <- predict(politer, list(x=newx)) # calcola il valore corrispondente del polinomio

plot(y~x, col="black", pch=1, main="Regressione polinomiale di terzo grado", xlab="Valori prefissati (x)", ylab="Risultati delle misure effettuate (y)") # predispone il grafico e riporta i dati

lines(newx, newy, col="black", lty=2,lwd=1) # traccia il polinomio

#

options(scipen=0) # ripristina la notazione scientifica

#

Le prime due righe di codice servono semplicemente ad assegnare (<-) i valori inclusi nella funzione c() al vettore x e al vettore y che in questo modo conterranno i nostri dati.

Poi mediante la funzione options() con l'argomento scipen=999 facciamo si che i risultati siano espressi in formato fisso, cosa che li renderà meglio leggibili (questo almeno a mio modo di vedere: se volete lasciare i risultati nella notazione scientifica eliminate questa riga di codice o, forse meglio, anteponete a questa riga di codice il simbolo #).

La quarta riga di codice assegna (<-) all'oggetto politer i risultati del calcolo della regressione polinomiale y~poly(x) di terzo grado (degree=3), l'argomento raw=TRUE è impiegato perché non ci interessano i polinomi ortogonali (l'opzione di default è raw=FALSE).

Se siete interessati ad approfondimenti, ricordo che potete visualizzare ad esempio la struttura dell'oggetto politer digitando str(politer) e potete avere informazioni sulle funzioni impiegate nello script digitando help(nomedellafunzione).

Alla quinta riga di codice la funzione summary(politer) consente infine di visualizzare le statistiche generate:

Call:

lm(formula = y ~ poly(x, degree = 3, raw = TRUE))

Residuals:

       Min         1Q     Median         3Q        Max 

-0.0032624 -0.0012957  0.0006167  0.0014060  0.0019740 

Coefficients:

                                         Estimate       Std. Error t value             Pr(>|t|)    

(Intercept)                       0.6966205128205  0.0025292032241 275.431 < 0.0000000000000002 ***

poly(x, degree = 3, raw = TRUE)1 -0.0003180913177  0.0000662417487  -4.802             0.000552 ***

poly(x, degree = 3, raw = TRUE)2  0.0000097653837  0.0000004730826  20.642       0.000000000381 ***

poly(x, degree = 3, raw = TRUE)3 -0.0000000200739  0.0000000009739 -20.612       0.000000000387 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.001865 on 11 degrees of freedom

Multiple R-squared:  0.9997,    Adjusted R-squared:  0.9996 

F-statistic: 1.081e+04 on 3 and 11 DF,  p-value: < 0.00000000000000022

La voce Residuals riporta le statistiche delle differenze che, dopo l'approssimazione dei punti con il polinomio di terzo grado, residuano tra i punti stessi e la curva calcolata: differenze piccole come quelle qui osservate indicano un'ottima approssimazione della curva ai dati.

La voce Coefficients riporta dall'alto verso il basso il valore stimato (Estimate) dell'intercetta a (0.6966205128205), del coefficiente di primo grado b (-0.0003180913177), del coefficiente di secondo grado c (0.0000097653837) e del coefficiente di terzo grado d (-0.0000000200739) e ci consente di scrivere l'equazione del polinomio che sarà quindi 

y = 0.6966205128205 - 0.0003180913177· x + 0.0000097653837· x² - 0.0000000200739· x³     

   

Per ciascun coefficiente viene calcolato il valore del test t di Student come rapporto tra il valore stimato (Estimate) e il suo errore standard (Std. Error) e viene infine riportato il valore di probabilità (Pr(>|t|)) di osservare per caso il valore di t riportato: se tale probabilità è sufficientemente bassa ci assumiamo il rischio di considerare significativo il coefficiente in questione. Nel nostro caso essendo i valori di p sempre molto piccoli (di gran lunga inferiori al valore p = 0.05 generalmente considerato come valore soglia per la significatività) possiamo affermare che tutti e tre i coefficienti del polinomio di terzo grado sono significativi e lo rendono un'equazione in grado di descrivere adeguatamente i dati forniti. 

Il coefficiente di correlazione multipla R² (Multiple R-squared: 0.9997) molto elevato  con una statistica F significativa (p-value: < 0.00000000000000022) conferma ulteriormente la bontà dell'approssimazione. E niente meglio del grafico dei punti e del polinomio che li approssima può confermarla.

Il grafico è tracciato in modo da garantire una risoluzione grafica elevata:

→ l'intervallo compreso tra il valore minimo min(x) e il valore massimo max(x) della variabile in ascisse viene suddiviso per ottenere 1000 valori (ma potete cambiarne il numero) salvati nell'oggetto newx e in corrispondenza dei quali calcolare il valore del polinomio;

→ con la funzione predict() sono calcolati e riportati nell'oggetto newy i valori delle ordinate del polinomio corrispondenti ai valori delle ascisse contenuti in newx ;

→ la funzione plot() viene impiegata per riportare sul grafico di dati contenuti in x e in y;

→ la funzione lines() viene impiegata per sovrapporre ai dati il polinomio.

Potete riportare i dati modificando il simbolo impiegato (pch=...) e il suo colore (col="..."), come pure tracciare il polinomio modificando l'aspetto o stile della linea tracciata (lty=...), il suo colore (col="...") e il suo spessore (lwd=...) [2].

Da notare che il grafico è stato tracciato esclusivamente all'interno dei dati, quindi in un'ottica di interpolazione: e questo ci ricorda anche che, trattandosi di un modello empirico, è da escludere l'impiego del polinomio per la estrapolazione, cioè per trarre conclusioni al di fuori dei dati impiegati per calcolarlo. 

A chiudere lo script è la funzione options() che con l'argomento scipen=0 ripristina l'espressione dei numeri nella notazione scientifica [3].

Lo script può essere riutilizzato semplicemente riportando manualmente in x e in y i nuovi dati, o meglio ancora importandoli da un file [4].

----------

[1] Il metodo dei minimi quadrati prevede - ed è quello che abbiamo assunto essere il nostro caso - che la variabile in ascisse (x) abbia valori prefissati, che la variabile in ordinate (y) sia rappresentata da misure (affette da un errore casuale distribuito normalmente) effettuate corrispondenza delle x e minimizza la somma dei quadrati delle differenze y - y' ovvero delle differenze tra ciascun valore misurato y e il valore y' = f(x) corrispondente calcolato mediante la regressione (lineare o polinomiale). Esempi possono essere la risposta y di un sensore misurata in corrispondenza di una serie di concentrazioni date x di una sostanza chimica, la misurazione di un evento y in corrispondenza di una serie data di tempi x, e così via.

[2] Trovate le relative indicazioni:

→ nel post Cambiare i simboli dei punti di R, 

→ nel post Visualizzare i colori disponibili in R 

→ nel post Cambiare gli stili delle linee di R

[3] La notazione scientifica è una notazione esponenziale, cioè una notazione che consente di scrivere un numero N assegnato come prodotto di un opportuno numero a per la potenza di un altro numero b^k essendo b la base della notazione esponenziale. In particolare la notazione scientifica è una notazione esponenziale in base 10 (b = 10). Tuttavia "notazione esponenziale" non è sinonimo di "notazione scientifica" in quanto esistono notazioni esponenziali che impiegano una base diversa da quella impiegata dalla notazione scientifica: ne sono un esempio la notazione esponenziale in base 2 (b = 2) impiegata in informatica, e la notazione esponenziale in base e (b = e) - essendo il numero di Nepero e = 2.718 281 828 459... - impiegata in campo matematico e in campo scientifico.

[Precisazione aggiuntiva: in R il separatore delle cifre decimali è il punto (.) e come già riportato altrove questa convenzione per ragioni di omogeneità viene adottata non solo negli script ma anche nei file di dati e in tutto il testo di questo sito, quindi anche nell'espressione del numero e di Nepero. Per le norme che regolano la punteggiatura all'interno dei numeri rimando al post Come separare i decimali e come raggruppare le cifre].

[4]  Alla pagina Indice nel capitolo R - gestione dei dati trovate i collegamenti ai post che contengono gli script per importare i dati da file esterni.