Capitolo 14 Appendix 3: Per chi vuole approfondire un po’…

[Intro da fare]

14.1 Capitolo 3: Progettare un esperimento

14.1.1 Organizzare un esperimento di diserbo chimico

Si suppone che gli erbicidi A, B e C siano più efficaci di D, E ed F verso Solanum nigrum, una comune pianta infestante delle colture di pomodoro. L’obiettivo generale della ricerca sarà quello di trovare un’efficace soluzione per l’eliminazione di Solanum nigrum dal pomodoro. Gli obiettivi specifici saranno:

  1. valutare l’efficacia erbicida di A, B e C, confrontandola con quella di D, E ed F
  2. valutare la selettività degli anzidetti erbicidi verso il pomodoro

Il fattore sperimentale oggetto di studio sarà il diserbo del pomodoro, con 5 livelli inseriti in prova (6 trattamenti sperimentali): A, B, C, D, E ed F. Inoltre, si ritiene opportuno inserire in prova un testimone non trattato (NT), che ci permetterà di quantificare la percentuale di malerbe controllate. In totale, avremo quindi sette tesi sperimentali.

Questo esperimento verrà eseguito in vaso e, di conseguenza, potremo realizzare sei repliche con un disegno sperimentale a randomizzazione completa. La variabile rilevata, tre settimane dopo il trattamento, sarà il peso della biomassa presente in ogni vasetto.

14.1.2 Organizzare un esperimento di confronto varietale

L’ipotesi è che le varietà di girasole A, B e C non abbiano la stessa base genetica e quindi non siano tutte ugualmente produttive. L’obiettivo generale è quello di capire quale tra A, B e C sia più adatta alle condizioni pedoclimatiche della collina Umbra. Gli obiettivi specifici sono quelli di valutare:

  1. produttività di A, B e C
  2. stabilità produttiva di A, B e C

Il fattore sperimentale in studio sarà la varietà di girasole con 3 livelli inclusi in prova (varietà A, B e C). Come testimone, inseriremo la varietà di riferimento per la zona (D). Dato che eseguiremo questa prova su un terreno nel quale vi sono due chiari gradienti di fertilità, disegneremo l’esperimento considerando due fattori di blocco: trasversale e longitudinale (spiego meglio tra poco…). Poiché dobbiamo valutare la stabilità produttiva, dovremo ripetere l’esperimento più volte (es. in tre anni diversi) e quindi avremo un secondo fattore sperimentale, incrociato con il primo.

Questo esperimento verrà realizzato in pieno campo, su parcelle di dimensioni 2 m x 8 m, seguendo uno schema sperimentale a quadrato latino con quattro repliche. Dovendo misurare la stabilità produttiva, cioè l’oscillazione di produzione da un ambiente all’altro, questa prova dovrà essere ripetuta in più anni (es. tre anni).

Per ognuno degli anni di prova, la mappa contiene una griglia 4 x 4, nella quale possiamo identificare quattro colonne e quattro righe. Dato che abbiamo presupposto l’esistenza di un gradiente trasversale e lungitudinale (tra righe e tra colonne), l’allocazione dei trattamenti dovrà esser fatta in modo che ognuno di essi si trovi su ogni riga e ogni colonna (Quadrato latino). Un’aspetto fondamentale è comunque quello di definire una diversa randomizzazione in ogni anno/località, per evitare che le stesse varietà siano sempre nelle stesse posizioni, che potrebbe dare origine a dubbi di confounding. La definizione delle randomizzazioni per il secondo e terzo anno è lasciata per esercizio.

Anche in questo caso potremo chiedere ad R di aiutarci a trovare la combinazione corretta (anche se questo potrebbe essere comodamente fatto a mano).

library(agricolae)
## Warning: package 'agricolae' was built under R version 3.5.2
## 
## Attaching package: 'agricolae'
## The following object is masked from 'package:aomisc':
## 
##     AMMI
trt <- c("A", "B", "C", "D")
designLS <- design.lsd(trt, seed=543, serie=2)
designLS$book
##    plots row col trt
## 1    101   1   1   C
## 2    102   1   2   A
## 3    103   1   3   B
## 4    104   1   4   D
## 5    201   2   1   D
## 6    202   2   2   B
## 7    203   2   3   C
## 8    204   2   4   A
## 9    301   3   1   B
## 10   302   3   2   D
## 11   303   3   3   A
## 12   304   3   4   C
## 13   401   4   1   A
## 14   402   4   2   C
## 15   403   4   3   D
## 16   404   4   4   B
Schema sperimentale a quadrato latino per l'Esempio 2 (un anno)

Figure 14.1: Schema sperimentale a quadrato latino per l’Esempio 2 (un anno)

Un’altro aspetto da considerare è la metodica impiegata per la determinazione del peso di 1000 semi. Abbiamo già visto che, per aumentare la precisione e la rappresentatività, da tutta la granella raccolta da una parcella preleviamo quattro lotti da 1000 semi, di cui determinare il peso. In questo modo, per ogni trattamento avremo 16 valori (quattro repliche x quattro lotti per replica). Ovviamente non possiamo affermare di avere 16 repliche, in quanto solo le parcelle sono da considerare repliche, in quanto ricevono il trattamento (varietà) in modo indipendente. I quattro lotti raccolti da ogni parcella sono unità osservazionali (perché ne viene rilevato il peso), ma non unità sperimentali, perché appartengono alla stessa parcella e non sono indipendenti. I quattro lotti si dicono sub-repliche, quindi il disegno ha quattro repliche e quattro sub-repliche per replica (disegno a quadrato latino con sottocampionamento). I due strati di errore (variabilità tra repliche e variabilità tra sub-repliche entro replica), devono essere mantenuti separati in fase di analisi, altrimenti l’analisi è invalida, perché è condotta come se avessimo un più alto grado di precisione (16 repliche) rispetto a quello che abbiamo effettivamente (una sorta di millantato credito!).

