#Import Excel
library(readxl)
df <- read_excel("C:/Users/micnl/Google Drive/Research/Bergman Clinics/ROBOT NEUROSPINE/rdata.xlsx")
str(df)
table(df$TLIF)

#reformat
cols <- c("TLIF", "listhesis", "ldh" ,"fbss", "stenosis", "pseudoarthrosis", "radiculopathy", "clbp", "smoker", "male", "caucasian", "asa3", "opioids", "asthma", "priorsurgery", "mcidback", "mcidleg", "mcidodi")
df[cols] <- lapply(df[cols], factor)
str(df)

#Matching
library(MatchIt)
m.df <- matchit(TLIF ~ age + male + prevasback + prevasleg + preodi, data = df, method = "nearest", ratio = 1)
summary(m.df)
df <- match.data(m.df)
str(df)
table(df$TLIF)

#Summary Statistics

#Numerical
a <- df$prevasleg
y <- df$TLIF
hist(a)
#TLIF
mean(a[y==1], na.rm = T)
sd(a[y==1], na.rm = T)
#PLIF
mean(a[y==0], na.rm = T)
sd(a[y==0], na.rm = T)
#Test
library(exactRankTests)
exactRankTests::wilcox.exact(a~y, paired = F)

#categorical
b <- df$mcid
y <- df$TLIF
tab <- table(b,y)
tab
round(prop.table(tab,2) * 100, 1)
chisq.test(tab, correct = F)


#bOXPLOT
A1 <- df$changevasback[y==1]
A2 <- df$changevasback[y==0]
B1 <- df$changevasleg[y==1]
B2 <- df$changevasleg[y==0]
C1 <- df$changeodi[y==1]
C2 <- df$changeodi[y==0]
DF <- data.frame(A1, A2, B1, B2, C1, C2)
jpeg("Figure1.jpeg", units = "in", height = 5, width = 7, res = 1200)
boxplot(DF, col = rainbow(2, s = 0.3, v = 1), at = c(1:2,4:5,7:8), xaxt = "n", ylim = c(-100,100), xlab = "Outcome Measure", ylab = "Change Score (%)")
axis(side = 1, at = c(1.5,4.5,7.5), labels = c("VAS Back Pain","VAS Leg Pain", "ODI"))
abline(a = c(-30,0), lwd = 2, lty = "dashed")
legend("topleft", fill = rainbow(2, s = 0.5), legend = c("TLIF","PLIF"), horiz = T)
dev.off()