Boxplot con media e deviazione standard

Quando si crea un boxplot in R, automaticamente vengono rappresentati graficamente la mediana, il primo e terzo quartile (i due "hinges", ossia gli estremi del rettangolo) e l'intervallo di confidenza al 95% della mediana (gli estremi dei "notches", ossia le linee).

Tuttavia possiamo voler rappresentare graficamente la media, la deviazione standard (media + DS, media - DS) e i valori massimo e minimo.

Di seguito il codice per rappresentare graficamente, usando ggplot2 tale versione del boxplot; si aggiunge al grafico anche la rappresentazione dei singoli valori come punti, spostati grazie alla funzione jitter.

	library(ggplot2)
	# create fictitious data
	a <- runif(10)
	b <- runif(12)
	c <- runif(7)
	d <- runif(15)
	# data groups
	group <- factor(rep(1:4, c(10, 12, 7, 15)))
	# dataframe
	mydata <- data.frame(c(a,b,c,d), group)
	names(mydata) <- c("value", "group")
	# function for computing mean, DS, max and min values
	min.mean.sd.max <- function(x) {
	r <- c(min(x), mean(x) - sd(x), mean(x), mean(x) + sd(x), max(x))
	names(r) <- c("ymin", "lower", "middle", "upper", "ymax")
	r
	}
	# ggplot code
	p1 <- ggplot(aes(y = value, x = factor(group)), data = mydata)
	p1 <- p1 + stat_summary(fun.data = min.mean.sd.max, geom = "boxplot") + geom_jitter(position=position_jitter(width=.2), size=3) + ggtitle("Boxplot con media, 95%CI, valore min. e max.") + xlab("Gruppi") + ylab("Valori")

view raw ggplot1 hosted with ❤ by GitHub

Cochran Q Test in R (Nonparametric Test for Categorical Data)

Il Cochran Q test è un test non parametrico sviluppato per l'analisi di valori nominali. L'ipotesi che si vuole testare è: nel dataset costituito da k (k maggiore o uguale a 2) campioni dipendenti, esistono almeno due campioni che provengono da differenti popolazioni?



Il test Q di Cochran rappresenta una estensione del test di McNemar, il quale può essere applicato con analoghi risultati quando k = 2.
Ma vediamo subito un esempio, per comprendere il campo di applicazione di questo test.

---

Una società svolge una ricerca di mercato, e chiede a 12 soggetti di compilare un questionario. Tra le domande, figura il quesito: "acquisteresti il prodotto X da una o più delle seguenti 3 marche (A, B, C)?". I 12 soggetti rispondono con delle crocette per marcare la loro preferenza di marca. Le possibili risposte variano da "nessuna crocetta" a "tre crocette" sulle 3 marche.
Risulta subito evidente che i dati di cui disponiamo sono binomiali (categorici): 0 indica che la marca in questione non è stata barrata, mentre 1 indica che la marca è stata crocettata.

La tabella seguente presenta i dati che vengono così raccolti:



Le ipotesi del problema sono: la proporzione di "sì" per ciascuna condizione sono uguali; l'ipotesi alternativa è che almeno una delle proporzioni differisce dalle altre.



La statistica test è:



I simboli che qui compaiono sono:

k = numero di gruppi (o condizioni);


La tabella seguente mostra i calcoli che vengono eseguiti:



Calcoliamo quindi il Q di Cochran:



Il Q di Cochran segue una distribuzione Chi-quadro, con k-1 gradi di libertà. Il valore critico tabulato è:



Poichè Q è maggiore di Chi-quadro, rifiutiamo l'ipotesi nulla H0.

Si pone adesso il problema di identificare quale/i gruppo/i differisce significativamente dagli altri, rendendo quindi significativo il test Q.
Per fare questo, possiamo eseguire il test di McNemar, a coppie. Del test ne parlerò in un prossimo post, ma è sufficiente una rapida ricerca su Google per ottenere pagine esaustive.
Calcolando il test di McNemar per le 3 coppie di paragone, dobbiamo correggere il livelli di significatività (in quanto stiamo eseguendo pur sempre dei confronti multipli). Poichè le coppie da analizzare sono 3, dividiamo il valore di alpha per 3, e quindi il livello di significatività deve essere fissato pari a 0.0167.

---

Soluzione in R