[Inserire immagine]

Al termine del ciclo colturale, si misurerà il peso di mille semi. Per questo, prenderemo dalla produzione di granella di ogni parcella, quattro sub-campioni da mille semi, da sottoporre a pesate.

14.1.3 Organizzare un esperimento fattoriale

Nella barbabietola da zucchero, il diserbo localizzato lungo la fila consente di diminuire l’impiego di erbicidi. Tuttavia, se la coltura precedente ha prodotto semi e se non abbiamo effettuato una lavorazione profonda per interrarli, la coltura sarà più infestata e quindi sarà più difficile ottenere una buona produttività con il diserbo parziale. Su questa ipotesi costruiamo un esperimento volto a valutare l’interazione tra lavorazione del terreno e diserbo chimico. Per raggiungere questo obiettivo generale, proveremo a valutare se:

  1. il diserbo parziale consente di ottenere produzioni comparabili a quelle del diserbo totale
  2. l’effetto erbicida è indipendente dalla lavorazione prescelta

In questo caso avremo due fattori sperimentali incrociati: il diserbo, con due livelli (totale o parziale, localizzato sulla fila) e la lavorazione, con tre livelli (aratura profonda, aratura superficiale e minimum tillage). Non vi è la necessità di un testimone, ma avremo la necessità di un fattore di blocco. In totale, avremo sei tesi sperimentali.

In questo caso abbiamo un disegno fattoriale con due livelli a blocchi randomizzati. Nel principio, questo disegno non ha nulla di diverso da quello relativo all’esempio 1, fatto salvo un minor numero di trattamenti (solo 6). Anche in questo caso, ci facciamo aiutare da R.

trt <- c(3,2) # factorial 3x2
design2way <-design.ab(trt, r=4, serie=2,
  design="rcbd", seed=777)
book <- design2way$book
levels(book$A) <- c("PROF", "SUP", "MIN")
levels(book$B) <- c("TOT", "PARZ")
book
##    plots block    A    B
## 1    101     1  SUP PARZ
## 2    102     1 PROF PARZ
## 3    103     1 PROF  TOT
## 4    104     1  MIN  TOT
## 5    105     1  SUP  TOT
## 6    106     1  MIN PARZ
## 7    107     2  MIN  TOT
## 8    108     2  SUP  TOT
## 9    109     2  MIN PARZ
## 10   110     2 PROF  TOT
## 11   111     2  SUP PARZ
## 12   112     2 PROF PARZ
## 13   113     3  MIN  TOT
## 14   114     3  SUP  TOT
## 15   115     3 PROF PARZ
## 16   116     3  MIN PARZ
## 17   117     3  SUP PARZ
## 18   118     3 PROF  TOT
## 19   119     4  MIN PARZ
## 20   120     4 PROF  TOT
## 21   121     4 PROF PARZ
## 22   122     4  MIN  TOT
## 23   123     4  SUP  TOT
## 24   124     4  SUP PARZ

La mappa risultante è visibile più sotto.

Schema sperimentale fattoriale a blocchi randomizzati per l'Esempio 3

Figure 14.2: Schema sperimentale fattoriale a blocchi randomizzati per l’Esempio 3

Questo disegno è totalmente appropriato, ma ci costringe a lasciare parecchio spazio tra una parcella e l’altra, per poter manovrare con la macchina per la lavorazione del terreno. Sarebbe utile raggruppare le parcelle caratterizzate dalla stessa lavorazione, in modo da poter lavorare su superfici più ampie. Ne guadagnerebbe l’uniformità dell’esperimento e l’accuratezza dei risultati. Possiamo quindi immaginare un disegno a un fattore, con parcelle di dimensione doppia (main-plots), sulle quali eseguire, in modo randomizzato le lavorazioni del terreno. Succesivamente, ogni main-plot può essere suddivisa in due e, su ognuna delle due metà, possono essere allocati in modo random i due trattamenti di diserbo. In questo modo ci troviamo ad operare con parcelle di due dimensioni diverse: le main-plots per le lavorazioni e le sub-plots per il diserbo. Questo tipo di schema prende il nome di parcella suddivisa (split-plot), ed è piuttosto comune nella sperimentazione di pieno campo.

Proviamo ad utilizzare R per redigere il piano sperimentale.

lavorazione <- c("PROF", "SUP", "MIN")
diserbo <- c("TOT", "PARZ")
designSPLIT <- design.split(lavorazione, diserbo,
  r=4, serie=2, seed=777)
