ANOVA - BLOOD PRESSURE

Abstract

Il nostro database, estrapolato dal sito “https://www.statslectures.com/r-scripts-datasets”, prende a confronto una terapia innovativa per ridurre la pressione arteriosa sistolica (PAS) contro placebo (assenza di farmaco). La pressione è definita sistolica quando è costantemente superiore a 120 mmHg ma inferiore a 140 mmHg. Il proposito di questo elaborato è valutare l’efficacia di un nuovo farmaco anti-ipertensivo prima dell’immissione in commercio. L’ipotesi zero è che non ci sia differenza tra i trattamenti. Utilizzeremo il test ANOVA per confrontare i due gruppi di dati ottenuti dall’esperimento, confrontando la variabilità interna a questi gruppi con la variabilità tra i gruppi.

INTRODUZIONE:

Prendendo a riferimento 25 osservazioni di pazienti, di cui si riconosce il limite del campionamento ridotto, si confronta la riduzione media della PAS (variabile dipendente - v_dip) in relazione all’utilizzo del tipo di terapia (variabile indipendente_v_indip)

Una volta inserito il database nel sistema di R, assegnati i vettori “PAS” e “TYPE”, e caratterizzati i livelli del fattore, si va a delineare il boxplot come segue:

plot(PAS~TYPE, col = c("blue", 602, "turquoise"))

Già da questo boxplot iniziale si può vedere una differenza tra le medie di PAS in chi ha assunto il farmaco e chi no; la media di P (cioè coloro che hanno assunto il placebo) è di pochissimo sopra lo zero. la media di A (cioè coloro che hanno assunto il farmaco sperimentale) si avvicina ma non raggiunge il valore 10. Una osservazione è fuori dall’intervallo di confidenza del grafico in colore blu. Il che è un primo segnale sull’andamento di questo studio.

Proseguo con la verifica delle assunzioni di ANOVA, cioè : 1) La variabile dipendente quantitativa deve avere una distribuzione Normale in tutte le popolazioni corrispondenti ai gruppi campionati; 2) Le varianze in tutte le popolazioni corrispondenti ai gruppi campionati devono essere uguali o non devono contenere valori anomali;

boxplot(PAS, main="Boxplot variabile effetto", col = "green")

Per compredere se la variabile dipendente quantitativa si distribuisce come una Normale possiamo avvalerci sia dell’utilizzo di metodi grafiCi e differenti test statistici. Il primo grafico necessita del calcolo dei residui standardizzati, che devono ricadere nell’intervallo +/- 1.96, ossia nell’area di probabilità +/- 0.025, in quanto il test è di tipo bidirezionale, con alfa=0.05. In altre parole, dovremmo avere tutte le osservazioni disposte casualmente ma nel intervallo delineato, invece abbiamo un valore fuori dall’intervallo il che metterebbe in discussione l’applicabilità del test.

plot(resid, ylim=c(-3,3), main = "Residui standardizzati",
     xlab = "Osservazioni", ylab = "Residui")
abline(h=-1.96, col="red")

abline(h=1.96, col="red")

La realizzazione del grafico quantile-quantile plot non chiarisce ancora il nostro dubbio. La retta, definita intercetta, ossia il valore che assume la y quando la x è uguale a 0 ed il coefficiente di regressione, che esprime la variazione della y data una variazione unitaria della x. il calcolo dell’intercetta e del suo coefficiente angolare effettuato da R ci restituisce un grafico che non si discosta dalle attese, la retta attraversa il più vicino possibile i punti osservati.

qqnorm(PAS, main="Normalità dei dati", pch=20, col="gray")      
qqline(PAS, col="darkblue")

Il grafico raffigura una bisettrice attorno a cui sono presenti i punti del nostro data, salvo eccezioni in 3 osservazioni inferiori al primo quantile e 3 osservazioni nel quarto quantile.

I test grafici non ci hanno restituito con sicurezza la prima assunzione di ANOVA. Per dimostrare la nostra ipotesi ci siamo avvalsi del test non parametrico di *Shpariro-Wilk**.

shapiro.test(PAS)

    Shapiro-Wilk normality test

data:  PAS
W = 0.96797, p-value = 0.5941

