Capitolo 13 Breve introduzione ai modelli misti

In questo libro abbiamo affrontato le questioni più rilevanti in relazione all’attività di ricerca e sperimentazione, almeno quelle che dovrebbero far parte del bagaglio culturale di ogni laureato in agraria o altre discipline biologiche. La trattazione non è affatto esaustiva e quindi, se si avesse l’intenzioni intraprendere una carriera professionale nel mondo dell’analisi dei dati, non si potrebbe prescindere dalla lettura di altri testi di approfondimento, che sono stati via via segnalati. Ci sembra comunque opportuno, in questo capitolo e quello successivo, dare almeno alcune informazioni di base su altri modelli un po’ più complessi, il cui impiego è abbastanza comune nella sperimentazione di pieno campo.

13.1 Raggruppamenti tra parcelle

Nel capitolo 2 abbiamo parlato dei disegni a split-plot o a strip-plot, caratterizzati da trattamenti che non sono allocati a caso, ma a gruppi di parcelle contigue.

Ad esempio, nella Figura 2.9 abbiamo mostrato come, in un esperimento a split-plot, i livello del fattore principale vengano allocati alle main-plots, ognuna delle quali è poi suddivisa in tante sub-plots, alle quali vengono allocati i livelli del trattamento secondario. Analogamente, in uno schema a strip-plot, il blocco viene diviso in righe e colonne, alle quali vengono rispettivamente allocati i livelli dei due fattori sperimentali.

Questi gruppi di parcelle trattate in modo non indipendente costituiscono dei veri e propri blocking factors, che si aggiungono agli altri fattori sperimentali di cui abbiamo parlato finora. Trascurare questi elementi di raggruppamento nel corso dell’analisi dei dati è sbagliato, in quanto si viene a rompere l’indipendenza dei residui. Ad esempio, in un disegno a split-plot, se una main-plot è più fertile di un’altra, tutte le misurazioni effettuate su di essa condivideranno lo stesso effetto positivo e saranno quindi più simili tra di loro che le misurazioni effettuate su main-plots diverse. Insomma, con uno schema a split-plot, i dati ottenuti sulla stessa main-plot non sono indipendenti, ma sono in qualche modo ‘apparentati’; più propriamente, si parla di correlazione intra-classe, un concetto simile, ma non totalmente coincidente con quello di correlazione di Pearson. Se l’effetto prodotto dalle main-plots o dagli altri blocking factor non viene incluso nel modello statistico descrittivo, esso rimane tra gli effetti non spiegati e, di conseguenza, viene incluso nei residui, che, pertanto, saranno correlati e quindi non indipendenti tra di loro.

In generale, dobbiamo tener presente che, ogni volta in cui l’allocazione dei trattamenti comporta una qualche forma di raggruppamento tra parcelle, l’elemento di raggruppamento, o meglio, l’effetto da esso provocato, deve essere introdotto nel modello, in modo che nei residui rimangano solo effetti puramente stocastici ed ignoti.

13.2 Esperimenti a split-plot

Consideriamo nuovamente l’esempio presentato nel capitolo 12, nel quale avevamo due fattori sperimentali a confronto: il diserbo, con due livelli (totale o localizzato sulla fila) e la lavorazione, con tre livelli (aratura profonda, aratura superficiale e minimum tillage). In precedenza abbiamo analizzato i dati come se l’esperimento fosse stato disegnato a blocchi randomizzati, secondo la mappa di campo riportata nel Capitolo 2, in Figura 2.8. In realtà, questo esperimento era stato disegnato seguendo uno schema a split-plot, come indicato in Figura 2.9; vediamo quindi come procedere all’analisi dei dati in modo corretto.

13.2.1 Definizione del modello lineare

Il modello ANOVA per un disegno a split-plot è simile a quello dell’ANOVA fattoriale, fatta salva la presenza di un effetto aggiuntivo, cioè quello delle parcelle principali, che costituiscono un elemento di raggruppamento delle osservazioni:

\[Y_{ijk} = \mu + \gamma_k + \alpha_i + \theta_{ik} + \beta_j + \alpha\beta_{ij} + \varepsilon_{ijk}\]

dove \(\gamma\) è l’effetto del blocco \(k\), \(\alpha\) è l’effetto della lavorazione \(i\), \(\beta\) è l’effetto del diserbo \(j\), \(\alpha\beta\) è l’effetto dell’interazione per la specifica combinazione della lavorazione \(i\) e del diserbo \(j\). Abbiamo incluso anche \(\theta\) che è l’effetto della main-plot; dato che ogni parcella principale è identificata univocamente dal blocco \(k\) a cui appartiene e dalla lavorazione \(i\) in essa eseguita, abbiamo utilizzato i pedici corrispondenti.

