Capitolo 3 Modelli matematici ed analisi dei dati

L’eredità galileiana ci porta ad immaginare che il funzionamento della natura sia basato su una serie di relazioni causa-effetto, descrivibili utilizzando il linguaggio universale della matematica. La conoscenza esatta di queste relazioni, nella teoria, ci potrebbe consentire di prevedere l’andamento di ogni fenomeno naturale, almeno quelli osservabili con sufficiente precisione.

In effetti era proprio questa l’ambizione più grande degli scienziati all’inizio dell’ottocento: conoscendo lo stato iniziale di un sistema, doveva essere possibile prevederne l’evoluzione futura. In realtà si è ben presto scoperto che si trattava di un’ambizione irrealistica, non tanto e non solo per la comparsa della meccanica quantistica un secolo dopo, ma anche per l’aumento di importanza degli studi in ambito psichiatrico e medico/biologico. Questi studi, infatti, venivano (e vengono) eseguiti su organismi viventi molto complessi, che, se sottoposti allo stesso stimolo, danno risposte altamente variabili e, spesso, anche difficilmente misurabili e controllabili. Immaginiamo quanto possa essere difficile quantificare uno stato legato ad una patologia mentale e individuare un pattern di risposta ad un certo stimolo, ad esempio farmacologico.

Queste difficoltà fecero prevalere, tra i biologi, la convinzione che la natura funzionasse in base a meccanismi deterministici ben definiti, anche se difficilmente osservabili nella pratica sperimentale, a causa dei numerosi elementi d’incertezza che si manifestavano nel corso delle osservazioni sperimentali. Insomma, la natura è perfetta, ma l’osservazione è fallace, perché influenzata dalla presenza di una forte componente stocastica e imprevedibile, che va sotto il nome generico di ’errore sperimentale’.

Dell’errore sperimentale abbiamo già parlato nei capitoli precedenti. Abbiamo anche visto che Ronald Fisher, nel suo famoso testo “Il disegno degli esperimenti” ha posto le basi per una corretta metodica sperimentale, volta a minimizzare l’impatto della componente stocastica e, soprattutto, ad impedire che essa possa confondersi con gli effetti degli stimoli sperimentali in studio. Minimizzare, tuttavia, non significa eliminare ed è evidente che, pur con tutti gli sforzi possibili, i risultati sperimentali saranno influenzati sempre e comunque da una certa quota di variabilità stocastica. Si vengono quindi a creare due contrapposte situazioni:

  1. la verità ‘vera’, immanente, che ‘agisce’ in base a meccanismi deterministici causa-effetto ben definiti.
  2. La ‘verità’ sperimentale, che si produce a partire dalla verità ‘vera’, per effetto dell’azione di elementi perturbativi casuali, che non ci permettono di osservare la verità ‘vera’.

Tenendo conto di questo, nella logica galileiana, possiamo provare ad utilizzare dei modelli matematici per descrivere le nostre osservazioni sperimentali e come esse si producono. Tuttavia, nessun modello è da considerare completo, se non è caratterizzato dalla capacità di descrivere sia la componente deterministica che quella stocastica, che sono alla base della formazione delle osservazioni sperimentali.

3.1 Verità ‘vera’ e modelli deterministici

In semplice linguaggio algebrico, possiamo porre un modello deterministico causa-effetto utilizzando una funzione di questo tipo:

\[ Y_E = f(X, \theta) \]

dove \(Y_E\) è l’effetto atteso dello stimolo \(X\), secondo la funzione \(f\), basata su una collezioni di parametri \(\theta\).

In questo modello vi sono una serie di componenti, che proviamo a guardare un po’ più nel dettaglio.

La risposta attesa (\(Y_E\)) è l’oggetto del nostro studio e può assumere le forme più disparate: spesso è numerica, ma a volte rappresenta una qualità. In questo libro consideremo solo la situazione in cui \(Y\) è rappresentato da una sola variabile (analisi univariata), ma esistono casi in cui viene osservata e analizzata la risposta di soggetti in relazione a molte variabili (analisi multivariata).

Lo stimolo sperimentale (\(X\)) è costituito da una o più variabili continue, discrete o categoriche, che rappresenta/ano il/i trattamento/i sperimentale/i (fattore/i sperimentale/i). Insieme ad \(Y\), è l’elemento noto di un esperimento, in quanto viene definito a priori con il disegno sperimentale.

La ‘funzione’ di risposta (\(f\)) è un’equazione, altrimenti detta ‘modello matematico’. L’equazione può essere lineare o non-lineare ed è selezionata o in base a considerazioni di carattere biologico, o con un approccio puramente empirico, nel quale osservo la risposta e scelgo un’equazione la cui forma si adatta bene ad essa.

I parametri di una funzione (\(\theta\)) sono un insieme di valori numerici che definiscono il modello. Nel prossimo paragrafo vedremo qualche esempio.

3.1.1 Qualche esempio di modello deterministico

Il modello più semplice è il cosiddetto modello della media:

\[ Y = \mu \]

