www = "http://www.biostatisticaumg.it/umgrisk/magnagraecia.csv"
magnagraecia = read.csv(www, header = TRUE)
attach(magnagraecia)



# The mean ages for patients with fibroids and sarcomas
t.test(AGE ~ OUTCOME)



# LDH isoenzymes levels between the benign and cancer patients
wilcox.test(LDH1 ~ OUTCOME)
tapply(LDH1, OUTCOME, mean)
tapply(LDH1, OUTCOME, median)
wilcox.test(LDH2 ~ OUTCOME)
tapply(LDH2, OUTCOME, mean)
tapply(LDH2, OUTCOME, median)
wilcox.test(LDH3 ~ OUTCOME)
tapply(LDH3, OUTCOME, mean)
tapply(LDH3, OUTCOME, median)
wilcox.test(LDH4 ~ OUTCOME)
tapply(LDH4, OUTCOME, mean)
tapply(LDH4, OUTCOME, median)
wilcox.test(LDH5 ~ OUTCOME)
tapply(LDH5, OUTCOME, mean)
tapply(LDH5, OUTCOME, median)



# unpublished Figure by Di Cello et al. depicting UMG
par(mar=c(5.1,4.1,4.1,2.1), mfrow = c(1,2))
plot(LDH1, LDH3, col = "white", main = "LDH3 versus LDH1 in fibroids and sarcomas")
x = seq(from = 10, to = 40, by = 0.1)
y = 29 - 24/x
lines(x,y, lty = 4, col = "blue", lwd = 3)
points(LDH1[OUTCOME == "benignant"], LDH3[OUTCOME == "benignant"], pch = 21, bg = "green", col = "darkgreen")
points(LDH1[OUTCOME == "malignant"], LDH3[OUTCOME == "malignant"], pch = 23, bg = "red", col = "black")
text(25, 13, "green bullets: fibroids", col = "darkgreen")
text(17, 31, "red diamonds: sarcomas", col = "red")
text(16, 26, "blue line:", col = "blue")
text(16, 24, "LDH3 = 29 - 24/LDH1", col = "blue")
plot(LDH1, LDH3, col = "white", main = "U.M.G. risk index in fibroids and sarcomas")
x = seq(from = 10, to = 40, by = 0.1)
y = 29 - 24/x
lines(x,y, lty = 4, col = "blue", lwd = 3)
text(LDH1[OUTCOME == "malignant"], LDH3[OUTCOME == "malignant"], UMG[OUTCOME == "malignant"], col = "red")
text(LDH1[OUTCOME == "benignant"], LDH3[OUTCOME == "benignant"], UMG[OUTCOME == "benignant"], col = "darkgreen")
text(20, 34, "high-risk group", col = "red")
text(20, 14, "low-risk group", col = "darkgreen")
text(20, 13, "LDH3 + 24/LDH1 < 29", col = "darkgreen")



# unpublished Figure by Di Cello et al. depicting conditional tree classification according to UMG
# install.packages("party")
# citation("party")
library(party)
model = ctree(OUTCOME ~ UMG)
plot(model)



# Table 1
# set.seed is here for reproducibility purpose
# install.packages("pROC")
# citation("pROC")
library(pROC)
set.seed(1234)
(roc1 = roc(OUTCOME ~ LDH1, auc=TRUE, ci=TRUE))
(roc2 = roc(OUTCOME ~ LDH2, auc=TRUE, ci=TRUE))
(roc3 = roc(OUTCOME ~ LDH3, auc=TRUE, ci=TRUE))
(roc4 = roc(OUTCOME ~ LDH4, auc=TRUE, ci=TRUE))
(roc5 = roc(OUTCOME ~ LDH5, auc=TRUE, ci=TRUE))
(rocumg = roc(OUTCOME ~ UMG, auc=TRUE, ci=TRUE))



# Table 2 accuracy of U.M.G. risk index
table(OUTCOME, UMG > 29)
(sensit = 43/43)
(specif = 2202/(2202+9))
(plr = sensit/(1-specif))
(nlr = (1-sensit)/specif)



# Table 3 ROC-AUC and confidence intervals
# set.seed is here for reproducibility purpose
set.seed(4321)
ci.se(roc1, specificities = c(0.98, 0.99, 0.995), conf.level = 0.99)
ci.se(roc2, specificities = c(0.98, 0.99, 0.995), conf.level = 0.99)
ci.se(roc3, specificities = c(0.98, 0.99, 0.995), conf.level = 0.99)
ci.se(roc4, specificities = c(0.98, 0.99, 0.995), conf.level = 0.99)
ci.se(roc5, specificities = c(0.98, 0.99, 0.995), conf.level = 0.99)
ci.se(rocumg, specificities = c(0.98, 0.99, 0.995), conf.level = 0.99)


# Table 3 comparisons of ROC-AUC
# set.seed is here for reproducibility purpose
set.seed(1234)
roc.test(roc1, rocumg)
roc.test(roc2, rocumg)
roc.test(roc3, rocumg, method = "delong")
roc.test(roc4, rocumg, method = "delong")
roc.test(roc5, rocumg, method = "delong")