#install.packages('Rcmdr', dependencies=TRUE)
#install.packages('Biodem',dependencies=TRUE )

#library(Rcmdr)
#library(Biodem)
library(MASS)

T- Test

BodyTemperature<- read.table('BodyTemperature.txt', 
           sep=' ',
           header=TRUE)

Il dataset BodyTemperature è composto da:

per il momento concentriamoci sulla variable Temperature che riguarda la temperatura corporea misurata per un dato paziente (ogni riga del dataset rappresenta un paziente). Supponiamo di voler testare l’ipotesi che la temperatura della popolazione sia diversa da \(98.6\). Questo si traduce nel seguente set di ipotesi. ( Ricordiamo che l’ipotesi alternativa è sempre un ipotesi di ricerca )

\[ H_0: \mu=98.6 \,\, \text{vs} \,\, H_A: \mu \ne 98.6 \]

per ora supponiamo che \(\sigma^2=1\) (ovvero che la varianza della popolazione sia nota e uguale a 1). Inolre supponiamo di aver selezionato un campione di \(100\) individui dalla popolazione. Sappiamo che \[\bar{X} \sim N(\mu, \sigma^2/n)\]

Dove \('\sim N'\) distribuito come una variabile casuale Normale. Supponendo che l’ipotesi nulla sia vera (ovvero che la popolazione abbia una media di \(98.6\)) la distribuzione diventa

\[ \bar{X}| H_o\sim N(98.6, 1/100) \]

Dove \(\sigma^2=1\) e \(n=100\). Con la notazione \(\bar{X}|H_0\) si è voluto stressare l’idea che questa è la distribuzione della media campionaria dato che ( così si legge il simbolo ‘\(|\)’) \(H_0\) è vera.

Sappiamo che la Statistica Test (ST) per il test sulla media sotto \(H_0\) è

\[ z=\frac{\bar{x}-\mu_0}{s/\sqrt{n}}\sim t_{n-1} \]

e che un intervallo \(100(1-\alpha)\%\) per la media ha la seguente forma

\[ \bar{x}\pm t_{\alpha/2}\frac{s}{\sqrt{n}} \]

La funzione t.test, se correttamente utilizzata, resituisce sia il test che un intervallo di confidenza.

t.test(BodyTemperature$Temperature, # variabile che vogliamo testare
       mu=98.6, # ipotesi nulla
       alternative = "two.sided") # ipotesi a due code
## 
##  One Sample t-test
## 
## data:  BodyTemperature$Temperature
## t = -2.8216, df = 99, p-value = 0.005775
## alternative hypothesis: true mean is not equal to 98.6
## 95 percent confidence interval:
##  98.14013 98.51987
## sample estimates:
## mean of x 
##     98.33

L’output del comando t.test per prima cosa ci dice qual’è l’potesi alternativa che stiamo specificando alternative hypothesis: true mean is not equal to 98.6, questo toglie ogni dubbio. Stiamo effettivamente testando in seguente set di ipotesi

\[ H_0: \mu=98.6 \,\, \text{vs} \,\, H_A: \mu \ne 98.6 \]

Inoltre, restituisce il t-score (\(t=-2.8216\)), i gradi di libertà (degree of freedom) (\(\nu=(n-1)=99\)) ed il p-value \(0.005775\). Ci da una stima della media \(\bar{x}=98.33\) e un intervallo di confidenza al \(95\%\) \([98.14013; 98.51987]\).

Con la funzione t.test possiamo settare diverse opzioni. Immaginiamo di voler testare il seguente set di ipotesi

\[ H_0: \mu\ = 98.6 \,\, \text{vs} \,\, H_A: \mu > 98.6 \]

Qui la nostra ipotesi di ricerca è la seguente: capire se, in base al nostro campione, è ragionevole suppore che la temperatura della popolazione sia maggiore di \(98.6\). In t.test questo è possibile farlo settando alternative = “greater”.

t.test(BodyTemperature$Temperature, #  variabile che vogliamo testare
       mu=98.6, # ipotesi nulla
       alternative = "greater") # ipotesi a una sola coda
## 
##  One Sample t-test
## 
## data:  BodyTemperature$Temperature
## t = -2.8216, df = 99, p-value = 0.9971
## alternative hypothesis: true mean is greater than 98.6
## 95 percent confidence interval:
##  98.17112      Inf
## sample estimates:
## mean of x 
##     98.33

La dicitura alternative hypothesis: true mean is greater than 98.6 conferma che quello che stiamo lanciando in R è esattamente il set di potesi che vogliamo testare.

Nel caso di un set di potesi del tipo

\[ H_0: \mu = 98.6 \,\, \text{vs} \,\, H_A: \mu < 98.6 \]

t.test(BodyTemperature$Temperature, # La variabile che vogliamo testare
       mu=98.6, # ipotesi nulla
       alternative = "less") # ipotesi a una sola coda
