Are my p-values reliable?

Several collegues of mine have an ‘unfriendly’ relationship with biometry. They just look for p-values. When I ask them: “Do you really need a p-value for this statement?”, they usually answer: “No, I don’t. But editors do”. I must admit that this is indeed true: the role of p-values has been extremely important on the development of rigorous scientific procedures and it is no wonder that adding such a measure of ‘strength’ to all main statements in a scientific paper would seem highly appropriate. However, not everybody agrees on this and there are now several dissenting voices (see for example this post here). I think that these voices should be taken as warning messages against the abuse of p-values as the sole and final arbiter of whether the null hypothesis is true or false. In particular, the importance of high sample sizes and confirmatory studies has been frequently highlighted.

To my opinion, there is nothing wrong with p-values, though we should all learn to understand their real meaning. In particular, it is necessary to reinforce the idea that p-values will never measure the probability of the hypothesis given a certain body of data! On the contrary, they measure the probability of the data given a certain null hypothesis. Such a difference may seem small for us biologists, but it is of extreme importance! In other words, we should never say:

1. there is less than 5% probability that the null is true;
2. there is less than 5% probability that I’m wrong, if I reject the null.

These statements are nonsense, in frequentist statistics. What we can say is:

1. If the null is right, there is less than 5% probability that I obtain such a dataset. Therefore I reject the null.

Obviously, it is fundamental that we make sure that our p-values make sense, which is a far too often neglected aspect! All statistical tests make assumptions and it is our sole responsibility to ensure that those assumptions were met, before reporting p-values (and editors should make sure that the authors made sure about this!)

With reference to F tests in the ANOVA, these rely on the basic assumptions for linear models, i.e. errors are normal, independent and identically distributed (i.i.d). Independendence should be ensured by an appropriate experimental design, but we should always check that:

1. there are no outliers;
2. errors are gaussian;
3. errors are homoscedastic.

How should we do this?

A case study

Let’s consider the following dataset:

sugarData <- structure(data.frame(
  Treat = structure(c(1L, 1L, 1L, 1L, 1L, 2L, 2L, 
2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 
5L, 5L, 6L, 6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 
8L, 9L, 9L, 9L, 9L, 9L, 10L, 10L, 10L, 10L, 10L, 11L, 11L, 11L, 
11L, 11L, 12L, 12L, 12L, 12L, 12L, 13L, 13L, 13L, 13L, 13L, 14L, 
14L, 14L, 14L, 14L, 15L, 15L, 15L, 15L, 15L, 16L, 16L, 16L, 16L, 
16L, 17L, 17L, 17L, 17L, 17L, 18L, 18L, 18L, 18L, 18L, 19L, 19L, 
19L, 19L, 19L, 20L, 20L, 20L, 20L, 20L, 21L, 21L, 21L, 21L, 21L, 
22L, 22L, 22L, 22L, 22L, 23L, 23L, 23L, 23L, 23L, 24L, 24L, 24L, 
24L, 24L), .Label = c("1", "2", "3", "4", "5", "6", "7", "8", 
"9", "10", "11", "12", "13", "14", "15", "16", "17", "18", "19", 
"20", "21", "22", "23", "24"), class = "factor"), 
  Block = structure(c(1L, 
2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 
3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 
4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 
5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 
1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 
2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 
3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 
4L, 5L, 1L, 2L, 3L, 4L, 5L), .Label = c("1", "2", "3", "4", "5"
), class = "factor"), 
  Sugar = c(4.812, 5.515, 5.546, 4.444, 5.895, 
3.138, 3.113, 2.007, 2.877, 3.306, 5.516, 8.151, 7.258, 9.246, 
6.212, 7.614, 9.283, 7.019, 8.015, 5.117, 11.477, 10.671, 10.567, 
9.8, 10.111, 10.019, 8.686, 8.904, 8.956, 7.975, 4.356, 5.689, 
7.35, 4.767, 5.122, 7.591, 8.071, 9.086, 9.285, 8.312, 5.806, 
6.524, 7.365, 7.728, 8.784, 6.063, 5.833, 3.623, 7.196, 6.464, 
8.117, 9.526, 9.901, 10.175, 8.911, 7.608, 8.442, 8.272, 7.738, 
9.684, 3.699, 4.14, 3.073, 3.979, 2.199, 3.618, 3.604, 2.917, 
5.287, 1.698, 2.593, 3.808, 4.93, 2.649, 4.607, 6.584, 6.437, 
7.934, 7.365, 5, 6.508, 5.74, 3.999, 7.025, 6.699, 5.699, 5.918, 
5.341, 5.119, 6.984, 4.945, 5.471, 6.563, 3.81, 4.531, 5.345, 
6.564, 6.81, 6.878, 8.586, 5.968, 6.306, 6.28, 4.634, 6.103, 
6.85, 5.914, 6.447, 7.146, 7.686, 9.014, 7.96, 9.765, 5.587, 
7.254, 6.712, 7.4, 8.494, 7.766, 7.497)), 
  .Names = c("Treat", "Block", "Sugar"))