mydata <- matrix(c(
                 1, 1, 0,
                 0, 1, 0,
                 1, 1, 1,
                 0, 1, 0,
                 0, 1, 0,
                 0, 1, 1,
                 0, 0, 0,
                 0, 1, 0,
                 1, 1, 0,
                 0, 1, 0,
                 0, 0, 0,
                 0, 0, 1),
              nrow = 12,
              byrow = T,
              dimnames = list(1 : 12,
                 c("MarcaA", "MarcaB", "MarcaC")))

Utilizziamo il seguente codice, che ho trovato a questo link:

cochranq.test <- function(mat)
{
  k <- ncol(mat)

  C <- sum(colSums(mat) ^ 2)
  R <- sum(rowSums(mat) ^ 2)
  T <- sum(rowSums(mat))

  num <- (k - 1) * ((k * C) - (T ^ 2))
  den <- (k * T) - R

  Q <- num / den

  df <- k - 1
  names(df) <- "df"
  names(Q) <- "Cochran's Q"

  p.val <- pchisq(Q, df, lower = FALSE)

  QVAL <- list(statistic = Q, parameter = df, p.value = p.val,
               method = "Cochran's Q Test for Dependent Samples",
               data.name = deparse(substitute(mat)))
  class(QVAL) <- "htest"
  return(QVAL)
}

Applichiamo la funzione:

cochranq.test(mydata)

        Cochran's Q Test for Dependent Samples

data:  mydata 
Cochran's Q = 8, df = 2, p-value = 0.01832

In alternativa, possiamo usare la funzione cochran.test della library outliers.

Adesso analizziamo le coppie di gruppi con il McNemar's Chi-squared Test for Count Data:

mydatadf <- as.data.frame(mydata)

AB <- with(mydatadf, table(MarcaA, MarcaB))
BC <- with(mydatadf, table(MarcaB, MarcaC))
AC <- with(mydatadf, table(MarcaA, MarcaC))

mcnemar.test(AB)

        McNemar's Chi-squared test with continuity correction

data:  AB 
McNemar's chi-squared = 4.1667, df = 1, p-value = 0.04123

mcnemar.test(BC)

        McNemar's Chi-squared test with continuity correction

data:  BC 
McNemar's chi-squared = 3.125, df = 1, p-value = 0.0771

mcnemar.test(AC, correct=T)

        McNemar's Chi-squared test

data:  AC 
McNemar's chi-squared = 0, df = 1, p-value = 1

Test di Friedman: confronti multipli non parametrici

Il test di Friedman viene usato per comparare tra loro più di due campioni dipendenti, in termini di popolazione. Il confronto che viene fatto si basa sulla mediana della popolazione.
La statistica test (Fr) viene calcolata con le seguenti due formule.


La prima si applica se non ci sono ripetizioni, mentre la seconda nel caso in cui ci sono ripetizioni nella assegnazione del rango.





Nella formula, n è il numero di righe nella tabella dei dati; k è il numero delle colonne (o condizioni); Ri è la somma dei ranghi per ciascuna colonna; Cf è la correzione apportata nel caso, appunto, in cui siano presenti delle ripetizioni dei ranghi, ed equivale a (1/4)nk(k+1)^2; rij è il rango assegnato a ciascuna cella della tabella.
I gradi di libertà sono pari a df = k-1.
La statistica test segue una distribuzione di Friedman; per campioni molto grandi (molti gruppi, o molti elementi per gruppo) allora si usa la distribuzione Chi-quadro, per la scelta del valore critico.

Consideriamo i seguenti esempi.

---

