1 Εύρεση διαστημάτων εμπιστοσύνης

1.1 Για την κανονική κατανομή

Η γενική προσέγγιση υπολογισμού του CI είναι: CI = σημειακή εκτίμηση (π.χ. μέσος)+/-(critical value)x(τυπικό σφάλμα). Η κρίσιμη τιμή για την κανονική κατανομή βρίσκεται με την εντολή qnorm() στην R. Ας χρησιμοποιήσουμε το σετ δεδομένων iris για να βρούμε τα διαστήματα εμπιστοσύνης για τον μέσο του είδους setosa για τη μεταβλητή Sepal.Length

# R
str(iris)
## 'data.frame':    150 obs. of  5 variables:
##  $ Sepal.Length: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
##  $ Sepal.Width : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
##  $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
##  $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
##  $ Species     : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...

Δημιουργούμε ένα νέο σετ δεδομένων, το υποσύνολο του iris μόνο για το είδος setosa και βρίσκουμε τον αριθμό των ατόμων του δείγματος, τον μέσο, την τυπική απόκλιση και το τυπικό σφάλμα.

# R
setosa <- subset(iris, Species == "setosa")
n <- length(setosa$Sepal.Length)
mean <- mean(setosa$Sepal.Length)
sd <- sd(setosa$Sepal.Length)
se <- sd/sqrt(n)

Βρίσκουμε το διάστημα εμπιστοσύνης για επίπεδο 95%.

# R
ci1<- qnorm(0.025, mean = mean, sd = sd, lower.tail = TRUE)
ci1
## [1] 4.315133
ci2<- qnorm(0.025, mean = mean, sd = sd, lower.tail = FALSE)
ci2
## [1] 5.696867

Το διάστημα εμπιστοσύνης είναι 4.31 ως 5.7 cm με μέση τιμή 5 cm. Σχεδιάζουμε ένα ιστόγραμμα της κατανομής των δεδομένων και ένα boxplot.

# R
hist(setosa$Sepal.Length,
main = "Sepal length of Iris setosa",
xlab = "Sepal length (cm)",
ylab = "Frequency"
)

boxplot(setosa$Sepal.Length)

1.2 Για την κατανομή t - Student

Επειδή δεν γνωρίζουμε την πραγματική τυπική απόκλιση του πληθυσμού και χρησιμοποιούμε το τυπικό σφάλμα από την τυπική απόκλιση του δείγματος, πιο σωστό είναι να χρησιμοποιήσουμε την κατανομή t - Student.

# R
alpha=0.05
df=n-1
t.score = qt((alpha/2), df,lower.tail=F)
margin.error <- t.score * se
lower.bound <- mean - margin.error
upper.bound <- mean + margin.error
print(c(lower.bound,upper.bound))
## [1] 4.905824 5.106176

1.3 Boxplot και violin plot

Κατασκευάζουμε ένα boxplot για τα τρία είδη χωριστά, προκειμένου να συγκρίνουμε τις κατανομές τους

# R
setosa <- subset(iris, Species == "setosa")
versicolor <- subset(iris, Species == "versicolor")
virginica <- subset(iris, Species == "virginica")
boxplot(setosa$Sepal.Length, versicolor$Sepal.Length, virginica$Sepal.Length, names=c("setosa", "versicolor", "virginica"))

Τα boxplot είναι πολύ χρήσιμα για να συγκρίνουν συνοπτικά διάφορα σετ δεδομένων, αλλά είναι υπερβολικά λιτά και απλοποιούν την παρεχόμενη πληροφορία. Για το λόγο αυτό, μπορούμε να κατασκευάσουμε violin plots. Ένα violin plot (γράφημα βιολιού) απεικονίζει κατανομές αριθμητικών δεδομένων για μία ή περισσότερες ομάδες χρησιμοποιώντας καμπύλες πυκνότητας (density plots). Το πλάτος κάθε καμπύλης αντιστοιχεί με την κατά προσέγγιση συχνότητα των σημείων δεδομένων σε κάθε περιοχή. Θα χρειαστούμε το πακέτο vioplot.