book <- designSPLIT$book
book
##    plots splots block lavorazione diserbo
## 1    101      1     1         SUP    PARZ
## 2    101      2     1         SUP     TOT
## 3    102      1     1        PROF     TOT
## 4    102      2     1        PROF    PARZ
## 5    103      1     1         MIN    PARZ
## 6    103      2     1         MIN     TOT
## 7    104      1     2         SUP    PARZ
## 8    104      2     2         SUP     TOT
## 9    105      1     2         MIN     TOT
## 10   105      2     2         MIN    PARZ
## 11   106      1     2        PROF     TOT
## 12   106      2     2        PROF    PARZ
## 13   107      1     3         MIN     TOT
## 14   107      2     3         MIN    PARZ
## 15   108      1     3         SUP     TOT
## 16   108      2     3         SUP    PARZ
## 17   109      1     3        PROF     TOT
## 18   109      2     3        PROF    PARZ
## 19   110      1     4        PROF    PARZ
## 20   110      2     4        PROF     TOT
## 21   111      1     4         MIN     TOT
## 22   111      2     4         MIN    PARZ
## 23   112      1     4         SUP    PARZ
## 24   112      2     4         SUP     TOT
Schema sperimentale split-plot a blocchi randomizzati per l'Esempio 3

Figure 14.3: Schema sperimentale split-plot a blocchi randomizzati per l’Esempio 3

In alcune circostanze, soprattutto nelle prove di diserbo chimico, potrebbe trovare applicazione un altro tipo di schema sperimentale, nel quale, in ogni blocco, un trattamento viene applicato a tutte le parcelle di una riga e l’altro trattamento a tutte le parcelle di una colonna. Ad esempio, il disegno sottostante mostra una prova nella quale il terreno è stato diserbato in una striscia nel senso della lunghezza e, dopo il diserbo, le colture sono state seminate in striscia, nel senso della larghezza. Questo disegno è detto strip-plot ed è molto comodo perché consente di lavorare velocemente.

Schema sperimentale a strip-plot

Figure 14.4: Schema sperimentale a strip-plot

14.1.4 Organizzare un esperimento con una coltura poliennale

Vogliamo porre a confronto tre varietà di erba medica (A, B e C) e, considerando che l’erba medica è una coltura poliennale, vogliamo capire se il giudizio di merito è indipendente dall’anno di coltivazione. I nostri obiettivi specifici saranno:

  1. valutare la produttività media delle varietà in prova
  2. valutare le oscillazione nei quattro anni di durata del cotico erboso

Il fattore sperimentale in studio sarà la varietà di erba medica con 3 livelli inclusi in prova (varietà A, B e C) ai quali aggiungiamo il riferimento di zona (D) come testimone. Come nel caso del girasole, dovremo valutare la stabilità produttiva negli anni, ma, dato che abbiamo una coltura poliennale, non avremo bisogno di ripetere la prova, ma potremo ripetere le osservazioni per quattro anni sulla stessa prova.

La prova di erba medica è fondamentalmente un esperimento a blocchi randomizzati, il cui piano è riportato più sotto. Tuttavia, si tratta di una coltura poiliennale nella quale ripeteremo le misurazioni ogni anno sulle stesse parcelle. le misure ripetute non sono randomizzate (non possono esserlo), ma seguono una metrica temporale. Proprio per questo sviluppo lungo la scala del tempo, i dati che si raccolgono in questi esperimenti a misure ripetute sono detti dati longitudinali. Guardando bene il disegno si capisce anche per si parla di split-plot nel tempo. Esempi affini sono relativi all’analisi di accrescimento con misure non distruttive (esempio l’altezza) oppure i prelievi di terreno a profondità diverse, anche se, in quest’ultimo caso, la metrica delle misure ripetute è spaziale, non temporale.

Si può notare una certa analogia con il sottocampionamento illustrato più sopra, nel senso che vengono prese più misure per parcella. Tuttavia, bisogna tener presente che nel sottocampionamento le diverse misure sono solo repliche e non vi è nessuna esigenza di distinguere tra quelle prese nella stessa parcella. Invece, nel caso delle misure ripetute ognuna di esse ha interesse individuale, in quanto espressione di un’anno particolare.

Schema sperimentale a blocchi randomizzati con misure ripetute (split-plot in time)

Figure 14.5: Schema sperimentale a blocchi randomizzati con misure ripetute (split-plot in time)

14.1.5 Utilizzare R per disegnare gli esperimenti

Negli esperimenti più semplici lo schema sperimentale può essere pianificato a mano. Per esperimenti complessi potremo invece utilizzare il computer; in R, potremo utilizzare, ad esempio, il package agricolae (de Mendiburu 2019), seguento il codice che troverete nei paragrafi seguenti.

[Spostare qui gli esempi, lasciando sopra gli schemi]

14.2 Capitolo 4: Modelli matematici a ‘due facce’

14.2.1 La distribuzione t di Student

La distribuzione t di Student è analoga per forma ad una distribuzione normale con media 0 e deviazione standard 1. Rispetto a questa, 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).

Per verificare l’entità della sovradispersione, proviamo a disegnare su un grafico una curva normale standardizzata ed una serie di curve di t, con 2, 6 e 24 gradi di libertà.

par(mfrow = c(1, 1))
curve(dnorm(x),-3, +3, col="Black", xlab="",
     ylab="Densità")
curve(dt(x, 2), add=TRUE, col = "blue")
curve(dt(x,6), add=TRUE, col = "red")
curve(dt(x,24), add=TRUE, col = "green")

14.2.2 La distribuzione 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.

Col codice che segue, possiamo disegnare la distribuzione di F con \(\nu_1 = \nu_2 = 3\) e possiamo calcolare la probabilità di estrarre da questa distribuzione un valore pari o superiore a 5. Inoltre, calcoliamo anche il 95° percentile, utilizzando le apposite funzioni in R.

