Capitolo 9 Combinazioni lineari/nonlineari dei parametri

Dopo aver dimostrato che un modello fornisca una buona descrizione dei dati sperimentali e che non vi siano violazioni degli assunti di base, è possibile testare diverse ipotesi rilevanti che siano ‘traducibili’ in combinazioni lineari o non lineari dei parametri del modello, cioè delle funzioni che coinvolgono i parametri ed un insieme di coefficienti; per comprendere questo concetto e la sua importanza, iniziamo con un esempio.

9.1 Esempio 9.1: Determinare l’effetto medio degli erbicidi

Torniamo al nostro esperimento in cui abbiamo confrontato due erbicidi e la loro miscela con il testimone non trattato, in un disegno sperimentale completamente randomizzato con quattro repliche. Nel capitolo precedente, abbiamo visto che non vi sono problemi con nessuna delle assunzioni di base per il fitting dei modelli lineari, mentre l’analisi della varianza ed il relativo test di F hanno mostrato come esistano differenze significative tra alcuni dei trattamenti. Nel riquadro sottostante riportimo nuovamente il codice R per l’adattamento del modello.

# Example 9.1
# Determining the average effect of herbicides in the 'mixture' data

# Loading the packages
library(statforbiology)
library(multcomp)

# Loading the data
dataset <- getAgroData("mixture")

# Fitting a one-way ANOVA model
model <- lm(Weight ~ Treat, data=dataset)
summary(model)
## 
## Call:
## lm(formula = Weight ~ Treat, data = dataset)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -6.360 -2.469  0.380  2.567  6.025 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            9.175      1.959   4.684 0.000529 ***
## TreatMixture_378      -4.047      2.770  -1.461 0.169679    
## TreatRimsulfuron_30    7.685      2.770   2.774 0.016832 *  
## TreatUnweeded         17.598      2.770   6.352 3.65e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.918 on 12 degrees of freedom
## Multiple R-squared:  0.8554, Adjusted R-squared:  0.8193 
## F-statistic: 23.66 on 3 and 12 DF,  p-value: 2.509e-05
# Variance partitioning
anova(model)
## Analysis of Variance Table
## 
## Response: Weight
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## Treat      3 1089.53  363.18  23.663 2.509e-05 ***
## Residuals 12  184.18   15.35                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

A questo punto, è rilevante chiedersi: “Il controllo delle piante infestanti è stato (in media) una pratica efficace?”. Questa domanda ‘biologica’ può essere tradotta in una combinazione lineare dei parametri stimati, cioè: \(\mu\) = 9.175, \(\alpha_2\) = -4,0475, \(\alpha_3\) = 7.685 e \(\alpha_4\) = 17.5975, dove \(\mu\) è la media del primo erbicida (in ordine alfabetico), \(\alpha_2\), \(\alpha_3\) e \(\alpha_4\) sono, rispettivamente, le differenze tra le medie del secondo, terzo e quarto trattamento e la media del primo (vedi Capitolo 4). Con queste stime, possiamo calcolare la media dei vasi non trattati come \(\mu_u = \mu + \alpha_4\), e la media dei vasi trattati come \(\mu_t = [\mu + (\mu + \alpha_2) + (\mu + \alpha_3)]/3\), ovvero \(\mu_t = \mu + 1/3 \, \alpha_2 + 1/3 \, \alpha_3\). Di conseguenza, la differenza tra vasi trattati e non trattati è uguale alla seguente combinazione lineare: \(\delta_{tu} = \mu_t - \mu_{ut} = 1/3 \, \alpha_2 + 1/3 \, \alpha_3 - \alpha_4\). Il risultato di questa combinazione (\(\delta_{tu}\)) risponde alla domanda precedente: se \(\delta_{tu} < 0\) possiamo concludere che il trattamento erbicida è, in media, efficace nel ridurre la biomassa di Solanum nigrum.

La combinazione lineare sopra definita può essere espressa come somma dei prodotti di un insieme di parametri per un insieme di coefficienti:

\[C = 0 \cdot \mu + \frac{1}{3} \cdot \alpha_2 + \frac{1}{3} \cdot \alpha_3 - 1 \cdot \alpha_4\]

In questa espressione, la matrice dei coefficienti è:

\[K = \left[0 \,\,\, 1/3 \,\,\, 1/3 \,\, -1 \right]\]

Vedremo che in R esistono alcune funzioni che ricevono in input l’elenco dei parametri stimati e la matrice dei coefficienti e restituiscono il valore della combinazione (in questo caso: \(C = \delta_{tu} = -16,385\)), insieme al suo errore standard ottenuto tramite la legge di propagazione degli errori, considerando le varianze e le covarianze dei parametri del modello. Maggiori dettagli sulla propagazione degli errori sono reperibili nel blog associato a questo libro4.

Una delle funzioni più importanti è flessibili è glht(), contenuta nel package multcomp, il cui impiego è mostrato nel riquadro sottostante, insieme a quello della funzione confint() che può essere utilizzata per ottenere l’intervallo di confidenza.

# Example 9.1 [continuation]

# Fitting linear combinations with the 'glht()' function in R
library(multcomp)

# Setting the matrix of coefficients
K <- matrix(c(0, 1/3, 1/3, -1),
            nrow = 1, ncol = 4, byrow = T)
row.names(K) <- c("delta_tu")
K
##          [,1]      [,2]      [,3] [,4]
## delta_tu    0 0.3333333 0.3333333   -1
# Fitting the contrast
con <- glht(model, linfct = K)
summary(con)
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = Weight ~ Treat, data = dataset)
## 
## Linear Hypotheses:
##               Estimate Std. Error t value Pr(>|t|)    
## delta_tu == 0  -16.385      2.262  -7.244 1.02e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
# Confidence intervals 
confint(con)
## 
##   Simultaneous Confidence Intervals
## 
## Fit: lm(formula = Weight ~ Treat, data = dataset)
## 
## Quantile = 2.1788
## 95% family-wise confidence level
##  
## 
## Linear Hypotheses:
##               Estimate lwr      upr     
## delta_tu == 0 -16.3850 -21.3132 -11.4568

È importante notare che la stima \(\delta_{tu}\) può essere divisa per il suo errore standard per ottenere un test t per l’ipotesi nulla che la stima osservata non sia significativamente diversa da zero (vedere Sezione 7.1). Il valore P corrispondente è molto basso, consentendo il rifiuto dell’ipotesi nulla. Pertanto, concludiamo che l’effetto medio dell’erbicida è statisticamente significativo.

9.2 Esempio 9.2: Determinare l’efficiacia erbicida media

I concetti precedenti possono essere estesi a qualsiasi altra funzione dei parametri del modello, anche se non può essere ricondotta alla somma dei prodotti dei coefficienti per i parametri (come nell’Equazione riportata nel precedente paragrafo) oppure se è intrinsecamente non lineare (cioè contiene logaritmi o altre funzioni che non siano semplici somme di parametri). L’unica differenza è che, con le combinazioni non lineari, gli errori standard non possono essere calcolati con il metodo di propagazione degli errori, ma debbono essere approssimati utilizzando il cosiddetto metodo delta (Weisberg, 2005)5.