# R
library("vioplot")
## Warning: package 'zoo' was built under R version 4.3.0
vioplot(setosa$Sepal.Length, versicolor$Sepal.Length, virginica$Sepal.Length, names=c("setosa", "versicolor", "virginica"))

Μπορούμε να βελτιώσουμε το γράφημα προσθέτοντας χρώμα και λεζάντα.

# R
vioplot(setosa$Sepal.Length, versicolor$Sepal.Length, virginica$Sepal.Length, names=c("setosa", "versicolor", "virginica"), main = "Sepal Length", col=c("lightgreen", "lightblue", "palevioletred"))
legend("topleft", legend=c("setosa", "versicolor", "virginica"), fill=c("lightgreen", "lightblue", "palevioletred"), cex = 1.5)

2 Έλεγχος κανονικότητας

2.1 Έλεγχος κανονικότητας με το σετ δεδομένων Arthritis

Παράδειγμα με το σετ δεδομένων “Arthritis” και το πακέτο vcd. Πρόκειται για κλινική μελέτη για 84 ασθενείς με ρευματοειδή αρθρίτιδα διαφορετικού φύλου και ηλικίας, που δέχτηκαν θεραπεία με ένα κανονικό φάρμακο ή με εικονική θεραπεία (placebo). Ο βαθμός βελτίωσης της υγείας τους αξιολογήθηκε με βάση μια κλίμακα κλινικών συμπτωμάτων. Κάνουμε μια αρχική οπτική εκτίμηση με ένα γράφημα της κατανομής των τιμών της ηλικίας.

# R
library(vcd)
str(Arthritis)
## 'data.frame':    84 obs. of  5 variables:
##  $ ID       : int  57 46 77 17 36 23 75 39 33 55 ...
##  $ Treatment: Factor w/ 2 levels "Placebo","Treated": 2 2 2 2 2 2 2 2 2 2 ...
##  $ Sex      : Factor w/ 2 levels "Female","Male": 2 2 2 2 2 2 2 2 2 2 ...
##  $ Age      : int  27 29 30 32 46 58 59 59 63 63 ...
##  $ Improved : Ord.factor w/ 3 levels "None"<"Some"<..: 2 1 1 3 3 3 1 3 1 1 ...
hist(Arthritis$Age, col="dark red", breaks=10, xlab="Age", ylab="No of patients", main="All patients", las=1)

qqnorm(Arthritis$Age, pch=19)

Κατασκευάζουμε ένα νέο Q-Q plot με την ευθεία της κανονικής κατανομής για σύγκριση.

# R
nAge<-rnorm(length(Arthritis$Age), mean=mean(Arthritis$Age), sd=sd(Arthritis$Age))
qqnorm(Arthritis$Age, pch=19) 
qqline(nAge)

Μπορούμε εύκολα να κατασκευάσουμε ένα Q-Q plot με ενσωματωμένο το διάστημα εμπιστοσύνης, μέσα από το πακέτο ggpubr.

# R
library(ggpubr)
## Warning: package 'ggpubr' was built under R version 4.3.0
ggqqplot(Arthritis$Age)

Ο αριθμητικός έλεγχος κανονικότητας γίνεται με το τεστ Shapiro-Wilk. Aποδίδει δύο τιμές, την τιμή W που είναι το στατιστικό και την p-value. Στην πιο κάτω ανάλυση βλέπουμε ότι τα δεδομένα δεν είναι κανονικά κατανεμημένα.

# R
shapiro.test(Arthritis$Age)
## 
##  Shapiro-Wilk normality test
## 
## data:  Arthritis$Age
## W = 0.91913, p-value = 5.813e-05

3 Σύγκριση δύο μέσων τιμών

3.1 Σύγκριση δύο μέσων τιμών με το σετ δεδομένων chickwts