Concentriamoci per un attimo sulle main-plots: esiste una differenza sostanziale tra questo effetto e gli altri effetti in gioco. In particolare, mentre i livelli delle lavorazioni e del diserbo li abbiamo attentamente prescelti, perché eravamo specificatamente interessati ad essi, nel caso delle main-plots non abbiamo nessun interesse specifico, abbiamo operato una selezione casuale da un universo più grosso, quello di tutte le main-plots possibili nel nostro appezzamento. In altre parole, possiamo dire che le main-plots potrebbero essere cambiate a piacimento senza snaturare le finalità del nostro esperimento, cosa che non possiamo altrettanto dire dei livelli delle lavorazioni e del diserbo. In termini statistici diciamo che lavorazioni e diserbo sono effetti fissi mentre le main-plots sono un effetto random, che si aggiunge ad \(\varepsilon\).

Abbiamo quindi definito un modello che, a differenza di tutti quelli che abbiamo considerato finora, ha due effetti random, entrambi assunti come gaussiani, con media zero e deviazioni standard rispettivamente pari a \(\sigma_{\theta}\) e \(\sigma\). I modelli che hanno più di un effetto random si chiamano modelli misti e, di conseguenza, i dati provenienti da esperimenti a split-plot e, vedremo tra poco, a strip-plot debbono essere analizzati utilizzando un modello misto.

La piattaforma dei modelli misti presenta molte differenze concettuali rispetto alla piattaforma dei modelli lineari, che abbiamo utilizzato finora. Data la complessità, l’argomento richiede un corso di livello più avanzato e non può essere introdotto qui, anche se ci sembra opportuno mostrare almeno come trattare i dati provenienti da questi esperimenti così comuni in ambito agrario.

13.2.2 Model fitting con R

Come abbiamo visto, i risultati di questo esperimento sono disponibili nel file ‘beet.csv’, che può essere aperto direttamente dalla solita repository online, utilizzando il codice sottostante. Subito dopo aver caricato il file, trasformiamo le variabili indipendenti in ‘factors’, come al solito.

fileName <- "https://www.casaonofri.it/_datasets/beet.csv"
dataset <- read.csv(fileName, header=T)
dataset$Tillage <- as.factor(dataset$Tillage)
dataset$WeedControl <- as.factor(dataset$WeedControl)
dataset$Block <- as.factor(dataset$Block)
head(dataset)
##   Tillage WeedControl Block  Yield
## 1     MIN         TOT     1 11.614
## 2     MIN         TOT     2  9.283
## 3     MIN         TOT     3  7.019
## 4     MIN         TOT     4  8.015
## 5     MIN        PART     1  5.117
## 6     MIN        PART     2  4.306

Prima di tutto, dobbiamo creare una variabile che identifichi in modo univoco le main-plots. Per questo potremmo utilizzare una codifica numerica, oppure, più semplicemente, costruire un nuovo fattore combinando i livelli del blocco e della lavorazione, dato che ogni main-plot è identificata univocamente dal blocco di cui fa parte e dalla lavorazione che ha subito.

dataset$mainPlot <- with(dataset, factor(Block:Tillage))

In secondo luogo, il model fitting non può essere effettuato con la funzione lm(), che ci permette di definire un solo effetto random. Dobbiamo invece utilizzare una delle funzioni di R per i modelli misti, tra le quali proponiamo lmer(), che richiede due packages aggiuntivi: ‘lme4’ e ‘lmerTest’, che dobbiamo installare, se non lo abbiamo già fatto, e caricare in memoria.

install.packages("lme4") #only at first time
install.packages("lmerTest") #only at first time

La sintassi di lmer() è simile a quella di lm() e gli effetti random possono essere inseriti tra parentesi, utilizzando l’operatore ‘1|’, come mostrato nel box sottostante.

library(lme4)
library(lmerTest)
mod.split <- lmer(Yield ~ Block + Tillage * WeedControl +
                  (1|mainPlot), data=dataset)

Per quanto riguarda le assunzioni di base, i modelli misti sono sempre e comunque modelli linear e quindi fanno le stesse assunzioni di moscedasticità e normalità dei residui. L’ispezione grafica di questi utlimi può essere condotta nel modo usuale, se non che la funzione plot() restituisce solo il grafico dei residui verso i valori attesi (quello corrispondente all’argomento ‘which = 1’ in lm()). D’altro canto, non c’è un modo veloce per ottenere un QQ-plot e, pertanto, per la valutazione della normalità dei residui, possiamo utilizzare il test di Shapiro-Wilks, già indicato nel capitolo 8.

Analisi grafica dei residui per un modello misto

Figura 13.1: Analisi grafica dei residui per un modello misto 

shapiro.test(residuals(mod.split))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mod.split)
## W = 0.93838, p-value = 0.1501

