Farmacevtsko podjetje raziskuje učinek dveh novih zdravil (A in B) za zdravljenje visokega krvnega tlaka. 27 naključno izbranih bolnikov z visokim krvnim pritiskom je en mesec prejemalo zdravilo A in jim vsako jutro merili krvni tlak, iz katerega so izračunali povprečni jutranji krvni pritisk (Pritisk_A). V naslednjem mesecu smo zdravilo A zamenjali z zdravilom B in postopek merjenja ponovili (Pritisk_B). Lahko iz rezultatov sklepamo, da je učinkovitost zdravil različna?

####Enota: bolnik, velikost vzorca 27. 2 spremenljivki.

podatki <- read.table("./Krvni pritisk.csv", header = TRUE, sep=";")

head(podatki)
##   Pritisk_A Pritisk_B
## 1        82        90
## 2        80        90
## 3        88        99
## 4       108       118
## 5       116       123
## 6       136       137

Opis spremenljivk:

Zanima nas, če je razlika med učinkovitostjo zdravil. Gre za preizkus dvojic (vsako osebo testiramo 2x), imamo širok format podatkov (2 meritvi za eno osebo v isti vrstici).

H0: povprečjeA = povprečjeB; povprečje(d) = 0 -> povprečje diference je enako 0. H1: povprečjeA ni enako povp. B

PREDPOSTAVKE ZA PREIZKUS DVOJIC: - diference se normalno porazdeljujejo.

library(pastecs)
round(stat.desc(podatki), 2)
##              Pritisk_A Pritisk_B
## nbr.val          27.00     27.00
## nbr.null          0.00      0.00
## nbr.na            0.00      0.00
## min              63.00     76.00
## max             136.00    139.00
## range            73.00     63.00
## sum            2608.00   2754.00
## median           94.00     99.00
## mean             96.59    102.00
## SE.mean           3.46      3.06
## CI.mean.0.95      7.12      6.30
## var             323.87    253.62
## std.dev          18.00     15.93
## coef.var          0.19      0.16

Naredimo opisno statistiko. Učinkovitejše je bilo zdravilo A (nižji pristik). Zanima nas ali je statistično značilna razlika med povprečje A in B. Mediana in povprečje kažeta, da je zdravilo A boljše od B.

Če zelimo gledati, kje je večja variabilnost gledaš coef.var. V tem primeru so manjše razlike pri B. (0,16) Manjši kot je koeficient manjše so razlike.

podatki$Diferenca <- podatki$Pritisk_A - podatki$Pritisk_B

Imamo 2 načina preko katerih ugotavljamo ali je spremenljivka normalno porazdeljena: 1. SHAPIRO WILKOV TEST, 2. histogram.

library(ggplot2)
ggplot(podatki, aes(x = Diferenca)) +
  geom_histogram(position = "identity", binwidth = 3, colour = "black") +
  ylab("Frequency") +
  xlab("Diferenca")

## Porazdelitev asimetrična v desno. koeficient asimetrije bi bil pozitiven (večji od 0).

shapiro.test(podatki$Diferenca)
## 
##  Shapiro-Wilk normality test
## 
## data:  podatki$Diferenca
## W = 0.95383, p-value = 0.2652

Pri Shapiro-Wilkovem preizkusu je ničelna domneva(H0): Diference so porazdeljene normalno. H1: Diference niso porazdeljene normalno. Gledamo p-vrednost (0,2652). Glede na p-vrednost ne moremo zavrniti H0 -> iz tega sledi, da so diference normalno porazdeljene. Ta porazdelitev je glede na test normalna, zato uporabimo PARAMETRIČEN PREIZKUS.

t.test(podatki$Pritisk_A, podatki$Pritisk_B,
       paired = TRUE,
       alternative = "two.sided")
## 
##  Paired t-test
## 
## data:  podatki$Pritisk_A and podatki$Pritisk_B
## t = -4.6867, df = 26, p-value = 7.686e-05
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -7.779049 -3.035766
## sample estimates:
## mean difference 
##       -5.407407