In R, possiamo utilizzare la funzione gnlht() nel package statforbiology, che generalizza l’approccio della glht(). Come esempio, riconsideriamo il dataset “mixture” ed utilizziamo i parametri del modello per calcolare la riduzione percentuale della crescita delle infestanti rispetto al controllo non trattato (efficacia erbicida; HE). Per il primo erbicida, è:

\[\textrm{HE} = \frac{\mu_4 - \mu_1}{\mu_4} \times 100 = \frac{\mu_1 + \alpha_4 - \mu_1}{\mu_1 + \alpha_4} \times 100 = \frac{\alpha_4}{\mu_1 + \alpha_4} \times 100\]

In questa espressione, i parametri \(\mu_1\) e \(\alpha_4\) sono al denominatore e, pertanto, questa combinazione è non lineare. Le espressioni per ottenere gli HE per tutti gli altri erbicidi possono essere ottenute con semplici calcoli algebrici, come riportato nel succesivo riquadro.

Per utilizzare la funzione gnlht() bisogna codificare la lista delle funzioni da calcolare e passarla come argomento. Inoltre, il legame tra il nome dei parametri utilizzati in tali combinazioni e la posizione delle relative stime come risulta nell’output della funzione coef() viene creato con l’argomento parameterNames. Le efficacie dei tre erbicidi sono rispettivamente del 65%, 80% e 37% e sono tutte significativamente diverse da 0 (riquadro sottostante)

# Esempio 9.2
# Determinazione dell'efficacia media degli erbicidi nei dati della miscela

# Caricamento dei pacchetti e dei dati
library(statforbiology)
dataset <- getAgroData("mixture")

# Adattamento di un modello ANOVA ad una via
model <- lm(Weight ~ Treat, data=dataset)

# Adattamento di combinazioni non lineari con 'gnlht()'
library(statforbiology)
funList <- list(~ a4/(mu1 + a4) * 100,
                ~ (a4 - a2)/(mu1 + a4) * 100, 
                ~ (a4 - a3)/(mu1 + a4) * 100)
gnlht(model, funList, parameterNames = c("mu1", "a2", "a3", "a4"))
##                         Form Estimate       SE   t-value      p-value
## 1        a4/(mu1 + a4) * 100 65.72976 7.734313  8.498461 2.014359e-06
## 2 (a4 - a2)/(mu1 + a4) * 100 80.84788 7.449568 10.852694 1.469226e-07
## 3 (a4 - a3)/(mu1 + a4) * 100 37.02493 8.646543  4.282050 1.065143e-03

9.3 Combinazioni di interesse generale

Le funzioni lineari e non lineari dei parametri del modello sono molto utili per rispondere a specifiche domande di ricerca, anche se, soprattutto quando sono molte, la codifica ‘manuale’ della matrice dei coefficienti può rivelarsi un compito molto noioso. Fortunatamente, R ci ci mette a disposizione diverse matrici pre-definite per le combinazioni lineari o non lineari di maggior interesse generale, cioè:

  1. Medie Marginali attese (EMMs)
  2. Contrasti ortogonali (e non)
  3. Confronti a coppie
  4. Retrotrasformazioni
  5. Previsioni e previsioni inverse

Nelle situazioni in cui abbiamo uno o più predittori nominali, il nostro interesse principale, di solito, è quello di stimare le medie della popolazione per ogni livello di trattamento (ad esempio, \(\mu_1\), \(\mu_2\), \(\mu_3\), \(\ldots{}\) e \(\mu_n\), dove \(n\) è il numero di livelli di trattamento). Sebbene le medie aritmetiche dei gruppi di trattamento possano essere utilizzate come stimatori delle medie delle relative popolazioni (vedere la funzione tapply() nel Capitolo 3), è necessario tener conto che si tratta di stimatori che, in alcune circostanze, possono essere distorti, ad esempio quando i dati sono sbilanciati (numero diverso di repliche per ciascun gruppo di trattamento). Le medie marginali attese (Expected Marginal Means; EMMs) rappresentano un’importante alternativa: sono uguali alle medie aritmetiche quando i dati sono bilanciati, mentre forniscono stimatori migliori delle medie della popolazione quando i dati sono sbilanciati. Le EMMs possono essere facilmente ottenute come combinazioni lineari dei parametri del modello; ad esempio, tornando al dataset “mixture”, nell’esempio 9.1, dove i parametri stimati sono \(\mu = 9.175\), \(\alpha_2 = -4.0475\), \(\alpha_3 = 7.685\) e \(\alpha_4 = 17.5975\) ed hanno il significato biologico illustrato in precedenza (media del primo gruppo e differenze tra le medie degli altri gruppi rispetto al primo), le quattro EMMs possono essere ottenuti dalle seguenti combinazioni lineari:

\[\left\{ {\begin{array}{l} C_1 = \mu_1 = 1 \cdot \mu + 0 \cdot \alpha_2 + 0 \cdot\alpha_3 + 0 \cdot\alpha_4 = 9,18\\ C_2 = \mu_2 = 1 \cdot \mu + 1 \cdot \alpha_2 + 0 \cdot\alpha_3 + 0 \cdot\alpha_4 = 5,13\\ C_3 = \mu_3 = 1 \cdot \mu + 0 \cdot \alpha_2 + 1 \cdot\alpha_3 + 0 \cdot\alpha_4 = 16:86\\ C_4 = \mu_4 = 1 \cdot \mu + 0 \cdot \alpha_2 + 0 \cdot\alpha_3 + 1 \cdot\alpha_4 = 26,77\\ \end{array}} \right.\]

La matrice dei coefficienti è:

\[ K = \left[ {\begin{array}{llll} 1 & 0 & 0 & 0\\ 1 & 1 & 0 & 0\\ 1 & 0 & 1 & 0\\ 1 & 0 & 0 & 1 \end{array}} \right]\]

In linea di principio, per calcolare queste combinazioni si può utilizzare la funzione glht() nel package multcomp, come mostrato nel riquadro sottostante.

# Esempio 9.3
# Medie dei minimi quadrati per i dati 'mixture'

# Caricamento dei pacchetti
library(statforbiology)
library(multcomp)

# Caricamento dei dati
dataset <- getAgroData("mixture")

# Adattamento di un modello ANOVA ad una via
model <- lm(Weight ~ Treat, data=dataset)

# Stima degli EMM con la funzione 'glht()'

# Impostazione della matrice dei coefficienti
K <- matrix(c(1,0,0,0,1,1,0,0,1,0,1,0,1,0,0,1),
nrow = 4, ncol = 4, byrow = T)
row.names(K) <- c("mu1", "mu2", "mu3", "mu4")
K
##     [,1] [,2] [,3] [,4]
## mu1    1    0    0    0
## mu2    1    1    0    0
## mu3    1    0    1    0
## mu4    1    0    0    1
# Calcolo della combinazione lineare
con <- glht(model, linfct = K)
summary(con)
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = Weight ~ Treat, data = dataset)
## 
## Linear Hypotheses:
##          Estimate Std. Error t value Pr(>|t|)    
## mu1 == 0    9.175      1.959   4.684  0.00199 ** 
## mu2 == 0    5.127      1.959   2.618  0.08143 .  
## mu3 == 0   16.860      1.959   8.607  < 0.001 ***
## mu4 == 0   26.773      1.959  13.668  < 0.001 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
# Intervalli di confidenza
confint(con)
## 
##   Simultaneous Confidence Intervals
## 
## Fit: lm(formula = Weight ~ Treat, data = dataset)
## 
## Quantile = 2.8936
## 95% family-wise confidence level
##  
## 
## Linear Hypotheses:
##          Estimate lwr     upr    
## mu1 == 0  9.1750   3.5069 14.8431
## mu2 == 0  5.1275  -0.5406 10.7956
## mu3 == 0 16.8600  11.1919 22.5281
## mu4 == 0 26.7725  21.1044 32.4406