curve(df(x, 3, 3), 0, +3,col="Black",
     xlab="", ylab="Densità")

pf(5, 3, 3, lower.tail = F)
## [1] 0.890449
qf(0.95, 3, 3)
## [1] 9.276628

14.2.3 La distribuzione binomiale

Ogni esperimento per il quale ci sono 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

14.2.3.1 Esercizio

Da una popolazione di insetti che ha un rapporto tra maschi e femmine pari a 0.5, qual è la probabilità di campionare casualmente 2 maschi e 8 femmine?

dbinom(2, 10, 0.5)
## [1] 0.04394531

14.2.3.2 Esercizio

Riportare su un grafico la funzione di ripartizione binomiale, per p=0.5 e n=5. Costruire anche la densità di frequenza, utilizzando le opportune funzioni R.

prob <- 0.5
n <- 5
barplot(dbinom(seq(0, n, by=1), size=n, prob=prob),
          xlab="Successi", ylab="Probabilità",
          names.arg=seq(0,5))
Distribuzione di probabilità binomiale (sinistra) e probabilità binomiale cumulata (destra)

Figure 14.6: Distribuzione di probabilità binomiale (sinistra) e probabilità binomiale cumulata (destra)

barplot(pbinom(seq(0, n, by=1), size=n, prob=prob),
          xlab="Successi", ylab="Probabilità",
          names.arg=seq(0,5))
Distribuzione di probabilità binomiale (sinistra) e probabilità binomiale cumulata (destra)

Figure 14.6: Distribuzione di probabilità binomiale (sinistra) e probabilità binomiale cumulata (destra)

Allo stesso modo possiamo immaginare di estrarre 20 insetti a caso da una popolazione in cui il rapporto tra i sessi è 1:1. Questo esperimento può essere simulato con:

Y <- rbinom(1, size = 20, prob = 0.5)
Y
## [1] 10

Assumendo che il ‘successo’ sia ottenere una femmina, il computer ci restituisce il numero delle femmine.

14.3 Capitolo 5: Esperimenti stime ed incertezza

14.3.1 E’ realistico l’intervallo di confidenza?

Abbiamo visto che un metodo semplice per costruire un intervallo di confidenza è utilizzare il doppio dell’errore standard. Questo intervallo, se viene utilizzato come misura di precisione/incertezza, è sempre accettabile. Tuttavia, da un punto di vista strettamente probabilistico, è lecito chiedersi: ma è proprio vero che se io ripeto l’esperimento molte volte e calcolo sempre l’intervallo di confidenza, riesco a centrare la media \(\mu\) nel 95% dei casi?

Proviamo a rispondere a questa domanda con una simulazione Monte Carlo. Prendiamo la nostra popolazione (\(\mu = 120\) e \(\sigma = 12\)) ed estraiamo centomila campioni. Per ogni campione calcoliamo l’intervallo di confidenza della media (P = 0.95) considerando il doppio dell’errore standard. Verifichiamo poi se questo intervallo contiene il valore 120: se si, assegniamo al campionamento il valore 1 (successo), altrimenti assegniamo il valore 0.

result <- rep(0, 100000)
set.seed(1234)
for (i in 1:100000){
  sample <- rnorm(3, 120, 12)
  limInf<- mean(sample) - sd(sample)/sqrt(3) * 2 
  limSup<- mean(sample) + sd(sample)/sqrt(3) * 2
  if (limInf<= 120 & limSup>= 120) result[i] = 1
}
sum(result)/100000
## [1] 0.81656

La simulazione mostra che la risposta alla domanda precedente è no: il nostro intervallo di confidenza non è riuscito a centrare la media nel 95% dei casi; ciò è avvenuto in poco più dell’80% dei casi. In realtà, possiamo facilmente verificare, con altre simulazioni di Monte Carlo, che la copertura effettiva dell’intervallo di confidenza si avvicina al 95% solo se abbiamo un numero di repliche superiori a 15-20 circa.

result <- rep(0, 100000)
set.seed(1234)
for (i in 1:100000){
  n <- 15
  sample <- rnorm(n, 120, 12)
  limInf<- mean(sample) - sd(sample)/sqrt(n) * 2 
  limSup<- mean(sample) + sd(sample)/sqrt(n) * 2
  if (limInf<= 120 & limSup>= 120) result[i] = 1
}
sum(result)/100000
## [1] 0.93591

Insomma, quando gli esperimenti sono piccoli, con poche repliche, dovremmo trovare un metodo di calcolo un po’ più affidabile, se veramente vogliamo ottenere un grado di copertura pari a quello nominale (P = 0.95).

Il problema nasce dal fatto che, nella statistica T che abbiamo introdotto nel capitolo 5:

\[T = \frac{m - \mu}{\sigma_m}\] \(\sigma_m\) viene sostituito con \(s_m\), cioè il valore di deviazione standard stimato nel campione. Come tutte le stime, anche \(s\) è ’soggetto ad incertezza, il che aggiunge un elemento ulteriore di imprecisione nella sampling distribution di T. Insomma ci chiediamo, la sampling distribution di T, calcolata con \(s\) invece che \(\sigma\) è ancora normale? Verifichiamo questo aspetto empiricamente, con una nuova simulazione Monte Carlo. Questa volta facciamo la seguente operazione:

  1. campioniamo tre individui
  2. Calcoliamo il valore di T con la statistica precedente, utilizzando la deviazione standard del campione e lo salviamo
  3. Con un po’ di pazienza, ripetiamo il tutto 100’000 volte.