head(sugarData)
##   Treat Block Sugar
## 1     1     1 4.812
## 2     1     2 5.515
## 3     1     3 5.546
## 4     1     4 4.444
## 5     1     5 5.895
## 6     2     1 3.138

This dataset comes from an experiment based on one experimental treatment (Treat) with 24 levels, which is investigated in relation to its effect on the sugar yield of sugarbeet. This experiment was performed in a randomised complete blocks design and the ANOVA can be performed as:

mod.aov <- lm(Sugar ~ Block + Treat, data=sugarData)
anova(mod.aov)
## Analysis of Variance Table
## 
## Response: Sugar
##           Df Sum Sq Mean Sq F value Pr(>F)    
## Block      4   2.65  0.6627  0.6125 0.6547    
## Treat     23 443.57 19.2855 17.8262 <2e-16 ***
## Residuals 92  99.53  1.0819                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

In this case the F test yields highly significant results, but we should still be concerned about the meaning of such a p-value: is it reliable? We need to check for the basic assumptions.

Graphical analysis of residuals

A graphical inspection of residuals is enough to spot a main part of the problems with basic assumptions. This can be obtained by:

par(mfrow=c(2,2))
plot(mod.aov)

Which is quite reassuring: the ‘residual vs. fitted’ plot shows a random distribution of residuals (there are no visible patterns), while the ‘QQ plot’ shows that data lie along the bisector line. Very good, then!

Some possible ‘suspect’ patterns for the ‘residual vs. fitted’ plot are like this:

Testing for homoscedasticity

Sometimes a more formal test is required. One possibility is to use the Bartlett’s test for homoscedasticity:

bartlett.test(Sugar ~ Treat, data=sugarData)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  Sugar by Treat
## Bartlett's K-squared = 18.1506, df = 23, p-value = 0.7493

which confirms that there are no problems. Be warned: this test rely on the normality assumption! A more robust alternative is the Levene’s test (be careful to remove the Block effect!):

library(car)
## Warning: package 'car' was built under R version 3.0.3
leveneTest(Sugar ~ Treat, data=sugarData)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group 23  0.6047 0.9159
##       96

This is also reassuring, in our case.

Testing for normality

Normality is rarely an issue, as the ANOVA is pretty robust against this type of violation. If we need a formal test, we can use the Shapiro-Wilk test, that makes use of the standardised residuals from model fit (the same as used for the QQ-plot):

standardRes <- residuals(mod.aov)/summary(mod.aov)$sigma
shapiro.test(standardRes)
## 
##  Shapiro-Wilk normality test
## 
## data:  standardRes
## W = 0.9876, p-value = 0.3452

which confirms what we have already observed on the QQ-plot. Please, do not give too much importance to this test: the graphical inspection of the QQ-plot is much more useful!

Help me: I have problems!

Only after this fundamental checks we can be sure that the p-values in ANOVA are reliable. In our case, everything was OK, but what if some violation is discovered?

The simplest answer is that you need to transform your data. Instead of selecting a transformation by common sense (which is by the way a correct behaviour in most cases), we can select a family of transformations, such as:

\[ W = Y^\lambda;\,\, if \lambda \ne 0 \]

\[ W = log (Y); \,\, if \lambda = 0 \]

where W is the transformed variable, Y is the original variable and \(\lambda\) is the trasformation parameter. By changing this latter we obtain almost all the main transformations, i.e. 1 = untransformed data, 0 = logarithmic transformation, 0.5 = root transformation …

We can select \(\lambda\) by maximum likelihood, with the help of a graphical tool, that is available in R as:

library(MASS)
## Warning: package 'MASS' was built under R version 3.0.3
boxcox(mod.aov)

Obviously, in this case the maximum likelihood value for \(\lambda\) is close to 1, which confirms that this dataset does not require any correcting measure.

Should a transformation be required, do not be scared: transform your data and repeat the analyses. Unfortunately, the discussion and presentation of results becomes slightly more complex… You may take a look at my other post on this subject.

Further readings

• Ahrens, W. H., D. J. Cox, and G. Budwar. 1990, Use of the arcsin and square root trasformation for subjectively determined percentage data. Weed science 452-458.
• Anscombe, F. J. and J. W. Tukey. 1963, The examination and analysis of residuals. Technometrics 5: 141-160.
• Babbini, M., B. Chiandotto, G. Chisci, R. d. Cristofaro, G. A. Maccararo, N. Montanaro, F. Nicolis, E. Ottaviano, F. Salvi, and M. Turri. 1978. Biometria: principi e metodi per studenti e ricercatori biologi. Padova: P. 552.
• Box, G. E. P. and D. R. Cox. 1964, An analysis of transformations. Journal of the Royal Statistical Society, B-26, 211-243, discussion 244-252.
• Camussi , A., F. Moller , E. Ottaviano , and M. Sarli Gorla . 1986, Metodi statistici per la sperimentazione biologica. Ed. Zanichelli.
• Chatfield, C. 1985, The initial examination of data. J. R. Statist. Soc. A-148, 3, 214-253 A-148: 214-253.
• D’Elia, A. 2001, Un metodo grafico per la trasformazione di Box-Cox: aspetti esplorativi ed inferenziali. STATISTICA LXI: 630-648.
• Draper, N. R. and H. Smith. 1981, Applied regression. John Wiley & Sons, Inc. , New York, 2nd ed.
• Saskia, R. M. 1992, The Box-Cox transformation technique: a review. Statistician 41: 169-178.
• Snedecor, G. W. and W. G. Cochran. 1991. Statistical methods. AMES (Iowa): IOWA State University Press, VIII Edition. P. 503.