Dopo esserci accertati della validità delle assunzioni di base, possiamo procedere con la scomposizione della varianza, utilizzando l’usuale funzione anova(), che fornisce però un output leggermente diverso dal solito. Come secondo argomento è necessario indicare il metodo prescelto per la stima dei gradi di libertà, che non è così semplice come nel caso dei modelli lineari generali. Consigliamo il metodo ‘Kenward-Roger’ che, pur fornendo comunque un’approssimazione, ha dimostrato ottima affidabilità.

anova(mod.split, ddf="Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
##                      Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)  
## Block                3.6596  1.2199     3     6  0.6521 0.61016  
## Tillage             23.6565 11.8282     2     6  6.3228 0.03332 *
## WeedControl          3.3205  3.3205     1     9  1.7750 0.21552  
## Tillage:WeedControl 19.4641  9.7321     2     9  5.2023 0.03152 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Da qui in avanti l’analisi procede con il metodo usuale, già spiegato nel capitolo 12, utilizzando la funzione emmeans() per il calcolo delle medie marginali e per i confronti multipli. Questa parte di lavoro viene lasciata per esercizio.

13.3 Esperimenti a strip-plot

In alcune circostanze, soprattutto nelle prove di diserbo chimico, potrebbe trovare applicazione un altro tipo di schema sperimentale, nel quale ogni blocco viene diviso in tante righe quanti sono i livelli di un fattore sperimentale e tante colonne quanti sono i livelli dell’altro. In questo modo, il primo trattamento sperimentale viene applicato a tutte le parcelle di una riga e l’altro trattamento a tutte le parcelle di una colonna. Ovviamente, l’allocazione alle righe e alle colonne è casuale e cambia in ogni blocco.

Questo disegno è detto strip-plot ed è molto comodo perché consente di lavorare velocemente. Nel capitolo 2 abbiamo già presentato un esempio di come l’esperimento di cui stiamo parlando (interazione tra lavorazioni e diserbo chimico) avrebbe potuto essere pianificato, con uno schema a strip-plot (Figura 2.10).

13.3.1 Definizione del modello lineare

Dalla mappa si può osservare che, entro ogni blocco, le osservazioni sono raggruppate per riga e per colonna. Di conseguenza, anche le righe e le colonne debbono essere incluse nel modello lineare, per ripristinare l’indipendenza dei residui. Il modello è quindi:

\[Y_{ijk} = \mu + \gamma_k + \alpha_i + \theta_{ik} + \beta_j + \zeta_{jk} + \alpha\beta_{ij} + \varepsilon_{ijk}\]

Analogamente a quanto abbiamo visto per lo split-plot, le 12 righe sono univocamente definite da un livello di lavorazione e un blocco e, per questo, l’effetto \(\theta\) porta i pedici \(i\) e \(k\). D’altra parte, le 8 colonne sono definite da un livello di diserbo e un blocco, per cui l’effetto \(\zeta\) porta i pedici \(j\) e \(k\). Entrambi gli effetti sono da considerare random, per cui anche questo è un modello ‘misto’, con tre effetti random.

13.3.2 Model fitting con R

Il codice da utilizzare è analogo a quello precedentemente esposto per il disegno a split-plot e prevede l’introduzione degli effetti random ‘riga’ e ‘colonna’, che vengono inseriti tra parentesi e con l’operatore ‘1|’, come indicato nel box sottostante.

dataset$column <- with(dataset, factor(Block:Tillage))
dataset$row <- with(dataset, factor(Block:WeedControl))
mod.strip <- lmer(Yield ~ Block + Tillage * WeedControl +
                  (1|column) + (1|row), data=dataset)
anova(mod.strip, ddf = "Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
##                      Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
## Block                1.6323  0.5441     3 2.5022  0.3631 0.78795  
## Tillage             23.6565 11.8282     2 6.0000  7.8934 0.02089 *
## WeedControl          1.4810  1.4810     1 3.0000  0.9883 0.39343  
## Tillage:WeedControl 19.4641  9.7321     2 6.0000  6.4946 0.03155 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anche in questo caso, il calcolo delle medie e i confronti multipli vengono effettuati con la tecnica usuale, indicata nei capitoli precedenti.

13.4 Altre letture

  1. Bates, D., Mächler, M., Bolker, B., Walker, S., 2015. Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software 67. https://doi.org/10.18637/jss.v067.i01
  2. Gałecki, A., Burzykowski, T., 2013. Linear mixed-effects models using R: a step-by-step approach. Springer, Berlin.
  3. Kuznetsova, A., Brockhoff, P.B., Christensen, H.B., 2017. lmerTest Package: Tests in Linear Mixed Effects Models. Journal of Statistical Software 82, 1–26.