Il gestore di un locale ha 7 dipedenti, i quali sono sempre ritardatari. Per cercare di risolvere tale problema, adotta due strategie. Durante il primo mese (di strategia) decide di premiare i dipendenti che non fanno ritardo; il mese successivo invece punisce i dipedenti che ritardano. Raccoglie gli orari di arrivo (come differenza tra l'ora d'inizio del lavoro, e l'ora di arrivo), ne calcola le medie, e ottiene una tabella come la seguente:



Il gestore vuole sapere se uno dei due metodi applicati apporti un cambiamento nei ritardi dei suoi dipedenti. Le ipotesi del problema sono quindi:



Per prima cosa, assegnamo il rango a ciascuna cella, disponendo in ordine crescente le celle in base alla riga. Otteniamo la seguente tabella:



Calcoliamo ora la somma di ciascuna colonna della riga, ottenendo:



Applicando la prima equazione, con n=7 e k=3; otteniamo il valore:



Il valore critico, per n=7 e k=3 è:



Poichè quindi il valore calcolato è inferiore al valore tabulato, non possiamo rifiutare l'ipotesi nulla.

---

Un allenatore di pallavolo di 12 ragazzi vuole sperimentare due nuovi tipi di allenamento. Per valutarne l'efficacia, usa come parametro il numero di piegamenti che ciascun ragazzo riesce a fare in un minuto. Raccoglie i dati in una tabella: la prima colonna contiene il numero (medio) di piegamenti che i ragazzi riescono a fare di base (prima degli allenamenti), mentre le altre due colonne fanno riferimento ai risultati ottenuti dopo i due nuovi tipi di allenamento, durante il primo e il secondo mese.



A questi valori vengono assegnati i seguenti ranghi. Poichè sono presenti delle ripetizioni, le celle con i valori ripetuti conterranno la media dei ranghi:



Le ipotesi del problema sono:



Applichiamo la formula per calcolare il valore di F, con la correzione per le ripetizioni:



Essendo i dati numerosi, possiamo confrontare tale valore con la distribuzione Chi-quadro, a 2 gradi di libertà. Poichè il valore tabulato è stavolta inferiore al valore calcolato, concludiamo che la mediana di uno (o più) dei gruppi considerati è significativamente diversa (p < 0.05).
Per capire adesso quali confronti risultano significativi, possiamo applicare il Wilcoxon signed rank test.

Con questo test confrontiamo i 3 gruppi, a coppie. Poichè stiamo facendo un confronto multiplo, applichiamo la correzione di Bonferroni, impostando quindi alpha = 0.05 / 3 = 0.0167.
Lascio a voi il prosieguo dei calcoli. Qui riporto soltanto i risultati:

• Baseline - mese1: non c'è differenza (p > 0.0167)


• Baseline - mese2: i due gruppi sono differenti (p < 0.0167)


• Mese1 - mese2: i due gruppi sono differenti (p < 0.0167)

Qual è quindi il tipo di allenamento migliore? A voi la risposta.

---

---

Risolviamo i due problemi con R.

Problema 1:

mydata <- matrix(c(
     16, 13, 12,
     10,  5, 2,
     7,  8, 9,
     13, 11, 5,
     17,  2, 6,
     10,  7, 9,
     11,  6, 7),
     nrow = 7,
     byrow = T,
     dimnames = list(1 : 7,
           c("Baseline", "Mese 1", "Mese 2")))

friedman.test(mydata)

        Friedman rank sum test

data:  mydata 
Friedman chi-squared = 5.4286, df = 2, 
   p-value = 0.06625

---

Problema 2:

mydata <- matrix(c(
     66, 67, 69,
     49, 50, 56,
     51, 52, 49,
     65, 65, 69,
     42, 43, 46,
     38, 39, 40,
     33, 31, 39,
     41, 41, 44,
     46, 47, 48,
     45, 46, 46,
     36, 33, 34,
     51, 55, 67),
     nrow = 12,
     byrow = T,
     dimnames = list(1 : 12,
           c("Baseline", "Mese 1", "Mese 2")))

friedman.test(mydata)

        Friedman rank sum test

data:  mydata 
Friedman chi-squared = 10.9778, df = 2, 
 p-value = 0.004132

CDC Growth Charts implementate in R

Con gran fatica e sudore, sono riuscito a implementare l'intero set di Growth Charts messo a disposizione dalla CDC (Centers for Disease Control and Prevention).

Volendo evitare di riscrivere (e cercare nuovamente i links perduti dei vari documenti/codici), rinvio alla lettura del mio post inglese.
...

Implementation of the CDC Growth Charts in R

Per qualsiasi commento, proposta, potete scrivere su una qualsiasi delle due pagine, essendo equivalenti.

R è anche un ottimo sound editor!

Ebbene sì, le capacità di R! Sono sempre più convinto che possa essere un ottimo concorrente del ben più noto MatLab.

Per divertimento, ho creato una funzione, che riceve in input una sequenza di numeri, e come output produre la melodia che potreste udire, digitando la stessa sequenza con il vostro telefon di casa! =)
Supporta anche le lettere A, B, C, D, e i simboli * e #.