Varianza del campione e della popolazione: perchè le formule sono diverse?

Una delle domande più frequenti quando si inizia un corso di metodologia sperimentale è: “perchè la varianza del campione si calcola utilizzando come divisore n-1 invece di n”?

La risposta è non banale, soprattutto per gli studenti delle Facoltà biologiche, le cui basi di matematica sono spesso elementari. Occorre ricordare che la varianza di una popolazione è il valore atteso del quadrato degli scarti tra ogni dato e la media \( \mu \). Se abbiamo a che fare con una popolazione piccola, per esempio composta dai tre valori 99, 100 e 101, la media \( \mu \) è uguale a 100 e gli scostamenti sono rispettivamente pari a -1, 0 ed 1. La somma dei loro quadrati è pari a 2 e il valore atteso dello scarto quadratico (appunto lo scarto quadratico medio, cioè la varianza) è pari 2/3 = 0.667. In generale:

\[ {\sigma ^2} = \frac{{\sum\limits_{i = 1}^n {{{({X_i} - \mu )}^2}} }}{n} \]

Tuttavia, quando i tre valori anzidetti non costituiscono l'intera popolazione, ma rappresentano un campione estratto da essa e vogliamo stimare \( \sigma^2 \) (cioè la varianza della popolazione, non quella del campione), dobbiamo considerare che in mano non abbiamo \( \mu \), ma abbiamo la media del campione \( \bar{X} \). Quest'ultima non coincide con \( \mu \), anche se ne rappresenta una stima non distorta.

Se ammettiamo che la media del campione è diversa da \( \mu \) e \( \bar{X} - \mu = \varepsilon \), possiamo scrivere che:

\[ \sum\limits_{i = 1}^n {{{({X_i} - \mu )}^2}} = \sum\limits_{i = 1}^n {{{({X_i} - \bar X + \varepsilon )}^2} = \sum\limits_{i = 1}^n {{{({X_i} - \bar X)}^2} + n{\varepsilon ^2}} } \]

Da qui comprendiamo che la devianza del campione (somma dei quadrati degli scarti) è maggiore dalla devianza della popolazione se \( \varepsilon > 0 \). Possiamo quindi scrivere:

\[ \sum\limits_{i = 1}^n {{{({X_i} - \bar X)}^2} = } \sum\limits_{i = 1}^n {{{({X_i} - \mu )}^2}} - n{\left( {\bar X - \mu } \right)^2} \]

Ricordando la formula della varianza della popolazione, abbiamo:

\[ \sum\limits_{i = 1}^n {{{({X_i} - \bar X)}^2} = } n{\sigma ^2} - n{\left( {\bar X - \mu } \right)^2} \]

A questo punto possiamo vedere che il termine a destra costituisce lo scarto quadratico della media del campione rispetto alla media della popolazione, cioè il quadrato dell'errore standard della media. Quindi:

\[ \sum\limits_{i = 1}^n {{{({X_i} - \bar X)}^2} = } n{\sigma ^2} - n\left( {\frac{{{\sigma ^2}}}{n}} \right) \]

cioè:

\[ \sum\limits_{i = 1}^n {{{({X_i} - \bar X)}^2} = } (n - 1){\sigma ^2} \]

A questo punto il gioco è fatto. Infatti:

\[ \frac{{\sum\limits_{i = 1}^n {{{({X_i} - \bar X)}^2}} }}{{n - 1}} = {\sigma ^2} \]

Insomma, se abbiamo un insieme di dati e vogliamo sapere la sua varianza, dividiamo la somma dei quadrati degli scarti per n. Invece nel caso in cui l'insieme di dati rappresenta un campione estratto da una popolazione più ampia e vogliamo stimare la varianza della popolazione (non quella del campione) allora dobbiamo dividere per n-1. Non dimenticate che questo secondo caso è di gran lunga più frequente nella sperimentazione biologica!

Reporting results for designed experiments in applied biology

A great majority of designed experiments in applied biology (especially field experiments) is aimed at comparing several treatment levels on a randomised block design (or similar), with few replicates (normally 3 to 5). In the end, biologists are usually interested in reporting the results by way of a table/graph of observed means, together with some appropriate measure of uncertainty. The presence of this latter measure is to be regarded as mandatory for publication in all scientific journals, though its selection may pose a few problems and become an obstacle during the reviewing process.