Dall’output ottenuto il valore di W è molto vicino ad 1. Inoltre, partendo dall’assunto che l’ipotesi nulla sia che i dati seguono una distribuzione normale, con p-value > 0.05, posso accettare l’ipotesi nulla di normalità dei dati.

Per verificare il secondo assunto, Le varianze in tutte le popolazioni corrispondenti ai gruppi campionati devono essere uguali o non devono contenere valori anomali proseguo con i test di Levene e Bartlett.

leveneTest(PAS ~ TYPE, data = ANOVA.FARMACO)

Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  1  0.0456 0.8328
      23               

I gradi di libertà indicati sono 1 e 23, il valore di P è 0.8328. Quindi l’ipotesi di partenza è confermata.
Di seguito invece ho utilizzato il test di Bartlett per verificare se k campioni provengono da popolazioni con varianze uguali. Le varianze uguali tra le popolazioni sono chiamate omoscedasticità o omogeneità delle varianze.In questo caso, otteniamo p-value > 0.05. Quindi accettiamo (o per meglio dire, non rifiutiamo) l’ipotesi nulla: le varianze sono omogenee. Il test di ANOVA può essere utilizzato.

Inseriamo su R il comando di calcolo di ANOVA sul nostro dataset.

aov_1w

Call:
   aov(formula = PAS ~ TYPE, data = ANOVA.FARMACO)

Terms:
                     TYPE Residuals
Sum of Squares   190.4067 1487.4333
Deg. of Freedom         1        23

Residual standard error: 8.041829
Estimated effects may be unbalanced

summary(aov_1w)

            Df Sum Sq Mean Sq F value Pr(>F)  
TYPE         1  190.4  190.41   2.944 0.0996 .
Residuals   23 1487.4   64.67                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Il risultato della summary fornisce un dato di p-value di 0.0996, confermando con margine minimo l’ipotesi nulla cioè che non ci sia differenza tra i trattamenti e che il farmaco non abbia effetto rispetto al placebo.

Effettuiamo il test HSD di Tukey, per capire su quale gruppo esista o meno la differenza tra le medie:

hsd

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = PAS ~ TYPE, data = ANOVA.FARMACO)

$TYPE
         diff       lwr      upr     p adj
P-A -5.633333 -12.42487 1.158199 0.0996264

L’output evidenzia i due intervalli di confidenza dei valori tra le differenze delle medie. Stabilendo un valore di alpha = 0.05 è possibile affermare che la sicurezza di ottenere un valore di differenza tra le medie compreso negli intervalli di confidenza pari a 1-alpha.

plot(hsd)

Il grafico realizzato consente di visualizzare l’intervallo di confidenza all’interno del quale si pone lo 0 come variabile auspicabile; nel nostro caso, lo 0 rientra nell’intervallo plausibile. Non permette di sottolineare alcuna differenza tra le medie. Considerando l’ANOVA come generalizzazione del modello di regressione, esplicitiamo dunque la regressione (reg.lm)

summary(reg.lm)

Call:
lm(formula = PAS ~ TYPE, data = ANOVA.FARMACO)

Residuals:
     Min       1Q   Median       3Q      Max 
-17.3333  -5.7000  -0.3333   5.3000  12.6667 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)    9.333      2.076   4.495 0.000164 ***
TYPEP         -5.633      3.283  -1.716 0.099626 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 8.042 on 23 degrees of freedom
Multiple R-squared:  0.1135,    Adjusted R-squared:  0.07494 
F-statistic: 2.944 on 1 and 23 DF,  p-value: 0.09963

calcoliamo l’intercetta del campione:

aov_1w$coefficients

(Intercept)       TYPEP 
   9.333333   -5.633333 

CONCLUSIONI

Il farmaco ha una tendenza non statisticamente significativa. l’“ipotesi zero”, cioè che il farmaco e il placebo non avessero influenza sulle condizioni dei pazienti è confermata. Il valore p-value è di poco superiore a 0.05 il chè darebbe un margine alla considerazione della ipotesi che il farmaco abbia effetto positivo, ipotesi che potrebbe essere confutata allargando a più di 25 pazienti il campione statistico su cui testarlo.

Pierluigi Mancarella