#SIMULAZIONE MONTE CARLO - t di Student
set.seed(435)
result <- c()
for (i in 1:100000){
  sample3 <- rnorm(3, 120, 12)
  T <- (mean(sample3) - 120) / (sd(sample3)/sqrt(3))
  result[i] <- T
  }

Se riportiamo i valori ottenuti su una distribuzione di frequenze otteniamo il grafico sottostante.

#Plot sampling distribution
b <- seq(-600, 600, by=0.2)
hist(result, breaks = b, freq=F, xlab = expression(paste(m)), ylab="Density", xlim=c(-10,10), ylim=c(0,0.4), main="")
curve(dnorm(x, 0, 1), add=TRUE, col="blue")
curve(dt(x, 2), add=TRUE, col="red")

Vediamo che la sampling distribution di T calcolato utilizzando \(s\) invece che \(\sigma\) è solo approssimativamente normale. E’ facile vedere che questa approssimazione è sufficientemente buona solo se la numerosità del campione diviene abbastanza grande (es. \(n > 30)\), ma non certamente quando \(n\) = 3 (ve lo lascio per esercizio). In questo caso, la sampling distribution che osserviamo è più ‘dispersa’ di quella normale, con un maggior numero di valori sulle code.

Neyman scoprì che la sampling distribution di T poteva essere perfettamente descritta utilizzando la distribuzione t di Student, con un numero di gradi di libertà pari a quelli del campione (in questo caso 2), come vediamo nella figura sovrastante. In realtà questa conclusione era stata già raggiunta da William Sealy Gosset (1876 - 1937), uno statistico impiegato presso la fabbrica londinese della famosa birra Guinness, dove elaborava i dati relativi all’andamento del processo di maltazione. Egli, avendo definito questa nuova funzione di densità, per aggirare il divieto di pubblicazione imposto dal suo datore di lavoro, pubblicò i risultati sotto lo pseudonimo Student, da cui deriva il nome della distribuzione di densità.

Quindi, quando i campioni sono piccoli, il modo giusto di calcolare l’intervallo di confidenza è quello di utilizzare l’espressione seguente:

\[P \left( m + \textrm{qt}(0.025,n - 1) \cdot s_m \le \mu \le m + \textrm{qt}(0.975,n - 1) \cdot s_m \right) = 0.95\]

dove \(\textrm{qt}(0.025,n - 1)\) e \(\textrm{qt}(0.975,n - 1)\) sono rispettivamente il 2.5-esimo e il 97.5-esimo percentile della distribuzione t di Student, con n-1 gradi di libertà.

Nel capitolo 5 abbiamo utilizzato un esempio in cui abbiamo eseguito tre analisi chimiche da una soluzione erbicida di concentrazione pari a 120 mg/l, con uno strumento caratterizzato da un coefficiente di variabilità del 10%, che quindi, in assenza di errori sistematici, produce misure distribuite normalmente con media uguale a 120 e deviazione standard uguale a 12. Il campione osservato era

set.seed(1234)
Y <- rnorm(3, 120, 12)
Y
## [1] 125.1584 114.7349 105.6998

le statistiche descrittive sono:

m <- mean(Y)
s <- sd(Y)
m; s
## [1] 115.1977
## [1] 9.737554

I valori della distribuzione t di Student che lasciano al loro esterno il 5% delle varianti (2.5% per coda) sono:

qt(0.025, 2)
## [1] -4.302653
qt(0.975, 2)
## [1] 4.302653

Gli intervalli di confidenza sono pertanto:

m + qt(0.025, 2) * s/sqrt(3)
## [1] 91.00824
m + qt(0.975, 2) * s/sqrt(3)
## [1] 139.3871

E’ facile osservare che, se l’intervallo di confidenza è calcolato in questo modo, il suo coverage3 Operiamo con una simulazione Monte Carlo analoga a quella utilizzata nel capitolo 5.

result <- rep(0, 100000)
set.seed(1234)
for (i in 1:100000){
  sample <- rnorm(3, 120, 12)
  limInf<- mean(sample) + sd(sample)/sqrt(3) * qt(0.025, 2) 
  limSup<- mean(sample) + sd(sample)/sqrt(3) * qt(0.975, 2) 
  if (limInf<= 120 & limSup>= 120) result[i] = 1
}
sum(result)/100000
## [1] 0.94992

Ovviamente possiamo calcolare anche gli intervalli di confidenza al 99% di proababilità o qualunque altro intervallo di confidenza rilevante per il nostro studio.

14.3.2 Che cosa NON significa l’intervallo di confidenza?

Abbiamo già detto che l’intervallo di confidenza, calcolato su una serie di campionamenti ripetuti, contiene al suo interno la media vera e ignota della popolazione (\(\mu\)) con una probabilità pari a 0.95.

Tuttavia, la formula di Neyman si presta a cattive letture, che sono insensate da un punto di vista probabilistico, ma tuttavia molto frequenti nella pratica operativa. Ad esempio:

  1. NON E’ VERO CHE: c’è il 95% di probabilità che la media ‘vera’ della popolazione si trovi tra 91.0082383 e 139.3870891. La media vera della popolazione è sempre fissa e pari a 120 e non cambia affatto tra un campionamento e l’altro.
  2. NON E’ VERO CHE: ripetendo l’esperimento, il 95% delle stime che otteniamo cadono nell’intervallo 91.0082383 e 139.3870891. Una semplice simulazione mostra che quasi tutte le medie campionate cadono in quell’intervallo:
result <- rep(0, 100000)
set.seed(1234)
for (i in 1:100000){
  sample <- rnorm(3, 120, 12)
  if (mean(sample) <= 156.15 & mean(sample) >= 92.13) result[i] = 1
}
sum(result)/100000
## [1] 0.99996
  1. NON E’ VERO CHE: c’è il 95% di probabilità che l’affermazione ’la media vera è compresa tra 91.0082383 e 139.3870891 sia vera. Nelle normali condizioni sperimentali la media vera è ignota e non sapremo mai nulla su di essa: il nostro intervallo di confidenza può catturarla o no. Nel nostro esempio lo ha fatto, ed è tutto quello che possiamo dire.

Insomma, l’intervallo di confidenza vale per la sampling distribution e non vale per ogni singolo campionamento (esperimento). Pertanto, affermazioni del tipo: ”c’è il 95% di probabilità che \(\mu\) è compreso nell’intervallo di confidenza” oppure ”il valor più probabile di \(\mu\) è…” non sono corrette e anzi non hanno senso nella statistica tradizionale.

In altre parole, l’intervallo di confidenza è una sorta di polizza assicurativa che ci garantisce che, se operiamo continuativamente con le procedure indicate, al termine della nostra carriera avremo sbagliato in non più del 5% dei casi.

14.3.3 Popolazioni non gaussiane

Nel capitolo 5 abbiamo presentato un esempio in cui avevamo campionato da una distribuzione normale, riscontrando una sampling distribution per la media campionaria anch’essa normale. Ma che succede se la distribuzione di partenza è non-normale? La sampling distribution di uno stimatore è ancora normale? Vediamo un nuovo esempio.

Immaginiamo di avere 4’000’000 di semi ben mischiati (in modo che non ci siano raggruppamenti non casuali di qualche tipo), che costituiscono la nostra popolazione di partenza. Vogliamo appurare la frequenza relativa (p) dei semi dormienti. Questa informazione, nella realtà, esiste (\(\pi\) = 0.25), ma non è nota.

Dato l’elevato numero di ‘soggetti’, non possiamo testare la germinabilità di tutti i semi, ma dobbiamo necessariamente prelevare un campione casuale di 40 soggetti; ogni seme viene saggiato e, dato che la popolazione è molto numerosa, l’estrazione di un seme non modifica sensibilmente la proporzione di quelli dormienti nella popolazione (esperimenti indipendenti).

14.3.3.1 Il modello dei dati

Dopo aver descritto la popolazione e l’esperimento, ci chiediamo quale sia il modello matematico che genera i nostri dati (numero di successi su 40 semi estratti). Il disegno sperimentale ci assicura che ogni estrazione è totalmente indipendente dalla precedente e dalla successiva ed ha due soli risultati possibili, cioè successo (seme dormiente), o insuccesso (seme germinabile). Di conseguenza, ogni singola estrazione si configura come un esperimento Bernoulliano, con probabilità di successo pari a \(\pi\), il cui valore ‘vero’ esiste, è fisso, pre-determinato (esiste ancor prima di organizzare l’esperimento), anche se incognito e inconoscibile, a meno di non voler/poter esaminare tutti i semi disponibili. L’insieme delle 40 estrazioni (40 esperimenti Bernoulliani) può produrre un ventaglio di risultati possibili, da 40 successi a 40 insuccessi, per un totale di 41 possibili ‘outcomes’.

E’ evidente che i 41 possibili risultati non sono ugualmente probabili e si può dimostrare che la probabilità di ottenere k successi (con k che va da 0 ad n; n è al numero delle estrazioni) dipende da \(\pi\) ed è descrivibile matematicamente con la distribuzione binomiale \(\phi\):

\[\phi(k, n, p) = \frac{n!}{(n-k)!k!} p^k (1 - p)^{(n-k)}\]

Abbiamo quindi definito il modello matematico che descrive la probabilità di tutti i possibili risultati del nostro esperimento e quindi può in qualche modo essere considerato il ‘meccanismo’ che ‘genera’ i dati sperimentali osservati. Si tratta di un meccanismo puramente ‘stocastico’ nel quale è solo il caso che, attraverso il campionamento, determina il risultato dell’esperimento.

Con queste informazioni, possiamo simulare un esperimento con R, ottenendo i seguenti risultati:

set.seed(23456789)
rbinom(1, 40, 0.25)
## [1] 10

Abbiamo ottenuto 9 successi su 40, cioè 9 semi dormienti su 40 saggiati.

14.3.3.2 Stima dei parametri

Dovendo stimare la quantità \(\pi\), la statistica tradizionale trascura totalmente le nostre aspettative sul fenomeno e utilizza soltanto i risultati dell’esperimento. Chiamiamo p la quantità stimata e, dato che abbiamo contato nove semi dormienti, concludiamo che p = 0.225, in quanto questa, con le informazioni che abbiamo, è la cosa più verosimile. Anche in questo caso vi è chiara discrasia tra la verità ‘vera’ e l’osservazione sperimentale (tra \(\pi\) e \(p\)).

14.3.3.3 Sampling distribution

