Hmsc.R

#' @title Hmsc
#'
#' @description Creates an \code{Hmsc}-class object
#'
#' @param Y a matrix of species occurences or abundances
#' @param XFormula a \code{\link{formula}}-class object for fixed effects
#'    (linear regression)
#' @param XData a data frame of measured covariates for fixed effects with
#'    \code{\link{formula}}-based specification
#' @param X a matrix of measured covariates for fixed effects with direct specification
#' @param XScale a boolean flag indicating whether to scale covariates for the fixed effects
#' @param XSelect a list describing how variable selection is to be applied
#' @param XRRRData a data frame of covariates for reduced-rank regression
#' @param XRRRFormula \code{\link{formula}} for reduced-rank regression
#' @param XRRR a matrix of covariates for reduced-rank regression
#' @param ncRRR number of covariates (linear combinations) for reduced-rank regression
#' @param XRRRScale a boolean flag indicating whether to scale covariates for reduced-rank regression
#' @param YScale a boolean flag whether to scale responses for which normal distribution is assumed
#' @param studyDesign a data frame of correspondence between sampling units and units on different levels of latent
#'   factors
#' @param ranLevels a named list of \code{HmscRandomLevel}-class objects, specifying the structure and data for random
#'   levels
#' @param ranLevelsUsed a vector with names of levels of latent factors that are used in the analysis
#' @param TrFormula a \code{\link{formula}}-class object for regression
#'    dependence of \eqn{\beta_{kj}} coefficients on species traits
#' @param TrData a data frame of measured species traits for
#'    \code{\link{formula}}-based specification
#' @param Tr a matrix of measured traits for direct specification
#' @param TrScale a boolean flag whether to scale values of species traits
#' @param phyloTree a phylogenetic tree (object of class \code{phylo} or \code{corPhyl}) for species in \code{Y}
#' @param C a phylogenic correlation matrix for species in \code{Y}
#' @param distr a string shortcut or \eqn{n_s \times 2} matrix specifying the observation models
#' @param truncateNumberOfFactors logical, reduces the maximal number of latent factor to be at most the number of species
#'
#' @return An object of \code{Hmsc} class without any posterior samples.
#'
#' @details Matrix \eqn{Y} may contain missing values, but it is not recommended to add a
#'   species/sampling unit with fully missing data, since those do not bring any new additional information.
#'
#'   Only one of \code{XFormula}-\code{XData} and \code{X} arguments can be specified. Similar requirement applies to
#'   \code{TrFormula}-\code{TrData} and \code{Tr}. It is recommended to use the specification with \code{\link{formula}},
#'   since that information enables additional features for postprocessing of the fitted model.
#'
#'   As default, scaling is applied for \code{X} and \code{Tr} matrices, but not for \code{Y} matrix. If the \code{X} and/or \code{Tr} matrices are
#'   scaled, the estimated parameters are back-transformed so that the estimated parameters correspond to the original
#'   \code{X} and \code{Tr} matrices, not the scaled ones. In contrast, if \code{Y} is scaled, the estimated parameters are not
#'   back-transformed because doing so is not possible for all model parameters. Thus, the estimated parameters
#'   correspond to the scaled \code{Y} matrix, not the original one. If the \code{Y} matrix is scaled, the predictions
#'   generated by \code{predict} are back-transformed, so that the predicted \code{Y} matrices are directly comparable
#'   to the original \code{Y} matrix. If default priors are assumed, it is recommended that all matrices (\code{X}, \code{Tr} and \code{Y}) are scaled.
#'
#'   The object \code{XSelect} is a list. Each object of the list \code{Xsel = XSelect[[i]]} is a named list with objects
#'   \code{Xsel$covGroup}, \code{Xsel$spGroup} and \code{Xsel$q}. The parameter \code{covGroup} is a vector containing
#'   the columns of the matrix \code{X} for which variable selection is applied. The parameter \code{spGroup}
#'   is a vector of length equal to the number of species \code{ns}, with values \code{1,...,ng},
#'   where \code{ng} is the number of
#'   groups of species for which variable selection is applied simultanously. The parameter
#'   \code{q} is a vector of length \code{ng}, containing the prior probabilities by which the variables
#'   are to be included.
#'   For example, choosing \code{covGroup = c(2,3)}, \code{spGroup = rep(1,ns)} and \code{q=0.1} either includes
#'   or excludes both of the covariates 2 and 3 simultaneously for all species. For another example, choosing \code{covGroup = c(2,3)},
#'   \code{spGroup = 1:ns} and \code{q=rep(0.1,ns)} either includes
#'   or excludes both of the covariates 2 and 3 separately for each species.
#'
#'   The included random levels are specified by the \code{ranLevels} and \code{ranLevelsUsed} arguments. The
#'   correspondence between units of each random level and rows of \code{Y} must be specified by a column of
#'   \code{studyDesign}, which corresponds to the name of a list item in \code{ranLevels}. It is
#'   possible to provide an arbitrary number of columns in \code{studyDesign} that are not listed in \code{ranLevels}.
#'   These do not affect the model formulation or fitting scheme, but can be utilized during certain functions
#'   postprocessing the results of statistical model fit.
#'
#'   The \code{distr} argument may be either a matrix, a string literal, or a vector of string literals. In the case of
#'   a matrix, the dimension must be \eqn{n_s \times 2}, where the first column defines the family of the observation
#'   model and the second argument defines the dispersion property. The elements of the first column must take values
#'   1-normal, 2-probit and 3-Poisson with log link function. The second argument stands for the dispersion parameter
#'   being fixed (0) or estimated (1). The default fixed values of the dispersion parameters are 1 for normal and probit,
#'   and 0.01 for Poisson (implemented as a limiting case of lognormally-overdispersed Poisson). Alternatively, a string
#'   literal shortcut can be given as a value to the \code{distr} argument, simultaniously specifying similar class of
#'   observation models for all species. The available shortcuts are \code{"normal"}, \code{"probit"}, \code{"poisson"},
#'   \code{"lognormal poisson"}. If \code{distr} is a vector of string literals, each element corresponds to one species,
#'   should be either \code{"normal"}, \code{"probit"}, \code{"poisson"}, \code{"lognormal poisson"},
#' and these can be abbreviated as long as they are unique strings.
#'   The matrix argument and the vector of string literals allows specifying different observation
#'   models for different species.
#'
#'   By default this constructor assigns default priors to the latent factors. Those priors are designed to be
#'   reasonably flat assuming that the covariates, species traits and normally distributed responses are scaled.
#'   In case when other priors needed to be specified, a call of \code{setPriors.Hmsc} methods should be made,
#'   where the particular priors may be specified.
#'
#' @seealso \code{\link{HmscRandomLevel}}, \code{\link{sampleMcmc}}, \code{\link{setPriors.Hmsc}}
#'
#' @examples
#' # Creating a Hmsc object without phylogeny, trait data or random levels
#' m = Hmsc(Y=TD$Y, XData=TD$X, XFormula=~x1+x2)
#'
#' # Creating a Hmsc object with phylogeny and traits
#' m = Hmsc(Y=TD$Y, XData=TD$X, XFormula=~x1+x2,
#' TrData=TD$Tr, TrFormula=~T1+T2, phyloTree=TD$phylo)
#'
#' # Creating a Hmsc object with 2 nested random levels (50 sampling units in 20 plots)
#' studyDesign = data.frame(sample = as.factor(1:50), plot = as.factor(sample(1:20,50,replace=TRUE)))
#' rL1 = HmscRandomLevel(units=TD$studyDesign$plot)
#' rL2 = HmscRandomLevel(units=TD$studyDesign$sample)
#' m = Hmsc(Y=TD$Y, XData=TD$X, XFormula=~x1+x2,
#' studyDesign=studyDesign,ranLevels=list("sample"=rL1,"plot"=rL2))
#'
#' @importFrom stats model.matrix
#' @importFrom ape vcv.phylo
#'
#' @export