Tuttavia, soprattutto quando il numero di medie è molto elevato, codificare la matrice dei coefficienti da zero può essere un compito arduo; pertanto, per il calcolo degli EMMs, si preferisce la funzione emmeans() nel pacchetto emmeans (Lenth, 2016), perché è molto più semplice da usare. Richiede due argomenti di base: l’oggetto modello e il predittore per il quale si desidera calcolare le EMMs (riquadro sottostante).

# Esempio 9.3 [continuazione]

# EMMs con la funzione 'emmeans()'

library(emmeans) # Carica il pacchetto
muj <- emmeans(model, ~Treat)
muj
##  Treat           emmean   SE df lower.CL upper.CL
##  Metribuzin__348   9.18 1.96 12     4.91     13.4
##  Mixture_378       5.13 1.96 12     0.86      9.4
##  Rimsulfuron_30   16.86 1.96 12    12.59     21.1
##  Unweeded         26.77 1.96 12    22.50     31.0
## 
## Confidence level used: 0.95

Oltre a essere stime affidabili delle medie della popolazione, le EMMs possono anche servire come punto di partenza per costruire ulteriori combinazioni lineari o non lineari, che saranno esaminate in seguito.

9.4 Example 9.4: Contrasti ortogonali (e non)

Nell’Esempio 9.1, abbiamo visto che una combinazione lineare dei parametri del modello può aiutare a rispondere ad una domanda ‘biologica’ rilevante, basata sulla differenza tra la media dei vasi non trattati e la media dei vasi trattati. Considerando sempre il dataset “mixture”, è possibile pianificare tre combinazioni che reppresentano differenze tra medie, cioè:

  1. trattato vs. non trattato (in media)
  2. miscela vs. erbicidi singoli (in media)
  3. metribuzin vs. rimsulfuron

Invece che combinare i parametri del modello come nell’Esempio 9.1, le tre combinazioni di cui sopra possono anche essere codificate come combinazioni tra le medie. Se consideriamo l’ordine con cui compaiono le quattro medie nel dataset muj in un riquadro precedente, in Code box 9.5, le tre combinazioni lineari possono essere scritte come:

\[\left\{ {\begin{array}{l} C_1 = \frac{1}{3} \, \mu_1 + \frac{1}{3} \, \mu_2 + \frac{1}{3} \,\mu_3 - \mu_4\\ C_2 = \frac{1}{2} \, \mu_1 - \mu_2 + \frac{1}{2} \,\mu_3 - 0 \cdot \mu_4\\ C_3 = \mu_1 - 0 \cdot \mu_2 - \mu_3 - 0 \cdot \mu_4 \end{array}} \right.\]

Le tre combinazioni lineari C1, C2 e C3 hanno la tipica proprietà che, per ognuna, la somma dei coefficienti è zero e, per questo, sono chiamate contrasti. Ne consegue che il risultato di un contrasto deve essere uguale a 0, quando la differenza tra medie o gruppi di medie che esso sottende è uguale a zero, di modo che sia possibile utilizzare un test di t per saggiare questa ipotesi nulla.

Inoltre, i tre contrasti sopra riportati hanno l’ulteriore caratteristica che i loro coefficienti sono reciprocamente non correlati e, pertanto, sono chiamati contrasti ortogonali. Normalmente, in un dataset esistono tanti contrasti ortogonali quante sono le medie meno una.

L’ortogonalità era considerata una proprietà desiderabile perché la somma delle devianze associate ad ogni contrasto corrisponde alla devianza del trattamento (vedere il riquadro sottostante). Pertanto, l’applicazione dei contrasti ortogonali era considerata un metodo elegante ed efficiente per sfruttare l’intera informazione contenuta (Chew, 1976). Oggigiorno, l’enfasi sull’ortogonalità è minore, ma può essere utile fornire un esempio di come i contrasti ortogonali possano essere calcolati con R, utilizzando la funzione contrast() nel package emmeans.

I tre passaggi tipici sono:

  1. preparare una lista con i coefficienti, assegnando un nome a ciascun contrasto (si noti che questa è una differenza rispetto alla funzione glht(), che richiede una matrice di coefficienti, non una lista);
  2. utilizzare la funzione contrast() nel package emmeans, passando i seguenti argomenti: (i) l’oggetto risultante dalla chiamata alla funzione emmeans(); (ii) la lista dei contrasti \(K\), come argomento method; (iii) il tipo di aggiustamento (ne parleremo più avanti, per ora utilizzate il comando così come mostrato nel riquadro sottostante)
# Example 9.4
# Orthogonal contrasts with the 'emmeans()' function

# Loading the packages
library(statforbiology)
library(emmeans)
library(multcomp)

# Loading the data
dataset <- getAgroData("mixture")

# Fitting a one-way ANOVA model
model <- lm(Weight ~ Treat, data=dataset)

# Creating a vector of coefficients for each contrast
C1 <- c(1/3, 1/3, 1/3, -1)
C2 <- c(1/2, -1, 1/2, 0)
C3 <- c(1, 0, -1, 0)

# Creating a list of contrasts (NOT a matrix)
K <- list("U vs T" = C1, 
          "mix vs single" = C2,
          "met vs rim" = C3)

# Fitting the contrasts, using the 'muj' object
muj <- emmeans(model, ~Treat)
con <- contrast(muj, method = K, adjust="none")
con
##  contrast      estimate   SE df t.ratio p.value
##  U vs T          -16.39 2.26 12  -7.244  <.0001
##  mix vs single     7.89 2.40 12   3.289  0.0065
##  met vs rim       -7.68 2.77 12  -2.774  0.0168
# Confidence intervals
confint(con)
##  contrast      estimate   SE df lower.CL upper.CL
##  U vs T          -16.39 2.26 12   -21.31   -11.46
##  mix vs single     7.89 2.40 12     2.66    13.12
##  met vs rim       -7.68 2.77 12   -13.72    -1.65
## 
## Confidence level used: 0.95

I contrasti ortogonali possono essere calcolati anche utilizzando la funzione aov(). In questo caso, l’elenco dei contrasti deve essere trasformato in una matrice e associato al fattore sperimentale tramite la funzione C() (notare la lettera maiuscola; vedere il riquadro sottostante). Utilizzando questa codifica, è possibile vedere che la somma dei quadrati per l’effetto del trattamento è completamente suddivisa in tre parti, una per ciascun contrasto, producendo valori P identici a quelli mostrati nel riquadro precedente.

# Example 9.4 [continuation]

# Fitting orthogonal contrasts and partitioning
# the sum of squares for the treatment effect

# Creating a vector of coefficients for each contrast (as in Code box 9.6)
C1 <- c(1/3, 1/3, 1/3, -1)
C2 <- c(1/2, -1, 1/2, 0)
C3 <- c(1, 0, -1, 0)

# Creating a list of contrasts (as in Code box 9.6)
K <- list("U vs T" = C1, 
          "mix vs single" = C2,
          "met vs rim" = C3)

# Transform the list of contrasts into a matrix
Kmat <- do.call(rbind, K)

# Associate the matrix of coefficients to the treatment factor
dataset$Treat2 <- C(factor(dataset$Treat), t(Kmat))

# Fit the model with the 'aov()' function
mod <- aov(Weight ~ Treat2, dataset)
summary(mod, split = list(Treat2 = list(1,2,3)))
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## Treat2        3 1089.5   363.2  23.663 2.51e-05 ***
##   Treat2: C1  1  805.4   805.4  52.476 1.02e-05 ***
##   Treat2: C2  1  166.0   166.0  10.816  0.00647 ** 
##   Treat2: C3  1  118.1   118.1   7.696  0.01683 *  
## Residuals    12  184.2    15.3                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

9.5 Esempio 9.5: I confronti multipli a coppie

Non sempre le prove sperimentali consentono di saggiare pochi contrasti pre-stabiliti, ma spesso è necessario confrontare tutte le possibili coppie di trattamenti (pairwise comparisons). In questo caso dovremmo definire un contrasto per ogni coppia di medie, anche se l’impiego del package ‘emmeans’ agevola, non di poco, il lavoro.

In particolare, possiamo immaginare due situazioni di riferimento: tutti contro tutti (confronti tipo “Tukey”) e tutti verso uno (confronti tipo “Dunnett”)(Hsu, 1996). Questo secondo tipo di contrasto può essere interessante, quando sia importante confrantare tutti i trattamenti verso un riferimento prescelto, ad esempio la miscela tra i due erbicidi. Con il metodo di Tukey, il numero di contrasti è molto elevato, ed è pari a \(n \times (n - 1) / 2\), dove \(n\) è il numero di medie da confrontare. Con il metodo di Dunnett, il numero di contrasti si riduce a \(n - 1\), e tale riduzione potrebbe essere importante per ragioni che verranno descritte più avanti.

Bisogna tener presente che l’adozione di confronti a coppie non è particolarmente elegante perché implica molti test non pianificati e privi di un obiettivo chiaro (una sorta di “spray and pray”; Cousens, 1988). Inoltre, secondo lo scienziato inglese John Tukey: “tutto ciò che sappiamo del mondo ci insegna che gli effetti di A e B sono sempre diversi, purchè noi consideriamo un numero di decimali sufficientemente elevato. Pertanto, chiedersi se gli effetti siano diversi è sciocco”. Ciononostante, i confronti a coppie possono essere molto utili; sono di moda da quasi un secolo e quindi vale la pena conoscerli.

Per il dataset “mixture”, ci sono quattro medie da confrontare e, quindi, un totale di sei confronti a coppie che possono essere definiti come segue:

\[\left\{ {\begin{array}{*{20}{c}} C_1 = \mu_1 - \mu_2 = 1 \cdot \mu_1 - 1 \cdot \mu_2 + 0 \cdot \mu_3 + 0 \cdot \mu_4 \\ C_2 = \mu_1 - \mu_3 = 1 \cdot \mu_1 - 0 \cdot \mu_2 - 1 \cdot \mu_3 + 0 \cdot \mu_4 \\ C_3 = \mu_1 - \mu_4 = 1 \cdot \mu_1 - 0 \cdot \mu_2 + 0 \cdot \mu_3 - 1 \cdot \mu_4 \\ C_4 = \mu_2 - \mu_3 = 0 \cdot \mu_1 + 1 \cdot \mu_2 - 1 \cdot \mu_3 + 0 \cdot \mu_4 \\ C_5 = \mu_2 - \mu_4 = 0 \cdot \mu_1 + 1 \cdot \mu_2 + 0 \cdot \mu_3 - 1 \cdot \mu_4 \\ C_6 = \mu_3 - \mu_4 = 0 \cdot \mu_1 + 0 \cdot \mu_2 + 1 \cdot \mu_3 - 1 \cdot \mu_4 \end{array}} \right.\]

La matrice dei coefficienti è:

\[K = \left[ {\begin{array}{rrrr} 1 & - 1 & 0 & 0\\ 1 & 0 & - 1 & 0\\ 1 & 0 & 0 & - 1\\ 0 & 1 & - 1 & 0\\ 0 & 1 & 0 & - 1\\ 0 & 0 & 1 & - 1 \end{array}} \right]\]

Fortunatamente, non è necessario codificare la matrice \(K\) da zero; nel quadro sottostante mostriamo come si possa ottenere il confronto multiplo con la funzione contrast() (come sopra) e passando il valore ‘tukey’ o ‘dunnett’ all’argomento ‘method’. In questo secondo caso, R confronta tutte le tesi con metribuzin, che è il primo livello in ordine alfabetico, mentre noi avremmo preferito confrontare tutte le tesi con la miscela. Per ottenere questo risultato basta aggiungere l’argomento ‘ref’ ed assegnare il valore ‘2’, considerando che la miscela è la seconda tesi in ordine alfabetico:

# Example 9.5
# Pairwise comparisons for the 'mixture' data

# Loading the packages
library(statforbiology)
library(multcomp)

# Loading the data
dataset <- getAgroData("mixture")

# Fitting a one-way ANOVA model
model <- lm(Weight ~ Treat, data=dataset)

# Determining the least square means
muj <- emmeans(model, ~Treat)

# All-pairwise comparisons (Tukey's method)
contrast(muj, adjust="none", method="tukey")
##  contrast                         estimate   SE df t.ratio p.value
##  Metribuzin__348 - Mixture_378        4.05 2.77 12   1.461  0.1697
##  Metribuzin__348 - Rimsulfuron_30    -7.68 2.77 12  -2.774  0.0168
##  Metribuzin__348 - Unweeded         -17.60 2.77 12  -6.352  <.0001
##  Mixture_378 - Rimsulfuron_30       -11.73 2.77 12  -4.235  0.0012
##  Mixture_378 - Unweeded             -21.64 2.77 12  -7.813  <.0001
##  Rimsulfuron_30 - Unweeded           -9.91 2.77 12  -3.578  0.0038
## Comparisons all-against-one (Dunnett's method)
con <- contrast(muj, adjust="none", method="dunnett", ref = 2)
con
##  contrast                      estimate   SE df t.ratio p.value
##  Metribuzin__348 - Mixture_378     4.05 2.77 12   1.461  0.1697
##  Rimsulfuron_30 - Mixture_378     11.73 2.77 12   4.235  0.0012
##  Unweeded - Mixture_378           21.64 2.77 12   7.813  <.0001
## Confidence intervals
confint(con)
##  contrast                      estimate   SE df lower.CL upper.CL
##  Metribuzin__348 - Mixture_378     4.05 2.77 12    -1.99     10.1
##  Rimsulfuron_30 - Mixture_378     11.73 2.77 12     5.70     17.8
##  Unweeded - Mixture_378           21.64 2.77 12    15.61     27.7
## 
## Confidence level used: 0.95

Il risultato delle elaborazioni sovrastanti mostra i contrasti con il loro errore standard e potrebbe essere interessante calcolare anche l’intervallo di confidenza per le differenze stimate. Ciò può esser fatto facilmente assegnando il risultato della funzione contrast() ad una variabile ed esplorando quest’ultima con il metodo confint().

9.5.1 Display a lettere

I risultati di un confronto multiplo a coppie possono essere presentati anche con un display a lettere, nel quale le medie seguite da lettere diverse sono significativamente diverse per un livello di probabilità di errore minore di quello dato.

Realizzare un display a lettere manualmente è piuttosto facile, utilizzando la seguente procedura:

  1. ordinare le medie in senso crescente/decrescente
  2. partire dalla prima media e aggiungere la lettera A a tutte quelle che non sono significativamente diverse
  3. passare alla seconda media e aggiungere la lettera B a tutte quelle che non sono significativamente diverse
  4. procedere analogamente con tutte le medie successive, finche non vengono più individuate differenze significative.

Con R si può sfruttare il package ‘emmeans’, utilizzando i comandi sottostanti.

multcomp::cld(muj, adjust="none", Letters=LETTERS)
##  Treat           emmean   SE df lower.CL upper.CL .group
##  Mixture_378       5.13 1.96 12     0.86      9.4  A    
##  Metribuzin__348   9.18 1.96 12     4.91     13.4  A    
##  Rimsulfuron_30   16.86 1.96 12    12.59     21.1   B   
##  Unweeded         26.77 1.96 12    22.50     31.0    C  
## 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping symbol,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.

9.5.2 Correzione per la molteplicità

Operando nel modo anzidetto, ogni contrasto/confronto ha una probabilità di errore del 5% (\(\alpha_C = 0.05\)). Se i contrasti/confronti sono più di uno (‘famiglia’ di n contrasti), la probabilità di sbagliarne almeno uno (maximum experimentwise error rate: \(\alpha_E\)) è data da:

\[\alpha_E = 1 - (1 - \alpha_C)^n\]

Bisogna premettere che l’anzidetta formula vale quando i contrasti sono totalmente indipendenti tra loro, cosa che quasi mai avviene, dato che, anche in un semplice modello ANOVA, i contrasti condividono la stessa varianza d’errore e sono quindi più o meno correlati tra di loro. Con contrasti non indipendenti la formula anzidetta fornisce una sovrastima di \(\alpha_E\) (per questo si parla di maximum experimentwise error rate).

Il numero di confronti a coppie per esperimento può essere anche molto elevato: se ho k medie il numero dei confronti possibili è pari a \(k \cdot (k-1)/2\). Di conseguenza, la probabilità di errore per esperimento (\(\alpha_E\)) può essere molto più alta del valore \(\alpha_C\) prefissato per confronto.

Ad esempio, se ho 15 medie, ho \((15 \cdot 14)/2 = 105\) confronti possibili. Se uso \(\alpha_C = 0.05\) per ogni confronto, la probabilità di sbagliarne almeno uno è pari (in caso di confronti indipendenti) a \(1 - (1 - 0.05)^105 = 0.995\). Sostanzialmente, vi è pressoché la certezza che in questo esperimento qualcosa sia sbagliato!

Per questo motivo, quando si elaborano i dati di un esperimento nel quale è necessario fare molti contrasti, o confronti, o, più in generale, molti test d’ipotesi simultanei, si potrebbe voler esprimere un giudizio globale (simultaneo) sull’intera famiglia di contrasti/confronti, minimizzando la possibilità che anche solo uno o pochi di essi siano sbagliati. In particolare, ciò potrebbe capitare quando:

  1. vogliamo trovare i migliori di k trattamenti, senza correre rischi di escluderne erroneamente qualcuno. In questa situazione, facendo ogni confronto con il 5% di probabilità di errore, la probabilità di escludere erroneamente anche solo un trattamento dal lotto dei migliori è molto più alta di quella prefissata, perché basta sbagliare anche uno solo dei k - 1 confronti con il migliore.
  2. Abbiamo utilizzato un display a lettere e intendiamo affermare che ‘i trattamenti seguiti da lettere diverse sono significativamente diversi’. In questo caso, stiamo tirando una conclusione basata sull’intera famiglia di confronti e non possiamo lecitamente riportare la probabilità di errore di un singolo confronto.

In tutte le condizioni analoghe a quelle più sopra accennate si pone il problema di applicare un aggiustamento per la molteplicità, in modo da rispettare un certo livello prestabilito di protezione per esperimento (e non per confronto). La prima possibilità che abbiamo è quella di aggiustare il P-level per ogni confronto, in modo da diminuire la probabilità di errore per l’intera famiglia di sei confronti. Utilizzando la formula che lega la probabilità d’errore per esperimento (\(\alpha_E\)) alla probabilità d’errore per confronto (\(\alpha_C\); vedi sopra), possiamo prendere la sesta colonna del dataframe ‘con’, quella che contiene i P-levels non corretti, e trasformarla come segue:

alphaC <- as.data.frame(con)[,6]
1 - (1 - alphaC)^6
## [1] 6.722991e-01 6.923077e-03 2.869757e-05

Più facilmente, possiamo arrivare allo stesso risultato con il package ‘emmeans’:

contrast(muj, method = "tukey", adjust = "sidak")
##  contrast                         estimate   SE df t.ratio p.value
##  Metribuzin__348 - Mixture_378        4.05 2.77 12   1.461  0.6723
##  Metribuzin__348 - Rimsulfuron_30    -7.68 2.77 12  -2.774  0.0968
##  Metribuzin__348 - Unweeded         -17.60 2.77 12  -6.352  0.0002
##  Mixture_378 - Rimsulfuron_30       -11.73 2.77 12  -4.235  0.0069
##  Mixture_378 - Unweeded             -21.64 2.77 12  -7.813  <.0001
##  Rimsulfuron_30 - Unweeded           -9.91 2.77 12  -3.578  0.0226
## 
## P value adjustment: sidak method for 6 tests

Vediamo che il secondo confronto, che era significativo, non è più tale adottando la correzione di Sidak.

Un’alternativa più nota (e semplice) è quella di utilizzare la diseguaglianza di Bonferroni:

\[\alpha_E = \alpha_C \cdot k\]

Quest’ultima è un po’ più conservativa della precedente, nel senso che fornisce un P-level aggiustato leggermente più alto dell’altra.

alphaC * 6
## [1] 1.018071e+00 6.943132e-03 2.869792e-05

Oppure possiamo utilizzare la funzione contrast():

contrast(muj, method = "pairwise", adjust = "bonferroni")
##  contrast                         estimate   SE df t.ratio p.value
##  Metribuzin__348 - Mixture_378        4.05 2.77 12   1.461  1.0000
##  Metribuzin__348 - Rimsulfuron_30    -7.68 2.77 12  -2.774  0.1010
##  Metribuzin__348 - Unweeded         -17.60 2.77 12  -6.352  0.0002
##  Mixture_378 - Rimsulfuron_30       -11.73 2.77 12  -4.235  0.0069
##  Mixture_378 - Unweeded             -21.64 2.77 12  -7.813  <.0001
##  Rimsulfuron_30 - Unweeded           -9.91 2.77 12  -3.578  0.0228
## 
## P value adjustment: bonferroni method for 6 tests

Esistono altre procedure di aggiustamento del P-level (metodi di Holm, Hochberg, Hommel), ma nessuna di queste tiene conto della correlazione eventualmente esistente tra i contrasti e tutte quindi sono da definirsi più o meno ‘conservative’.

Oltre che aggiustare il P-level, possiamo anche utilizzare altre procedure di aggiustamento, basate sulla distribuzione multivariata di t e in grado di tener conto dell’eventuale correlazione esistente tra i contrasti stessi. Una di queste procedure viene impiegata di default nella funzione contrast() nel package ‘emmeans’:

#Confronti multipli a coppie, basati sul t multivariato
contrast(muj, method="tukey")
##  contrast                         estimate   SE df t.ratio p.value
##  Metribuzin__348 - Mixture_378        4.05 2.77 12   1.461  0.4885
##  Metribuzin__348 - Rimsulfuron_30    -7.68 2.77 12  -2.774  0.0698
##  Metribuzin__348 - Unweeded         -17.60 2.77 12  -6.352  0.0002
##  Mixture_378 - Rimsulfuron_30       -11.73 2.77 12  -4.235  0.0055
##  Mixture_378 - Unweeded             -21.64 2.77 12  -7.813  <.0001
##  Rimsulfuron_30 - Unweeded           -9.91 2.77 12  -3.578  0.0173
## 
## P value adjustment: tukey method for comparing a family of 4 estimates

Possiamo notare che i P-levels sono leggermente più bassi di quelli ottenuti con Bonferroni, che conferma quindi di essere una procedura molto conservativa, mentre l’impiego del t multivariato consente di rispettare esattamente il tasso di errore ‘per esperimento’.

Ovviamente la correzione per la molteplicità ed il conseguente innalzamento del P-level sono fortemente dipendenti dal numero di contrasti effettuati; se utilizziamo un set di confronti tipo ‘dunnett’ invece che tipo ‘pairwise’, avremo meno confronti e quindi una correzione più ‘leggera’.

#Confronti multipli a coppie, basati sul t multivariato
contrast(muj, method="dunnett")
##  contrast                         estimate   SE df t.ratio p.value
##  Mixture_378 - Metribuzin__348       -4.05 2.77 12  -1.461  0.3711
##  Rimsulfuron_30 - Metribuzin__348     7.68 2.77 12   2.774  0.0442
##  Unweeded - Metribuzin__348          17.60 2.77 12   6.352  0.0001
## 
## P value adjustment: dunnettx method for 3 tests

9.5.3 Procedure ‘classiche’ di confronto multiplo

In questo libro, i confronti a coppie sono inseriti nel contesto dei contrasti lineari, come, ad esempio, in Hsu (1996). In passato, i confronti multipli si basavano sulla determinazione di una differenza critica, in base alla quale due medie venivano considerate significativamente diverse quando la loro differenza superava la differenza critica. Nella letteratura scientifica si trovano molte procedure di questo tipo, tra le quali segnaliamo:

  1. Minima Differenza Significativa (MDS o LSD)
  2. Honest Significant Difference di Tukey
  3. Test di Dunnett
  4. Test di Duncan
  5. Test di Newman-Keuls
  6. Test di confronto multiplo di Tukey

Il primo metodo corrisponde esattamente a quello che abbiamo utilizzato all’inizio, per fare confronti multipli ‘tutti contro tutti’, senza correzione per la molteplicità. Il secondo e il terzo metodo corrispondono rispettivamente al test di confronto ‘tutti verso tutti’ e ‘uno verso tutti’ indicati in precedenza, con correzione per la molteplicità.

Non è necessario dettagliare gli altri test, in quanto, seppur siano ancora molto utilizzati, vengono ormai ritenuti obsoleti e sconsigliabili, da parecchi ricercatori. Chi vuole, trova altre informazioni nella letteratura indicata in fondo al capitolo; mostriamo solo come la HSD di Tukey HSD possa essere calcolata ed utilizzata per in confronto multiplo, utilizzando la funzione aov() per la stima del modello e la funzione TukeyHSD() per il confronto multiplo con correzione per la molteplicità (Code Box 9.13).

# Esempio 9.5 [continuazione]
# Tukey HSD per un modello lineare adattato con la funzione 'aov()' 
mod2 <- aov(Weight ~ Treat, data = dataset)
TukeyHSD(mod2)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Weight ~ Treat, data = dataset)
## 
## $Treat
##                                   diff        lwr       upr     p adj
## Mixture_378-Metribuzin__348    -4.0475 -12.271978  4.176978 0.4884620
## Rimsulfuron_30-Metribuzin__348  7.6850  -0.539478 15.909478 0.0698118
## Unweeded-Metribuzin__348       17.5975   9.373022 25.821978 0.0001854
## Rimsulfuron_30-Mixture_378     11.7325   3.508022 19.956978 0.0055015
## Unweeded-Mixture_378           21.6450  13.420522 29.869478 0.0000247
## Unweeded-Rimsulfuron_30         9.9125   1.688022 18.136978 0.0172572

Per concludere la sezione sui confronti multipli vogliamo dare alcuni consigli pratici di buon senso.

  1. Quando è possibile, pianificare gli esperimenti in modo da ottenere le risposte cercate con pochi contrasti di interesse. In questo modo il problema della molteplicità è minimizzato.
  2. Non usare mai contrasti con serie di dati quantitative. In questo caso la regressione è l’approccio corretto e ne parleremo in un prossimo capitolo. In generale, utilizzare i contrasti solo se sono coerenti con la logica dell’esperimento.
  3. Evitare di utilizzare i confronti multipli nei disegno fattoriale per confrontare tutte le combinazioni dei fattori sperimentali e preferire invece il confronto dei lvelli di un fattore all’interno dei livelli dell’altri.
  4. Pianificare esattamente il numero di contrasti necessari ed eseguirli, fornendo il valore del contrasto e il suo errore standard.
  5. Decidere è necessario aggiustare il p-level (e gli intervalli di confidenza) per la molteplicità (tasso di errore comparisonwise o experimentwise).
  6. Se si decide di aggiustare il p-level, considerare che le procedure di Bonferroni o Sidak possono essere eccessivamente protette. Preferire quindi le procedure di aggiustamento basate sulla distribuzione t multivariata, il che, a livello di confronto multiplo con dati bilanciati, è equivalente ad utilizzate la Tukey HSD o il test di Dunnett.
  7. Evitare le procedure di Duncan e Newmann-Keuls: non danno il livello di protezione cercato e, inoltre, non sono basate su una differenza critica costante (quindi sono difficili da discutere).

9.6 Esempio 9.6: confronti multipli con dati trasformati

Un altro aspetto di cui tener conto è relativo al fatto che, nel caso in cui si fosse reso necessario qualche tipo di trasformazione stabilizzante, il test di confronto multiplo deve opportunamente adeguato. Per capire questo aspetto, possiamo riprendere in mano il dataset ‘insects.csv’ già esplorato nel capitolo precedente e relativo alle quindici piante trattate con tre insetticidi (cinque piante per insetticida), su ciascuna delle quali, alcune settimane dopo il trattamento, sono stati contati gli insetti presenti e vitali. Ricarichiamo il dataset e, tenendo conto dell’esigenza emerse nel capitolo precedente, analizziamo i dati previa trasformazione logaritmica. Quest’ultima viene eseguita direttamente nel corso del ‘model fitting’ e non nel dataset, per un motivo che sarà chiaro in seguito.

fileName <- "https://www.casaonofri.it/_datasets/insects.csv"
dataset <- read.csv(fileName, header = T)
dataset$Insecticide <- factor(dataset$Insecticide)
mod <- lm(log(Count) ~ Insecticide, data = dataset)

A questo punto, possiamo calcolare le medie marginali attese, che, tuttavia, sono su scala logaritmica, come indicato nel codice sottostante.

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

Presentare le medie su scala logaritmica, in molti casi, potrebbe non essere di immediata o facile lettura. Per questo, potremmo pensare di retro-trasformare le medie, utilizzando la trasformazione inversa di quella logaritmica, cioè l’antilogaritmo. Ad esempio, per la prima media:

\[e^{6.34} = 566.7963\]

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

mean(dataset[dataset$Insecticide == "T1","Count"])
## [1] 589.8

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

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

retroMedie <- emmeans(mod, ~Insecticide, type = "response")
retroMedie
##  Insecticide response     SE df lower.CL upper.CL
##  T1             568.6 101.01 12    386.1    837.3
##  T2             335.1  59.54 12    227.6    493.5
##  T3              51.9   9.22 12     35.2     76.4
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale

Con questo abbiamo tutto quello che ci serve: stime ed errori standard, che, ovviamente, sono diversi per le diverse medie retro-trasformate, coerentemente con la mancanza di omoscedasticità. Il test di confronto multiplo è, pertanto:

library(multcomp)
cld(retroMedie, Letters = LETTERS)
##  Insecticide response     SE df lower.CL upper.CL .group
##  T3              51.9   9.22 12     35.2     76.4  A    
##  T2             335.1  59.54 12    227.6    493.5   B   
##  T1             568.6 101.01 12    386.1    837.3   B   
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## P value adjustment: tukey method for comparing a family of 3 estimates 
## Tests are performed on the log scale 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping symbol,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.

La retrotrasformazione è immediata in quanto la funzione emmeans() è stata in grado di riconoscere la trasformazione stabilizzante applicata all’interno del ‘model fitting’. La situazione diviene più complessa se la trasformazione è applicata prima del ‘model fitting’ oppure se non viene automaticamente riconosciuta. Ad esempio, una trasformazione inversa non viene automaticamente riconosciuta e, di conseguenza, la retrotrasformazione non viene effettuata.

mod2 <- lm(1/Count ~ Insecticide, data = dataset)
emmeans(mod2, ~Insecticide, type = "response")
##  Insecticide  emmean      SE df lower.CL upper.CL
##  T1          0.00182 0.00279 12 -0.00426  0.00789
##  T2          0.00317 0.00279 12 -0.00291  0.00924
##  T3          0.02120 0.00279 12  0.01513  0.02727
## 
## Confidence level used: 0.95

In questo caso, dobbiamo eseguire la trasformazione prima del ‘model fitting’ e poi cambiare la griglia di riferimento per il modello (update(ref_grid(mod2), ....)) specificando la trasformazione effettuata (tran = "inverse"). La griglia così modificata viene passata al posto del modello alla funzione emmeans(), come indicato di seguito.

dataset$InvCount <- 1/dataset$Count
mod3 <- lm(InvCount ~ Insecticide, data = dataset)
updGrid <- update(ref_grid(mod3), tran = "inverse")
emmeans(updGrid, ~Insecticide, type = "response")
##  Insecticide response    SE df lower.CL upper.CL
##  T1             550.9 845.9 12    126.8      Inf
##  T2             315.8 278.0 12    108.2      Inf
##  T3              47.2   6.2 12     36.7       66
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the inverse scale

Questo metodo si può utilizzare con numerose funzioni, come, ad esempio “identity”, “1/mu^2”, “inverse”, “reciprocal”, “log10”, “log2”, “asin.sqrt” e “asinh.sqrt”.

Un approccio analogo può essere utilizzato se vogliamo effettuare la trasformazione di Box-Cox nella sua forma semplificata (\(W = Y^{\lambda}\)), utilizzando un valore \(\lambda\) diverso da quelli ‘semplici’ (cioè diverso da 1, 0, 0.5 o -1). Ad esempio, se vogliamo utilizzare \(\lambda = 0.25\), dobbiamo operare in modo simile a quello descritto in precedenza, utilizzando la funzione make.tran(), come evidenziato nel box seguente.

dataset$lCount <- dataset$Count^0.25
mod4 <- lm(lCount ~ Insecticide, data = dataset)
updGrid <- update(ref_grid(mod4), tran = make.tran("power", 0.25))
emmeans(updGrid, ~Insecticide, type = "response")
##  Insecticide response   SE df lower.CL upper.CL
##  T1             573.5 78.7 12    420.3    765.3
##  T2             340.5 53.2 12    238.5    472.1
##  T3              53.1 13.2 12     29.6     88.3
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the mu^0.25 scale

9.7 Esempio 9.7: Previsioni e previsioni inverse

Per i modelli di regressione lineare o non lineare, lavorare con le medie è solitamente irrilevante, mentre può essere necessario utilizzare il modello per fare previsioni, ovvero per calcolare la risposta attesa con un certo valore della variabile indipendente oppure per calcolare quale valore di una variabile indipendente consente di ottenere la risposta voluta (previsione inversa). Anche le previsioni possono essere considerate combinazioni lineari o non lineari dei parametri del modello e possono essere calcolate con R, utilizzando una delle numerose funzioni disponibili, come predict(), emmeans(), glht() e gnlht().

A titolo di esempio, possiamo considerare i dati ‘Ammi94’, relativi a un esperimento per definire la relazione tra la densità di Ammi majus (un’importante infestante di diverse colture primaverili nell’Italia centrale) espressa come numero di piante per metro quadrato e la resa del girasole in tonnellate per ettaro. Le diverse densità delle infestanti sono state ottenute con semina mauale, mentre i dati di resa rappresentano la media di quattro repliche (Onofri et al., 1994). Il dataset è disponibile nel pacchetto statforbiology e nel riquadro successivo, si vede come venga adattato un semplice modello di regressione lineare e le due stime (intercetta e pendenza, ovvero \(\beta_0\) e \(\beta_1\)) vengono recuperate e visualizzate.

# Example 9.7
# Predictions and inverse predictions for linear regression
# Weed density/Yield relationship for Ammi majus

# Loading the packages
library(statforbiology)
library(multcomp)
library(emmeans)

# Loading the data
dataset <- getAgroData("Ammi94")
dataset
##   Density  Yield
## 1       0 2.4866
## 2      16 2.1806
## 3      21 1.8456
## 4      36 1.8412
## 5      49 1.6486
# Fitting a regression model
mod <- lm(Yield ~ Density, data = dataset)
summary(mod)
## 
## Call:
## lm(formula = Yield ~ Density, data = dataset)
## 
## Residuals:
##        1        2        3        4        5 
##  0.08403  0.04167 -0.21094  0.03182  0.05342 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.402566   0.108697  22.103 0.000203 ***
## Density     -0.016477   0.003667  -4.494 0.020567 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.138 on 3 degrees of freedom
## Multiple R-squared:  0.8707, Adjusted R-squared:  0.8275 
## F-statistic: 20.19 on 1 and 3 DF,  p-value: 0.02057

In un esperimento del genere, vi sono due interessi fondamentali:

  1. Prevedere la resa con misurazioni precoci della densità delle infestanti (previsione);
  2. Determinare la densità delle infestanti che provoca un certo livello di perdita di resa (previsione inversa).

Queste stime possono essere ottenute come combinazioni lineari o non lineari. Ad esempio, immaginiamo di voler determinare il livello di resa previsto con 10 piante di Ammi majus per metro quadrato. Ci sono tre possibilità:

  1. Utilizzare la funzione glht() per adattare un contrasto \(C = 1 \times \beta_0 + 10 \times \beta_1\), dove i parametri \(\beta_0\) e \(\beta_1\) sono, rispettivamente, l’intercetta e la pendenza del modello lineare;
  2. Utilizzare la funzione predict() e passare un dataframe (non un vettore!) contenente il livello di densità per la previsione (argomento ‘newdata’);
  3. Utilizzare la funzione emmeans(), passando il livello di densità come dataframe all’argomento ‘at’.

Tutte e tre le possibilità sono illustrate nel riquadro sottostante; i risultati sono, ovviamente, del tutto equivalenti.

# Example 9.7 [continuation]

# Three methods for predictions

# 1.Fitting with 'glht()'##################
# Setting a contrast matrix
K <- matrix(c(1, 10), nrow = 1, ncol = 2)
summary(glht(mod, linfct = K))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = Yield ~ Density, data = dataset)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value Pr(>|t|)    
## 1 == 0  2.23779    0.08123   27.55 0.000105 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
# 2.Fitting with 'predict()' ################
predict(mod, newdata = data.frame(Density= c(10)), se.fit = T)
## $fit
##        1 
## 2.237793 
## 
## $se.fit
## [1] 0.08123191
## 
## $df
## [1] 3
## 
## $residual.scale
## [1] 0.1380352
# 3.Fitting with 'emmeans()' ####################
emmeans(mod, ~Density, at = data.frame(Density= c(10)))
##  Density emmean     SE df lower.CL upper.CL
##       10   2.24 0.0812  3     1.98      2.5
## 
## Confidence level used: 0.95

