R/Drda1.R In rda: Shrunken Centroids Regularized Discriminant Analysis

Defines functions Drda1

```Drda1 <- function(x, y, xnew=NULL, ynew=NULL, prior=table(y)/length(y),
alpha=seq(0.01, 1, len=10), delta=seq(0, 3, len=10),
genelist=FALSE, trace=FALSE){
this.call <- match.call()

testdata.x <- (is.null(xnew) || missing(xnew))
testdata.y <- (is.null(ynew) || missing(ynew))

mu <- apply(x, 1, mean)
x <- x-mu

yhat <- array(0, dim=c(length(alpha), length(delta), length(y)))

if (!(testdata.x)){
xnew <- xnew-mu
yhat.new <- array(0, dim=c(length(alpha), length(delta),
ncol(xnew)))
posterior <- array(0, dim=c(length(alpha), length(delta),
ncol(xnew), length(prior)))
}

if (genelist){
gene.list <- array(0, dim=c(length(alpha), length(delta), nrow(x)))
dn2 <- list(round(alpha, 3), round(delta, 3), seq(nrow(x)))
names(dn2) <- c("alpha", "delta", "gene #")
dimnames(gene.list) <- dn2
}

Y <- model.matrix( ~ factor(y)-1)
dimnames(Y) <- list(NULL, NULL)

xbar <- scale(x%*%Y, FALSE, table(y))

xstar <- x-xbar[, unclass(factor(y))]

DD <- sqrt(apply(xstar^2, 1, sum)/length(y))

xstar1 <- xstar/DD
xout <- t(xstar1)%*%xstar1
sout <- svd(xout)

svalue <- sout\$d
svaluePos <- seq(svalue)[svalue > 0]
svalue <- sqrt(svalue[svaluePos])

V <- scale(xstar1%*%sout\$u[, svaluePos], FALSE, svalue)
svalue <- svalue/sqrt(length(y))

xbar1 <- xbar/DD
Vtxbar <- t(V)%*%xbar1

error <- ngene <- testerror <- matrix(0, length(alpha), length(delta))
if (!(testdata.y)){
testerror <- error
}

for(i in seq(along=alpha)){
dd <- 1/(svalue^2*alpha[i]+(1-alpha[i]))-1/(1-alpha[i])
temp <- Vtxbar*dd

coefmat <- (V%*%temp+xbar1/(1-alpha[i]))/DD
for(j in seq(along=delta)){
coefmat1 <- I(abs(coefmat)>delta[j])*(abs(coefmat)-
delta[j])*sign(coefmat)

if(genelist){
gene.list[i, j, ] <- I(apply(I(abs(coefmat1) > 0), 1, sum) > 0)
}

ngene[i, j] <- sum(apply(I(abs(coefmat1) > 0), 1, sum) > 0)

coefmat2 <- coefmat1*DD
ystar <- t(V)%*%coefmat2
intvec <- -(alpha[i]*apply(ystar^2*svalue^2, 2, sum)+
(1-alpha[i])*apply(coefmat2^2, 2, sum))/2+log(prior)

dhat <- scale(t(x) %*% coefmat1, -intvec, FALSE)
dimnames(dhat) <- list(NULL, names(prior))
yhat[i, j, ] <- softmax(dhat)

error[i, j] <- sum(y != yhat[i, j, ])

if(!(testdata.x)){
dhat <- scale(t(xnew)%*%coefmat1, -intvec, FALSE)
posterior[i, j, , ] <- dhat
dimnames(dhat) <- list(NULL, names(prior))
yhat.new[i, j, ] <- softmax(dhat)
if(!(testdata.y)){
testerror[i, j] <- sum(ynew != yhat.new[i, j, ])
}
}

if(trace){
if(!(testdata.x || testdata.y)){
cat("alpha=", alpha[i], "delta=", delta[j], "test error=",
format(round(testerror[i, j], 3)), "\n")
}
else{
cat("alpha=", alpha[i], "delta=", delta[j], "error=",
format(round(error[i, j], 3)), "\n")
}
}
}
}
dn <- list(round(alpha, 3), round(delta, 3))
names(dn) <- c("alpha", "delta")
dimnames(error) <- dimnames(ngene) <- dn

dn3 <- list(round(alpha, 3), round(delta, 3), seq(ncol(x)))
names(dn3) <- c("alpha", "delta", "sample")
dimnames(yhat) <- dn3

obj <- list(alpha = alpha, delta = delta, prior=prior,
error=error, ngene=ngene, centroids=xbar,
centroid.overall=mu, call=this.call, yhat=yhat)
if(!(testdata.x)) {
dn4 <- list(round(alpha, 3), round(delta, 3), seq(ncol(xnew)))
names(dn4) <- c("alpha", "delta", "sample")
dimnames(yhat.new) <- dn4
obj\$yhat.new <- yhat.new

dn1 <- list(round(alpha, 3), round(delta, 3), seq(ncol(xnew)),
names(prior))
names(dn1) <- c("alpha", "delta", "sample #", "class")
dimnames(posterior) <- dn1
obj\$posterior <- posterior

if(!(testdata.y)) {
dimnames(testerror) <- dn
obj\$testerror <- testerror
}
}
if(genelist) {
obj\$gene.list <- gene.list
}
obj
}
```

Try the rda package in your browser

Any scripts or data that you put into this service are public.

rda documentation built on July 22, 2018, 1:01 a.m.