Θα δούμε ένα παράδειγμα για το βάρος που αποκτούν τα κοτόπουλα μιας κτηνοτροφικής μονάδας με βάση την κύρια πηγή πρωτεΐνης στην τροφή τους. Tο σετ δεδομένων ονομάζεται chickwts και περιέχει τιμές βάρους για έξι διαφορετικά είδη τροφών.

# R
str(chickwts)
## 'data.frame':    71 obs. of  2 variables:
##  $ weight: num  179 160 136 227 217 168 108 124 143 140 ...
##  $ feed  : Factor w/ 6 levels "casein","horsebean",..: 2 2 2 2 2 2 2 2 2 2 ...
levels(chickwts$feed)
## [1] "casein"    "horsebean" "linseed"   "meatmeal"  "soybean"   "sunflower"

To ερώτημά είναι αν η μέση τιμή του βάρους των ζώων που έχουν τραφεί με μια συγκεκριμένη τροφή είναι μεγαλύτερη από αυτήν των ζώων που έχουν τραφεί με κάποια άλλη. Θα αναπαραστήσουμε γραφικά τα δεδομένα. Χρησιμοποιούμε τη δομή του dataframe για να δημιουργήσουμε θηκογράμματα (boxplots) με την χρήση της qplot() από το πακέτο ggplot2.

# R
library(ggplot2)
qplot(feed, weight, data=chickwts, geom="boxplot", col=feed)
## Warning: `qplot()` was deprecated in ggplot2 3.4.0.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Θέλουμε να συγκρίνουμε το βάρος των ζώων που τρέφονται με λιναρόσπορο με αυτά που τρέφονται με κουκιά. Αρχικά δημιουργούμε σετ δεδομένων τια τις δύο αυτές τροφές και ελέγχουμε την κανονικότητα.

# R
linseed <- subset(chickwts, feed == "linseed")
horsebean <- subset(chickwts, feed == "horsebean")
shapiro.test(linseed$weight)
## 
##  Shapiro-Wilk normality test
## 
## data:  linseed$weight
## W = 0.96931, p-value = 0.9035
shapiro.test(horsebean$weight)
## 
##  Shapiro-Wilk normality test
## 
## data:  horsebean$weight
## W = 0.93758, p-value = 0.5264

Βλέπουμε ότι δεν μπορούμε να απορρίψουμε την μηδενική υπόθεση, ότι οι δύο κατανομές αυτές είναι κανονικές. Συγκρίνουμε του δύο μέσους με το t-test και με την συνάρτηση t.test() για 95% επιπεδο εμπιστοσύνης (default).

# R
t.test(linseed$weight,horsebean$weight)
## 
##  Welch Two Sample t-test
## 
## data:  linseed$weight and horsebean$weight
## t = 3.0172, df = 19.769, p-value = 0.006869
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  18.0403 99.0597
## sample estimates:
## mean of x mean of y 
##    218.75    160.20

Επαναλαμβάνουμε για 99% επίπεδο εμπιστοσύνης.

# R
t.test(linseed$weight,horsebean$weight, conf.level = 0.99)
## 
##  Welch Two Sample t-test
## 
## data:  linseed$weight and horsebean$weight
## t = 3.0172, df = 19.769, p-value = 0.006869
## alternative hypothesis: true difference in means is not equal to 0
## 99 percent confidence interval:
##    3.267538 113.832462
## sample estimates:
## mean of x mean of y 
##    218.75    160.20

Μποορύμε να ορίσουμε το t-test ως object (π.χ. student) και να καλέσουμε τα στοιχεία του (π.χ. p.value) για ευκολία.

# R
student<-t.test(linseed$weight,horsebean$weight)
student$p.value
## [1] 0.006869064
student$conf.int
## [1] 18.0403 99.0597
## attr(,"conf.level")
## [1] 0.95

3.2 Σύγκριση δύο μέσων τιμών με το σετ δεδομένων Arthritis (μη παραμετρικός έλεγχος)

