#' Take a series of Poisson count signals x, with data on different samples
#' in each row, and compute estimate for logit(p) and its standard error
#' optionally estimate the "effect" of a covariate g if g is provided.
MLE.se.logitp <- function(x, g=NULL, reflect=FALSE, minobs=1, pseudocounts=0.5, all=FALSE, center=FALSE, repara=TRUE, forcebin=FALSE, lm.approx=TRUE, disp="add"){
if(!is.numeric(x)) stop("Error: invalid parameter 'x': 'x' must be numeric")
if(!is.null(g)) if(!(is.numeric(g)|is.factor(g))) stop("Error: invalid parameter 'g', 'g' must be numeric or factor or NULL")
if(!is.logical(center)) stop("Error: invalid parameter 'center', 'center' must be bool")#need this
if(!is.logical(reflect)) stop("Error: invalid parameter 'reflect', 'reflect' must be bool")
if(!(((minobs%%1)==0)&minobs>0)) stop("Error: invalid parameter 'minobs', 'minobs' must be positive integer")
if(!(is.numeric(pseudocounts)&pseudocounts>0)) stop("Error: invalid parameter 'pseudocounts', 'pseudocounts' must be a positive number")
if(!is.logical(all)) stop("Error: invalid parameter 'all', 'all' must be bool")
if(!is.logical(repara)) stop("Error: invalid parameter 'repara', 'repara' must be bool")
if(!is.logical(forcebin)) stop("Error: invalid parameter 'forcebin', 'forcebin' must be bool")
if(!is.logical(lm.approx)) stop("Error: invalid parameter 'lm.approx', 'lm.approx' must be bool")
if(!is.element(disp,c("add","mult"))) stop("Error: invalid parameter 'disp', 'disp' must be either 'add' or 'mult' ")
if(is.vector(x)){dim(x)<- c(1,length(x))} #change x to matrix
nsig = nrow(x)
if(!is.null(g)) if(length(g)!=nsig) stop("Error: covariate g for all samples are not provided")
if(nsig==1){forcebin=TRUE} #if only one observation, don't allow overdispersion
J = log2(ncol(x)); if((J%%1)!=0){reflect=TRUE} #if ncol(x) is not a power of 2, reflect x
if(reflect==TRUE) reflect.indices=reflectSignal(x) #reflect signal; this function is pseudo calling x by reference
if(!is.null(g)){
if(is.factor(g))
g.num = as.numeric(levels(g))[g]
else{
g.num=g
if(length(unique(g))==2)
g=factor(g)
}
#weights for quantitative covariate
if(center==TRUE)
w=unique(sort(g.num-mean(g.num)))
else
w=c(0,1)
}
y = haar.aggregate(x)
zdat=glm.approx(y, g, minobs=minobs, pseudocounts=pseudocounts, center=center, all=all, forcebin=forcebin, repara=repara, lm.approx=lm.approx, disp=disp)
if(is.null(g)){
return(list(alpha = zdat[1,], alpha.se = zdat[2,], beta = NULL, beta.se = NULL, cov.val = NULL, center=center, repara=repara, w=NULL))
}else{
return(list(alpha = zdat[1,], alpha.se = zdat[2,], beta = zdat[3,], beta.se = zdat[4,], cov.val = zdat[5,], center=center, repara=repara, w=w))
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.