## 
##  One Sample t-test
## 
## data:  BodyTemperature$Temperature
## t = -2.8216, df = 99, p-value = 0.002887
## alternative hypothesis: true mean is less than 98.6
## 95 percent confidence interval:
##      -Inf 98.48888
## sample estimates:
## mean of x 
##     98.33

Possiamo ottenere diverse ampiezze per il nostro intervallo di confidenza ( invece del \(95\%\) usato di default). Basta modificare l’opzione conf.level:

t.test(BodyTemperature$Temperature, # variabile che vogliamo testare
       mu=98.6, # ipotesi nulla
       alternative = "two.sided", # tipo di test
       conf.level = 0.9) # intervallo di confidenza
## 
##  One Sample t-test
## 
## data:  BodyTemperature$Temperature
## t = -2.8216, df = 99, p-value = 0.005775
## alternative hypothesis: true mean is not equal to 98.6
## 90 percent confidence interval:
##  98.17112 98.48888
## sample estimates:
## mean of x 
##     98.33

Assunzione di Normalità

Nell’esempio precedente abbiamo considerato un campione molto grande \(n=100\), quindi anche qualora i dati non fossero distribuiti come una v.c. Normale possiamo usare il Teorema del Limite Centrale. Tuttavia, vale la pena verificare la Normalità (inteso come distribuzione Normale).

{qqnorm(BodyTemperature$Temperature)
qqline(BodyTemperature$Temperature, col='green')}

Generalmente il qqplot è uno strumento utile per capire se i nostri dati sono distribuiti come una v.c. Normale. Senza entrare nel dettaglio, per essere Normali, dobbiamo osservare che i punti siano distribuiti sulla linea \(x=y\) ovvero una retta a \(45\) gradi. Vediamo alcuni esempi di qqplot basati su dati non distribuiti come una v.c. Normale.

par(mfrow=c(1,2))
n <- 200

for(i in 1:2) {x <- rnorm(n); 
qqnorm(x, 
       main='Normal Data'); 
qqline(x,, col='green')}

for(i in 1:2) {x <- exp(rnorm(n)); 
qqnorm(x, 
       main='Dati NON Normali');
qqline(x,, col='green')} # DISTRIBUZIONE NON NORMALE 

for(i in 1:2) {x <- rcauchy(n); 
qqnorm(x, 
       main='Dati NON Normali'); 
qqline(x,, col='green')} # DISTRIBUZIONE NON NORMALE 

for(i in 1:2) {x <- runif(n); 
qqnorm(x, 
       main='Dati NON Normali'); 
qqline(x,, col='green')} # DISTRIBUZIONE NON NORMALE

Non vogliamo entrare nel dettaglio qui, lo scopo è semplicemente studiare la forma dei qqplot quando i dati non sono distribuiti come una v.c. Normale.

Quando il sample size è piccolo (< 30) bisogna prestare molta attenzione. Proviamo a ripercorrere lo stesso esempio con un sample size di \(n=25\)

par(mfrow=c(1,3))
n <- 25

for(i in 1:6) {x <- rnorm(n); qqnorm(x); qqline(x,, col='green')}

par(mfrow=c(1,3))
for(i in 1:6) {x <- exp(rnorm(n)); 
qqnorm(x, 
       main='Dati NON Normali'); 
qqline(x,, col='green')} # DISTRIBUZIONE NON NORMALE 

par(mfrow=c(1,3))
for(i in 1:6) {x <- rcauchy(n); 
qqnorm(x,  
       main='Dati NON Normali'); 
qqline(x,, col='green')} # DISTRIBUZIONE NON NORMALE 

par(mfrow=c(1,3))
for(i in 1:6) {x <- runif(n);
qqnorm(x,  
       main='Dati NON Normali'); 
qqline(x,, col='green')} # DISTRIBUZIONE NON NORMALE

Take Home Message: Prestate sempre attenzione, mai basare l’analisi su un singolo strumento. In statistica bisogna ragionare sulla natura dei dati. Mai applicare una procedura senza pensare a cosa si sta facendo.

Test di Normalità: Shapiro Wilk

Esiste un metodo molto utile che consente di testare l’ipotesi di normalità dei dati, lo Shapiro-Wilk (S-W) test. Il set di ipotesi da testare è il seguente

\[ H_0: \text{Dati sono Normali}\,\,\, \text{vs}\,\,\, H_A: \text{Dati NON Normali} \] in R è tutto molto semplice

shapiro.test(BodyTemperature$Temperature)
## 
##  Shapiro-Wilk normality test
## 
## data:  BodyTemperature$Temperature
## W = 0.98171, p-value = 0.1803

