### R code from vignette source 'vignettes/lmdme/inst/doc/lmdme-vignette.Rnw'
### Encoding: UTF-8

###################################################
### code chunk number 1: Loading libraries
###################################################
library(lmdme)


###################################################
### code chunk number 2: Dataset loading
###################################################
library(stemHypoxia)
data(stemHypoxia) #This will load M and design objects in memory
timeIndex   <- design$time %in% c("0.5","1","5")  #time levels
oxygenIndex <- design$oxygen %in% c("1","5","21") #oxygen levels
design<-design[ timeIndex & oxygenIndex,]# Both time & oxygen
design$time  <-as.factor(design$time)  
design$oxygen<-as.factor(design$oxygen)
rownames(M)<-M[,1] #Gene ID as row.names of M
M <- M[,colnames(M) %in% design$samplename] #Just what is needed


###################################################
### code chunk number 3: Exploring M and design objects
###################################################
head(design)
head(M)[,1:3]


###################################################
### code chunk number 4: Exploring M and design objects
###################################################
library(lmdme)
fit <- lmdme(model=~time*oxygen,data=M,design=design)
fit


###################################################
### code chunk number 5: ASCA decomposition
###################################################
id<-F.p.values(fit,term="time:oxygen")[[1]]<0.001
sum(id) #The amount of genes for further exploration
decomposition(fit,decomposition="pca", type="coefficient", 
              term="time:oxygen", subset=id, scale="row")
biplot(fit,xlabs=rep("o",sum(id)), mfcol=NULL)


###################################################
### code chunk number 6: ASCA_biplot
###################################################
biplot(fit,xlabs=rep("o",sum(id)), mfcol=NULL)


###################################################
### code chunk number 7: PLS decomposition
###################################################
fit.plsr<-fit
decomposition(fit.plsr,decomposition="plsr", type="coefficient", 
              term="time:oxygen", subset=id,scale="row")
biplot(fit.plsr, which = "loadings", xlabs=rep("o",sum(id)), 
      ylabs=colnames(coefficients(fit.plsr,term="time:oxygen")[[1]]),
       var.axes=TRUE, mfcol=NULL)


###################################################
### code chunk number 8: PLS_biplot
###################################################
biplot(fit.plsr, which = "loadings", xlabs=rep("o",sum(id)), 
      ylabs=colnames(coefficients(fit.plsr,term="time:oxygen")[[1]]),var.axes=TRUE, mfcol=NULL)


###################################################
### code chunk number 9: Loadingplot
###################################################
loadingplot(fit,term.x="time",term.y="oxygen")


###################################################
### code chunk number 10: Session Info
###################################################
sessionInfo()