6.2 Alcune funzioni R per la statistica inferenziale

In questa sezione mostriamo come effettuare con R le analisi inferenziali presentate nel corso.

6.2.1 Inferenza sulla media di una popolazione normale

E’ possibile fare inferenza su una media \(\mu\) di una popolazione normale assumendo che la varianza \(\sigma^2\) della popolazione sia nota oppure no. In questa sezione ci occupiamo solo del secondo caso perché si tratta della situazione più frequente e rilevante nella pratica dell’analisi dei dati. Inoltre, né R né Radiant mettono a disposizione funzioni o comandi che consentono di ottenere rapidamente i risultati nel caso in cui la varianza fosse nota.

Oltre alle funzioni mean e sd già discusse nella Sezione 6.1.1.2.5 per la stima rispettivamente della media e della deviazione standard di un popolazione, R mette a disposizione la funzione t.test() che consente di gestire tutti i casi presentati nel corso per quanto riguarda l’inferenza su medie. La funzione t.test() restituisce contemporaneamente l’intervallo di confidenza e il test relativo al caso considerato.

Nel caso di inferenza su una singola media per una popolazione normale, gli argomenti da specificare sono:

• x, variabile su cui effettuare le analisi
• alternative, tipo di ipotesi alternativa che vogliamo usare, ovvero bilaterale (two.sided) o unilaterale (less o greater, a seconda del segno previsto da \(H_1\))
• mu, valore della media che vogliamo testare sotto \(H_0\) (la quantità che nel corso è stata chiamata \(\mu_0\))
• conf.level, livello di confidenza che si desidera usare

Vediamo subito un esempio usando ancora una volta i dati contenuti nel data frame forbes94. In particolare, definiamo una nuova variabile che chiamiamo ROS, ovvero return on sales, calcolata come il rapporto percentuale dei profitti e i ricavi di vendita:

> forbes94 <- transform(forbes94, ROS = Profits/Sales*100)

Procediamo ora a calcolare l’intervallo di confidenza al 90% per la media della variabile ROS. Inoltre, insieme all’intervallo effettueremo un test per verificare se i dati supportano la conclusione che la media del ROS per la popolazione di aziende da cui questo campione è stato estratto sia diversa dal 6%. In altri termini, indicata con \(\mu\) la media non nota del ROS nella popolazione, con il test confronteremo le seguenti ipotesi

\[\begin{equation*} H_0: \mu = 6\% \qquad \mbox{vs.} \qquad H_1: \mu \ne 6\%. \end{equation*}\]

Come detto, tutti questi risultati possono essere recuperati in un colpo solo tramite la funzione t.test() come segue:

> t.test(x = forbes94$ROS, alternative = "two.sided", mu = 6, conf.level = 0.90)

    One Sample t-test

data:  forbes94$ROS
t = 4.2, df = 800, p-value = 3e-05
alternative hypothesis: true mean is not equal to 6
90 percent confidence interval:
 6.780 7.798
sample estimates:
mean of x 
    7.289 

I risultati indicano che l’intervallo di confidenza al 90% è pari a \((6.7799, 7.7978)\), mentre il test restituisce un p-value molto piccolo, riportato come 3.377e-05⁵¹, il quale permette di rifiutare ampiamente l’ipotesi anche usando un livello di significatività \(\alpha\) di 0.01.

6.2.2 Inferenza sulla proporzione di successi in una popolazione bernoulliana

In analogia a quanto visto per la media di una popolazione normale, la funzione prop.test() fornisce i medesimi calcoli per la proporzione di successi di una popolazione bernoulliana. Gli argomenti della funzione sono:

• x, numero di successi osservati nelle \(n\) osservazioni campionarie
• n, dimensione del campione
• p, valore della proporzione \(p_0\) che vogliamo testare sotto l’ipotesi nulla
• alternative, tipo di test che vogliamo effettuare, ovvero se bilaterale (two.sided) o unilaterale (less o greater, a seconda della direzione richiesta in \(H_1\))
• conf.level, livello di confidenza che si desidera usare
• correct, con cui è possibile indicare ad R se usare o meno una correzione, detta di Yates, che non è stata presentata nel corso e che quindi noi fisseremo sempre al valore FALSE

Anche in questo caso, la funzione prop.test() calcola automaticamente sia l’intervallo di confidenza sia il test sulla proporzione.

Supponiamo di voler stimare la proporzione di CEO con il titolo di MBA nella popolazione di riferimento usando il campione disponibile nel data frame forbes94. In aggiunta, vogliamo anche calcolarne l’intervallo di confidenza al 99% e ci interessa testare l’ipotesi nulla che la proporzione di CEO in possesso del titolo di MBA nella popolazione sia uguale al 30%. Il codice seguente consente di effettuare congiuntamente queste analisi:

