R e il mondo aleatorio
- Parte II Statistiche in R
Paola Lecca, CIBIO - UNITN
Corso di Matematica e Stastistica 2
Testo di riferimento
S. M. Iacus, G. Masarotto, Laboratorio di statistica con R, McGraw-Hill.
Verifica di ipotesi sulle proporzioni
Esempio: supponiamo che un sedicente indovino affermi di leggere nel pensiero
delle persone. Lo volgiamo mettere alla prova con il seguente esperimento.
Prendiamo un mazzo di 40 carte, 20 rosse 20 nere e ci sistemiamo in una stanza ben
separata da quella dove si trova l’indovino.
Estraiamo in successione in successione una carta, guardiamo di colore è, è, prima
di rimetterla nel mazzo, chiediamo all’indovino di indovinare il colore della carta
estratta.
Ci chiediamo poi se alla luce di questi esperimenti possiamo affermare che
l’indivino sia realmente dotato di poteri paranormali.
Verifica di ipotesi sulle proporzioni
Se l’indovino fosse realmente dotato di questi potere, la propabilità che indivini il
colore della carta sarebbe uguale a 1.
Appelliamoci a un principio di “buona fede” per cui ammettiamo ch eanche nelle
migliori condizioni, a causa di possibili disturbi/interferenze con i suoi presunti poteri,
l’indovino non abbia probabilità uguale a 1 di indovinare il colore della carta estratta.
Poniamoci allora la seguente domanda:
l’indovino indovina per puro caso?
Impostiamo dunque un test statistico come segue.
L’ipotesi nulla è che l’indovino indovini il colore della carta a caso, quindi che la
probabilità di indovinare sia 50%, così come quella di sbagliare il colore della carta
estratta.
Verifica delle ipotesi sulle proporzioni
Supponimao che vengano eseguite 50 estrazioni e che l’indovino abbia individuato 32
estrazioni su 50. Siamo quindi difronte a uno schema binomiale con n=50 e p uguale
alla probabilità che l’indovino indovini una carta.
In R il test si esegue così:
binom.test(32, 50, alternative=“greater”)
Inputs della funzione binom.test
binom.test {stats}
R Documentation
Exact Binomial Test
Description
Performs an exact test of a simple null hypothesis about the probability of success in a Bernoulli
experiment.
Usage
binom.test(x, n, p = 0.5, alternative = c("two.sided", "less", "greater"), conf.level = 0.95)
Arguments
x
number of successes, or a vector of length 2
giving the numbers of successes and failures,
respectively.
n
number of trials; ignored if x has length 2.
p
hypothesized probability of success.
alternative
indicates the alternative hypothesis and must
be one of "two.sided", "greater" or "less". You
can specify just the initial letter.
conf.level
confidence level for the returned confidence
interval.
Outputs della funzione binom.text
Value
A list with class "htest" containing the following components:
statistic
the number of successes.
parameter
the number of trials.
p.value
the p-value of the test.
conf.int
a confidence interval for the
probability of success.
estimate
the estimated probability of
success.
null.value
the probability of success
under the null, p.
alternative
a character string describing
the alternative hypothesis.
method
the character string "Exact
binomial test".
data.name
a character string giving the
names of the data.
binom.test in R: esempi
Circa uguale al valore che otteniamo con
> binom.test(32, 50, alternative="greater")
Z <- (32/50, -0.5/sqrt(0.5*0.5/50)
pnorm(Z, lower.tail=FALSE)s
Exact binomial test
data: 32 and 50
number of successes = 32, number of trials = 50, p-value = 0.03245
alternative hypothesis: true probability of success is greater than 0.5
95 percent confidence interval:
0.5142308 1.0000000
sample estimates:
probability of success
0.64
Nota: Il livello di confideva è posto uguale al 95%
E con prop.test?
> prop.test(32, 50, alternative="greater")
1-sample proportions test with continuity correction
data: 32 out of 50, null probability 0.5
X-squared = 3.38, df = 1, p-value = 0.033
alternative hypothesis: true p is greater than 0.5
95 percent confidence interval:
0.5137352 1.0000000
sample estimates:
p
0.64
binom.test in R: esempi
> binom.test(32, 50, alternative="greater", conf.level=0.01)
Exact binomial test
data: 32 and 50
number of successes = 32, number of trials = 50, p-value = 0.03245
alternative hypothesis: true probability of success is greater than 0.5
1 percent confidence interval:
0.7742634 1.0000000
sample estimates:
probability of success
0.64
In questo caso il test non ha ragione di rifiutare l’ipotesi nulla. Spiegate perchè.
L’enigma dell’indovino
Come possiamo risolvere l’enigma sulle presunte capacità paranormali
dell’indovino?
Provate a rispondere a questa domanda giustificando la vostra risposta:
avere un campione maggiore di estrazioni (prove ripetute) potrebbefornirci una
statistica più robusta sulla base della quale poter ragionare?
Verifica di ipotesi sulla varianza(1/5)
Se abbiamo un campione di n dati di tipo gaussiano con media
e varianza
allora
CASO 1. Consideriamo le seguenti ipotesi
Se pensiamo che il valore di
sia maggiore di un certo valore
costruire un test che rifiuta H0 se
allora possiamo
Verifica di ipotesi sulla varianza (2/5)
Il valore c sidetermina in relazione all’errore del I tipo che è definito come segue
ovvero
il che implica
dove
è il quantile di ordine
gradi di libertà.
della distribuzione Chi-quadro con n - 1
Verifica di ipotesi sulla varianza (3/5)
Quindi il test di livello a è del tipo:
Rifiuto
in favore di
se
cioè quando
(F1)
Verifica di ipotesi sulla varianza (4/5)
CASO 2. Consideriamo le seguenti ipotesi
L’errore di I tipo si calcola come segue
Allora il test è del tipo: rifiutiamo H0 in favore di H1 se
(F2)
Verifica di ipotesi sulla varianza (5/5)
Eserizio: scrivere un codice in R per la verifica di ipotesi sulla varianza.
Suggerimenti:
a) si implementino le formule (F1) ed (F2) viste nelle slide precedenti
b) utilizzare la funzione qchisq (eseguite help(qchisq) per accedere all’aiuto).
RISOLVETE L’ESERCIZIO SCRIVENDO UN CODICE IN R.
INVIATE ENTRO IL 21 MAGGIO 2013 IL CODICE ALL’INSEGNANTE PER E-MAIL
ALLEGATO A UN BREVE TESTO DI SPIEGAZIONE.
Verifica di ipotesi su due
campioni
Verifica di ipotesi su due proporzioni
Farmaco
Infartuati
Non
infartuati
Totali
Placebo
239
10795
11034
Aspirina
139
10686
11037
Si vuole accertare se l’aspirina ha effetti positivi sulla prevenzione dell’infarto.
In R possiamo usare la funzione prop.test, come segue.
# Per un test a due code
prop.test( c(239, 139), c(11034, 11037) )
# Per un test a una coda
prop.test( c(239, 139), c(11034, 11037) , alternative=“greater”)
Confronto tra le medie di gruppi (1/3)
Si construisce un test t per verificare l’uguaglianza delle medie come segue
(F3)
dove n1 ed n2 sono le dimensioni dei due campioni ed
con
ed
uguali alle varianze campionarie dei due campioni.
Questa statistica test t si distribuisce come una t di Studnet con n1 + n2 -2 gradi di libertà.
Confronto tra le medie di gruppi (2/3)
Confronto tra le medie di gruppi (3/3)
Esercizio da svolgere in R.
Su due campioni di autovetture guidate da un gruppo di uomini e da un gruppo di
donne, sono state calcolate le seguenti spese annuali:
- per i gruppo di n1=5 uomini si è avuto 540 Euro in media
- per il gruppo di n2=7 donne si è avuto 300 Euro in media.
Supponiamo a = 5% e
Suggerimento: implementate la formula (F3). Vi sarà poi utile usare la funzione qt
(che significa quantile) come segue
qt(0.95, df=10)
Confronto tra le medie di gruppi
Esercizio da svolgere in R.
Su due campioni di autovetture guidate da un gruppo di uomini e da un gruppo di
donne, sono state calcolate le seguenti spese annuali:
- per i gruppo di n1=5 uomini si è avuto 540 Euro in media
- per il gruppo di n2=7 donne si è avuto 300 Euro in media.
Le varianze campionarie sono
Euro ripettivamente.
Supponiamo a = 5%.
Suggerimento: implementate la formula (F3). Vi sarà poi utile usare la funzione qt
(che significa quantile) come segue
qt(0.95, df=10)
SI tratta di un test t. Quindi calcoliamo le quantità in gioco.
Quindi:
Se vogliamo testare l’ipotesi alternativa H1: m1 > m2, dobbiamo calcolare il valore di
soglia
Se a = 5%, qt(0.95, df=19) in R dà 1.81. Siccome t < 1.81 non rifiutimao l’ipotesi nulla
in favore del fatto che i guidatori uominispendono di più in media ripestto alla donne
per la manutenzione delle automobili.
Per utilizzare l a funzione t.test di R in questo caso dobbiamo disporre di due insiemi di
dati su cui effettuare il test.
Per semplicità simuliamo due campioni casuali estratti dalla normale ipotizzando come
media e varianza quelle del nostro esempio, come segue
> m <- rnorm(5, mean=540, sd=299)
> f <- rnomr(7, mean=300, sd=238)
> t.test(m, f, alternative="greater")
Welch Two Sample t-test
data: m and f
t = 5.7483, df = 9.941, p-value = 9.493e-05
alternative hypothesis: true difference in means is greater than 0
95 percent confidence interval:
437.463 Inf
sample estimates:
mean of x mean of y
801.2664 162.1743