Πώς προχωράμε στη σύγκριση δειγμάτων που δεν είναι κανονικά κατανεμημένα; Θα καταφύγουμε σε κάποιο μη-παραμετρικό έλεγχο όπως αυτός που πραγματοποιείται με την μέθοδο Mann-Whitney ή Wilcoxon Rank-Sum (που είναι πρακτικά ισοδύναμες). Στο σετ δεδομένων Arthritis δεν είχαμε κανονικά κατανεμημένες τιμές ηλικίας. Θέλουμε να δούμε αν το προφίλ των ασθενών που αποκρίνονται καλύτερα στη θεραπεία έχει μια συγκεκριμένη ηλικιακή τάση. Χωρίζουμε τους ασθενείς με βάση το αποτέλεσμα της θεραπείας.

# R
None <- subset(Arthritis, Improved == "None")
Marked <- subset(Arthritis, Improved == "Marked")
shapiro.test(None$Age)
## 
##  Shapiro-Wilk normality test
## 
## data:  None$Age
## W = 0.94019, p-value = 0.029
shapiro.test(Marked$Age)
## 
##  Shapiro-Wilk normality test
## 
## data:  Marked$Age
## W = 0.91064, p-value = 0.02049

Καμία από τις δύο νέες ομάδες δεδομένων δεν ακολουθούν την κανονική κατανομή. Εφαρμόζουμε τον μη-παραμετρικό έλεγχο Wilcoxon Rank Sum.

# R
wilcox.test(None$Age,Marked$Age)
## Warning in wilcox.test.default(None$Age, Marked$Age): cannot compute exact
## p-value with ties
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  None$Age and Marked$Age
## W = 394.5, p-value = 0.02056
## alternative hypothesis: true location shift is not equal to 0

Όπως βλέπουμε από το p-value που είναι μικρότερο του 0.05, υπάρχουν διαφορές μεταξύ των δύο ομάδων. Η μηδενική υπόθεση (H0: δεν υπάρχουν διαφορές ανάμεσα στον μέσο των δύο διαφορετικών ομάδων) απορρίπτεται. Άρα η διατροφή με διαφορετική πηγή πρωτεϊνών μαίζει ρόλο στο βάρος των πτηνών.

3.3 Σύγκριση δύο μέσων τιμών σε δεδομένα κλινικών ερευνών

