Supervised Learning Continued - What Could Go Wrong?¶

## Supervised

# Simulate noisy data matrix (EXPRS)
set.seed(123)
# We'll use 2 groups of 20 subjects - think 20 cases and 20 controls
n=20
# Simulate 1000 genes
m=1000

# randomly generate a matrix of 'expression levels' -- any continuous variable we may be interested in
EXPRS=matrix(rnorm(2*n*m),2*n,m)

# Just naming rows and columns
rownames(EXPRS)=paste("patient",1:(2*n),sep="")
colnames(EXPRS)=paste("gene exp",1:m,sep="")

# The group labels are assigned arbitrarily - i.e. we are just randomly assigning
#                case/control status with no reference to gene expression
grp=rep(0:1,c(n,n))

#Pick the top 10 features based on the
#two-sample $t$-test

# load the library 'genefilter' - part of the Bioconductor suite
library(genefilter)

# rowttests is a genefilter function that performs a student t test on rows.

ttest.data=rowttests(t(EXPRS), factor(grp))
head(ttest.data)

# Extract the absolute value of the statistic
stats=abs(ttest.data$statistic)

# 'order' will return the indices of ascending values of it's argument
ii=order(-stats)

#Filter out all genes (i.e. columns of the matrix) except the top 10 rated by stats

TOPEXPRS=EXPRS[, ii[1:10]]
dim(TOPEXPRS)

plot(TOPEXPRS, col=grp+1,
      main="Gene Expression")
legend(-2, 2, c("Case", "Control"), pch=1, col=1:2)

# Fit 3-NN
library(class)
mod0=knn(train=TOPEXPRS,test=TOPEXPRS,cl=grp,k=3)
print(mod0)

# Error Resubstituion
table(mod0,grp)

# Naive CV
mod1=knn.cv(TOPEXPRS,grp,k=3)
table(mod1,grp)

# Oh! Super fancy code! I'll rewrite if there is time. Otherwise, we need to skip through.

# Proper CV
top.features=function(EXP,resp,test,fsnum)
  {
    top.features.i=function(i,EXP,resp,test,fsnum)
      {
        stats=abs(mt.teststat(EXP[,-i],resp[-i],test=test))
        ii=order(-stats)[1:fsnum]
        rownames(EXP)[ii]
      }
    sapply(1:ncol(EXP),top.features.i,EXP=EXP,resp=resp,test=test,fsnum=fsnum)
  }

# This function evaluates the knn

knn.loocv=function(EXP,resp,test,k,fsnum,tabulate=FALSE,permute=FALSE)
  {
    if(permute)
      resp=sample(resp)
    topfeat=top.features(EXP,resp,test,fsnum)
    pids=rownames(EXP)
    EXP=t(EXP)
    colnames(EXP)=as.character(pids)
    knn.loocv.i=function(i,EXP,resp,k,topfeat)
      {
        ii=topfeat[,i]
        mod=knn(train=EXP[-i,ii],test=EXP[i,ii],cl=resp[-i],k=k)[1]
      }
    out=sapply(1:nrow(EXP),knn.loocv.i,EXP=EXP,resp=resp,k=k,topfeat=topfeat)
    if(tabulate)
      out=ftable(pred=out,obs=resp)
    return(out)
}

library(multtest)
knn.loocv(t(EXPRS),as.integer(grp),"t.equalvar",3,10,TRUE)