Cosa succede se ripetiamo l’esperimento? Come abbiamo imparato a fare, possiamo cercare una risposta attraverso la simulazione Monte Carlo, ricorrendo ad un generatore di numeri casuali da una distribuzione binomiale con n = 40 e \(\pi\) = 0.25 (in R si usa la funzione ‘rbinom(numeroDatiCasuali, n, p)’). Il codice è più semplice, in quanto non è necessario impostare un ciclo iterativo:

set.seed(1234)
result <- rbinom(10000000, 40, 0.25)

Esploriamo i risultati ottenuti:

result_p <- result/40
mean(result_p)
## [1] 0.2500129
sd(result_p)
## [1] 0.0684611

Osserviamo subito che, anche se i singoli esperimenti portano a stime diverse da \(\pi\), la media di \(p\) tende ad essere uguale a \(\pi\). L’errore standard (deviazione standard della sampling distribution) è 0.0685. Fino a qui, non vie è nulla di diverso dall’esempio precedente, se teniamo presente che la deviazione standard della popolazione originale (che è binomiale) è pari a \(\sqrt{p \times (1 - p)}\), quindi l’errore standard è \(\sqrt{0.25 \times 0.75 / 40} = 0.0685\).

Rimane da stabilire se la sampling distribution di \(p\) è normale. Possiamo utilizzare i 10’000’000 di valori ottenuti per costruire una distribuzione empirica di frequenze, come nel codice sottostante.

breaks <- seq(0, 0.7, by=0.025)
freqAss <- as.numeric( table(cut(result_p, breaks) ) ) 
freqRel <- freqAss/length(result_p)
density <- freqRel/0.025
p_oss <- breaks[2:length(breaks)]

plot(density ~ p_oss, type = "h",
     xlab = expression(paste(bar(p))),
     ylab="Density", 
    main="Sampling distribution per p", 
    xlim=c(0,0.6) )

curve(dnorm(x, 0.25, 0.0685), add=TRUE, col="red")

Vediamo che sampling distribution è approssimativamente normale con media pari a 0.25 e deviazione standard pari a 0.0685. Lo percepiamo chiaramente dal grafico soprastante, ma c’è una spiegazione scientifica per questo, basata sul TEOREMA DEL LIMITE CENTRALE:

  1. La sampling distribution di una statistica ottenuta da campioni casuali e indipendenti è approssimativamente normale, indipendentemente dalla distribuzione della popolazione da cui i campioni sono stati estratti.
  2. La media della sampling distribution è uguale al valore della statistica calcolata sulla popolazione originale, la deviazione standard della sampling distribution (errore standard) è pari alla deviazione standard della popolazione originale divisa per la radice quadrata della numerosità di un campione.

14.4 Capitolo 6. Introduzione al test d’ipotesi

14.4.1 Simulazione Monte Carlo di un test t di Student

La sampling distribution per T potrebbe essere ottenuta empiricamente, utilizzando una simulazione MONTE CARLO ed immaginando di estrarre numerose coppie di campioni, dalla stessa distribuzione normale, analogamente a quanto abbiamo fatto nell’esempio precedente. Se l’ipotesi nulla è vera, possiamo immaginare che questa distribuzione gaussiana abbia una media pari a (70.2 + 85.4)/2 = 77.8 e una deviazione standard pari alla deviazione standard delle dieci osservazioni (tutte insieme, senza distinzioni di trattamento), cioè 5.71.

Il codice da utilizzare in R per le simulazioni è il seguente:

A <- c(65, 68, 69, 71, 78)
P <- c(80, 81, 84, 88, 94)
media <- mean(c(A, P))
devSt <- sd(c(A, P))
set.seed(1234)
result <- rep(0, 100000)
for (i in 1:100000){
  sample1 <- rnorm(5, media, devSt)
  sample2 <- rnorm(5, media, devSt)
  SED <- sqrt( (sd(sample1)/sqrt(5))^2 +
                 (sd(sample2)/sqrt(5))^2 )
  result[i] <- (mean(sample1) - mean(sample2)) / SED
}

I risultati delle 100’000 simulazioni sono riportati nel grafico sottostante. Possiamo notare che, dei 100’000 valori di T osservati assumendo vera l’ipotesi nulla, solo l’un per mille sono superiori a quello da noi osservato e altrettanti sono inferiori a -4.5217. In totale, la probabilità di osservare un valore di T così alto in valore assoluto e dello 0.21 %.

SED_obs <- sqrt( (sd(A)/sqrt(5))^2 +
                   (sd(P)/sqrt(5))^2 )
T_obs <- (mean(A) - mean(P))/SED_obs
(length(result[result < T_obs]) + 
    length(result[result > - T_obs])) /100000
## [1] 0.00164
#Codice Grafico 
b <- seq(-12, 12, by=0.25)
hist(result, breaks = b, freq=F, xlab = expression(paste(m)), ylab="Density", xlim=c(-10,10), ylim=c(0,0.45), main="")

curve(dnorm(x), add=T, col="blue")

curve(dt(x, 8), add=T, col="red")
abline(v = 4.52, lty = 2)
abline(v = -4.52, lty = 2)
text(5, 0.4, label="4.52", adj=0, col = "red")
text(-5, 0.4, label="-4.52", adj=1, col = "red")

14.4.2 Tipologie alternative di test t di Student

Il test t può essere di tre tipi:

  1. Appaiato. In questo caso le misure sono prese a coppia sullo stesso soggetto e non sono quindi indipendenti.
  2. Omoscedastico. Le misure sono prese su soggetti diversi (indipendenti) e possiamo suppore che i due campioni provengano da due popolazioni con la stessa varianza.
  3. Eteroscedastico. Le misure sono prese su soggetti diversi, ma le varianze non sono omogenee.