Dove \(W\) è la statistica test, mentre p-value è la probabilità di osservare un valore più estremo della statistica test calcolata dato che \(H_0\) è vera. In questo caso notiamo che il p-value è grande, quindi non rigettiamo la nulla. Il risultato del test insieme al qq-plot ci pone nella posizione di concludere che è ragionevole considerare, i dati testati, come provenienti da una v.c. Normale.

Notiamo che in genere il Test di S-W viene utilizzato quando il nostro campione è basso (<50).

set.seed(23)
x <- rcauchy(25) # generiamo 25 osservazioni non normali
shapiro.test(x)
## 
##  Shapiro-Wilk normality test
## 
## data:  x
## W = 0.67341, p-value = 3.263e-06

Infatti, rifiutiamo la nulla. FUNZIONA !!!!

Test sulla media per distribuzioni non normali: Wilcoxon Test

Quando non si lavora con dati NON normali, entra in gioco la statistica non parametrica. Non entreremo nei dettagli tecnici, in \(R\) tutto è molto semplice.

dta <- read.csv("insurance.csv", header=T) # carichiamo alcuni dati
str(dta) # verifichiamo la struttura 
## 'data.frame':    1338 obs. of  7 variables:
##  $ age     : int  19 18 28 33 32 31 46 37 37 60 ...
##  $ sex     : chr  "female" "male" "male" "male" ...
##  $ bmi     : num  27.9 33.8 33 22.7 28.9 ...
##  $ children: int  0 1 3 0 0 0 1 3 2 0 ...
##  $ smoker  : chr  "yes" "no" "no" "no" ...
##  $ region  : chr  "southwest" "southeast" "southeast" "northwest" ...
##  $ charges : num  16885 1726 4449 21984 3867 ...

In questo caso lavoriamo con la variable charges nel dataset insurance. Controlliamo velocemente l’ assunzione di normalità tramite il qqplot.

{qqnorm(dta$charges)
qqline(dta$charges, col='green')}

Qui i dati non sembrano distribuiti come una v.c. Normale, non serve andare oltre ( effettuare un test ). A questo punto, estraiamo un campione di \(n=200\) e verifichiamo se la media di changes è uguale a \(13270.42\). Ovviamente NON possiamo usare il comando t.test ( i dati non sono nornali ). Dobbaimo effettuare il test tramite tecniche non parametriche. utilizzeremo il Wilcoxon Test.

N <- nrow(dta) # Numerosità della popolazione
n <- 200        # Numerosità campionaria

lista <- sample(1:N, n, replace = F) 
              # Stiamo campionando n elementi 
              # da una popolazione di ampiezza 
              # N senza reinserire il capione estratto.

campione <- dta[lista,]

mean(campione$charges) 
## [1] 13323.05
              # media del campione
              #Consiglio: Lanciate questa 
              #chunck di codice più volte... 
              #notate come il valore della media 
              #campionaria cambia. Sapreste spiegare il perchè ?

A questo punto siamo pronti per il test

wilcox.test(dta$charges, # Variabile da testare
            mu=13270.42, # Ipotesi nulla
            alternative = "two.sided") # Ipotesi alternativa
## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  dta$charges
## V = 346290, p-value = 6.599e-13
## alternative hypothesis: true location is not equal to 13270.42
wilcox.test(dta$charges, mu=13270.42, alternative = "less")
## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  dta$charges
## V = 346290, p-value = 3.299e-13
## alternative hypothesis: true location is less than 13270.42
wilcox.test(dta$charges, mu=13270.42, alternative = "two.sided")
## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  dta$charges
## V = 346290, p-value = 6.599e-13
## alternative hypothesis: true location is not equal to 13270.42

Consiglio di esplorare la funzione ?wilcox.test.

Ricordiamo che il test di Wilcoxon è ideato per lavorare con solo dati continui. Inoltre, in caso di campioni molto grandi (\(n\) molto grande) abbiamo sempre il supporto del Teorema del Limite Centrale. Dunque, il test di Wilcoxon potrebbe non essere necessario.

Test Confronto Medie

Quando si parla di test per il confronto tra medie, bisogna fare tre distinzioni.

Test confronto medie: Dati Normali e varianze note

Nel caso (1) la statistica test risulta

\[ ST=\frac{(\bar{x}_1-\bar{x}_2)- (\mu_1-\mu_2)}{\sqrt{\frac{\sigma^2_1}{n_1}+\frac{\sigma^2_2}{n_2}}} \]

dove \(\bar{x}_1\) e \(\sigma^2_1\) è la media e varianza del primo gruppo, \(\bar{x}_2\) è a media e varianza del secondo gruppo.

Test confronto medie: Dati Normali e varianze NON note e uguali

  1. E’ il caso che più frequente ( molto improbabile conoscere la varianza di una popolazione ). La statistica test risulta essere

\[ ST=\frac{(\bar{x}_1-\bar{x}_2)- (\mu_1-\mu_2)}{\sqrt{\frac{s^2_p}{n_1}+\frac{s^2_p}{n_2}}} \] dove \(s^2_p\) è una media ponderata delle due varianze campionarie

