set.seed(24)

library(MASS)
nsamp <- 100
SMat <- matrix(c(1.0, 0.6, 0.3, 0.6, 1.0, -0.3, 0.3, -0.3, 1.0), nrow=3, ncol=3)
datmat1 <- mvrnorm(n = nsamp, mu = rep(-1.0, 3), Sigma = SMat)
datmat2 <- mvrnorm(n = nsamp, mu = rep(1.0, 3), Sigma = SMat)
datmat3 <- mvrnorm(n = nsamp, mu = c(0.0,0.5,-1.0), Sigma = SMat)
datmat <- rbind(datmat1, datmat2, datmat3)
groups <- factor(c(rep("gene1", nsamp), rep("gene2", nsamp), rep("gene3", nsamp)))
dat <- data.frame(groups, datmat)
names(dat) <- c("groups","marker1","marker2","marker3")
pairs(dat)
pairs(dat, col=1+as.numeric(dat$groups))


# 3D Scatterplot
install.packages("scatterplot3d")
library(scatterplot3d)
s3d <- scatterplot3d(dat$marker1,dat$marker2,dat$marker3, xlab="marker 1", ylab = "marker 2", zlab = "marker 3", pch=19, color = 1+as.numeric(dat$groups), main="3D Scatterplot")
legend(s3d$xyz.convert(3.5, 2, 3), levels(dat$groups), col = 2:4, pch = 19)

#### Logistic regression

dat2 <- subset(dat, subset=groups!="gene3")
dat2 <- droplevels(dat2)

glm1 <- glm(groups ~ marker1 + marker2 + marker3, family="binomial", data=dat2)

summary(glm1)
confint(glm1)

exp(coefficients(glm1))
exp(confint(glm1))

#########  LDA ####

library(MASS)

ldafit <- lda(groups ~ marker1 + marker2 + marker3, data=dat)
ldafit

ldapred <- predict(ldafit)
names(ldapred)
head(ldapred$class)
head(ldapred$posterior)
head(ldapred$x)

table(ldapred$class,dat$groups)
mean(ldapred$class==dat$groups)


set.seed(15)
datmat1b <- mvrnorm(n = nsamp, mu = rep(-1.0, 3), Sigma = SMat)
datmat2b <- mvrnorm(n = nsamp, mu = rep(1.0, 3), Sigma = SMat)
datmat3b <- mvrnorm(n = nsamp, mu = c(0.0,0.5,-1.0), Sigma = SMat)
datmatb <- rbind(datmat1b, datmat2b, datmat3b)

dat3 <- data.frame(groups, datmatb)
names(dat3) <- c("groups","marker1","marker2","marker3")


ldapred3 <- predict(ldafit,newdata=dat3)
names(ldapred3)
head(ldapred3$class)
head(ldapred3$posterior)
head(ldapred3$x)

table(ldapred3$class,dat3$groups)
mean(ldapred3$class==dat3$groups)