Con questo modello si vuole indicare che un’osservazione dovrebbe conformarsi ad un certo valore atteso, come effetto di un unico stimolo sperimentale noto. Ha un solo parametro, cioè \(\mu\).

Un modello analogo, ma più complesso è il cosidetto modello ANOVA:

\[ Y = \left\{ {\begin{array}{ll} \mu_1 & se \quad X = A \\ \mu_2 & se \quad X = B \end{array}} \right. \]

In questo caso la risposta attesa dipende dallo stimolo sperimentale X, che è composto di due trattamenti: se il soggetto è trattato con A, fornisce la risposta \(\mu_1\), se è trattato con B fornisce \(\mu_2\). Il trattamento sperimentale è costituito da una variabile categorica con due modalità, ma l’estensione a più modalità è immediata. In questo caso, abbiamo due parametri (\(\mu_1\) e \(\mu_2\)).

Un ulteriore esempio di modello che vedremo in questo libro è la regressione lineare semplice, dove la relazione di risposta è descrivibile con una retta, la cui equazione generale è:

\[ Y = a + b \times X \]

In questo caso sia la \(Y\) che la \(X\) sono variabili quantitative e vi sono due parametri \(a\) e \(b\).

I modelli finora descritti sono lineari, ma esistono numerosi funzioni che descrivono relazioni curvilinee. Come, ad esempio, la parabola:

\[ Y = a + b \, X + c \, X^2\] caratterizzata da tre parametri, oppure la funzione esponenziale:

\[ Y = a \, e^{b X} \]

caratterizzata da due parametri (\(a\) e \(b\)), mentre \(e\) è l’operatore di Nepero. Delle funzioni curvilinea parleremo al termine di questo libro.

3.2 Genesi deterministica delle osservazioni sperimentali

Le equazioni esposte più sopra sono tutte nella loro forma generale e non sono utilizzabili se ai parametri non viene sostituito un valore numerico. E’proprio questo il nostro punto di partenza: un fenomeno scientifico può essere descritto attraverso un’equazione e dei parametri.

Ad esempio, potremmo immaginare che la risposta produttiva (Y) di una coltura (es. frumento) dipenda dalla dose della concimazione azotata (X), seguendo un andamento lineare, con \(a = 20\) e \(b = 0.3\). Rifacendosi alla notazione esposta sopra diremo che \(\theta\) è l’insieme dei valori 20 e 0.3.

Se questo assunto è vero, un eventuale esperimento di concimazione azotata, in assenza di qualunque fattore perturbativo, dovrebbe dare risultati assolutamento prevedili. Ad esempio, se concimiamo il frumento con 0, 30, 60, 90, 120 e 150 kg/ha di azoto, dovremmo osservare le risposte riportate qui di seguito.

X <- c(0, 30, 60, 90, 120, 150)
Ye <- 20 + 0.3 * X
Ye
## [1] 20 29 38 47 56 65

Qualcuno potrebbe obiettare che si tratta di una situazione irrealistica, perché la relazione tra concimazione azotata e produzione del frumento non è lineare. Non importa; in questo momento vogliamo semplicemente illustrare qual è la genesi delle osservazioni sperimentali. Fondamentalmente postuliamo l’esistenza di un modello matematico deterministico (causa-effetto) che, in assenza di errore sperimentale, è in grado di descrivere il comportamento della natura in una data situazione pedo-climatica.

3.3 Errore sperimentale e modelli stocastici

Tuttavia esiste un problema: l’errore sperimentale, puramente stocastico, confonde le nostre osservazioni e le rende diverse da quanto previsto dal modello deterministico. Come fare per incorporare questi effetti stocastici in un modello matematico? Ne parliamo immediatamente attraverso un esempio, semplice, ma concreto.

Immaginiamo un campo di frumento, di notevoli dimensioni. Un campo nella valle del Tevere, con milioni di piante geneticamente simili, in un ambiente abbastanza uniforme da un punto di vista del microclima. Immaginiamo di voler determinare l’altezza delle piante di questo campo. Immaginiamo di non avere limitazioni e di poter determinare l’altezza di tutte le piante dell’appezzamento. La nostra popolazione di soggetti diviene una popolazione di misure. Come è fatta questa popolazione? Quali sono le sue caratteristiche? Basiamoci sulla nostra esperienza professionale, cercando di utilizzare argomenti che possano essere condivisibili per l’intera comunità scientifica.

In primo luogo, possiamo dire che l’altezza delle piante dovrebbe essere quella dettata dal loro patrimonio genetico. Poniamo che questa altezza sia \(\mu = 100\); di conseguenza, il valore atteso di altezza è:

\[Y_E = \mu = 100\]

In realtà, questo valore \(\mu\) non è osservabile, perché, pur nella loro relativa omogeneità, l’ambiente e il patrimonio genetico non sono totalmente uguali, il che fa si che ogni pianta abbia la sua altezza, diversa da quella delle altre piante circostanti. Quindi non osserveremo esattamente \(Y_E\), ma, al suo posto osserveremo \(Y_O \neq Y_E\), con:

\[ Y_O = \mu + \varepsilon \]

dove \(\varepsilon\) è appunto una misura della componente stocastica individuale che distorce i risultati. Che valori potrebbe assumere \(\varepsilon\)? I valori individuali non possiamo conoscerli, ma possiamo fare alcune considerazioni probabilistiche: se il valore atteso è 100, trovare altezze comprese tra 99 e 101 cm dovrebbe essere molto più frequente che non trovare altezze pari a 40 cm o 180 cm. Quindi trovare valori di \(\varepsilon\) vicini allo 0 dovrebbe essere molto più frequente che non trovarli lontani.

Ci chiediamo, esiste una funzione che ci permetta di assegnare valori di probabilità alle diverse altezze che troviamo in un campo di frumento? La risposta è si: funzioni di questo tipo si chiamano funzioni di probabilità.

3.3.1 Funzioni di probabilità

Se avessimo rilevato una qualità del soggetto, come il sesso (M/F), la mortalità (vivo/morto), la germinabilità (germinato/non germinato), avremmo una variabile categorica nominale e potremmo calcolare le probabilità definita come rapporto tra il numero degli eventi favorevoli e il numero totale di eventi possibili (probabilità ‘frequentista’).

Ad esempio, immaginiamo di aver rilevato il numero di germogli ‘laterali’ di accestimento di 20 piante di frumento e di averne trovate 4 con 0 germigli, 6 con 1 germoglio, 8 con due germogli e 2 con tre germogli. La funzione che assegna la probabilità P ad ognuno dei quattro eventi X possibili (funzione di probabilità) è:

\[ P(X) = \left\{ \begin{array}{l} 4/20 = 0.2 \,\,\,\,\,\,se\,\,\,\,\,\,X = 0 \\ 6/20 = 0.3 \,\,\,\,\,\,se\,\,\,\,\,\,X = 1 \\ 8/20 = 0.4\,\,\,\,\,\,se\,\,\,\,\,\, X = 2 \\ 2/20 = 0.1 \,\,\,\,\,\,se\,\,\,\,\,\,X = 3 \\ \end{array} \right. \]

Viene ad essere definita una distribuzione di probabilità, che ha due caratteristiche importanti:

  1. P(X) è sempre non-negativo (ovvio! le probabilità sono solo positive o uguali a 0);
  2. la somma delle probabilità di tutti gli eventi è sempre pari ad 1 (ovvio anche questo: la probabilità che capiti uno qualunque degli eventi è sempre 1).

Se gli eventi possibili sono ordinabili (come nel caso precedente), oltre alla funzione di probabilità, si può definire anche la funzione di probabilità cumulata, detta anche funzione di ripartizione con la quale si assegna ad ogni evento la sua probabilità più quella di tutti gli eventi ‘inferiori’. Nell’esempio precedente:

\[ P(X) = \left\{ \begin{array}{l} 0.2\,\,\,\,\,\,se\,\,\,\,\,\,X \leq 0 \\ 0.5\,\,\,\,\,\,se\,\,\,\,\,\,X \leq 1 \\ 0.9\,\,\,\,\,\,se\,\,\,\,\,\,X \leq 2 \\ 1.0\,\,\,\,\,\,se\,\,\,\,\,\,X \leq 3 \\ \end{array} \right. \]

Per una distribuzione di probabilità come questa (classi numeriche ordinate), considerando il valore centrale di ogni classe, possiamo calcolare la media (valore atteso) come:

\[ \mu = E(X) = \sum{\left[ x_i \cdot P(X = x_i ) \right]} \]

e la varianza come:

\[\sigma ^2 = Var(X) = E\left[ {X - E(X)} \right]^2 = \sum{ \left[ {\left( {x_i - \mu } \right)^2 \cdot P(X = x_i )} \right]}\]

In questo caso specifico, la media è pari a:

mu <- 0 * 0.2 + 1 * 0.3 + 2 * 0.4 + 3 * 0.1
mu
## [1] 1.4

e la varianza è pari a:

(0 - mu)^2 * 0.2 + (1 - mu)^2 * 0.3 + (2 - mu)^2 * 0.3 +
  (3 - mu)^2 * 0.2
## [1] 1.06

Mediamente, le nostre piante hanno 1.4 germogli con una varianza pari a 1.06.

3.3.2 Funzioni di densità

Ciò che abbiamo detto finora non si applica al nostro caso, in quanto abbiamo rilevato una variabile continua (altezza) e non abbiamo intenzione di discretizzarla in classi. In teoria, le altezze che possiamo riscontrare sono pressoché infinite e non ha molto senso chiedersi, ad esempio, qual è la probabilità di trovare un individuo alto esattamente 96 cm (e non 96.0000001 cm, o 95.99999 cm). Capiamo da soli che questa probabilità è un infinitesimo.

Al contrario, come abbiamo visto più sopra, possiamo calcolare la probabilità di ottenere un valore compreso in un intervallo, per esempio da 80 a 90 cm. Tuttavia, abbiamo detto di non voler discretizzare, anche perché la probabilità dipenderebbe dall’ampiezza dell’intervallo prescelto, il che introdurrebbe un elemento arbitrario. Possiamo tuttavia pensare di calcolare la densità di probabilità, vale a dire il rapporto tra la probabilità di un intervallo e la sua ampiezza (cioè la probabilità per unità di ampiezza dell’intervallo; per questo si parla di densità). E’ evidente che se un intervallo diventa infinitamente piccolo anche la probabilità tende a zero con la stessa ‘velocità’, in modo che la densità di probabilità tende ad un numero finito (ricordate il limite del rapporto di polinomi?).

Insomma, con i fenomeni continui non possiamo lavorare con la probabilità dei singoli eventi, ma possiamo lavorare con la densità di probabilità e definire quindi apposite funzioni di densità. Analogamente alle funzioni di probabilità, le funzioni di densità debbono avere due caratteristiche:

  1. assumere solo valori non-negativi;
  2. la somma delle probabilità di tutti gli eventi possibili, calcolabile come integrale della funzione di densità, deve essere unitaria (anche in questo caso la densità di probabilità di tutti gli eventi possibili è pari ad 1).

Data una funzione di densità, possiamo costruire la corrispondente funzione di probabilità cumulata, facendo l’integrale per ogni evento pari o inferiore a quello dato. Più in generale, per variabili continue sia la funzione di ripartizione (probabilità cumulata), che la media o la devianza sono definite ricorrendo agli integrali:

\[ \begin{array}{l} P(X) = f(x) \\ P(X \le x) = \int\limits_{ - \infty }^x {f(x)} dx \\ \mu = E(X) = \int\limits_{ - \infty }^{ + \infty } {xf(x)} dx \\ \sigma ^2 = Var(X) = \int\limits_{ - \infty }^{ + \infty } {\left( {x - \mu } \right)^2 f(x)} dx \\ \end{array} \]

In pratica, vedremo che, a seconda della funzione di densità, è possibile adottare formule semplificate per le diverse statistiche descrittive.

3.4 La distribuzione normale (curva di Gauss)

Torniamo ancora alla nostra popolazione di misure, relative alle altezze del frumento nella media Valle del Tevere. E’ ragionevole pensare che, effettuando le misurazioni con uno strumento sufficientemente preciso e in presenza delle sole variazioni casuali (visto che abbiamo idealmente rimosso ogni differenza sistematica spiegabile), i risultati tendono a differire tra di loro, muovendosi intorno ad un valor medio, rispetto al quale le misure superiori ed inferiori sono equiprobabili e tendono ad essere più rare, via via che ci si allontana dal valore medio. Questo ragionamento ci porta verso una densità di probabilità (parliamo di variabili continue) a forma di campana, che potrebbe essere descritta con una funzione continua detta curva di Gauss.

La curva è descritta dalla seguente funzione di densità:

\[P(X) = \frac{1}{{\sigma \sqrt {2\pi } }}\exp \left[{\frac{\left( {X - \mu } \right)^2 }{2\sigma ^2 }} \right]\]

ove \(P(X)\) è la densità di probabilità di una certa misura \(X\), mentre \(\mu\) e \(\sigma\) sono rispettivamente la media e la deviazione standard della popolazione (per la dimostrazione si rimanda a testi specializzati). Le variabili casuali che possono essere descritte con la curva di Gauss, prendono il nome di variabili normali o normalmente distribuite.

Studiare le principali proprietà matematiche della curva di Gauss è estremamente utile. Ad esempio, senza voler entrare troppo nel dettaglio, guardando la curva di Gauss possiamo notare che:

  1. la forma della curva dipende da solo da \(\mu\) e \(\sigma\) (figura @ref{fig:figName51}). Ciò significa che, se prendiamo un gruppo di individui e partiamo dal presupposto (assunzione parametrica) che in relazione ad un determinato carattere quantitativo (es. produzione) la distribuzione di frequenza è normale (e quindi può essere descritta con una curva di Gauss), allora basta conoscere la media e la deviazione standard degli individui e immediatamente conosciamo l’intera distribuzione di frequenza (cioè l’intera popolazione di misure);
  2. la curva ha due asintoti e tende a 0 quando x tende a infinito. Questo ci dice che se assumiamo che un fenomeno è descrivibile con una curva di Gauss, allora assumiamo che tutte le misure sono possibili, anche se la loro frequenza decresce man mano che ci si allontana dalla media;
  3. la probabilità che la x assuma valori compresi in un certo intervallo è data dall’integrale della curva di Gauss in quell’intervallo. Ad esempio, la figura @ref{fig:figName52} mostra l’80° percentile, cioè la misura più alta dell’80% delle misure possibili;
  4. Se la curva di Gauss è stata costruita utilizzando le frequenze relative, l’integrale della funzione è uguale ad 1. Infatti la somma delle frequenze relative di tutte le varianti possibili non può che essere uguale ad 1;
  5. la curva è simmetrica. Questo indica che la frequenza dei valori superiori alla media è esattamente uguale alla frequenza dei valori inferiori alla media.
  6. dato \(\sigma\), possiamo dire che la frequenza dei valori superiori a \(\mu + \sigma\) è pari al 15.87% ed è uguale alla frequenza dei valori inferiori a \(\mu - \sigma\).
Distribuzioni normali con diversa media e deviazione standard (rispettivamente 5 e 1 a sinistra, 6 e 3 a destra

Figure 3.1: Distribuzioni normali con diversa media e deviazione standard (rispettivamente 5 e 1 a sinistra, 6 e 3 a destra

Integrale della curva di densità normale (80° percentile; sinistra) e curva di probabilità cumulata (destra)

Figure 3.2: Integrale della curva di densità normale (80° percentile; sinistra) e curva di probabilità cumulata (destra)

3.5 Modelli ‘a due facce’

A questo punto, sempre in relazione al nostro frumento, possiamo scrivere che l’altezza della pianta \(i\) è:

\[Y_i = \mu + \varepsilon_i\] dove abbiamo introdotto l’elemento stocastico \(\varepsilon\), per ogni soggetto \(i\), tale che:

\[ \varepsilon \sim N(0, \sigma) \]

cioè la componente stocastica \(\varepsilon\) è normalmente distribuita, con media 0 e deviazione standard \(\sigma\). E’abbastanza evidente che è possibile scrivere:

\[Y_i \sim N(\mu, \varepsilon)\]

cioè che l’altezza del frumento è normalmente distribuita con media \(\mu\) e deviazione standard \(\sigma\). Dato che si tratta di una semplice traslazione di una distribuzione normale lungo l’asse delle ascisse (come in figura @ref{figName51}), le due espressioni sono totalmente equivalenti.

Ora si può dire che conosciamo perfettamente la popolazione di partenza se conosciamo \(\mu\) e \(\sigma\), cioè la parte (‘faccia’) deterministica e la parte (‘faccia’) stocastica del modello. Se quindi immaginiamo che \(\mu = 100\) (come abbiamo detto in precedenza) e \(\sigma = 8\), possiamo risolvere alcuni semplici esercizi, utilizzando le funzioni di calcolo di probabilità di R.

ESERCIZIO 1. Calcolare la densità di un’altezza pari a 120 cm.

dnorm(120, mean = 100, sd = 8)
## [1] 0.002191038

ESERCIZIO 2. Qual è la probabilità di ottenere piante con altezza inferiore a 80 cm?

pnorm(80, mean = 100, sd = 8)
## [1] 0.006209665

ESERCIZIO 3. Qual è la probabilità di ottenere piante con altezza superiore a 110 cm?

pnorm(110, mean = 100, sd = 8, lower.tail = F)
## [1] 0.1056498

Si utilizza l’argomento lower.tail=FALSE, in quanto stiamo cercando la probabilità di un’a concentrazione pari o superiore ad 1.1, e non pari od inferiore.’altezza pari o superiore a 110 cm (upper-tail) e non quella pari o inferiore a 110 cm (lower-tail), che sarebbe fornita di default. E’ totalmente equivalente utilizzare la funzione sottostante.

1 - pnorm(110, mean = 100, sd = 8)
## [1] 0.1056498

ESERCIZIO 4. Qual è la probabilità di ottenere piante con altezza compresa tra 80 e 110 cm?

pnorm(110, mean = 100, sd = 8) - pnorm(80, mean = 100, sd = 8)
## [1] 0.8881406

ESERCIZIO 5. Qual è quella misura che è superiore al 90% di tutte le piante del campo (90° percentile?

qnorm(0.9, 100, 8)
## [1] 110.2524

ESERCIZIO 6. Qual è quella misura che è inferiore al 20% di tutte le piante del campo (80° percentile?

qnorm(0.8, 100, 8)
## [1] 106.733
qnorm(0.2, 100, 8, lower.tail=F)
## [1] 106.733

ESERCIZIO 7. Quali sono quei due valori, simmetrici rispetto alla media e tali da formare un intervallo all’interno del quale cadono il 95% delle piante?

qnorm(0.025, 100, 8)
## [1] 84.32029
qnorm(0.975, 100, 8)
## [1] 115.6797

3.6 E allora?

Cerchiamo di ricapitolare. Le popolazioni di soggetti sperimentali e delle loro misure sono un oggetto largamente ignoto e inconoscibile. Infatti, le caratteristiche dei soggetti della popolazione sono, in parte, determinate in base a relazioni causa-effetto, ma, in altra parte, esse sono puramente stocastiche. Tuttavia è ragionevole supporre che esse seguano una qualche funzione di probabilità/densità (assunzione parametrica). Se questo è vero, allora possiamo utilizzare queste funzioni e i loro integrali per calcolare la probabilità di ottenere una certa misura o un certo insieme di misure.

3.7 Le simulazioni Monte Carlo

Se quanto abbiamo detto è vero, ogni esperimento scientifico non è altro che un’operazione di campionamento da una certa distribuzione di probabilità. Questo campionamento può essere simulato impiegando un generatore di numeri casuali. Immaginiamo di avere disegnato un esperimento con otto parcelle (repliche) per determinare la produzione del mais in un certo appezzamento e immaginiamo che queste parcelle costituiscano un campione rappresentativo della popolazione di parcelle di quell’appezzamento. Se la popolazione è distribuita normalmente, con media pari a 12 t/ha e deviazione standard pari ad 1.2 t/ha, allora possiamo simulare i risultati del nostro esperimento come segue:

set.seed(1234)
Y <- rnorm(8, 15, 1.2)
Y
## [1] 13.55152 15.33292 16.30133 12.18516 15.51495 15.60727 14.31031 14.34404
mean(Y)
## [1] 14.64344
sd(Y)
## [1] 1.328183

La generazione di numeri casuali con il computer viene fatta attraverso algoritmi che, a partire da un seed iniziale, forniscono sequenze che obbediscono a certe proprietà fondamentali (numeri pseudo-casuali). Il comando ‘set.seed(1234)’ ci permette di partire dallo stesso valore e quindi di ottenere lo stesso campione. Un’altra cosa da notare è che il nome della funzione che genara numeri casuali è formato con il nome della distribuzione (‘norm’) più il prefisso ‘r’. Questo è vero per tutte le altre distribuzioni in R (‘rbinom’, ‘rt’ e così via)

I valori campionati non riflettono le caratteristiche della popolazione, nel senso che la media e la deviazione standard del campione differiscono da quelle della popolazione. E’ esattamente ciò che capita durante un esperimento!

3.8 Analisi dei dati e ‘model fitting’

Secondo il principio illustrato in questo capitolo, un ricercatore arriverebbe a conoscere perfettamente la realtà se riuscisse ad individuare l’equazione e i parametri che governano il fenomeno in studio. Di conseguenza, l’ipotesi scientifica che sta alla base di un esperimento può essere posta sotto forma di modello matematico. Ad esempio, potremmo ipotizzare che la degradazione di un erbicida nel terreno segua una legge di decadimento esponenziale, rappresentabile, in genere, con l’equazione:

\[ C= a \, e^{-k T} + \varepsilon\]

dove \(C\) è la concentrazione dell’erbicida in un dato momento \(T\) ed \(a\) e \(k\) sono i parametri. Per quanto riguarda l’elemento stocastico, possiamo assumere che:

\[ \varepsilon \sim N(0, \sigma)\]

In modo equivalente, potremmo scrivere:

\[ C \sim N(a \, e^{-k T}, \sigma)\]

Questo modello è assolutamente generale; se vogliamo specificarlo per una certa situazione sperimentale, ad esempio per la degradazione di imazethapyr a 20°C, possiamo realizzare un esperimento nel quale contaminamo un terreno con questo erbicida, lo mettiamo a 20°C e, in tempi diversi, preleviamo aliquote di terreno da sottoporre a determinazione analitica. L’analisi dei dati raccolti consisterà nell’individuare \(a\), \(k\) e \(\sigma\), con una tecnica definita di model fitting.

Le diverse tecniche di analisi dei dati che descriveremo nei capitoli successivi sono accomunate dall’essere appunto tecniche di model fitting. Vedremo anche come queste tecniche possono essere utilizzate per verificare che le osservazioni sperimentali si conformino ad un dato modello (goodness of fit) oppure per confrontare due ipotesi alternative poste sotto forma di modelli diversi (model comparison).


3.9 Per approfondire un po’…

3.9.1 Generazione dei dati sperimentali: un esempio più realistico

Immaginiamo di conoscere la legge deterministica che definisce la risposta del frumento alla concimazione azotata. In particolare, immaginiamo che questa risposta produttiva sia fondamentalmente lineare:

\[Y_E = b_0 + b_1 X\]

ed immaginiamo che, senza concimazione (X = 0), la produzione sia pari a 25 q/ha (\(b_0\) = 25). Immaginiamo che l’incremente produttivo per kg di azoto somministrato sia pari a 0.15 q/ha (\(b_1\) = 0.15).

Per individuare questa legge naturale organizziamo un esperimento, con quattro dosi di azoto e quattro repliche. In questo esperimento, come in tutti gli esperimenti, agirà anche una componente stocastica, che in qualche modo sposterà la risposta osservata dalla risposta attesa:

\[Y_o = b_0 + b_1 X + \varepsilon \quad \textrm{con} \quad \varepsilon \sim N(0, \sigma)\]

Assumiamo che la componente stocastica \(\varepsilon\) sia distribuita normalmente, come media 0 e deviazione standard \(\sigma\), che immaginiamo essere pari a 2.5 q/ha.

Su questa base, generiamo i dati osservati, nelle seguenti fasi:

  1. generiamo i valori attesi (‘Yield_E’), utilizzando il modello deterministico di risposta. Le quattro repliche di una dose hanno, ovviamente, lo stesso valore;
  2. aggiungiamo l’elemento stocastico, specifico per ogni osservazione, campionando da una distribuzione normale;
  3. creiamo un dataframe con le osservazioni sperimentali
set.seed(1234)
Dose <- rep(c(0, 60, 120, 180), each=4) 
Yield_E <- 25 + 0.15 * Dose
epsilon <- rnorm(16, 0, 2.5)
Yield <- Yield_E + epsilon
dataset <- data.frame(Dose, Yield)
dataset
##    Dose    Yield
## 1     0 21.98234
## 2     0 25.69357
## 3     0 27.71110
## 4     0 19.13576
## 5    60 35.07281
## 6    60 35.26514
## 7    60 32.56315
## 8    60 32.63342
## 9   120 41.58887
## 10  120 40.77491
## 11  120 41.80702
## 12  120 40.50403
## 13  180 50.05937
## 14  180 52.16115
## 15  180 54.39874
## 16  180 51.72429
# write.csv(dataset, file = "Nwheat.csv", row.names = F)

Il dataset che otteniamo è realistico e lo utilizzeremo per l’analisi (regressione lineare semplice) nel capitolo 12.

3.9.2 Modelli stocastici non-normali

Oltre alla distribuzione gaussiana, che è largamente la più importante, esistono molti altri modelli stocastici, sia per eventi continui che discreti. Chi volesse approfondire queste distribuzioni trova informazioni in seguito. Per gli altri vogliamo solo far notare che le funzioni di R per il calcolo di probabilità hanno sempre la stessa sintassi, che, dato il nome della distribuzione (es. ‘norm’ per la distribuzione normale), assegna il prefisso ‘d’ per la funzione di probabilità/densità, ‘p’ per la funzione di probabilità cumulata, ‘q’ per la funzione quantile. Ad esempio, per la distribuzione binomiale, il nome è ‘binom’ e, di conseguenza, le funzioni in R sono: ‘dbinom()’ (probabilità), ‘pbinom()’ (probabilità cumulata) e ‘qbinom()’ (funzione ‘quantile’).

3.9.2.1 t di Student

La distribuzione t di Student è analoga, per forma, ad una distribuzione normale con media 0 e deviazione standard 1. Rispetto a quest’ultima, la dispersione è un po’ più ampia, nel senso la probabilità di avere valori lontani dalla media è più alta. In realtà, non esiste una sola distribuzione t di Student, ma ne esistono molte, caratterizzate da un diverso numero di gradi di libertà (\(\nu\)); maggiore è \(\nu\), minore la sovradispersione; se il numero di gradi di libertà è infinito, la distribuzione t di Student è identica alla normale standardizzata (distribuzione normale con media 0 e deviazione standard uguale ad 1).

La Figura 3.3 riporta un esempio di come diminuisce la sovradispersione all’aumentare del numero di gradi di libertà.

Distribuzione t di Student, con diversi gradi di libertà

Figure 3.3: Distribuzione t di Student, con diversi gradi di libertà

3.9.2.2 F di Fisher

La distribuzione F di Fisher è definita solo per valori positivi ed è fortemente asimmetrica. Anche in questo caso, abbiamo una famiglia di distribuzioni, che differiscono tra di loro per due parametri (gradi di libertà) \(\nu_1\) e \(\nu_2\). Solitamente questa distribuzione viene utilizzata per descrivere il rapporto tra le varianze di coppie di campioni estratti da un distribuzione normale standardizzata, per cui \(\nu_1\) e \(\nu_2\) sono i gradi di libertà del numeratore e del denominatore.

La Figura 3.4 mostra diverse funzioni di densità: possiamo notare come, all’aumentare dei gradi di libertà, la funzione F di Fisher tende ad assomigliare sempre di più ad una normale.

Distribuzioni F di Fisher, con $
u_1 = 
u_2 = 
u$ e $
u = 1$ (linea nera), $
u = 3$ (linea rossa), $
u = 10$ (linea blues), $
u = 30$ (linea violetta) e $
u = 60$ (linea marrone)

Figure 3.4: Distribuzioni F di Fisher, con $ u_1 = u_2 = u$ e $ u = 1$ (linea nera), $ u = 3$ (linea rossa), $ u = 10$ (linea blues), $ u = 30$ (linea violetta) e $ u = 60$ (linea marrone)

3.9.2.3 La distribuzione binomiale

Ogni esperimento per il quale ci siano solo due esiti possibili (successo ed insuccesso) e una certa probabilità di successo, viene detto esperimento Bernoulliano. Il tipico esempio è il lancio della moneta, nel quale possiamo ottenere solo testa o croce, con una probabilità di 0.5 (se la moneta non è truccata). In alcuni casi, potremmo avere una serie di esperimenti Bernoulliani indipendenti, con probabilità di successo costante (ad esempio, lanciare la moneta 10 volte) e potremmo essere interessati a conoscere la probabilità di ottenere k successi su n prove. Questa probabilità può essere descritta attraverso la funzione di probabilità binomiale.

Poniamo di sapere che in una Facoltà di Agraria con un numero molto elevato di studenti, il rapporto tra maschi e femmine sia pari a 0.7 e quindi che la probabilità di incontrare un maschio sia pari a \(p = 0.7\) (evento semplice). Deve essere estratto a sorte un viaggio studio per quattro studenti e, per una questione di pari opportunità, si preferirebbe che fossero premiati in ugual misura maschi e femmine (cioè si vogliono premiare due femmine). Qual è la probabilità che un simile evento si realizzi?

La probabilità cercata si può ottenere pensando che abbiamo un evento “estrazione” che può dare due risultati possibili (maschio o femmina) e che deve essere ripetuto quattro volte. Se consideriamo “successo” estrarre una femmina, allora la probabilità di successo in ogni estrazione è \(p = 0.3\) mentre quella di insuccesso (evento complementare) è pari a \(1 - p = q = 0.7\). Facciamo attenzione! Quanto abbiamo detto è vero solo se la popolazione è sufficientemente numerosa da pensare che la singola estrazione non cambia la probabilità degli eventi nelle successive (eventi indipendenti). La probabilità che su quattro estrazioni si abbiano 2 successi (due femmine) e due insuccessi (due maschi) è data da (teorema della probabilità composta):

\[0.3 \cdot 0.3 \cdot 0.7 \cdot 0.7 = 0.3^2 \cdot 0.7^2\]

In generale, data una popolazione molto numerosa, nella quale gli individui si presentano con due modalità possibili (in questo caso maschio e femmina) e posto di sapere che la frequenza con cui si presenta la prima modalità è pari a \(p\) (in questo caso la frequenza delle femmine è pari a 0.3), mentre la frequenza della seconda modalità è pari a \(q = 1 - p\), se vogliamo estrarre da questa popolazione \(n\) elementi, la probabilità che \(k\) di questi presentino la prima modalità (successo) è data da:

\[p^k \cdot q^{(n-k)}\]

La formula di cui sopra, tuttavia, non risolve il nostro problema, in quanto noi vogliamo che vengano estratte due femmine, indipendentemente dall’ordine con cui esse vengono estratte (prima, seconda, terza o quarta estrazione), mentre la probabilità che abbiamo appena calcolato è quella relativa all’evento in cui le due femmine sono estratte al primo e secondo posto.

Di conseguenza (teorema della probabilità totale), alla probabilità dell’evento indicato in precedenza (estrazione di due femmine in prima e seconda posizione) dobbiamo sommare la probabilità di tutti gli altri eventi utili (due femmine in seconda e terza posizione, oppure in terza e seconda, oppure in terza e quarta e così via). Il numero delle combinazioni possibili per 2 femmine in quattro estrazioni (combinazione di 4 elementi di classe 2) è dato dal coefficiente binomiale:

\[\left( {\begin{array}{*{20}c} n \\ k \\ \end{array}} \right) = \frac{n!}{(n - k)!k!} \]

Moltiplicando le due equazioni date in precedenza otteniamo la funzione di probabilità binomiale:

\[P(X = x_i ) = \frac{{n!}}{{(n - k)!k!}} \cdot p^k \cdot q^{(n - k)} \]

Nel caso specifico otteniamo il risultato:

\[P(X = 2) = \frac{4!}{(4 - 2)!2!} \cdot 0.3^2 \cdot 0.7^{(4 - 2)} = 0.2646 \]

che è appunto la probabilità cercata.

In R, utilizziamo la funzione ‘dbinom(successi, prove, probabilità semplice)’ per calcolare la probabilità di ottenere \(k\) successi in \(n\) prove:

dbinom(2, 4, 0.3)
## [1] 0.2646

La funzione binomiale è un modello stocastico e si può dimostrare che il valore atteso (media) è uguale ad \(n\cdot p\), mentre la varianza è pari a \(n\cdot p \cdot q\):

La funzione di ripartizione (probabilità cumulata) si calcola in R con la funzione ‘pbinom(successi, prove, probabilità semplice)’. Nell’esempio, se vogliamo sapere la probabilità totale di estrarre meno di tre femmine (2 femmine o meno), possiamo operare in questo modo:

pbinom(2,4,0.3)
## [1] 0.9163

Che risulta anche dalla somma della probabilità di estrarre 0, 1, 2 femmine:

zero <- dbinom(0,4,0.3)
uno <- dbinom(1,4,0.3)
due <- dbinom(2,4,0.3)
zero + uno + due
## [1] 0.9163

La funzione di ripartizione può anche essere utilizzata al contrario, per determinare i quantili, cioè il numero di successi che corrispondono ad una probabilità cumulata pari ad alfa:

qbinom(0.9163,4,0.3)
## [1] 2

3.9.3 Altre letture

  1. Bolker, B.M., 2008. Ecological models and data in R. Princeton University Press, Books.
  2. Schabenberger, O., Pierce, F.J., 2002. Contemporary statistical models for the plant and soil sciences. Taylor & Francis, CRC Press, Books.