\[ s_p^2=\frac{(n_1-1)s^2_1+(n_2-1)s^2_2}{n_1+n_2-2} \]

\[ s^2_1=\sum_{i=1}^{n_1}(x_i-\bar{x}_1)^2/(n_1-1) \] \[ s^2_2=\sum_{i=1}^{n_2}(x_i-\bar{x}_2)^2/(n_2-1) \] ## Test confronto medie: Dati Normali e varianze NON note e NON uguali

Nel caso di Varianze non uguali, la statistica test è

\[ ST=\frac{(\bar{x}_1-\bar{x}_2)- (\mu_1-\mu_2)}{\sqrt{\frac{s^2_1}{n_1}+\frac{s^2_2}{n_2}}} \] dove

\[ s^2_1=\sum_{i=1}^{n_1}(x_i-\bar{x}_1)^2/(n_1-1) \\ s^2_2=\sum_{i=1}^{n_2}(x_i-\bar{x}_2)^2/(n_2-1) \]

qui il problema è calcolare il giusto valore della t-student. Il valore critico sarà

\[ t_{\alpha/2}=\frac{\omega_1 t_1+ \omega_2 t_2}{\omega_1+ \omega_2} \]

dove \(\omega_1=s_1^2/n_1\), \(\omega_2=s_2^2/n_2\), \(t_1=t_{\alpha/2}\) per \(n_1-1\) gradi di libertà e \(t_2=t_{\alpha/2}\)per \(n_2-1\) gradi di libertà.

NOTA: Questa opzione è integrata nella funzione t.test tramite l’opzione var.equal. Di default è impostata come FALSE.

Esempio

Immaginiamo di confrontare la temperatura corporea in due diversi gruppi (Maschi vs Donne). Senza entrare nel dettaglio, è possibile testare \[ H_O: \mu_{Maschi}=\mu_{Donne} \,\,\, \text{vs}\,\,\, \mu_{Maschi}\ne\mu_{Donne} \]

dove \(\mu_{Maschi}\) è la media della temperatura del gruppo dei maschi e \(\mu_{Donne}\) è la media della popolazione della temperatura del gruppo delle donne. Qui siamo in un contesto in cui, la varianza non si conosce e non c’è ragione per considerarla differente nei due gruppi (in base all’ informazione a priori). Ma come possiamo verificare euristicamente che le varianze sono uguali ?

BodyTemperature$Gender <- as.factor(BodyTemperature$Gender)

plot(y=BodyTemperature$Temperature, 
     x=BodyTemperature$Gender, 
     main='Temperatura per genere4343',
     xlab='Genere',
     ylab='Temperatura')

var(BodyTemperature$Temperature[BodyTemperature$Gender=='M'])
## [1] 1.103121
var(BodyTemperature$Temperature[BodyTemperature$Gender=='F'])
## [1] 0.720502
t.test(BodyTemperature$Temperature~BodyTemperature$Gender, alternative = 'two.sided', conf.level = .95)
## 
##  Welch Two Sample t-test
## 
## data:  BodyTemperature$Temperature by BodyTemperature$Gender
## t = 1.3526, df = 92.265, p-value = 0.1795
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -0.1212510  0.6390581
## sample estimates:
## mean in group F mean in group M 
##        98.45686        98.19796

Con la funzione t.test possiamo effettuare diversi test semplicemente cambiando alcuni comandi. Immaginiamo di voler testare l’ipotesi che la differenza sia maggiore di zero ( alternative = ‘greater’) oppure che sia minore di zero ( alternative = ‘less’).

t.test(BodyTemperature$Temperature~BodyTemperature$Gender, 
       alternative = 'greater', 
       conf.level = .95)
## 
##  Welch Two Sample t-test
## 
## data:  BodyTemperature$Temperature by BodyTemperature$Gender
## t = 1.3526, df = 92.265, p-value = 0.08975
## alternative hypothesis: true difference in means between group F and group M is greater than 0
## 95 percent confidence interval:
##  -0.05914149         Inf
## sample estimates:
## mean in group F mean in group M 
##        98.45686        98.19796
t.test(BodyTemperature$Temperature~BodyTemperature$Gender, 
       alternative = 'less', 
       conf.level = .95)
## 
##  Welch Two Sample t-test
## 
## data:  BodyTemperature$Temperature by BodyTemperature$Gender
## t = 1.3526, df = 92.265, p-value = 0.9103
## alternative hypothesis: true difference in means between group F and group M is less than 0
## 95 percent confidence interval:
##       -Inf 0.5769486
## sample estimates:
## mean in group F mean in group M 
##        98.45686        98.19796
#t.test(BodyTemperature$Temperature~BodyTemperature$Gender, 
# alternative =  'less', 
# paired=TRUE,
# conf.level = .95)