Questa funzione richiede il pacchetto sound , ed ecco di seguito il codice (commentato in inglese, per comodità mia, così da postarlo anche sulla versione inglese di questo blog).

	PlayTel <- function(x=character){
	# load package
	require(sound)
	# control if the input is a character
	if (!is.character(x)) stop("x must be a character")
	# split the input
	num <- strsplit(as.character(x), "")
	len <- length(num[[1]])
	# create the sequence of time
	t <- seq(0, 0.2, length=44100 * 0.2)
	# control if all characters are in the group of those supported
	poss <- c("1", "2", "3", "4", "5", "6", "7", "8", "9", "0", "*", "#", "A", "a", "B", "b", "C", "c", "D", "d")
	if(prod(num[[1]] %in% poss) == 0) stop("Character not supported")
	# create a silence of 0.1 sec
	ts <- as.Sample(seq(0, 0.1, length=44100 * 0.1), 44100, 16)
	# create an empty list
	ll <- vector("list", length=len)
	# here are the frequencies of the tones =)
	t1 <- appendSample(as.Sample((sin(4379*t)+sin(7596*t))/2, 44100, 16), ts)
	t2 <- appendSample(as.Sample((sin(4379*t)+sin(8394*t))/2, 44100, 16), ts)
	t3 <- appendSample(as.Sample((sin(4379*t)+sin(9280*t))/2, 44100, 16), ts)
	t4 <- appendSample(as.Sample((sin(4838*t)+sin(7596*t))/2, 44100, 16), ts)
	t5 <- appendSample(as.Sample((sin(4838*t)+sin(8394*t))/2, 44100, 16), ts)
	t6 <- appendSample(as.Sample((sin(4838*t)+sin(9280*t))/2, 44100, 16), ts)
	t7 <- appendSample(as.Sample((sin(5353*t)+sin(7596*t))/2, 44100, 16), ts)
	t8 <- appendSample(as.Sample((sin(5353*t)+sin(8394*t))/2, 44100, 16), ts)
	t9 <- appendSample(as.Sample((sin(5353*t)+sin(9280*t))/2, 44100, 16), ts)
	t0 <- appendSample(as.Sample((sin(5912*t)+sin(8394*t))/2, 44100, 16), ts)
	tas <- appendSample(as.Sample((sin(5912*t)+sin(7596*t))/2, 44100, 16), ts)
	tca <- appendSample(as.Sample((sin(5912*t)+sin(9280*t))/2, 44100, 16), ts)
	ta <- appendSample(as.Sample((sin(4379*t)+sin(10260*t))/2, 44100, 16), ts)
	tb <- appendSample(as.Sample((sin(4838*t)+sin(10260*t))/2, 44100, 16), ts)
	tc <- appendSample(as.Sample((sin(5353*t)+sin(10260*t))/2, 44100, 16), ts)
	td <- appendSample(as.Sample((sin(5912*t)+sin(10260*t))/2, 44100, 16), ts)
	# full the list with the tone sequence
	for(i in 1:len){
	if(num[[1]][i] == "1") ll[[i]] = t1
	if(num[[1]][i] == "2") ll[[i]] = t2
	if(num[[1]][i] == "3") ll[[i]] = t3
	if(num[[1]][i] == "4") ll[[i]] = t4
	if(num[[1]][i] == "5") ll[[i]] = t5
	if(num[[1]][i] == "6") ll[[i]] = t6
	if(num[[1]][i] == "7") ll[[i]] = t7
	if(num[[1]][i] == "8") ll[[i]] = t8
	if(num[[1]][i] == "9") ll[[i]] = t9
	if(num[[1]][i] == "0") ll[[i]] = t0
	if(num[[1]][i] == "*") ll[[i]] = tas
	if(num[[1]][i] == "#") ll[[i]] = tca
	if(num[[1]][i] == "A" | num[[1]][i] == "a") ll[[i]] = ta
	if(num[[1]][i] == "B" | num[[1]][i] == "b") ll[[i]] = tb
	if(num[[1]][i] == "C" | num[[1]][i] == "c") ll[[i]] = tc
	if(num[[1]][i] == "D" | num[[1]][i] == "d") ll[[i]] = td
	}
	# make sound
	s <- as.Sample(unlist(ll))
	# make stereo sound
	sst <- stereo(s, s)
	sst
	}

view raw PlayTel.R hosted with ❤ by GitHub

Adesso potete comporre la vostra melodia telefonica =)

s2 <- PlayTel("556c885a4623#")

E' possibile ascoltare il file wave appena creato con:

play(s2)

(NOTA: in Windows 7 non esiste alcun programma (come il vecchio mplay32.exe), che possa essere usato dal prompt dei comandi. Per questo motivo, possedendo io stesso Windows 7, non sono stato capace di settare correttamente il pacchetto per ascoltare i suoni prodotti. Con le altre versione di Windows è invece possibile)

E' infine possibile salvare il file audio:

saveSample(s2, "tel.wav")

(Il salvataggio del file wave funziona perfettamente anche con Windows 7)

E questo è un esempio:



Have fun!! =)