Nel nostro esempio vediamo che le varianze dei campioni sono piuttosto simili e quindi adottiamo un test t omoscedastico (‘var.equal = T’).

Se dovessimo supporre che i due campioni provengono da popolazioni con varianze diverse, allora si porrebbe il problema di stabilire il numero di gradi di libertà del SEM. Abbiamo visto che se le varianze dei due campioni sono uguali (o meglio, sono due stime della stessa varianza), la varianza della somma/differenza ha un ha un numero di gradi di libertà pari alla somma dei gradi di libertà delle due varianze. Se le varianze fossero diverse, il numero di gradi di libertà della loro combinazione lineare (somma o differenza) si dovrebbe approssimare con la formula di Satterthwaite:

\[DF_s \simeq \frac{ \left( s^2_1 + s^2_2 \right)^2 }{ \frac{(s^2_1)^2}{DF_1} + \frac{(s^2_2)^2}{DF_2} }\]

Vediamo che se le varianze e i gradi di libertà sono uguali, la formula precedente riduce a:

\[DF_s = 2 \times DF\]

Nel nostro caso, se fosse \(s^2_1 \neq s^2_2\) avremmo un numero frazionario di gradi di libertà:

dfS <- (var(A) + var(P))^2 / 
  ((var(A)^2)/4 + (var(P)^2)/4)
dfS
## [1] 7.79772

Il risultato può essere riscontrato con:

t.test(A, P, var.equal=F)
## 
##  Two Sample t-test
## 
## data:  A and P
## t = -4.5217, df = 8, p-value = 0.001945
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -22.951742  -7.448258
## sample estimates:
## mean of x mean of y 
##      70.2      85.4

Se invece avessimo rilevato le misure accoppiate su quattro individui avremmo solo 4 gradi di libertà:

t.test(A, P, var.equal=T, paired=T)
## 
##  Paired t-test
## 
## data:  A and P
## t = -22.915, df = 4, p-value = 2.149e-05
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -17.04169 -13.35831
## sample estimates:
## mean of the differences 
##                   -15.2

14.4.3 Simulazione di un test di chi quadro

La simulazione di un test di \(\chi^2\) può esser fatta utilizzando la funzione ‘r2dtable()’ che produce il numero voluto di tabelle di contingenza, con righe e colonne indipendenti r rispettando i totali marginali voluti. Le tabelle prodotte (nel nostro caso 10’000) sono restituite come lista, quindi possiamo utilizzare la funzione ‘lapply()’ per applicare ad ogni elemento della lista la funzione che restituisce il \(\chi^2\) (‘chiSim’).

chiSim <- function(x) summary(as.table(x))$stat
set.seed(1234)
tabs <- r2dtable(10000, apply(tab, 1, sum), apply(tab, 2, sum))
chiVals <- as.numeric( lapply( tabs, chiSim) )
length(chiVals[chiVals > 9.768])
## [1] 435

Vediamo che vi sono 19 valori più alti di quello da noi osservato (p = 0.0019).

14.4.4 Errori di prima e di seconda specie

[da fare]

Bates, D. M., and D. G. Watts. 1988. Nonlinear Regression Analysis & Its Applications. Books: John Wiley & Sons, Inc.

Box, G. E. P., and D. R. Cox. 1964. “An Analysis of Transformations.” Journal of the Royal Statistical Society B-26: 211–52.

Carroll, R. J., and D. Ruppert. 1988. Transformation and Weighting in Regression. Books: Chapman and Hall.

Daniel, Johnnie. 2011. Sampling Essentials: Practical Guidelines for Making Sampling Choices. USA: SAGE.

de Mendiburu, Felipe. 2019. Agricolae: Statistical Procedures for Agricultural Research. https://CRAN.R-project.org/package=agricolae.

Draper, N. R., and H. Smith. 1998. Applied Regression Analysis. III. Books: John Wiley & Sons, Inc.

LeClerg, E. L., W. H. Leonard, and A. G. Clark. 1962. Field Plot Technique. Books: Burgess Publishing Company.

Pannacci, E., D. Pettorossi, and F. Tei. 2013. “Phytotoxic Effects of Aqueous Extracts of Sunflower on Seed Germination and Growth of Sinapis Alba L., Triticum Aestivum L. and Lolium Multiflorum Lam.” Allelopathy Journal 32 (1): 23.

Ratkowsky, David A. 1990. Handbook of Nonlinear Regression Models. Books: Marcel Dekker Inc.

Ritz, C., and J. C. Streibig. 2008. Nonlinear Regression with R. Books: Springer-Verlag New York Inc.

Ritz, Christian, Florent Baty, Jens C. Streibig, and Daniel Gerhard. 2015. “Dose-Response Analysis Using R.” Edited by Yinglin Xia. PLOS ONE 10 (12): e0146021. doi:10.1371/journal.pone.0146021.

References

de Mendiburu, Felipe. 2019. Agricolae: Statistical Procedures for Agricultural Research. https://CRAN.R-project.org/package=agricolae.


  1. Con il termine inglese coverage si intende, in un esperimento ripetuto un elevatissimo numero di volte, l’effettiva percentuale di campioni, per i quali l’intervallo di confidenza, calcolato per un certo P nominale (es. P = 0.95), contiene effettivamente la media \(\mu\) della popolazione.