add plrs package

git-svn-id: https://hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/plrs@71140 bc3139a8-67e5-0310-9ffc-ced21a209358

DESCRIPTION

@@ -0,0 +1,15 @@
+Package: plrs
+Version: 0.99.0
+Type: Package
+Title: Piecewise Linear Regression Splines (PLRS) for the association
+        between DNA copy number and gene expression
+Author: Gwenael G.R. Leday
+Maintainer: Gwenael G.R. Leday to <gleday@few.vu.nl>
+Depends: R (>= 2.10), Biobase
+Imports: BiocGenerics, CGHbase, graphics, grDevices, ic.infer, marray,
+         methods, quadprog, Rcsdp, stats, stats4, utils
+Suggests: mvtnorm, methods
+Description: The present package implements a flexible framework for
+    modeling the relationship between DNA copy number and gene
+    expression data using Piecewise Linear Regression Splines (PLRS).
+License: GPL (>=2.0)

NAMESPACE

@@ -0,0 +1,43 @@
+import(methods)
+
+importMethodsFrom(marray, lines, points)
+
+importFrom(BiocGenerics, residuals)
+importFrom(Biobase, exprs)
+importFrom(CGHbase, segmented, probloss, probnorm, probgain, probamp)
+importFrom(graphics, abline, box, polygon)
+importFrom(grDevices, dev.off, jpeg, pdf)
+importFrom(ic.infer, ic.weights, pbetabar)
+importFrom(quadprog, solve.QP)
+importFrom(Rcsdp, csdp)
+importFrom(stats, lm, p.adjust, pf, quantile, rnorm, runif, vcov)
+importFrom(stats4, coef, summary, plot)
+importFrom(utils, combn, head)
+
+export("plrs.sim","plrs.cb","plrs.series")
+exportClasses(
+     "plrs",
+     "plrs.select",
+     "plrs.series" 
+)
+exportMethods(
+    "coef",
+    "fitted",
+    "effects",
+    "knots",
+    "model.matrix",
+    "plot",
+    "predict",
+    "print",
+    "residuals",
+    "show",
+    "summary"
+)
+export(
+	"criteria",
+	"modify.conf",
+	"plrs",
+	"plrs.select",
+	"plrs.test"
+)
+

R/AllClasses.r

@@ -0,0 +1,33 @@
+################
+## All Classes
+################
+
+setClass("plrs",
+         representation(coefficients = "numeric",
+                        fitted.values = "numeric",
+                        residuals = "numeric",
+                        X = "matrix",
+                        data = "list",
+                        mdata = "list",
+                        QP = "list",
+                        test = "list",
+                        cb = "list",
+                        selected = "logical",
+                        type = "character",
+                        call.arg = "list")
+)
+setClass("plrs.select",
+         representation(table = "matrix",
+                        model = "plrs",
+                        crit = "character")
+)
+setClass("plrs.series",
+         representation(coefficients = "matrix",
+                        effects = "list",
+                        test = "matrix",
+                        general = "matrix",
+                        modelsType = "list",
+                        call.arg = "list")
+)
+
+

R/AllGenerics.r

@@ -0,0 +1,6 @@
+################
+## All Generics
+################
+
+setGeneric("fitted",    function(object,...) standardGeneric("fitted"))
+setGeneric("predict",   function(object,...) standardGeneric("predict"))

R/Bmat.r