Hmsc = function(Y, XFormula=~., XData=NULL, X=NULL, XScale=TRUE,
                XSelect=NULL,
                XRRRData=NULL, XRRRFormula=~.-1, XRRR=NULL, ncRRR=2, XRRRScale=TRUE,
                YScale = FALSE,
                studyDesign=NULL, ranLevels=NULL, ranLevelsUsed=names(ranLevels),
                TrFormula=NULL, TrData=NULL, Tr=NULL, TrScale=TRUE,
                phyloTree=NULL, C=NULL,
                distr="normal", truncateNumberOfFactors=TRUE){

   hM = structure(list(
      Y = NULL,
      XData=NULL, XFormula=NULL, X=NULL, XScaled=NULL, XSelect=NULL,
      XRRRData=NULL, XRRRFormula=NULL, XRRR=NULL, XRRRScaled=NULL,
      YScaled=NULL, XInterceptInd=NULL,
      studyDesign=NULL, ranLevels=NULL, ranLevelsUsed=NULL,
      dfPi = NULL, rL=NULL, Pi=NULL,
      TrData=NULL, TrFormula=NULL, Tr=NULL, TrScaled=NULL, TrInterceptInd=NULL,
      C=NULL, phyloTree=NULL,
      distr = NULL,

      # dimensions
      ny = NULL,
      ns = NULL,
      nc = NULL,
      ncNRRR = NULL,
      ncRRR = NULL,
      ncORRR = NULL,
      ncsel = NULL,
      nr = NULL,
      nt = NULL,
      nf = NULL,
      ncr = NULL,
      ncs = NULL,
      np = NULL,

      # names
      spNames = NULL,
      covNames = NULL,
      trNames = NULL,
      rLNames = NULL,

      # scaling
      XScalePar=NULL, XRRRScalePar=NULL, YScalePar=NULL, TrScalePar=NULL,

      # priors
      V0=NULL, f0=NULL,
      mGamma=NULL, UGamma=NULL,
      aSigma=NULL, bSigma=NULL,
      nu=NULL, a1=NULL, b1=NULL, a2=NULL, b2=NULL,
      rhopw=NULL,
      nuRRR=NULL, a1RRR=NULL, b1RRR=NULL, a2RRR=NULL, b2RRR=NULL,

      # sampling parameters
      samples=NULL, transient=NULL, thin=NULL, verbose=NULL, adaptNf=NULL,
      initPar=NULL, repN=NULL,
      randSeed=NULL,

      # posterior
      postList=NULL), class="Hmsc")

   ## take care that Y is a matrix (for data.frame, tibble, numeric vector)
   Y <- as.matrix(Y)
   if(!(is.numeric(Y) || is.logical(Y))) {  # allow binary data as TRUE/FALSE
      stop("Hmsc.setData: Y must be a numeric matrix of sampling units times species")
   }
   hM$Y = Y
   hM$ny = nrow(Y)
   hM$ns = ncol(Y)
   if(is.null(colnames(hM$Y))){
      colnames(hM$Y) = sprintf(sprintf("sp%%.%dd",ceiling(log10(hM$ns))), 1:hM$ns)
   }
   hM$spNames = colnames(hM$Y)

   # linear regression covariates
   if(!is.null(XData) && !is.null(X)){
      stop("Hmsc.setData: only single of XData and X arguments must be specified")
   }
    if(!is.null(XData)) {
        if (inherits(XData, "list")) {
            if(length(XData) != hM$ns){
                stop("Hmsc.setData: the length of XData list argument must be equal to the number of species")
            }
            if(any(!unlist(lapply(XData, is.data.frame)))){
                stop("Hmsc.setData: each element of X list must be a data.frame")
            }
            ## check agasint derived classes of data.frame (such as tibble)
            if (any(sapply(XData, function(a) class(a)[1L] != "data.frame"))) {
                for(i in seq_len(length(XData)))
                    XData[[i]] <- as.data.frame(XData, stringsAsFactors = TRUE)
            }
            if(any(unlist(lapply(XData, function(a) nrow(a) != hM$ny)))){
                stop("Hmsc.setData: for each element of XData list the number of rows must be equal to the number of sampling units")
            }
            if(any(unlist(lapply(XData, function(a) any(is.na(a)))))){
                stop("Hmsc.setData: all elements of XData list must contain no NA values")
            }
            hM$XData = XData
            hM$XFormula = XFormula
            hM$X = lapply(XData, function(a) model.matrix(XFormula, a))
            hM$nc = ncol(hM$X[[1]])
        }
        else if (is.data.frame(XData)) {
            if(nrow(XData) != hM$ny){
                stop("Hmsc.setData: the number of rows in XData must be equal to the number of sampling units")
            }
            if(any(is.na(XData))){
                stop("Hmsc.setData: XData must contain no NA values")
            }
            ## check against derived classes of data.frame (such as tibble)
            if (class(XData)[1L] != "data.frame")
                XData <- as.data.frame(XData, stringsAsFactors = TRUE)
            hM$XData = XData
            hM$XFormula = XFormula
            hM$X = model.matrix(XFormula, XData)
            hM$nc = ncol(hM$X)
        } else {
            stop("Hmsc.setData: XData must either a data.frame or a list of data.frame objects")
        }
    }
   if(!is.null(X)){
      switch(class(X)[1L],
             list={
                if(length(X) != hM$ns){
                   stop("Hmsc.setData: the length of X list argument must be equal to the number of species")
                }
                if(any(!unlist(lapply(X, is.matrix)))){
                   stop("Hmsc.setData: each element of X list must be a matrix")
                }
                if(any(unlist(lapply(X, function(a) nrow(a) != hM$ny)))){
                   stop("Hmsc.setData: for each element of X list the number of rows must be equal to the number of sampling units")
                }
                if(any(unlist(lapply(X, function(a) any(is.na(a)))))){
                   stop("Hmsc.setData: all elements of X list must contain no NA values")
                }
                hM$XData = NULL
                hM$XFormula = NULL
                hM$X = X
                hM$nc = ncol(hM$X[[1]])
             },
             matrix={
                if(nrow(X) != hM$ny){
                   stop("Hmsc.setData: the number of rows in X must be equal to the number of sampling units")
                }
                if(any(is.na(X))){
                   stop("Hmsc.setData: X must contain no NA values")
                }
                hM$XData = NULL
                hM$XFormula = NULL
                hM$X = X
                hM$nc = ncol(hM$X)
             },
             {
                stop("Hmsc.setData: X must either a matrix or a list of matrix objects")
             }
      )
   }
   if(is.null(XData) && is.null(X)){
      X = matrix(NA,hM$ny,0)
      hM$nc = 0
   }

   switch(class(hM$X)[1L],
          matrix = {
             if(is.null(colnames(hM$X))){
                colnames(hM$X) = sprintf(sprintf("cov%%.%dd",ceiling(log10(hM$nc))), 1:hM$nc)
             }
          },
          if(is.null(colnames(hM$X[[1]]))){
             list = {
                for(j in 1:hM$ns)
                   colnames(hM$X[[j]]) = sprintf(sprintf("cov%%.%dd",ceiling(log10(hM$nc))), 1:hM$nc)
             }
          }
   )

   switch (class(hM$X)[1L],
           "matrix" = {hM$covNames = colnames(hM$X)},
           "list" = {hM$covNames = colnames(hM$X[[1]])}
   )

   if(identical(XScale,FALSE)){
      hM$XScalePar = rbind(rep(0,hM$nc), rep(1,hM$nc))
      hM$XScaled = hM$X
      hM$XInterceptInd = NULL
   } else{
      switch(class(hM$X)[1L],
             matrix = {
                XStack = hM$X
             },
             list = {
                XStack = Reduce(rbind, hM$X)
             }
      )
      XInterceptInd = which(colnames(XStack) %in% c("Intercept","(Intercept)"))
      if(length(XInterceptInd)>1){
         stop("Hmsc.setData: only one column of X matrix could be named Intercept or (Intercept)")
      }
      if(!all(XStack[,XInterceptInd] == 1)){
         stop("Hmsc.setData: intercept column in X matrix must be a column of ones")
      }
      if(length(XInterceptInd)==1){
         hM$XInterceptInd = XInterceptInd
      } else
         hM$XInterceptInd = NULL
      XScalePar = rbind(rep(0,hM$nc), rep(1,hM$nc))
      XScaled = XStack
      if(identical(XScale,TRUE)){
         scaleInd = apply(XStack, 2, function(a) !all(a %in% c(0,1)))
      } else{
         scaleInd = XScale
      }
      scaleInd[XInterceptInd] = FALSE
      if(length(XInterceptInd)>0){
         sc = scale(XStack)
         XScalePar[,scaleInd] = rbind(attr(sc,"scaled:center"), attr(sc,"scaled:scale"))[,scaleInd]
      } else{
         sc = scale(XStack, center=FALSE)
         XScalePar[2,scaleInd] = attr(sc,"scaled:scale")[scaleInd]
      }
      XScaled[,scaleInd] = sc[,scaleInd]
      hM$XScalePar = XScalePar
      switch(class(hM$X)[1L],
             matrix = {
                hM$XScaled = XScaled
             },
             list = {
                hM$XScaled = lapply(split(XScaled,rep(1:hM$ns,each=hM$ny)), function(a) matrix(a,hM$ny,hM$nc))
             }
      )
   }


   hM$ncsel = length(XSelect)
   hM$XSelect = XSelect
   for (i in seq_len(hM$ncsel)){
      XSel = hM$XSelect[[i]]
      if(max(XSel$covGroup)>hM$nc){
         stop("Hmsc.setData: covGroup for XSelect cannot have values greater than number of columns in X")
      }
   }


   #covariates for reduced-rank regression
   hM$ncNRRR = hM$nc
   if(!is.null(XRRRData)){
      if(!is.data.frame(XRRRData))
      {
         stop("Hmsc.setData: XRRRData must be a data.frame")
      }
      if(nrow(XRRRData) != hM$ny){
         stop("Hmsc.setData: the number of rows in XRRRData must be equal to the number of sampling units")
      }
      if(any(is.na(XRRRData))){
         stop("Hmsc.setData: XRRRData must contain no NA values")
      }
      hM$XRRRData = XRRRData
      hM$XRRRFormula = XRRRFormula
      hM$XRRR = model.matrix(XRRRFormula, XRRRData)
      hM$ncORRR = ncol(hM$XRRR)
      hM$ncRRR = ncRRR
   }
   if(!is.null(XRRR)){
      if(!is.matrix(XRRR))
      {
         stop("Hmsc.setData: XRRR must be a matrix")
      }
      if(nrow(XRRR) != hM$ny){
         stop("Hmsc.setData: the number of rows in XRRR must be equal to the number of sampling units")
      }
      if(any(is.na(XRRR))){
         stop("Hmsc.setData: XRRR must contain no NA values")
      }
      hM$XRRRData = NULL
      hM$XRRRFormula = NULL
      hM$XRRR = XRRR
      hM$ncORRR = ncol(hM$XRRR)
      hM$ncRRR = ncRRR
   }
   if(is.null(XRRRData) && is.null(XRRR)){
      X = matrix(NA,hM$ny,0)
      hM$XRRR = NULL
      hM$ncORRR = 0
      hM$ncRRR = 0
   }
   if(hM$ncRRR>0){
      if(is.null(colnames(hM$XRRR))){
         colnames(hM$XRRR) = sprintf(sprintf("covRRR%%.%dd",ceiling(log10(hM$ncORRR))), 1:hM$ncORRR)
      }
      for(k in seq_len(hM$ncRRR)){
         hM$covNames=c(hM$covNames,paste0("XRRR_",as.character(k)))
      }
      hM$nc = hM$ncNRRR + hM$ncRRR

      if(identical(XRRRScale,FALSE)){
         hM$XRRRScalePar = rbind(rep(0,hM$ncORRR), rep(1,hM$ncORRR))
         hM$XRRRScaled = hM$XRRR
      } else {
         if(identical(XScale,FALSE)){
            stop("Hmsc.setData: XRRR can't be scaled if X is not scaled")
         }
         XRRRStack = hM$XRRR
         if(identical(XRRRScale,TRUE)){
            XRRRscaleInd = apply(XRRRStack, 2, function(a) !all(a %in% c(0,1)))
         } else{
            XRRRscaleInd = XRRRScale
         }
         XRRRScalePar = rbind(rep(0,hM$ncORRR), rep(1,hM$ncORRR))
         XRRRScaled=XRRRStack
         if(length(XInterceptInd)>0){
            XRRRsc = scale(XRRRStack)
            XRRRScalePar[,XRRRscaleInd] = rbind(attr(XRRRsc,"scaled:center"), attr(XRRRsc,"scaled:scale"))[,XRRRscaleInd]
         } else{
            XRRRsc = scale(XRRRStack, center=FALSE)
            XRRRScalePar[2,XRRRscaleInd] = attr(sc,"scaled:scale")[XRRRscaleInd]
         }
         XRRRScaled[,XRRRscaleInd] = XRRRsc[,XRRRscaleInd]
         hM$XRRRScalePar = XRRRScalePar
         hM$XRRRScaled = XRRRScaled
      }
   }

   # traits
   if(!is.null(TrData)){
      if(!is.null(Tr)){
         stop("Hmsc.setData: at maximum one of TrData and Tr arguments can be specified")
      }
      if(is.null(TrFormula)){
         stop("Hmsc.setData: TrFormula argument must be specified if TrData is provided")
      }
   }
   if(!is.null(TrData)){
      if(nrow(TrData) != hM$ns){
         stop("Hmsc.setData: the number of rows in TrData should be equal to number of columns in Y")
      }
      if(any(is.na(TrData))){
         stop("Hmsc.setData: TrData parameter must not contain any NA values")
      }
      hM$TrData = TrData
      hM$TrFormula = TrFormula
      hM$Tr = model.matrix(TrFormula, TrData)
   }
   if(!is.null(Tr)){
      if(!is.matrix(Tr)){
         stop("Hmsc.setData: Tr must be a matrix")
      }
      if(nrow(Tr) != hM$ns){
         stop("Hmsc.setData: the number of rows in Tr should be equal to number of columns in Y")
      }
      if(any(is.na(Tr))){
         stop("Hmsc.setData: Tr parameter must not contain any NA values")
      }
      hM$TrData = NULL
      hM$Tr = Tr
   }
   if(is.null(hM$Tr)){
      hM$Tr = matrix(1,hM$ns,1)
   }
   hM$nt = ncol(hM$Tr)
   if(is.null(colnames(hM$Tr))){
      colnames(hM$Tr) = sprintf(sprintf("tr%%.%dd",ceiling(log10(hM$nt))), 1:hM$nt)
   }
   hM$trNames = colnames(hM$Tr)

   if(identical(TrScale,FALSE)){
      hM$TrScalePar = rbind(rep(0,hM$nt), rep(1,hM$nt))
      hM$TrScaled = hM$Tr
      hM$TrInterceptInd = NULL
   } else{
      TrInterceptInd = which(colnames(hM$Tr) %in% c("Intercept","(Intercept)"))
      if(length(TrInterceptInd)>1){
         stop("Hmsc.setData: only one column of Tr matrix could be named Intercept or (Intercept)")
      }
      if(!all(hM$Tr[,TrInterceptInd]==1)){
         stop("Hmsc.setData: intercept column in Tr matrix must be a column of ones")
      }
      if(length(TrInterceptInd)==1){
         hM$TrInterceptInd = TrInterceptInd
      } else
         hM$TrInterceptInd = NULL
      TrScalePar = rbind(rep(0,hM$nt), rep(1,hM$nt))
      TrScaled = hM$Tr
      if(identical(TrScale,TRUE)){
         scaleInd = apply(hM$Tr, 2, function(a) !all(a %in% c(0,1)))
      } else{
         scaleInd = TrScale
      }
      scaleInd[TrInterceptInd] = FALSE
      if(length(TrInterceptInd)>0){
         sc = scale(hM$Tr)
         TrScalePar[,scaleInd] = rbind(attr(sc,"scaled:center"), attr(sc,"scaled:scale"))[,scaleInd]
      } else{
         sc = scale(hM$Tr, center=FALSE)
         TrScalePar[2,scaleInd] = attr(sc,"scaled:scale")[scaleInd]
      }
      TrScaled[,scaleInd] = sc[,scaleInd]
      hM$TrScalePar = TrScalePar
      hM$TrScaled = TrScaled
   }

   # phylogeny
   if(!is.null(C) && !is.null(phyloTree)){
      stop("Hmsc.setData: at maximum one of phyloTree and C arguments can be specified")
   }
   if(!is.null(phyloTree)){
      corM = vcv.phylo(phyloTree, model="Brownian", corr=TRUE)
      corM = corM[hM$spNames,hM$spNames]
      hM$phyloTree = phyloTree
      hM$C = corM
   }
   if(!is.null(C)){
      if(any(dim(C) != hM$ns)){
         stop("Hmsc.setData: the size of square matrix C must be equal to number of species")
      }
      hM$C = C
   }

   # latent factors
   if(is.null(studyDesign)){
      hM$dfPi = NULL
      hM$Pi = matrix(NA,hM$ny,0)
      hM$np = integer(0)
      hM$nr = 0
      hM$rLNames = character(0)
      if(!is.null(ranLevels)){
         if(length(ranLevels) > 0){
            stop("Hmsc.setData: studyDesign is empty, but ranLevels is not")
         }
      }
   } else {
      if(nrow(studyDesign) != hM$ny){
         stop("Hmsc.setData: the number of rows in studyDesign must be equal to number of rows in Y")
      }
      if(!all(ranLevelsUsed %in% names(ranLevels))){
         stop("Hmsc.setData: ranLevels must contain named elements corresponding to all levels listed in ranLevelsUsed")
      }
      if(!all(ranLevelsUsed %in% colnames(studyDesign))){
         stop("Hmsc.setData: studyDesign must contain named columns corresponding to all levels listed in ranLevelsUsed")
      }
      hM$studyDesign = studyDesign
      hM$ranLevels = ranLevels
      hM$ranLevelsUsed = ranLevelsUsed

      hM$dfPi = studyDesign[,ranLevelsUsed,drop=FALSE]
      hM$rL = ranLevels[ranLevelsUsed]
      hM$rLNames = colnames(hM$dfPi)

      hM$Pi = matrix(NA,hM$ny,ncol(hM$dfPi),dimnames=list(NULL,hM$rLNames))
      for(r in seq_len(ncol(hM$dfPi)))
         hM$Pi[,r] = as.numeric(as.factor(hM$dfPi[,r]))
      hM$np = apply(hM$Pi, 2, function(a) return(length(unique(a))))
      hM$nr = ncol(hM$Pi)
      if (truncateNumberOfFactors){
         for (r in  seq_len(hM$nr)){
            hM$rL[[r]]$nfMax = min(hM$rL[[r]]$nfMax,hM$ns)
            hM$rL[[r]]$nfMin = min(hM$rL[[r]]$nfMin,hM$rL[[r]]$nfMax)
         }
      }
   }

   ## distr can be given as a matrix: check its dims
   if (is.matrix(distr)) {
      if (NROW(distr) != hM$ns)
         stop("No. of rows in distr matrix must be equal to the no. of species")
      if (NCOL(distr) < 2) # we later warn on unused extra columns
         stop("distr matrix should have 2 columns")
   }

   ## allow abbreviation of 'distr' and check that it is one of the
   ## following known ones
   knownDistributions <- c("normal", "probit", "poisson", "lognormal poisson")

   if(length(distr)==1){
      switch (match.arg(distr, knownDistributions),
              "normal" = {
                 distr = matrix(0,hM$ns,2)
                 distr[,1] = 1
                 distr[,2] = 1
              },
              "probit" = {
                 distr = matrix(0,hM$ns,2)
                 distr[,1] = 2
                 distr[,2] = 0
              },
              "poisson" = {
                 distr = matrix(0,hM$ns,2)
                 distr[,1] = 3
                 distr[,2] = 0
              },
              "lognormal poisson" = {
                 distr = matrix(0,hM$ns,2)
                 distr[,1] = 3
                 distr[,2] = 1
              }
      )
   }
   if(length(distr) > 1 && !is.matrix(distr)){
      if (length(distr) != hM$ns)
         stop("length of distr should be 1 or equal to the number of species")
      distr2 = matrix(0,hM$ns,2)
      for (i in 1:hM$ns){
         switch (match.arg(distr[i], knownDistributions),
                 "normal" = {
                    distr2[i,1] = 1
                    distr2[i,2] = 1
                 },
                 "probit" = {
                    distr2[i,1] = 2
                    distr2[i,2] = 0
                 },
                 "poisson" = {
                    distr2[i,1] = 3
                    distr2[i,2] = 0
                 },
                 "lognormal poisson" = {
                    distr2[i,1] = 3
                    distr2[i,2] = 1
                 }
         )
      }
      distr=distr2
   }
   ## we had 4-column distr matrix in ancient versions
   if (NCOL(distr) > 2) {
      warning("Keeping only two first columns of 'distr' matrix")
      distr <- distr[, 1:2, drop=FALSE]
   }
   colnames(distr) = c("family","variance")
   if(any(distr[,1]==0)){
      stop("Hmsc.setData: some of the distributions ill defined")
   }
   hM$distr = distr

   #scaling of response
   if(identical(YScale,FALSE)){
      hM$YScalePar = rbind(rep(0,hM$ns), rep(1,hM$ns))
      hM$YScaled = hM$Y
   } else{
      scaleInd = which(hM$distr[,1]==1)
      YScalePar = rbind(rep(0,hM$ns), rep(1,hM$ns))
      YScaled = hM$Y
      if(length(scaleInd)>0){
         sc = scale(hM$Y)
         YScalePar[,scaleInd] = rbind(attr(sc,"scaled:center"), attr(sc,"scaled:scale"))[,scaleInd]
         YScaled[,scaleInd] = sc[,scaleInd]
      }
      hM$YScalePar = YScalePar
      hM$YScaled = YScaled
   }

   hM = setPriors(hM, setDefault=TRUE)
   hM$call <- match.call()
   hM
}