Test confronto medie (gruppi > 2) dati distribuiti come una v.c. Normale

Nel caso volessimo testare più di due gruppo dobbiamo spostaci sul comando anova, le opzioni sono simili a quelle viste con t.test ( niete paura !!!). Vediamo un esempio, useremo il dataset birtwt incorporato nel paccheto MASS.

head(birthwt, n=10) 
##    low age lwt race smoke ptl ht ui ftv  bwt
## 85   0  19 182    2     0   0  0  1   0 2523
## 86   0  33 155    3     0   0  0  0   3 2551
## 87   0  20 105    1     1   0  0  0   1 2557
## 88   0  21 108    1     1   0  0  1   2 2594
## 89   0  18 107    1     1   0  0  1   0 2600
## 91   0  21 124    3     0   0  0  0   0 2622
## 92   0  22 118    1     0   0  0  0   1 2637
## 93   0  17 103    3     0   0  0  0   1 2637
## 94   0  29 123    1     1   0  0  0   1 2663
## 95   0  26 113    1     1   0  0  0   0 2665
# dataset estratto dalla libreria MASS, assicuratevi di lanciare    # library(MASS) prima.

Questo dataset include il peso alla nascita (in grammi) di 189 neonati insieme ad alcune caratteristiche (ad esempio, età, stato di fumatore) delle loro madri. I dati sono stati raccolti presso il Baystate Medical Center, Springfield, MA, durante il 1986.

Le variabili incluse sono:

Notiamo però qualcosa di importante… low, race, smoke, ht, e ui sono categoriche (e.s., il sesso è categorico). Quindi prima di procedere dobbiamo formattare i dati. Formattare i dati è un passo fondamentale.

birthwt$low<- as.factor(birthwt$low)
birthwt$race<- as.factor(birthwt$race)
birthwt$smoke<- as.factor(birthwt$smoke)
birthwt$ht<- as.factor(birthwt$ht)
birthwt$ui<- as.factor(birthwt$ui)

Adesso siamo pronti per effettuare un test su medie, nel caso di più di due gruppi. In questo caso useremo race come gruppi

\[ H_0: \mu_1=\mu_2=\mu_3 \,\, \text{vs}\,\, H_A:\mu_1\ne \mu_2 \ne \mu_3 \]

ANOVA_test<- aov(birthwt$bwt~birthwt$race)
summary(ANOVA_test)
##               Df   Sum Sq Mean Sq F value  Pr(>F)   
## birthwt$race   2  5015725 2507863   4.913 0.00834 **
## Residuals    186 94953931  510505                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Non accettiamo l’ipotesi nulla ad un livello \(\alpha=0.01\). Nel caso volessimo vedere graficamente ( tramite boxplots ) le tre distribuzioni del peso ?

plot(birthwt$bwt~birthwt$race)

Test su confronto tra le medie per dati Non distribuiti come una v.c. Normale

Anche in questo caso, esiste un modo per effettuare il test per confronto di medie con dati non Normali. Visualizziamo la variabile charges per maschi e donne.

boxplot(campione$charges ~  campione$sex,
        main='Boxplot di charges per genere', 
        xlab='Genere',
        ylab='charges')

wilcox.test(campione$charges ~  campione$sex, alternative = "greater")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  campione$charges by campione$sex
## W = 4487, p-value = 0.8811
## alternative hypothesis: true location shift is greater than 0
wilcox.test(campione$charges ~  campione$sex, alternative = "less")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  campione$charges by campione$sex
## W = 4487, p-value = 0.1194
## alternative hypothesis: true location shift is less than 0
wilcox.test(campione$charges ~  campione$sex, alternative = "two.sided")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  campione$charges by campione$sex
## W = 4487, p-value = 0.2389
## alternative hypothesis: true location shift is not equal to 0

Test su confronto tra le medie per distribuzioni Normali e Non Normali su Campioni Dipendendi

Potrebbe capitare di avere a che fare con campioni dipendenti (es. misuriamo due caratteristiche dello stesso individuo in due istanti diversi). In questo caso dobbiamo modificare leggermente la statistica test

\[ ST=\frac{\bar{d}-\mu_{d_0}}{s_{\bar{d}}} \]

dove \(\bar{d}=n^{-1}\sum^n d_j\) è la media campionaria delle differenze e dove \(\{d_j=x_{1j}-x_{2j}\}\). \(\mu_{d_0}\) è la media delle differenze della popolazione, \(s_{\bar{d}}=s_{d}/\sqrt{n}\), \(n\) è la numerosità campionaria delle differenze ( stiamo assumendo che \(n_1=n_2=n\) !!!), \(s_{\bar{d}}\) è la deviazione standard delle differenze campionarie.

Quando \(H_0\) è vera \(ST\) si distribuisce come una v.c. T-student con \(n-1\) gradi di libertà.