@@ -0,0 +1,47 @@
+# Internal function
+# Build design matrix with truncated power basis functions
+
+# Author: Gwenael G.R. Leday
+
+.Bmat <- function(cghseg, knots, continuous=TRUE, general.intercept=TRUE){
+
+	nKnots <- length(knots)
+	nval <- nKnots + 1
+	xx <- cghseg
+	mat <- matrix(cghseg)
+	if(general.intercept) mat <- cbind(rep(1, length(xx)) ,mat)
+
+	# Truncated p-th power function
+	tpower <- function(x, t, p) (x - t) ^ p * (x > t)
+
+	if(continuous){
+		if(nKnots>=1) mat <- cbind(mat, tpower(xx, knots[1], 1))
+		if(nKnots>=2) mat <- cbind(mat, tpower(xx, knots[2], 1))
+		if(nKnots>=3) mat <- cbind(mat, tpower(xx, knots[3], 1))
+		nval2 <- 0
+	}
+	else{
+		if(nKnots>=1){
+			mat <- cbind(mat, tpower(xx, knots[1], 0))
+			mat <- cbind(mat, tpower(xx, knots[1], 1))
+		}
+		if(nKnots>=2){
+			mat <- cbind(mat, tpower(xx, knots[2], 0))
+			mat <- cbind(mat, tpower(xx, knots[2], 1))
+		}
+		if(nKnots>=3){
+			mat <- cbind(mat, tpower(xx, knots[3], 0))
+			mat <- cbind(mat, tpower(xx, knots[3], 1))
+		}
+		nval2 <- nval - 1
+	}
+
+	# Labels
+	theta.1 <- paste(paste(paste("theta",0:(nval-1),sep=""),".",sep=""),1,sep="")
+	if(!general.intercept & nval2==0) theta.0 <- NULL
+	else{ theta.0 <- paste(paste(paste("theta",0:nval2,sep=""),".",sep=""),0,sep="")}
+	colnames(mat) <- sort(c(theta.0, theta.1))
+
+	return(mat)
+}
+

R/convertToTime.r

@@ -0,0 +1,13 @@
+# Internal function
+# Convert seconds into a "HH:MM:SS" format
+
+# Author: Gwenael G.R. Leday
+
+.convertToTime <- function(x){
+	h <- as.character(x%/%3600)
+	m <- as.character((x%%3600)%/%60)
+	s <- as.character(round((x%%3600)%%60))
+	if(nchar(m)==1) m <- paste(0,m,sep="")
+	if(nchar(s)==1) s <- paste(0,s,sep="")
+	return(paste(h,m,s,sep=":"))
+}

R/criteria.r

@@ -0,0 +1,77 @@
+# Extract AIC, AICC, BIC and OSAIC from a 'plrs' object
+
+# Author: Gwenael G.R. Leday
+
+criteria <- function(obj, crit = "all"){
+
+	# Error handling
+	if(class(obj)!="plrs") stop("An object of class \"plrs\" is required")
+	if(length(crit)==0) stop("No criterion")
+	if(length(crit)==1) if(crit=="all") crit <- c("aic","aicc","bic","osaic")
+	if(!is.character(crit)) stop("'crit' has to be a character")
+	if(!"all"%in%crit & !"aic"%in%crit & !"aicc"%in%crit &
+         !"bic"%in%crit & !"osaic"%in%crit) stop("Wrong criterion")
+	bool <- TRUE
+	if("osaic"%in%crit){
+		if(!obj@call.arg$constr) bool <- FALSE
+		else{
+			if(length(coef(obj))!=1){
+				if(obj@QP$meq == nrow(t(obj@QP$Amat))) bool <- FALSE
+			}
+			else{
+				bool <- FALSE
+			}
+		}
+	}
+
+	rss <- sum(residuals(obj)^2)
+	n <- length(obj@data$expr)
+
+	# Compute criteria for the segmented model
+	if(rss != 0){
+		if(obj@call.arg$constr & length(coef(obj))>1) p <- length(coef(obj)) - obj@QP$meq
+		else p <- length(coef(obj))
+
+		if("aic" %in% crit || "aicc" %in% crit || !bool) aic <- n*log(rss/n) + 2*p
+		if("aicc" %in% crit) aicc <- aic + ((2*p)*(p+1))/(n-p-1)
+		if("bic" %in% crit) bic <- n*log(rss/n) + log(n)*p
+		if("osaic" %in% crit){ # As defined in Hughes et al (2003)
+			if(bool){
+				library(mvtnorm)
+				# Calculate weights
+				obj.unconstr <- .plrs.fit(X=obj@X, matconstr=NULL, mytype=obj@type, modelFull=obj, callArg=obj@call.arg)
+				cov.mat <- obj.unconstr@QP$vcov/obj.unconstr@QP$S2
+				uiw <- t(obj@QP$Amat)
+
+				# Aknowledgement: Taken from package ic.infer
+				if(obj@QP$meq==0){
+					#wt.bar <- ic.weights(uiw %*% cov.mat %*% t(uiw))
+					wt.bar <- ic.weights(crossprod(t(crossprod(t(uiw),cov.mat)),t(uiw)))
+				}
+          	 	 	else{
+					wt.bar <- ic.weights(solve(solve(crossprod(t(crossprod(t(uiw),cov.mat)),t(uiw)))[-(1:obj@QP$meq), -(1:obj@QP$meq)]))
+				}
+				wt.bar <- wt.bar[sort(names(wt.bar))]
+				r <- (length(wt.bar)-1)
+				osaic.npar <- sum(wt.bar*(p-r+0:r))
+				osaic <- n*log(rss/n) + 2*osaic.npar
+			}
+			else{
+				osaic <- aic
+			}
+		}
+	}
+	else{
+		if("aic" %in% crit) aic <- -Inf
+		if("aicc" %in% crit) aicc <- -Inf
+		if("bic" %in% crit) bic <- -Inf
+		if("osaic" %in% crit) osaic <- -Inf
+	}
+
+	out <- NULL
+	if("aic" %in% crit) out <- c(out, list("aic" = aic))
+	if("aicc" %in% crit) out <- c(out, list("aicc" = aicc))
+	if("bic" %in% crit) out <- c(out, list("bic" = bic))
+	if("osaic" %in% crit) out <- c(out, list("osaic" = osaic))
+	return(out)
+}