> prop.test(x = sum(forbes94$MBA == "1"), n = nrow(forbes94), p = 0.3,
+   alternative = "two.sided", conf.level = 0.99, correct = FALSE)

    1-sample proportions test without continuity
    correction

data:  sum(forbes94$MBA == "1") out of nrow(forbes94), null probability 0.3
X-squared = 5, df = 1, p-value = 0.03
alternative hypothesis: true p is not equal to 0.3
99 percent confidence interval:
 0.2257 0.3057
sample estimates:
     p 
0.2637 

Notate che per calcolare il numero di successi abbiamo usato il codice sum(forbes94$MBA == "1"), ovvero abbiamo prima verificato con l’operatore logico == per quali osservazioni abbiamo osservato un “successo” (cioè un valore pari a "1") e abbiamo calcolato la somma (sum()) del vettore logico risultante⁵².

I risultati indicano che l’intervallo di confidenza al 99% è pari a \((0.2257, 0.3057)\)⁵³, mentre il p-value del test, pari a 0.02526, ci suggerisce che esiste sufficiente evidenza empirica proveniente dai dati per rifiutare l’ipotesi nulla \(H_0: p = 0.30\) ad un livello di significatività \(\alpha\) di 0.05, ma non di 0.01.

6.2.3 Inferenza sul confronto tra le medie di due popolazioni normali

Il caso della differenza tra due medie è uno di quelli più utilizzati nelle applicazioni e nel corso ne sono state presentate alcune varianti, in particolare:

• confronto tra due medie di due popolazioni normali con varianze note, campioni indipendenti
• confronto tra due medie di due popolazioni normali con varianze non note ma assunte uguali, campioni indipendenti
• confronto tra due medie di due popolazioni normali con varianze non note, campioni dipendenti

In R non sono disponibili strumenti per il primo caso, che quindi non presenteremo qui.

Per confrontare due medie di due popolazioni normali in R possiamo utilizzare ancora la funzione t.test(), ma rispetto al caso di una singola media dovremo ora fornire i seguenti argomenti:

• x, valori osservati per il primo campione
• y, valori osservati per il secondo campione
• alternative, tipo di test che vogliamo effettuare, ovvero se bilaterale (two.sided) o unilaterale (less o greater, a seconda della direzione richiesta in \(H_1\))
• mu, valore della differenza delle medie che vogliamo testare sotto l’ipotesi nulla
• paired, che indica se i campioni sono indipendenti (FALSE) o dipendenti (TRUE)
• var.equal, che indica se le varianze (non note) delle due popolazioni nel caso di campioni indipendenti si assumono uguali (TRUE) o diverse (FALSE)
• conf.level, livello di confidenza che ci interessa usare

6.2.3.1 Campioni indipendenti

Vediamo subito un esempio di confronto tra due medie di due popolazioni normali con varianze non note ma uguali e campioni indipendenti. In particolare, confrontiamo la media della variabile ROS nei due gruppi identificati dalla variabile MasterPhd, che indica quali dei CEO nel 1994 possedeva un titolo di Master o di PhD. Calcoliamo l’intervallo di confidenza al 95% e effettuiamo il test per verificare se le due medie nella popolazione sono uguali, ovvero

\[\begin{equation*} H_0: \mu_{(\mbox{MasterPhd}=0)} = \mu_{(\mbox{MasterPhd}=1)} \quad \mbox{vs.} \quad H_1: \mu_{(\mbox{MasterPhd}=0)} \ne \mu_{(\mbox{MasterPhd}=1)}. \end{equation*}\]

Il codice seguente consente di ottenere queste analisi:

> ROS_0 <- forbes94$ROS[forbes94$MasterPhd == "0"]
> ROS_1 <- forbes94$ROS[forbes94$MasterPhd == "1"]
> t.test(ROS_0, ROS_1, alternative = "two.sided", mu = 0, paired = FALSE,
+   var.equal = TRUE, conf.level = 0.95)

    Two Sample t-test

data:  ROS_0 and ROS_1
t = -1.8, df = 790, p-value = 0.08
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -2.2948  0.1289
sample estimates:
mean of x mean of y 
    6.741     7.824 

La media campionaria nel secondo campione risulta essere maggiore, ma i risultati del test suggeriscono che non sembrano esserci differenze significative tra le medie del ROS nelle due popolazioni di aziende, quelle il cui CEO non ha un Master o un Phd e quelle in cui il CEO ha un Master o un Phd (almeno per i valori standard 1% e 5% del livello di significatività).