How do we make a reasonable choice? Two things should be kept in mind: (1) for each treatment level we have a sample of measures taken from a wider population: it is the population that matters, not the sample! (2) Both the mean and the uncertainty measure should preferably be on the same measurement unit. Bearing this in mind, the choice mainly depends on what distributional assumptions we are willing to make about the population that generated our dataset. Here are some simple suggestions.

1. If we do not want to/cannot make any assumptions about the original population and we would only like to describe the variability of each sample around the mean, we should better report the standard deviation (SD) for each sample.

2. If we can assume that our samples come from normal populations and we would like to describe the variability of the means, that could be expected by repeatedly performing the same experiment, we'd better report the standard error for each mean (SE). This latter is a good aid to express the repeatability of results beyond the observed sample, though it is not a good measure of the variability of data around the mean (the SE is always lower than the SD). Therefore, we may consider adding to the table/graph the numerosity of samples (n), so that everybody can go back from the SE to the SD.

3. If we can assume that the sampled populations are normal and with similar variances (homoscedasticity) we will probably go for an ANOVA. In this case, we can describe the variability of the means (as in #2) by using the pooled standard error of a mean (SEM), derived from the residual mean square in ANOVA. One and the same value is ok for all means, if the design is balanced! Also in this case, we may consider adding the numerosity of samples, to help readers go back and determine the corresponding SD.

4. In the same situation as above (normality and homoscedasticity), we might like to guide the reader to making comparisons among treatment levels. To this aim, we might like to replace the SEM with the pooled standard error of a difference (SED). Reporting the degrees of freedom for SEDs may be recommended, to help readers calculate a critical difference for each comparison (Least Significant Difference or Honest Significant Difference). This latter display has some limitations, which should be carefully considered. First of all, the adoption of pairwise comparisons is inelegant whean dealing with a quantitative series of treatment levels (that should be preferably compared on a regression settings). Secondly, in case of unbalanced data, there will be different SED values for each mean, depending on which comparison we intend to make.

What should we do when we have adopted some stabilising transformation prior to ANOVA?

Linear models (ANOVA and regression) are routinely used in applied biology. Whenever normality and homoscedasticity of residuals cannot be assumed, we may go for the adoption of some sort of stabilising transformations (logarithm, square root, arcsin square root…), prior to ANOVA. In this case, the above suggestions still hold if we intend to show the means of transformed data. However, such a display may hinder the clarity of results, as the original measurement unit is lost. What should we do, then? Here are some further suggestions.

1. As usual, we can present the means of original data with SDs (as in #1 above). This is clearly less than optimal, if we want to suggest more than the bare variability of the observed sample. Furthermore, please remember that the means of original data may not be a good measure of central tendency, if the original population is strongly 'asymmetric' (skewed)!

2. The best option in this case is to show back-transformed means (i.e. if we have done, for example, a square root transformation, we can take the means of transformed data and square them). Back-transformed means 'estimate' the medians of the original populations, which may be regarded as better measures of central tendency for skewed data.

3. To complete our work, we should add to each back-transformed mean its back-transformed standard error, by using the so-called delta method. I'll clear this up with an example.

An example with counts

Consider the following dataset, that represents the counts of insects on 15 independent leaves, treated with the insecticides A, B and C (5 replicates):

  Level   1   2   3   4   5   Mean
      A 448 906 484 477 634  589.8
      B 211 276 415 587 298  357.4
      C  50  90  73  44  26   56.6 

This variable represent a count and, as expected, variances for the treatment levels A, B and C are sensibly different. A logarithmic transformation can do the trick here and produce a normal and homoscedastic new dataset. Therefore we take the log-transformed variable and submit it to ANOVA:

Analysis of Variance Table

              Df Sum Sq Mean Sq F value    Pr(>F)    
Level          3 448.19 149.396  946.66 1.624e-14 ***
Residuals     12   1.89   0.158                       

If we were to report the means of the log-transformed variable (log-Means), we might follow the above suggestion and show:

Level  log-Means
A          6.343
B          5.815
C          3.985
SEM        0.178

where the SEM is obtained by taking \( \sqrt{0.158/5} \). Unfortunately, we loose clarity: how many insects did we have on each leaf?

A possible way out of this is to back-transform the above log-Means (for example \( exp(6.343)=568.499 \)) and use the delta method to back-transform the standard error. This is straightforward: (1) take the first derivative of the back-transform function [in this case the first derivative of exp(X)=exp(X)] and (2) multiply it by the standard error of the transformed data. For the case of level A, it is: \( exp(6.343) \times 0.178 = 101.19 \). Therefore, data might be presented as follows:

Level       Back-transformed 
            Means   (SE)
A           568.5 (101.19)
B           335.1 (59.68)
C           51.88 (9.57)

Far clearer, isn't it? If we had done a square root transformation, the back-transform function would be \( X^2 \). The first derivative would be \( 2 \times X \) and this should be multiplied by the SEM to get a back-transformed standard error. If you want to know something more about the delta-method you might start from my post here. Some collegues and I have discussed this issues in our paper
'Current statistical issues in Weed Research' (Onofri A., Carbonell E., Piepho H.-P., Mortimer A.M. & Cousens R.D., 2010. Weed Research, 50, 5-24).

In all cases, whatever measure of uncertainty you might like to use, do not forget to state it clearly and be consistent throughout the paper!

Abbandonare Word: sogno o realtà?

Da almeno 25 anni il mio PC trabocca di files .doc, .xls e .ppt. Non è stata una vera e propria scelta: quando ho iniziato a lavorare Microsoft Office sembrava l'unica possibilità per fare tutto bene, all'interno di un sistema di lavoro razionale. A quel tempo, i miei colleghi biologi più anziani che lavoravano su PC scrivevano i loro lavori in WordStar, utilizzavano Lotus per la gestione dei dati e DB3 per i riferimenti bibliografici e non era affatto facile mettere a collaborare questi tre mondi separati. Office 3 fu un vero schock: ad esempio era finalmente possibile creare tabelle nel foglio elettronico, incollarle nell'elaboratore testi ed unirle al resto del lavoro!

Con il passare del tempo tuttavia, la mia frustrazione è andata crescendo. Infatti, ho iniziato ad accorgermi di come Excel e qualunque altro foglio elettronico confondessero la visualizzazione dei dati con la loro gestione, cosicché ero spesso indotto ad immettere i dati in un formato poco razionale per la loro analisi. D'altra parte, mi sono accorto di come Word e qualunque altro word processor confondessero la composizione del testo con la formattazione, il che era vero (all'ennesima potenza) anche per Power Point e qualunque altro programma per le presentazioni. Il risultato era che finivo per mischiare fasi di lavoro che avrebbero dovuto essere logicamente separate, con spreco di tempo, energie e, in ultima analisi, bassa qualità di lavoro.

Il massimo della frustrazione l'ho raggiunto qualche anno fa, quando ho tentato di aprire con WORD 2010 il file della mia tesi di dottorato, scritto nel 1993 con WORD 2.0. Allora è comparsa una minacciosa finestra di dialogo con qualcosa di simile a: “Questo file non è leggibile per motivi di sicurezza”. Anche se ho poi scoperto che il problema poteva essere aggirato con una semplice modifica del file di registro, per qualche istante sono rimasto senza fiato: “è mai possibile che non sono libero di leggere la mia tesi di dottorato, nel mio computer, con lo stesso programma di videoscrittura usato al tempo (magari in una versione precedente)?”. Chi come me tiene al suo lavoro, lo immagazzina in files di vario formato e lo riusa e rielabora di continuo, non può non essere preoccupato di bloccarlo dentro un formato che non è di pubblico dominio. Questo vale anche per chi, come me, usa molti ambienti di lavoro diversi: dal PC in ufficio al MacBook a casa, all' Ipad, all' Iphone. In questo caso la portabilità è essenziale!

Da allora ho scoperto il formato testo: portabile e immarcescibile, una sicurezza totale. E, col tempo, ho totalmente ripensato il mio flusso: i miei files ora sono quasi tutti in quel formato, magari con diverse estensioni (.r, .mmd, .tex …) a seconda dell'uso. Li raccolgo tutti su Dropbox, in modo che siano disponibili su tutti i miei dispositivi, purchè online e compongo tutti i miei lavori (appunti, relazioni, lavori scientifici) utilizzando un semplice editor di testi (prediligo EditPad, ma uso anche TexStudio o RStudio, per alcune applicazioni). Se debbo inserire della formattazione o altri elementi non testuali (link, immagini ed altro), utilizzo Markdown, un semplice (si impara in un pomeriggio) linguaggio di mark-up che descrive la struttura del documento e lo stile (cosa dovrà costituire il titolo del capitolo o del paragrafo, cosa dovrà essere formattato come elenco puntato o numerato, cosa dovrà andare in grassetto o corsivo). A seconda del formato richiesto per l'output, utilizzo Pandoc per tradurre il file di testo originale da Markdown a PDF (per la stampa) o ad HTML (per il web).

Per l'analisi dei dati, salvo tutti i miei datasets in formato .csv (comma delineated) e li elaboro con R e RStudio, conservando i relativi scripts come files di testo con estensione .r, in modo da essere in grado di ripetere le analisi quando voglio.

Da tempo, le mie presentazioni sono in formato PDF, preparate con Latex e Beamer. Non ho mai problemi di caratteri non disponibili, scritte che non entrano nello schermo, files illeggibili, dovunque vada a fare le mie lezioni e qualunque strumento di presentazione io abbia a disposizione.

La mia collezione di riferimenti bibliografici è anche essa salvata in formato testo e redatta secondo le specifiche di bibtex; sono in grado di citarla facilmente in modo automatizzato (cite while you write).

Infine, quasi tutti i programmi sopra descritti sono di pubblico dominio, spesso freeware e sempre multipiattaforma. Cosa cercare di meglio? Di recente ho letto un post analogo, ma più legato al mondo APPLE che al PC. Potete leggerlo qui e vi consiglio proprio di farlo, anche perché vi sono utili links a siti che trattano Markdown, Latex e Rstudio.

Sono riuscito ad eradicare Office e Word dalla mia vita professionale? Assolutamente no!

Il motivo è che, come dicevo, sono un biologo e, nel mio mondo, tutti, o la maggior parte dei miei colleghi compongono i loro testi in Word, salvano i loro dati in Excel e preparano le loro presentazioni in Power Point. Oppure tengono la loro bibliografia in End Note o altri strumenti simili. E condividono i files risultanti, aspettandosi di ricevere in cambio lo stesso tipo di file. Anche le case editrici spesso chiedono il formato Word per la pubblicazione dei lavori e quindi, alla fine, sono costretto ad investire parte del mio tempo per trasformare i miei files di testo in una versione più 'compatibile' con il 'resto del mondo'. Certo è che, quando il lavoro è collaborativo, utilizzare i miei strumenti preferiti diviene quasi impossibile.

Mi ha colpito (positivamente) un post che raccomanda di non utilizzare il formato .doc (o .docx) per gli attachment di post elettronica. Si può leggere qui. Riusciremo anche noi biologi ad avere un futuro non dico libero da Office (sarebbe esagerato), ma almeno nel quale sia possibile che ognuno, oltre ai suoi strumenti preferiti, sia in grado di operare con una pluralità di strumenti diversi e quindi di interagire anche con chi, per sua scelta, non usi altro che software libero? Forse la mia generazione non ce la farà, ma spero che una scuola e un'Università più responsabili educhino le future generazioni ad un uso più razionale del proprio tempo di computing.

Confidence intervals: am I uncousciously a Bayesian?

Last week I made a survey among the students attending my course in 'Experimental Methods in Agriculture'. One of the questions was:

QUESTION: “I sampled 100 seeds from a big population of Holy Clover (Onobrychis viciifoliae Scop) and found that their average weight was \( \bar{Y} \) = 15.5 mg. The confidence interval for the population mean was 13-18. What is the meaning of such a statement?”.

The possible answers were:

1. There is 95% probability that the population mean is in the interval 13 to 18;
2. If we sample repeatedly from our population of Holy Clover, the estimated confidence intervals will contain the true mean in 95% of cases;
3. The true population mean is certainly between 13 and 18;
4. The true population mean can take any value between 13 and 18.

I asked my students to select the correct answer without looking at textbooks or class notes, just using their memory and intuition. This survey came after the first half of the course, approximately one month later than the lecture about point and interval estimation.

In the end, 75% of my students choose answer (1), while none of them choose answer (2). This came out quite as a shock to me: in a frequentist setting, the correct answer is clearly (2). Indeed, it should be intuitively clear that there is a 'true' (fixed) average weight \( \mu \) for my seed population, but the problem is that I will never come to known it exactely, as the population is too big for me to measure its weight. Therefore, I am forced into taking a small sample and measuring its average weight. My intuition suggests that further samples will show different average weights, but the true \( \mu \) will always be there, unchanged and unknown. Recalling the frequentist definition of probability (from Wikipedia: “the limit of the relative frequency of an event in a large number of trials”) it would seem pretty clear that it makes no sense to attach any sort of probability to the true value of \( \mu \), as this is not going to change at all, during my sampling! This is why answer (1) is nonsensical in a frequentist setting. Furthermore, the confidence interval (13 to 18) that I built from my sample may either contain \( \mu \) or not, but I have no way to favour one of the two situations. And the extremes of the interval (13 and 18) are actually meaningless: when I repeat the sampling I'll very likely get different values and a different interval.

If the above reasoning is so clear, why does not answer (2) come out as a natural choice for students? Why are they intuitively embracing the Bayesian perspective of answer (1), even though I am pretty sure that they have never been exposed to Bayesian thinking at all (as all agriculture students, at least in Italy)? It is clear that I have not done a good job in conveying the correct message, during my lecture! To use the same words as Dennis (Discussion: Should Ecologists become Bayesians. Ecological applications, 6, 1095-1103), I was probably suggesting more than a frequentist confidence interval delivers.

I am not the only one in this position: I am sure than most of my collegues biologists see confidence intervals very much like it is described in answer (1). I do not think we are to blame. Indeed, I went through the books I used to study when I was a PhD student and found, for example, that Sokal and Rohlf (Biometry. W.H. Freeman and Company, 1981) at page 141, report this equation (7.4):

\[ P\{ \bar{Y} - 1.96 \sigma _{\bar{Y}} \leq \mu \leq \bar{Y} + 1.96 \sigma_{\bar{Y}} \} \]

and comment:“Thus the probability, P, is 0.95 that the term \( \bar{Y} - 1.96 \sigma _{\bar{Y}} \) is less then or equal to the parametric mean \( \mu \) and that the term \( \bar{Y} + 1.96 \sigma_{\bar{Y}} \) is greater than or equal to \( \mu \).” Similar statements can be seen in Snedecor and Cochran (Statistical methods. IOWA State University Press, 1991) and in almost all biometry books I have at hand.

The above equation is derived by using simple math and it is certainly correct. But I am wondering: does it make sense in a frequentist setting? I'll leave the answer to the statisticians. As a biologist, I have to admit that, like my students and my collegues, I feel somewhat disappointed by the crude meaning of frequentist confidence intervals. Though the correct answer is (2), I find myself wishing it were (1): that would really be satisfactory! This more or less unconscious feeling may have influenced my lecture about point and interval estimation.

In the end, considering the ironic claim of IJ Good (“People who do not know they area Bayesians are called non-Bayesians”; cited in Kery, 2010. Introduction to WinBugs for ecologists. Academic Press) I am asking myself: am I one of those who are Bayesians, but do not know, yet?

The delta method for dummies

Reporting variances (or standard deviations) for all estimates is normally regarded as a basic requirement for scientific publications. In my work as consultant I am very frequently asked such questions: “I have the measure Q with variance q and the measure W with variance w. What is the variance of exp(Q + W)?”

I recently wrote a post relating to linear transformations (see here), but I feel that being able to menage nonlinear transformations is much more useful. Frequentist statistics offer an approximate solution to this problem, that is normally known as the 'delta' method.

The basic idea is to create a linear approximation to our nonlinear function and use the former to approximate the variance, by using the same methods as given in my above mentioned post. Indeed, every nonlinear function G (close to linear, differentiable) of X can be locally approximated by the tangent line on a point \(\mu\) (in the close vicinity of X), with slope equal to the first derivative G'(mu) (Taylor series approximation to first order):

\[G(X) \simeq G(\mu ) + G'(\mu )(X - \mu )\]

Having transformed a nonlinear function into a linear function, we can approximate the variance by using:

\[Var[G(X)] \simeq [G'(\mu)]^2 Var(X) \]

In words: if I have a nonlinear (close to linear) function G(X), where X is a random variable with small variance and normal errors, the variance of G(X) is approximated by knowing the first derivative G' and the variance of X

Example 1: exp(X)

To calculate the derivatives, we can use an R function:

D(expression(G(x), “x”)

In this case:

D(expression(exp(b0)), "b0")

## exp(b0)

If we have:

X <- 3.5
VarX <- 0.15
exp(X)

## [1] 33.12

VarX * exp(X)^2

## [1] 164.5

In R there is a shortcut function to calculate the delta method, from the car package. The syntax is:

library(car)
deltaMethod(object = c(X = X), g = "exp(X)", vcov = matrix(VarX, 1, 1))

##        Estimate    SE
## exp(X)    33.12 12.83

By using this function we obtain the standard error, i.e. the standard deviation of the sampling distribution. By squaring the above SE (12.82556) we see that the results are obviously the same.

Example 2: 1/X

The first derivative is:

D(expression(1/x), "x")

## -(1/x^2)

Delta standard errors may be obtained as:

X <- 3.5
VarX <- 0.15
1/X

## [1] 0.2857

VarX * ((-1/X)^2)^2

## [1] 0.0009996

With R:

library(car)
deltaMethod(object = c(X = X), g = "1/X", vcov = matrix(VarX, 1, 1))

##     Estimate      SE
## 1/X   0.2857 0.03162

Example 3: log(X)

The first derivative is:

D(expression(log(x)), "x")

## 1/x

Delta standard errors may be obtained as:

X <- 3.5
VarX <- 0.15
log(X)

## [1] 1.253

VarX * (1/X)^2

## [1] 0.01224

Example 4: Inverse logit function

First derivative is:

D(expression(exp(X)/(1 + exp(X))), "X")

## exp(X)/(1 + exp(X)) - exp(X) * exp(X)/(1 + exp(X))^2

Delta standard errors may be obtained as:

X <- 3.5
VarX <- 1.5
exp(X)/(1 + exp(X))

## [1] 0.9707

VarX * (exp(X)/(1 + exp(X)) - exp(X) * exp(X)/(1 + exp(X))^2)^2

## [1] 0.001214

---

The extension to the case of two random variables is straightforward: we need the partial derivatives in each variable, plus all covariances

\[ var[g(x, y,..)] = D S D^T \]

where D is the matrix of partial derivatives with respect to all the variables, S is the variance-covariance matrix of all variables.

Example 5: X Y

First derivatives are:

D(expression(X * Y), "X")

## Y

D(expression(X * Y), "Y")

## X

Delta standard errors may be obtained as:

X <- 3.5
Y <- 4.2
VarX <- 1.5
VarY <- 1.2
CoVar <- 0.8
X * Y

## [1] 14.7

(Y^2) * VarX + (X^2) * VarY + 2 * X * Y * CoVar

## [1] 64.68

With R:

library(car)
sigma <- matrix(c(VarX, CoVar, CoVar, VarY), 2, 2, byrow = T)
deltaMethod(object = c(X = X, Y = Y), g = "X*Y", vcov = sigma)

##       Estimate    SE
## X * Y     14.7 8.042

Example 6: X/Y

First partial derivatives are:

D(expression(X/Y), "X")

## 1/Y

D(expression(X/Y), "Y")

## -(X/Y^2)

Delta standard errors may be calculated as:

X <- 3.5
Y <- 4.5
VarX <- 0.9
VarY <- 0.8
CoVar <- 0.5
X/Y

## [1] 0.7778

(1/Y)^2 * VarX + (-(X/(Y^2)))^2 * VarY + 2 * (1/Y) * (-X/(Y^2)) * CoVar

## [1] 0.02993

With R:

library(car)
sigma <- matrix(c(VarX, CoVar, CoVar, VarY), 2, 2, byrow = T)
deltaMethod(object = c(X = X, Y = Y), g = "X/Y", vcov = sigma)

##     Estimate    SE
## X/Y   0.7778 0.173

Example 7: exp(X + Y)

Partial derivatives are:

D(expression(exp(X + Y)), "X")

## exp(X + Y)

D(expression(exp(X + Y)), "Y")

## exp(X + Y)

Delta standard errors may be obtained as:

X <- 3.5
Y <- 0.5
VarX <- 0.15
VarY <- 0.08
CoVar <- 0.05
exp(X + Y)

## [1] 54.6

(exp(X + Y)^2) * VarX + (exp(X + Y)^2) * VarY + 2 * exp(X + Y) * exp(X + Y) * 
    CoVar

## [1] 983.7

With R:

library(car)
sigma <- matrix(c(VarX, CoVar, CoVar, VarY), 2, 2, byrow = T)
deltaMethod(object = c(X = X, Y = Y), g = "exp(X+Y)", vcov = sigma)

##            Estimate    SE
## exp(X + Y)     54.6 31.36

Two final notes relating to nonlinear transformations are very important:

1. The estimated variance is always approximate.
2. If the original variables are gaussian, the transformed variable is normally not gaussian.

The propagation of measurement errors

Measures are usually divided into two groups: direct and derived. In this latter case we have two separate direct measurements (Q and W) and need to calculate a linear combination Z = AQ + BW where A and B are two generic coefficients. We already know that Q and W are affected by random measurement errors (and thus they are random variables) and we may wonder how such errors propagate to Z.
We'll start by answering a simpler question, i.e. “I have the estimate Q with variance q, what is the variance (var) of A + B Q?” (A and B are two generic coefficients).
The answer is very simple:
\[ var(A + BQ) = A^2 var(Q) = A^2 q \]
and it is very easy to prove. I'll give you just an example: let's consider three values (e.g. 12, 15, 19) with mean equal to 15.33 and variance equal to 12.33. If we transform each of the three values by multiplying it by 3 and adding 5, we can easily verify that the variance of the transformed sample is 3² x 12.3333 = 111.
And what if we have two estimates Q and W, with variances respectively equal to q and w, and we want to estimate the variance of A Q + B W?
In this case we need to know whether the two quantities Q and W are independent or, on the contrary, what their covariance (COV) is. It can be easily proved that:
\[ var(AQ + BW) = {A^2}var(Q) + {B^2}var(W) + 2AB cov(XY) \]
A formal proof can be found in every statistics book, e.g. in Sokal and Rohlf (1981) at pag. 818. We'll show this by an R example:

Q <- c(12, 14, 11, 9)
W <- c(2, 4, 7, 8)
a <- 2
b <- 3
var(a * Q + b * W)

## [1] 35.58

a^2 * var(Q) + b^2 * var(W) + 2 * a * b * cov(Q, W)

## [1] 35.58

This implies that for two uncorrelated variables the variance of the sum equals to the sum of the variances.
Furthermore, if we remember the following equality:
\[r = \frac{cov(Q,W)}{\sqrt {var(Q)var(W)} }\]
where r is the correlation coefficient, then:
\[cov(Q,W) = r sd(Q)sd(W)\]
where sd is the standard deviation. We can modify the equation above accordingly, if we would like to use r, instead of the covariance.
And what if we have more than three quantities to combine?
The above equation can be generalised as (in words): the variance of the sum of n random variables (X1, X2, X3, … Xn) with n coefficient (B1, B2, B3, …, Bn) is equal to the product of the square of each coefficient by the variance of the corresponding variable, plus twice the product of each couple of coefficient by the covariance of the two corresponding variables.
This may be easily posed in matrix notation as:
var[A X] = ASA^T
where X is the vector of random variables to be combined, A is the vector of coefficients. For example, if we want to assess the variance of the combination 3 + 6 + 4, where all coefficients are equal to one, variances are all equal to 0.5 and covariances are all equal to 0.75, we can do this by the following code:

X <- matrix(c(3, 6, 4), 3, 1)
A <- matrix(c(1, 1, 1), 1, 3)
A %*% X

##      [,1]
## [1,]   13

(Sigma <- matrix(c(0.5, 0.75, 0.75, 0.75, 0.5, 0.75, 0.75, 0.75, 0.5), 3, 3, 
    byrow = T))

##      [,1] [,2] [,3]
## [1,] 0.50 0.75 0.75
## [2,] 0.75 0.50 0.75
## [3,] 0.75 0.75 0.50

A %*% Sigma %*% t(A)

##      [,1]
## [1,]    6

Two final notes relating to linear transformations are very important:

1. The estimated variance is always exact.
2. If the original variables are gaussian, also the transformed variable is gaussian.