R/estimateKnots.r

@@ -0,0 +1,46 @@
+# Estimate knots
+
+# Author: Gwenael G.R. Leday
+
+.estimateKnots <- function(cghseg, cghcall, probloss=NULL, probnorm=NULL, probgain=NULL, probamp=NULL){
+
+	val <- sort(unique(cghcall))
+	nval <- length(val)
+	knots <- NULL
+
+	# Midpoints of intermediate intervals
+	if(is.null(probloss) || is.null(probnorm) || is.null(probgain)){
+		if(nval>1){
+			for(k in 1:(nval-1)){
+				knots <- c(knots,(max(cghseg[cghcall==val[k]]) +  min(cghseg[cghcall==val[k+1]]))/2)
+			}
+		}
+	}
+	# Estimate using call probabilities
+	else{
+		p <- function(alpha, cghseg, leftState, rightState){
+			ind1 <- as.numeric(cghseg <= alpha)
+			ind2 <- 1-ind1
+			as.numeric(ind1%*%leftState + ind2%*%rightState)
+		}
+		for(j in 1:(length(val)-1)){
+			# Reduce search to adjacent states range
+			lb <- min(cghseg[cghcall==val[j]])
+			ub <- max(cghseg[cghcall==val[j+1]])
+
+			# Search maximum
+			xx <- seq(from=lb, to=ub, length.out=1000)
+			if(val[j]==-1) tomax <- sapply(xx, p, cghseg=cghseg, leftState=probloss, rightState=probnorm)
+			if(val[j]==0) tomax <- sapply(xx, p, cghseg=cghseg, leftState=probnorm, rightState=probgain)
+			if(val[j]==1) tomax <- sapply(xx, p, cghseg=cghseg, leftState=probgain, rightState=probamp)
+
+			# Take midpoints of interval
+			optvals <- xx[which(tomax==max(tomax))]
+			leftpts <- max(cghseg[cghseg<=min(optvals)])
+			rightpts <- min(cghseg[cghseg>min(optvals)])
+
+			knots <- c(knots, mean(c(leftpts, rightpts)))
+		}
+	}
+	return(knots)
+}