E’ tutto implementato nella funzione t.test e wilcox.test ( nel caso di dati non Normali).

# install.packages('PairedData') 

library("PairedData") # Assicuratevi di aver istallato PairedData
## Loading required package: gld
## Loading required package: mvtnorm
## Loading required package: lattice
## Loading required package: ggplot2
## 
## Attaching package: 'PairedData'
## The following object is masked from 'package:base':
## 
##     summary
data("anorexia")


head(anorexia)
##   Treat Prewt Postwt
## 1  Cont  80.7   80.2
## 2  Cont  89.4   80.1
## 3  Cont  91.8   86.4
## 4  Cont  74.0   86.3
## 5  Cont  78.1   76.1
## 6  Cont  88.3   78.1
t.test(anorexia$Prewt, anorexia$Postwt, paired = T, alternative = "greater")
## 
##  Paired t-test
## 
## data:  anorexia$Prewt and anorexia$Postwt
## t = -2.9376, df = 71, p-value = 0.9978
## alternative hypothesis: true mean difference is greater than 0
## 95 percent confidence interval:
##  -4.331953       Inf
## sample estimates:
## mean difference 
##       -2.763889
t.test(anorexia$Prewt, anorexia$Postwt, paired = T, alternative = "less")
## 
##  Paired t-test
## 
## data:  anorexia$Prewt and anorexia$Postwt
## t = -2.9376, df = 71, p-value = 0.002229
## alternative hypothesis: true mean difference is less than 0
## 95 percent confidence interval:
##       -Inf -1.195825
## sample estimates:
## mean difference 
##       -2.763889
t.test(anorexia$Prewt, anorexia$Postwt, paired = T, alternative = "two.sided")
## 
##  Paired t-test
## 
## data:  anorexia$Prewt and anorexia$Postwt
## t = -2.9376, df = 71, p-value = 0.004458
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -4.6399424 -0.8878354
## sample estimates:
## mean difference 
##       -2.763889
wilcox.test(anorexia$Prewt, anorexia$Postwt, paired = T, alternative = "greater")
## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  anorexia$Prewt and anorexia$Postwt
## V = 831.5, p-value = 0.9948
## alternative hypothesis: true location shift is greater than 0
wilcox.test(anorexia$Prewt, anorexia$Postwt, paired = T, alternative = "less")
## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  anorexia$Prewt and anorexia$Postwt
## V = 831.5, p-value = 0.005301
## alternative hypothesis: true location shift is less than 0
wilcox.test(anorexia$Prewt, anorexia$Postwt, paired = T, alternative = "two.sided")
## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  anorexia$Prewt and anorexia$Postwt
## V = 831.5, p-value = 0.0106
## alternative hypothesis: true location shift is not equal to 0

Dove Prewt è il peso del paziente prima dell’ inizio dello studio e Postwt è il peso del paziente dopo lo studio. Lanciate ?anorexia nel caso voleste sapere di più.

Test su Proporzioni

Nel caso di proporzioni la statistica test è

\[ ST=\frac{\hat{p}-p_0}{\sqrt{p_0(1-p_0)/n}} \]

dove \(\hat{p}\) è la stima della proporzione, \(p_0\) è la proporzione testata sotto \(H_0\) e \(n\) è la numerosità campionaria. Quando \(H_0\) è vera, la Statistica test è distribuita approssimativamente come una Normale Standard.

dta <- read.csv("insurance.csv", header=T)
#str(dati)
head(dta, n=5)
##   age    sex    bmi children smoker    region   charges
## 1  19 female 27.900        0    yes southwest 16884.924
## 2  18   male 33.770        1     no southeast  1725.552
## 3  28   male 33.000        3     no southeast  4449.462
## 4  33   male 22.705        0     no northwest 21984.471
## 5  32   male 28.880        0     no northwest  3866.855
n=length(dta$smoker)

Tot_smokers <- sum(dta$smoker=="yes")



cat(p.hat=Tot_smokers/n)
## 0.2047833
binom.test(Tot_smokers,# Numero di successi ( numero di fumatori )
           n, # Ampiezza campionaria 
           p=0.5, # Ipotesi nulla
           alternative = "two.sided")
## 
##  Exact binomial test
## 
## data:  Tot_smokers and n
## number of successes = 274, number of trials = 1338, p-value < 2.2e-16
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.1834540 0.2274173
## sample estimates:
## probability of success 
##              0.2047833
binom.test(Tot_smokers, 
           n, 
           p=0.5, 
           alternative = "less")
## 
##  Exact binomial test
## 
## data:  Tot_smokers and n
## number of successes = 274, number of trials = 1338, p-value < 2.2e-16
## alternative hypothesis: true probability of success is less than 0.5
## 95 percent confidence interval:
##  0.0000000 0.2237773
## sample estimates:
## probability of success 
##              0.2047833
binom.test(Tot_smokers, 
           n, 
           p=0.5,
           alternative = "two.sided")