Naredimo t-test (parametričen preizkus), pokličemo obe spremenljivki (pritisk A in pritisk B), Paired=TRUE je zato, ker so dvojice. “two.sided” je zato, ker v naprej nismo nič predpostavljali, da je zdravilo A ali B boljše (puščamo obe možnosti).

Na podlagi p-vrednosti zavrnemo ničelno domnevo(HO) pri p< 0,001. Ugotavljamo, da se povprečni pritisk pri zdravilu A razlikuje od povprečnega pritiska pri zdravilu A. Mean difference (-5,47…) je negativen. Od prve spremenljivke (pritisk_A) odštevamo drugo (pritisk_B). Mean difference (ocena povprečne diference) je negativna, zato je zdravilo A učinkovitejše. Če bi bilo v zapisu obratno (najprej B in potem A), se bi spremenil t (bil bi +), 95% interval zaupanja bi bil +, ter mean difference bi bil +. Vrstni red ne vpliva na ugotovitve

Kako preverim, kako velike so razlike? VELIKOST UČINKA. Za tak preizkus je cohens_d statistika

library(effectsize)
cohens_d(podatki$Diferenca)
## Cohen's d |         95% CI
## --------------------------
## -0.90     | [-1.34, -0.45]

-0.90. Pogledaš na slide 5,6 tabele. Vrednost gledamo absolutno (minus zanemarimo). 0,90 pade v velike razlike.

interpret_cohens_d(0.90, rules = "sawilowsky2009")
## [1] "large"
## (Rules: sawilowsky2009)

Wilcox test delamo takrat, ko imamo odvisne vzorce. UPORABIMO TAKRAT KO JE KRŠENA NORMALNOST (Če je p-vrednost pri Shapiro wilk testu manjša ali enaka od 0,05) -> direktno to uporabimo!

WILCOX TEST: Ničelna domneva(H0): lokacije porazdelitve spremenljivk sta ENAKI. (Če bi narisali porazdelitve lokacij pritisk A in pritisk B bi se prekrivale (bile bi enaki). H1: lokaciji porazdelitve sta različni. Glede na p-vrednost vidimo, da sta lokaciji različni.

wilcox.test(podatki$Pritisk_A, podatki$Pritisk_B, 
            paired = TRUE,
            correct = FALSE,
            exact = FALSE,
            alternative = "two.sided")
## 
##  Wilcoxon signed rank test
## 
## data:  podatki$Pritisk_A and podatki$Pritisk_B
## V = 33, p-value = 0.0004822
## alternative hypothesis: true location shift is not equal to 0
library(ggplot2)
Fig1 <- ggplot(podatki, aes(x = Pritisk_A)) +
  geom_histogram(binwidth = 5, colour="gray") +
  scale_x_continuous(breaks = seq(50, 150, 5), limits = c(50, 150)) +
  scale_y_continuous(limits = c(0, 5)) +
  ylab("Frequency")

Fig2 <- ggplot(podatki, aes(x = Pritisk_B)) +
  geom_histogram(binwidth = 5, colour="gray") +
  scale_x_continuous(breaks = seq(50, 150, 5), limits = c(50, 150)) +
  scale_y_continuous(limits = c(0, 5)) +
  ylab("Frequency")

library(ggpubr)
ggarrange(Fig1, Fig2,
          ncol = 1)

Ali sta ti dve lokaciji enaki? Gleda na sliko lahko vidimo, da A leži bolj levo (je bolj učinkovito zdravilo - manjši pritisk). Ugotavljamo, da sta lokaciji različni pri p<0,001.

Pri Wilcox testu (vsot rangov ali predznačnih rangov) za velikost učinka uporabimo rank biserialno koleracijo.

library(effectsize)
effectsize(wilcox.test(podatki$Pritisk_A, podatki$Pritisk_B, 
                       paired = TRUE,
                       correct = FALSE,
                       exact = FALSE,
                       alternative = "two.sided"))
## r (rank biserial) |         95% CI
## ----------------------------------
## -0.80             | [-0.91, -0.58]

Gledaš lestvico r=0,8. Zelo velike razlike, ali bo jemal zdravilo A ali B.

interpret_rank_biserial(0.80)
## [1] "very large"
## (Rules: funder2019)