6.2.3.2 Campioni dipendenti

In R possiamo affrontare il caso di campioni dipendenti ancora con la funzione t.test(), salvo che questa volta gli argomenti x e y devono includere lo stesso numero di osservazioni e l’argomento paired deve essere posto uguale a TRUE.

Vediamo un esempio usando i dati contenuti nel data frame supermarket già discusso nella Sezione 4.3.2. Dopo aver caricato i dati, creiamo due vettori che contengono rispettivamente i dati relativi ai due campioni, con l’accortezza che il primo campione riguarderà il giorno in cui la promozione era attiva e il secondo quello in non lo era. Successivamente, con la funzione t.test() calcoliamo l’intervallo di confidenza al 90% per la differenza delle medie ed eseguiamo il test

\[\begin{equation*} H_0: \mu_X - \mu_Y \le 0 \qquad \mbox{vs.} \qquad H_1: \mu_X - \mu_Y > 0, \end{equation*}\]

ovvero ci interessa verificare se i dati suggeriscono che la promozione ha permesso di aumentare il numero medio di visite ai negozi⁵⁴. Il codice seguente consente di effettuare queste analisi:

> load("supermarket.RData")
> program_yes <- supermarket$customers[supermarket$program == "yes"]
> program_no <- supermarket$customers[supermarket$program == "no"]
> t.test(x = program_yes, y = program_no, alternative = "greater",
+   paired = TRUE, conf.level = 0.90)

    Paired t-test

data:  program_yes and program_no
t = 2.1, df = 9, p-value = 0.03
alternative hypothesis: true difference in means is greater than 0
90 percent confidence interval:
 1.022   Inf
sample estimates:
mean of the differences 
                      3 

Notate che abbiamo fissato alternative = "greater" perché ora ci interessa un test unilaterale. I risultati indicano che, assumendo di usare un livello di significatività del 5%, la promozione sembra avere avuto un effetto significativo sulla media del numero di visite poiché il p-value del test (0.03266) è inferiore a 0.05.

6.2.4 Inferenza sull’indice di correlazione lineare

In analogia con altre funzioni presentate finora (ovvero t.test() e prop.test()), R mette a disposizione la funzione cor.test() per fare inferenza sull’indice di correlazione lineare \(\rho\) di una popolazione normale bivariata. Come per la funzione cor(), anche questa funzione richiede di specificare le due variabili di cui vogliamo valutare l’associazione, insieme agli argomenti alternative e conf.level comuni alle altre funzioni sopra citate. A titolo di esempio, riprendiamo il data frame forbes94 ed effettuiamo un test per valutare l’assenza di associazione lineare tra le variabili Salary e Age, ovvero

\[\begin{equation*} H_0: \rho_{(\mbox{Salary}, \mbox{Age})} = 0 \quad \mbox{vs.} \quad H_1: \rho_{(\mbox{Salary}, \mbox{Age})} \ne 0. \end{equation*}\]

Insieme al diagramma di dispersione delle due variabili, il codice seguente fornisce i risultati relativi a tale test:

> plot(x = forbes94$Age, y = forbes94$Salary, xlab = "Age", ylab = "Salary")

Figura 6.14: Diagramma di dispersione per la variabile Salary contro Age.

> cor.test(x = forbes94$Salary, y = forbes94$Age,
+   alternative = "two.sided", conf.level = 0.95)

    Pearson's product-moment correlation

data:  forbes94$Salary and forbes94$Age
t = 6.7, df = 790, p-value = 4e-11
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.1654 0.2975
sample estimates:
   cor 
0.2325 

L’indice di correlazione campionario tra le due variabili, pari a 0.2353, segnala un’associazione lineare debole, ma il basso valore del p-value (3.788e-11) evidenzia che tale associazione è altamente significativa.

---

51. Ricordiamo che tale notazione rappresenta un modo alternativo per indicare numeri molto vicini a zero, in questo caso \(3.377 \times 10^{-5} = 0.00003377\).↩

52. Ricordate che R converte internamente il valore logico FALSE in 0 e il valore logico TRUE in 1.↩

53. Questi valori non corrispondono esattamente al calcolo che abbiamo presentato nel corso, perché R internamente usa una formula leggermente diversa, ma più precisa. Ricordatevi infatti che le espressioni che abbiamo visto in questo caso sono delle approssimazioni basate sulla distribuzione normale.↩

54. Nel test \(X\) indica la popolazione di negozi in cui la promozione è attiva, mentre \(Y\) denota la popolazione di negozi in cui la promozione non è attiva.↩