## 
##  Exact binomial test
## 
## data:  Tot_smokers and n
## number of successes = 274, number of trials = 1338, p-value < 2.2e-16
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.1834540 0.2274173
## sample estimates:
## probability of success 
##              0.2047833
#se n<30, utilizzare prop.test

prop.test(Tot_smokers, n, p=0.5, alternative = "greater")
## 
##  1-sample proportions test with continuity correction
## 
## data:  Tot_smokers out of n, null probability 0.5
## X-squared = 465.26, df = 1, p-value = 1
## alternative hypothesis: true p is greater than 0.5
## 95 percent confidence interval:
##  0.1868805 1.0000000
## sample estimates:
##         p 
## 0.2047833
prop.test(Tot_smokers, n, p=0.5, alternative = "less")
## 
##  1-sample proportions test with continuity correction
## 
## data:  Tot_smokers out of n, null probability 0.5
## X-squared = 465.26, df = 1, p-value < 2.2e-16
## alternative hypothesis: true p is less than 0.5
## 95 percent confidence interval:
##  0.000000 0.223902
## sample estimates:
##         p 
## 0.2047833
prop.test(Tot_smokers, n, p=0.5, alternative = "two.sided")
## 
##  1-sample proportions test with continuity correction
## 
## data:  Tot_smokers out of n, null probability 0.5
## X-squared = 465.26, df = 1, p-value < 2.2e-16
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.1836621 0.2276239
## sample estimates:
##         p 
## 0.2047833

Test su confronto proporzioni

Nel caso di confronto di due proporzioni l’ipotesi nulla che si vuole testare è \[ H_0:p_1=p_2\,\,\,\text{vs}\,\,\,H_A:p_1\ne p_2 \] dove \(p_1\) è la proporzione ipotizzata per il primo gruppo e \(p_2\) quella per il secondo gruppo. Bisogna a questo punto calcolare una stima della varianza

\[ \hat{\sigma}_{\hat{p}_1-\hat{p}_2}=\sqrt{\frac{\bar{p}(1-\bar{p})}{n_1}+\frac{\bar{p}(1-\bar{p})}{n_2}} \] dove \(\bar{p}=\frac{x_1+x_2}{n_1+n_2}\), \(x_1\) e \(x_2\) sono il numero di successi nel primo e secondo gruppo mentre \(n_1\) e \(n_2\) sono la numerosità campionaria del primo e secondo gruppo. La statistica test è

\[ ST=\frac{(\hat{p}_1-\hat{p}_2)-(p_1-p_2)}{\hat{\sigma}_{\hat{p}_1-\hat{p}_2}} \]

la quale, sotto \(H_0\), è distruibuita asintoticamente come una Normale Standard. Come sempre in R le cose sono molto semplici.

N <- dim(dta)[1] # immagino che il dataste sia la popolazione
n <- 200 # setto un ampiezza campionaria
lista <- sample(1:N, n, replace = F) 
campione1 <- dta[lista,] # estraggo il primo campione

lista <- sample(1:N, n, replace = F)
campione2 <- dta[lista,] # estraggo il secondo campione


# numero di fumatori nel primo capione
Tot_Smokers1 <- sum(campione1$smoker=="yes") 


# numero di fumatori nel secondo capione
Tot_Smokers2 <- sum(campione2$smoker=="yes")


n1 <- nrow(campione1)
n2 <- nrow(campione2)

prop.test(c(Tot_Smokers1, Tot_Smokers2), c(n1,n2), 
          alternative = "greater")
## 
##  2-sample test for equality of proportions with continuity correction
## 
## data:  c(Tot_Smokers1, Tot_Smokers2) out of c(n1, n2)
## X-squared = 0, df = 1, p-value = 0.5
## alternative hypothesis: greater
## 95 percent confidence interval:
##  -0.06913186  1.00000000
## sample estimates:
## prop 1 prop 2 
##  0.150  0.155
prop.test(c(Tot_Smokers1, Tot_Smokers2), c(n1,n2), 
          alternative = "less")
## 
##  2-sample test for equality of proportions with continuity correction
## 
## data:  c(Tot_Smokers1, Tot_Smokers2) out of c(n1, n2)
## X-squared = 0, df = 1, p-value = 0.5
## alternative hypothesis: less
## 95 percent confidence interval:
##  -1.00000000  0.05913186
## sample estimates:
## prop 1 prop 2 
##  0.150  0.155
prop.test(c(Tot_Smokers1, Tot_Smokers2), c(n1,n2), 
          alternative = "two.sided")
## 
##  2-sample test for equality of proportions with continuity correction
## 
## data:  c(Tot_Smokers1, Tot_Smokers2) out of c(n1, n2)
## X-squared = 0, df = 1, p-value = 1
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.08045996  0.07045996
## sample estimates:
## prop 1 prop 2 
##  0.150  0.155

