R/fanova.RPm.R

In moviedo5/fda.usc: Functional Data Analysis and Utilities for Statistical Computing

#' @rdname fanova.RPm
#' @aliases fanova.RPm summary.fanova.RPm anova.RPm
#' @note anova.RPm deprecated.
#' @title Functional ANOVA with Random Project.
#' 
#' @description The procedure is based on the analysis of randomly chosen one-dimensional
#' projections. The function tests ANOVA models for functional data with
#' continuous covariates and perform special contrasts for the factors in the
#' formula.
#' 
#' @details \code{zproj} allows to change the generator process of the projections. This
#' can be done through the inclusion of a function or a collection of
#' projections generated outside the function. By default, for a functional
#' problem, the function \code{\link{rproc2fdata}} is used. For multivariate problems,
#' if no function is included, the projections are generated by a normalized
#' gaussian process of the same dimension as \code{object}. Any user function
#' can be included with the only limitation that the two first parameters are:
#' \itemize{
#' \item \code{n}: number of projections
#' \item \code{t}: discretization points for functional problems
#' \item \code{m}: number of columns for multivariate problems.
#' }
#' That functions must return a \code{fdata} or \code{matrix} object
#' respectively.  
#' 
#' The function allows user-defined contrasts.  The list of contrast to be used
#' for some of the factors in the formula.  Each contrast matrix in the list
#' has \code{r} rows, where \code{r} is the number of factor levels. The user
#' can also request special predetermined contrasts, for example using the
#' \code{\link{contr.helmert}}, \code{\link{contr.sum}} or
#' \code{\link{contr.treatment}} functions.
#'  
#' The function returns (by default) the significance of the variables using
#' the Bonferroni test and the False Discovery Rate test. Bootstrap procedure
#' provides more precision
#' 
#' @param object Functional response data.  Object with class fdata with
#' \code{n} curves discretizated in \code{m} points. For multivariate problems
#' \code{object} can be a \code{data.frame} or a \code{matrix}
#' @param formula as \link[stats]{formula} without response.
#' @param data.fac Explanatory variables.  Data frame with dimension (\code{n}
#' x \code{p}), where \code{p} are the number of factors or covariates
#' considered.
#' @param RP Vector of number of random projections.
#' @param alpha Confidence level, by default \code{alpha}=0.95.
#' @param zproj Function for generating the projections or an object that
#' contains that projections.
#' @param par.zproj List of parameters for \code{zproj} function.
#' @param hetero \code{logical}. If \code{TRUE} (by default) means heteroskedastic ANOVA.
#' @param pr \code{logical}. If \code{TRUE} prints intermediate results.
#' @param w Vector of weights (only for multivariate problems).
#' @param nboot Number of bootstrap samples, by default no bootstrap
#' computations, \code{nboot}=0.
#' @param contrast List of special contrast to be used ; by default no special
#' contrasts are used (\code{contrast}=\code{NULL}).
#' @param ndec Number of decimals.
#' @param \dots Further arguments passed to or from other methods.
#' 
#' @return An object with the following components:\cr 
#' \itemize{
#' \item \code{proj}: The projection value of each point on the curves. Matrix with dimensions
#' (\code{RP} x \code{m}), where \code{RP} is the number of projections and
#' \code{m} is the number of points observed in each projection curve. 
#' \item \code{mins}: Minimum number for each random projection. 
#' \item \code{result}: p-value for each random projection. 
#' \item \code{test.Bonf}: Significance (TRUE or FALSE) for
#' vector of random projections \code{RP} in columns and factors (and special
#' contrasts) in rows. 
#' \item \code{p.Bonf}: p-value for vector of random projections
#' \code{RP} in columns and factors (and special contrasts) in rows.
#' \item \code{test.fdr}: False Discovery Rate (TRUE or FALSE) for vector of random
#' projections \code{RP} in columns and factors (and special contrasts) in rows.
#' \item \code{p.fdr}: p-value of False Discovery Rate for vector of random
#' projections \code{RP} in columns and factors (and special contrasts) in rows.
#' \item \code{test.Boot}: False Discovery Rate (TRUE or FALSE) for vector of random
#' projections \code{RP} in columns and factors (and special contrasts) in rows.
#' \item \code{p.Boot}: p-value of Bootstrap sample for vector of random projections
#' \code{RP} in columns and factors (and special contrasts) in rows.
#' }
#' @note If \code{hetero=TRUE} then all factors must be categorical.
#' @author Juan A. Cuesta-Albertos, Manuel Febrero-Bande, Manuel Oviedo de la
#' Fuente\cr \email{manuel.oviedo@@udc.es}
#' @seealso See Also as: \code{\link{fanova.onefactor}}
#' @references Cuesta-Albertos, J.A., Febrero-Bande, M. \emph{A simple multiway
#' ANOVA for functional data.} TEST 2010, DOI \bold{10.1007/s11749-010-0185-3}.
#' @keywords anova
#' @examples
#' \dontrun{
#' # ex fanova.hetero
#' data(phoneme)
#' names(phoneme)
#' # A MV matrix obtained from functional data
#' data=as.data.frame(phoneme$learn$data[,c(1,seq(0,150,10)[-1])]) 
#' group=phoneme$classlearn
#' n=nrow(data)
#' group.rand=as.factor(sample(rep(1:3,len=n),n))
#' RP=c(2,5,15,30)
#' 
#' #ex 1: real factor and random factor
#' m03=data.frame(group,group.rand)
#' resul1=fanova.RPm(phoneme$learn,~group+group.rand,m03,RP=c(5,30))
#' summary(resul1)
#' 
#' #ex 2: real factor with special contrast
#' m0=data.frame(group)
#' cr5=contr.sum(5)   #each level vs last level
#' resul03c1=fanova.RPm(data,~group,m0,contrast=list(group=cr5))
#' summary(resul03c1)
#' 
#' #ex 3: random factor with special contrast. Same projs as ex 2.
#' m0=data.frame(group.rand)
#' zz=resul03c1$proj
#' cr3=contr.sum(3)   #each level vs last level
#' resul03c1=fanova.RPm(data,~group.rand,m0,contrast=list(group.rand=cr3),zproj=zz)
#' summary(resul03c1)
#' }
#'   
#' @export 
fanova.RPm=function(object,formula,data.fac,RP=min(30,ncol(object)),alpha=0.95,
                   zproj=NULL,par.zproj=list(norm=TRUE),hetero=TRUE,
                   pr=FALSE,w=rep(1,ncol(object)),nboot=0,contrast=NULL,...){
  if (class(object)[1] %in% c("matrix","data.frame")) dataM=data.matrix(object)
  else if (is.fdata(object)) dataM=object[["data"]]
  lfac=unlist(lapply(data.fac,is.factor))
  min.data.fac<-min(table(data.fac[,lfac]))
  if (min.data.fac==0)  warning("Contingency table of factor levels (data.fac argument) contains 0 counts  values")
  nrow=nrow(dataM);ncol=ncol(dataM)
  bonf=(1-alpha)/RP
  nprRP=max(RP)
  terms.fd=attr(terms(formula),"term.labels")
  fml=as.formula(paste("value ~ ", paste(terms.fd, collapse= "+")))
  if (is.null(nrow) || is.null(ncol)) stop("fdata must be a matrix")
  nterms=length(terms.fd)+1
#  
projMV=function(n,m,w=rep(1,m),norm=TRUE){
	modulo=function(z){sqrt(sum(z^2))}
	z=rnorm(n*m)
	z=matrix(z,nrow=n,ncol=m)
	z=t(t(z)*w)
	if (norm) {
	modu=apply(z,1,modulo)
	z=z/modu}
}
if ((class(object) %in% c("matrix","data.frame")) & is.null(zproj)) zproj=projMV  
if ((is.fdata(object)) & is.null(zproj)) zproj=rproc2fdata  
  if (is.fdata(zproj) & is.fdata(object)) { 
	if (nrow(zproj)>=nprRP) { 
			z=zproj[1:nprRP] 
			} else {
			stop(paste("Not enough functions in zproj",length(zproj)," to compute ",nprRP," projections"))  
			}
			} else if (is.matrix(zproj) & (class(object)[1] %in% c("matrix","data.frame"))){
	if (nrow(zproj)>=nprRP){
			z=zproj[1:nprRP,] } else {
			stop(paste("Not enough rows in zproj",nrow(zproj)," to compute ",nprRP," projections"))
			}
			}  else if (is.function(zproj)){
			if (is.fdata(object)) {z=do.call(zproj,modifyList(list(n=nprRP,t=object$argvals),par.zproj))}
			 else if (class(object)[1] %in% c("matrix","data.frame")) {z=do.call(zproj,modifyList(list(n=nprRP,m=ncol(object)),par.zproj))}
			} else {stop("Parameter zproj is neither an fdata object or a function")}

  ff=attr(terms.fd,"factors")
  ffcol=colnames(ff)
  if (is.null(contrast)){
    mat=matrix(NA,ncol=(nterms-1),nrow=nprRP)
    colnames(mat)=c(terms.fd)
    ncontrast=0
  }
  else {
     b=length(contrast)
     mf2=model.frame(formula,data.fac)
     uniq=apply(mf2,2,function(x){length(unique(x))})
     ncontrast=rep(0,len=b)
     name.contr=rep(NA,len=sum(ncontrast))
     bb=length(uniq)
     cnombres=ncontrast=rep(0,len=b)
     tnombres=rep(FALSE,len=length(ncontrast))
     tgroups=rep(FALSE,len=(length(ffcol)+1))
if (pr) {print(contrast);print("contrast     contrast  contrast")}
     for (i in 1:b)    {
       a=which(names(contrast)[i]==colnames(mf2))
        tgroups[a]=tnombres[a]=TRUE
         if (is.vector(contrast[[i]])) {
              if (pr) print("Is vector")
              print(a)
                          ncontrast[a]=1
                          contrast[[i]]=matrix(contrast[[i]],ncol=1)
                          }
       else ncontrast[a]=ncol(contrast[[i]])
       names(ncontrast)[a]=names(contrast[i])
                       }
            j=1;ji=1;jk=1
            for (i in 1:length(ncontrast))    {
            if (tgroups[ji])  {
                name.contr[j:(j+ncontrast[i]-1)]=paste("C",j:(j+ncontrast[i]-1),".",
                names(contrast[jk]),sep="")
               colnames(contrast[[jk]])=name.contr[j:(j+ncontrast[i]-1)]
              j=j+ncontrast[i];jk=jk+1
              }
              ji=ji+1
              }
     mat=matrix(NA,ncol=(nterms+sum(ncontrast)-1),nrow=nprRP)
     colnames(mat)=c(terms.fd,name.contr)
if (pr) {
   print(ncontrast);print("ncontrast")
   print(contrast);print("contrast")
   print(name.contr);print("name.contr")
   }
  }
  for (j in 1:nprRP){
#    value=data%*%z[j,]
	if (is(object,"fdata")) {
	value=inprod.fdata(object,z[j]) 
	} else if (is(object,"data.frame") | is(object,"matrix")) {
	value=dataM%*%z[j,]
	}
    mdata=as.data.frame(cbind(value,data.fac))
    colnames(mdata)=c("value",colnames(data.fac))
    result=aov(fml,data=mdata)
    out=summary(result)
    if (!hetero){
       result3=lm(fml,data=mdata,contrasts=contrast,...)
       if (pr) {print(result);print(result3)}
       if (!is.null(contrast)) {
          out4=summary(result3)
          if (pr) {print(out4);print("summary lm contrast")}
          if (length(out)==1) {
            mat[j,1:(nterms-1)]=out[[1]][1:(nterms-1),5]
            ind=nterms;ind2=2
            for (i in 1:bb)    {
                mat[j,ind:(ind+ncontrast[i]-1)]=out4$coefficients[ind2:(ind2+ncontrast[i]-1),4]
                ind=ind+ncontrast[i];ind2=ind2+ncontrast[i]-1
         }        }
        if (pr) {print(mat);print("mat")}}
    else {  out=summary(result)
           if (length(out)==1) mat[j,1:(nterms-1)]=out[[1]][1:(nterms-1),5] }
    }
  else {              ###############################
      out=fanova.hetero(object=mdata,fml,pr=pr,contrast=contrast)[[1]]
      mat[j,]=out[,4]
      if (pr) print(out)
      if (!is.null(contrast)) mat[j,nterms]=out[nterms,4]
    }
  }
  if (pr) {print(summary(mat));print(paste("Bonferroni:",round(bonf,6)))}
  tmat=matrix(NA,nrow=length(RP),ncol=ncol(mat));colnames(tmat)=colnames(mat);rownames(tmat)=paste("RP",RP,sep="")
  mins=matrix(NA,nrow=length(RP),ncol=ncol(mat));colnames(mins)=colnames(mat);rownames(mins)=rownames(tmat)
  tFDR=matrix(NA,nrow=length(RP),ncol=ncol(mat));colnames(tFDR)=colnames(mat);rownames(tFDR)=rownames(tmat)
  pFDR=matrix(NA,nrow=length(RP),ncol=ncol(mat));colnames(pFDR)=colnames(mat);rownames(pFDR)=rownames(tmat)
  for (l in 1:length(RP)){
      if (RP[l]==1)      tmat[l,]=mat[1,]
      else {
      tmat[l,1:(nterms-1+sum(ncontrast))]=
      apply(matrix(mat[1:RP[l],1:(nterms-1+sum(ncontrast))],ncol=(nterms-1+sum(ncontrast))),2,min)}
      if (RP[l]==1) mins[l,]=rep(1,len=ncol(mat))
      else {
      mins[l,]=apply(matrix(mat[1:RP[l],1:(nterms-1+sum(ncontrast))],ncol=(nterms-1+sum(ncontrast))),2,which.min)
      }
      if (RP[l]==1) tFDR[l,]=(mat[1,]<(1-alpha))
      else {
      tFDR[l,1:(nterms-1+sum(ncontrast))]=
      apply(matrix(mat[1:RP[l],1:(nterms-1+sum(ncontrast))],ncol=(nterms-1+sum(ncontrast))),2,FDR,alpha=alpha)}
      if (RP[l]==1) pFDR[l,]=mat[1,]
      else {
      pFDR[l,1:(nterms-1+sum(ncontrast))]=
      apply(matrix(mat[1:RP[l],1:(nterms-1+sum(ncontrast))],,ncol=(nterms-1+sum(ncontrast))),2,pvalue.FDR)}
  }
  tbonf = (tmat < matrix(bonf, nrow = length(RP), ncol = ncol(mat)))
  colnames(tbonf) = colnames(mat);rownames(tbonf) = rownames(tmat)
  pbonf=tmat*RP
  colnames(pbonf) = colnames(mat);rownames(pbonf) = rownames(tmat)
  if (is.null(contrast)) {
   resb=matrix(NA,nrow=length(RP),ncol=nterms-1)    }
  else { resb=matrix(NA,nrow=length(RP),ncol=(nterms-1+sum(ncontrast)))}
  if (pr) {print(tbonf);print(pbonf);print(tFDR);print(pFDR)}
  if (nboot>0){
     if (pr) print("bootstrap procedure") ############
     resboot=fanova.RPm.boot(dataM,formula,data.fac,RP=RP,alpha=alpha,
                            nboot=nboot,zproj=z,hetero=hetero,
                            contrast=contrast,pr=pr,...)
     if (pr) print(resboot)
     if (is.null(contrast)){nbuc=nterms-1}
     else {nbuc=nterms-1+sum(ncontrast)}
     for (l in 1:length(RP)){
         for (k in 1:nbuc){
          if ((nterms-1)==1)    Fn=ecdf(resboot[[l]][k])
          else   Fn=ecdf(resboot[[l]][,k])
          resb[l,k]=Fn(tmat[l,k])
          }}
     rownames(resb)=paste("RP",RP,sep="")
  }
pFDR[pFDR>1]=1
pbonf[pbonf>1]=1
res=list("proj"=z,"mins"=mins,"result"=mat,
"test.Bonf"=tbonf,"p.Bonf"=pbonf,"test.FDR"=tFDR,"p.FDR"=pFDR)
if (nboot>0) {
   resb[resb>1]=1
   res$test.Boot=(resb<(1-alpha));res$p.Boot=resb
  colnames(res$p.Boot)=colnames(mat);colnames(res$test.Boot)=colnames(mat)}
if (!is.null(contrast)) res$contrast <- contrast
class(res) <- "fanova.RPm"
return(res)
}