Immaginiamo ora di voler stimare il livello di densità delle malerbe che provoca una riduzione della resa del 10%. Considerando che la resa senza malerbe è di 2.403 tonnellate per ettaro (l’intercetta stimata più sopra), una riduzione della resa del 10% corrisponde all’ottenimento di un livello di resa pari a \(2.403 - 0.1 \times 2.403 = 2.163\), che è pari a \((2,163 - \beta_0)/\beta_1\). Quest’ultima equazione è intrinsecamente non lineare; pertanto, si deve utilizzare la funzione gnlht(), come mostrato qui sotto.

# Example 9.7 [continuation]

# Inverse prediction with gnlht()

gnlht(mod, func = list(~(21.53 - b0)/b1), 
      parameterNames = c("b0", "b1"))
##              Form  Estimate       SE  t-value   p-value
## 1 (21.53 - b0)/b1 -1160.835 263.7767 4.400825 0.0217523

9.8 Altre letture

  1. Bretz F, Hothorn T, Westfall P (2011) Multiple Comparisons Using R. CRC Press, Boca Raton, FL
  2. Chew V (1976) Comparing treatment means: a compendium. Hortscience 11(4):348–357
  3. Cousens R (1988) Misinterpretetion of results in weed research through inappropriate use of statistics. Weed Res 28:281–289
  4. Hsu JC (1996) Multiple Comparisons: Theory and Methods. Chapman & Hall, London
  5. Lenth RV (2016) Least-squares means: the R package lsmeans. J Stat Softw 69(1). https://doi.org/10.18637/jss.v069.i01.
  6. Onofri A, Tei F (1994) Competitive ability and threshold levels for three broadleaf weed species in sunflower. Weed Res 34:471–479
  7. Onofri A, Carbonell EA, Piepho HP, Mortimer AM, Cousens RD (2010) Current statistical issues in Weed Research. Weed Res 50(1):5–24
  8. Piepho HP (2004) An algorithm for a letter-based representation of all-pairwise comparisons. J Comput Graph Stat 13(2):456–466
  9. Weisberg S (2005) Applied Linear Regression, 3rd edn. Wiley, Hoboken

9.9 Esercizi

  1. Considerare gli esercizi da 5 a 7 nel Capitolo 4. Per tutti questi esercizi, dopo essersi accertati se è necessaria una trasformazione per soddisfare le ipotesi di base per i modelli lineari, eseguire il test di confronto multiplo, adottando una correzione appropriata per la molteplicità.

  1. vedi: ↩︎

  2. Chi fosse interessato può trovare informazioni sul nel blog associato, all’indirizzo \url{https://www.statforbiology.com/2024/stat_general_thedeltamethod_edit/↩︎