R/MLE.se.logitp.R

In shimlab/HMTree: Hidden Markov tree prior for analysis of multiple functional data

#' 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))
  }
}