θα ακολουθήσουμε την ίδια διαδικασία με ένα πραγματικό παράδειγμα με σετ δεδομένων από κλινικές έρευνες που σχετίζονται με τις ημέρες επιβίωσης καρκινοπαθών και διάφορους παράγοντες. Το σετ υπάρχει στο αρχείο “cancer1.csv” που πρέπει να κατεβάσετε από το eclass και που αποκτήθηκε από το kaggle (https://www.kaggle.com/datasets/imtkaggleteam/clinical-dataset/data). Προσοχή! Μεταφέρετε το αρχείο “cancer1.csv” στο working directory του R που δουλεύετε αυτήν την άσκηση.

# R
data<-read.csv("cancer1.csv")
str(data)
## 'data.frame':    95 obs. of  14 variables:
##  $ Patient.ID   : chr  "P109342" "P124450" "P131833" "P131888" ...
##  $ Survival_time: int  2329 2532 2271 2193 2387 2225 2240 2314 299 2295 ...
##  $ death        : int  1 0 0 0 0 0 0 0 1 0 ...
##  $ Tumor_size   : num  2.6 3.5 2 3 4 2.5 2.2 1.2 2.8 2 ...
##  $ Grade        : int  3 3 2 2 2 2 3 3 1 2 ...
##  $ Stage        : chr  "IB" "IB" "IA2" "IA3" ...
##  $ Age          : int  67 68 80 63 88 66 69 52 67 59 ...
##  $ Sex          : chr  "Male" "Female" "Female" "Male" ...
##  $ Cigarette    : chr  "Former" "Former" "Never" "Former" ...
##  $ Pack_per_year: num  50 52.5 0 47 0 ...
##  $ Therapy      : chr  "None" "None" "None" "None" ...
##  $ batch        : int  1 1 1 1 1 1 2 2 2 1 ...
##  $ EGFR         : chr  NA "Negative" "Negative" "Negative" ...
##  $ KRAS         : chr  NA "Negative" "Negative" "G12C" ...
hist(data$Survival_time,
main = "Survival time of patients",
xlab = "Survival time (days)",
ylab = "Frequency"
)

boxplot(data$Survival_time)

vioplot(data$Survival_time)

ggqqplot(data$Survival_time)

shapiro.test(data$Survival_time)
## 
##  Shapiro-Wilk normality test
## 
## data:  data$Survival_time
## W = 0.92043, p-value = 2.353e-05

Βλέπουμε ότι το σύνολο των τιμών του χρόνου επιβίωσης δεν ακολουθεί την κανονική κατανομή. Κάνουμε το ερώτημα, αν η επιβίωση των ασθενών σχετίζεται με το κάπνισμα (μεταβλητή Cigarette με τρεις κατηγορίες: “current”, “former” and “never”). Χωρίζουμε τα δεδομένα στις τρεις κατηγορίες και ελέγχουμε ξανά την κανονικότητα. Στη συνέχεια και ανά δύο κατηγορίες, θα κάνουμε t.test αν τα δεδομένα ακολουθούν κανονική κατανομή, ή τον μη-παραμετρικό έλεγχο Wilcoxon Rank Sum, αν δεν έχουμε κανονικές κατανομές.

# R
never <- subset(data, Cigarette == "Never")
current <- subset(data, Cigarette == "Current")
former <- subset(data, Cigarette == "Former")
shapiro.test(never$Survival_time)
## 
##  Shapiro-Wilk normality test
## 
## data:  never$Survival_time
## W = 0.8959, p-value = 0.02464
shapiro.test(current$Survival_time)
## 
##  Shapiro-Wilk normality test
## 
## data:  current$Survival_time
## W = 0.92154, p-value = 0.3317
shapiro.test(former$Survival_time)
## 
##  Shapiro-Wilk normality test
## 
## data:  former$Survival_time
## W = 0.92176, p-value = 0.0007286
wilcox.test(never$Survival_time,current$Survival_time)
## 
##  Wilcoxon rank sum exact test
## 
## data:  never$Survival_time and current$Survival_time
## W = 118, p-value = 0.9252
## alternative hypothesis: true location shift is not equal to 0
wilcox.test(never$Survival_time,former$Survival_time)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  never$Survival_time and former$Survival_time
## W = 610.5, p-value = 0.4701
## alternative hypothesis: true location shift is not equal to 0
wilcox.test(current$Survival_time,former$Survival_time)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  current$Survival_time and former$Survival_time
## W = 284, p-value = 0.3836
## alternative hypothesis: true location shift is not equal to 0

Η0: Οι μέσοι όροι για κάθε σύγκριση δεν διαφέρουν. Καμία από τις συγκρίσεις δεν έδωσε p-value < 0,05 άρα σε καμία από τις τρεις συγκρίσεις δεν μπορέσαμε να απορρίψουμε τη μηδενική υπόθεση. Τα επίπεδα καπνίσματος δεν φαίνονται να επηρεάζουν τον χρόνο επιβίωσης των ασθενών.

# R
boxplot(never$Survival_time, former$Survival_time, current$Survival_time, names=c("never", "former", "current"))

vioplot(never$Survival_time, former$Survival_time, current$Survival_time, names=c("never", "former", "current"), main = "Survival time vs. smoking history", col=c("lightgreen", "lightblue", "palevioletred"))
legend("topleft", legend=c("never", "former", "current"), fill=c("lightgreen", "lightblue", "palevioletred"), cex = 1.5)