Test di indipendenza Test Chi-quadrato e Test di Fisher

table(campione1$sex, campione1$smoker)
##         
##          no yes
##   female 82  15
##   male   88  15
(Tabella <- table(campione1$sex, campione1$smoker))
##         
##          no yes
##   female 82  15
##   male   88  15

Adesso assumiamo di voler testare l’indipendenza tra fumatori ( smoker ) e genere ( sex ), questo può essere fatto ricorrendo al test di chi quadrato.

chisq.test(Tabella)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  Tabella
## X-squared = 2.9725e-31, df = 1, p-value = 1

In teoria il test chi-quadrato è valido in contesti in cui il campione è molto grande…c’entra il teorema del limite centrale ovviamente. Quando in una tavola di contingenza notate che almeno una cella evidenzia una frequenza molto bassa (nella pratica, al di sotto di 5), si ricorre al test esatto di Fisher.

fisher.test(Tabella)
## 
##  Fisher's Exact Test for Count Data
## 
## data:  Tabella
## p-value = 1
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.396862 2.189358
## sample estimates:
## odds ratio 
##  0.9321649

Tests Omogeneità della Varianza

Il test di Levene è usato per testare se \(k\) campioni hanno varianza uguale. Alcuni test statistici, per esempio l’ANOVA, assumono che i diversi gruppi hanno varianze uguali. Il tets di Levene può essere utilizzato per verificare questa assunzione.

Il set di ipotesi che stiamo testando è

\[ H_0: \sigma^2_1=\sigma^2_2=\dots=\sigma^2_k\\ H_A: \sigma^2_i\ne \sigma_j^2 \,\,\text{ per almeno una coppia (i,j)} \]

library(car) # il test di Levene è implementato nel pacchetto car
## Loading required package: carData
campione1$region <- as.factor(campione1$region)
leveneTest(charges ~ region, data = campione1)
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   3  1.0429 0.3747
##       196

Il test di Levene nasce come alternativa al test di Barlett. Questo perchè il primo è meno sensibile del secondo nel caso in cui c’è violazione dell’ ipotesi di Normalità. Se si ha l’evidenza di star lavorando con dati Normali, o approssimativamente distribuiti come una v.c. Normale, allora il test di Barlett’s performa meglio. Il test non ha bisogno di nessuna libreria per essere usato, bartlett.test.

bartlett.test(charges ~ region, data = campione1)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  charges by region
## Bartlett's K-squared = 16.241, df = 3, p-value = 0.001012

Intervalli di confidenza con boostrap

I test che abbiamo discusso sino ad ora si basano sull’ ipotesi di normalità degli errori ( accertatevi sempre di star lavorando con dati normali). Il bootstrap consente di superare questa assunzione e costruire intervalli di confidenza anche se i dati non sono distribuiti normalmente. L’intuizione è di campionare ripetutamente da una distribuzione nota e poi calcolare le nostre stime. Non entreremo nei dettagli tecnici, fortunatamente i passaggi in R sono molto semplici

#install.packages(boot) assicuratevi di aver installato il pacchetto boot 
library(boot)
## 
## Attaching package: 'boot'
## The following object is masked from 'package:car':
## 
##     logit
## The following object is masked from 'package:lattice':
## 
##     melanoma
# immaginiamo di voler ottenere una stima bootrap di bwt.

L’unica cosa un pò fastidiosa è quella di definire una funzione in R. In pratica con la funzione samplemean stiamo calcolando la media campionaria \(\bar{x}\) per diverse replice boostrap.

samplemean <- function(x, d) {
  #Input
  # x= Osservazioni 
  # d= consente di effettuare estrazioni dal vettore x
  
  # Output
  # media campionara basata su un set di osservazioni definite in d
  
  
  return(mean(x[d]))
  
}

# Lanciamo il comando boot
b = boot(birthwt$bwt, # vettore su cui vogliamo fare Bootstrap
         samplemean,  # funzione iterare
         R=1000)      # numero di replice bootstrap


cat('Boostrap sample mean =',b$t0)
## Boostrap sample mean = 2944.587
cat("Standard deviation bootstrap Sample Mean = ", sd(b$t[,1]), "\n")
## Standard deviation bootstrap Sample Mean =  51.0165
ci = boot.ci(b, type="basic")
cat("Bootstrap sample mean 95% CI [", ci$basic[1,4], " - ", ci$basic[1,5], "]")
## Bootstrap sample mean 95% CI [ 2844.064  -  3045.568 ]

Il boostrap è uno strumento utilissimo, che consente di ottenere delle stime qualora si ha un sample size molto piccolo e una deviaizone molto forte dei dati dalla normalità. Non sottovalitatelo. Consiglio di dare uno sguardo al seguente articolo link