Offset imlemented to Hmsc models.

Fix for exproting GPP defined with Spatial object to Hmsc-HPC. Minor adjustment for LFix and LRan summation or subtraction. Minor adjustments to make the code for error and warning messages more dense.
hmsc-r · Mar 22, 2024 · f26855d · f26855d
1 parent f21adb8
commit f26855d
Show file tree

Hide file tree

Showing 17 changed files with 92 additions and 95 deletions.
diff --git a/R/Hmsc.R b/R/Hmsc.R
@@ -115,14 +115,14 @@
 Hmsc = function(Y, XFormula=~., XData=NULL, X=NULL, XScale=TRUE,
                 XSelect=NULL,
                 XRRRData=NULL, XRRRFormula=~.-1, XRRR=NULL, ncRRR=2, XRRRScale=TRUE,
-                YScale = FALSE,
+                YScale = FALSE, Loff=NULL,
                 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,
+      Y = NULL, Loff=NULL,
       XData=NULL, XFormula=NULL, X=NULL, XScaled=NULL, XSelect=NULL,
       XRRRData=NULL, XRRRFormula=NULL, XRRR=NULL, XRRRScaled=NULL,
       YScaled=NULL, XInterceptInd=NULL,
@@ -701,6 +701,16 @@ Hmsc = function(Y, XFormula=~., XData=NULL, X=NULL, XScale=TRUE,
       hM$YScaled = YScaled
    }
 
+   #offset to the latent predictor
+   if(!is.null(Loff)){
+      if(!is.matrix(Loff)) stop("Loff argument must be NULL or a numeric matrix")
+      if(nrow(Loff) != hM$ny) stop("number of rows in Loff argument must be equal to ny")
+      if(ncol(Loff) != hM$ns) stop("number of columns in Loff argument must be equal to ns")
+      hM$Loff = Loff
+   } else{
+      hM$Loff = NULL
+   }
+
    hM = setPriors(hM, setDefault=TRUE)
    hM$call <- match.call()
    hM$HmscVersion <- packageVersion("Hmsc")

diff --git a/R/biPlot.R b/R/biPlot.R
@@ -22,15 +22,15 @@
 #' # Construct an ordination biplot using two chosen latent factors from a previously fitted HMSC model
 #' etaPost = getPostEstimate(TD$m, "Eta")
 #' lambdaPost=getPostEstimate(TD$m, "Lambda")
-#' biPlot(TD$m, etaPost = etaPost, lambdaPost=lambdaPost, factors=c(1,2))
+#' biPlot(TD$m, etaPost=etaPost, lambdaPost=lambdaPost, factors=c(1,2))
 #'
 #'
 #' @importFrom graphics plot points text
 #' @importFrom grDevices colorRampPalette palette
 #'
 #' @export
 
-biPlot=function(hM, etaPost, lambdaPost, factors=c(1,2), colVar=NULL, colors = NULL, spNames=hM$spNames, ...){
+biPlot=function(hM, etaPost, lambdaPost, factors=c(1,2), colVar=NULL, colors=NULL, spNames=hM$spNames, ...){
    if(!is.null(colVar)){
       if (!is.null(hM$XData))
           col = hM$XData[,colVar]

diff --git a/R/c.Hmsc.R b/R/c.Hmsc.R
@@ -46,7 +46,7 @@
     ## phyloTree (but only C); call. What about *Names, data.frames
     ## (built to model.matrix)?
     checkItems <-
-        c("Y", "XData", "X", "XScaled", "XRRRData", "XRRRScaled",
+        c("Y", "Loff", "XData", "X", "XScaled", "XRRRData", "XRRRScaled",
           "YScaled", "XInterceptInd", "studyDesign", "ranLevels",
           "ranLevelsUsed", "dfPi", "rL", "Pi", "TrData","Tr",
           "TrScaled", "TrInterceptInd", "C", "distr", "ny", "ns",

diff --git a/R/computeInitialParameters.R b/R/computeInitialParameters.R
@@ -255,7 +255,7 @@ computeInitialParameters = function(hM, initPar, computeZ=TRUE){
       Z = LFix
 
    if(computeZ){
-      Z = updateZ(Y=hM$Y,Z=Z,Beta=Beta,iSigma=sigma^-1,Eta=Eta,Lambda=Lambda, X=XScaled,Pi=hM$Pi,dfPi=hM$dfPi,distr=hM$distr,rL=hM$rL)
+      Z = updateZ(Y=hM$Y,Z=Z,Beta=Beta,iSigma=sigma^-1,Eta=Eta,Lambda=Lambda, Loff=hM$Loff,X=XScaled,Pi=hM$Pi,dfPi=hM$dfPi,distr=hM$distr,rL=hM$rL)
    } else{
       Z = NULL
    }

diff --git a/R/computePredictedValues.R b/R/computePredictedValues.R
@@ -97,6 +97,7 @@ computePredictedValues = function(hM, partition=NULL, partition.sp=NULL, start=1
          for(r in seq_len(hM$nr)){
             dfPi[,r] = factor(hM$dfPi[train,r])
          }
+         LoffVal = hM$LoffVal[val,,drop=FALSE]
          switch(class(hM$X)[1L],
                 matrix = {
                    XTrain = hM$X[train,,drop=FALSE]
@@ -120,7 +121,7 @@ computePredictedValues = function(hM, partition=NULL, partition.sp=NULL, start=1
          ## below. If this disturbs, it should probably be fixed in
          ## Hmsc(), but currently ranLevels=character(0) seems to
          ## work.
-         hM1 = Hmsc(Y=hM$Y[train,,drop=FALSE], X=XTrain,
+         hM1 = Hmsc(Y=hM$Y[train,,drop=FALSE], Loff=hM$Loff[train,,drop=FALSE], X=XTrain,
                     XRRR=XRRRTrain, ncRRR = hM$ncRRR, XSelect = hM$XSelect,
                     distr=hM$distr, studyDesign=dfPi, Tr=hM$Tr, C=hM$C, ranLevels=hM$rL)
          setPriors(hM1, V0=hM$V0, f0=hM$f0, mGamma=hM$mGamma, UGamma=hM$UGamma, aSigma=hM$aSigma, bSigma=hM$bSigma, rhopw=hM$rhowp)
@@ -157,7 +158,7 @@ computePredictedValues = function(hM, partition=NULL, partition.sp=NULL, start=1
             dfPi[,r] = factor(hM$dfPi[val,r])
          }
          if(is.null(partition.sp)){
-            pred1 = predict(hM1, post=postList, X=XVal, XRRR=XRRRVal, studyDesign=dfPi, Yc=Yc[val,,drop=FALSE], mcmcStep=mcmcStep, expected=expected)
+            pred1 = predict(hM1, post=postList, Loff=LoffVal, X=XVal, XRRR=XRRRVal, studyDesign=dfPi, Yc=Yc[val,,drop=FALSE], mcmcStep=mcmcStep, expected=expected)
             pred1Array = abind(pred1,along=3)
          } else {
             hM2 = duplicate(hM)
@@ -182,7 +183,7 @@ computePredictedValues = function(hM, partition=NULL, partition.sp=NULL, start=1
                val.sp = (partition.sp==i)
                YcFull = hM$Y
                YcFull[val,val.sp] = NA
-               pred2 = predict(hM2, post=postList2, X=hM2$X, XRRR=hM2$XRRR,studyDesign=hM2$studyDesign, Yc=YcFull, mcmcStep=mcmcStep, expected=expected)
+               pred2 = predict(hM2, post=postList2, Loff=hM2$Loff, X=hM2$X, XRRR=hM2$XRRR, studyDesign=hM2$studyDesign, Yc=YcFull, mcmcStep=mcmcStep, expected=expected)
                pred2Array = abind(pred2,along=3)
                pred1Array[,val.sp,] = pred2Array[val,val.sp,]
             }

diff --git a/R/computePredictor.HmscRandomLevel.R b/R/computePredictor.HmscRandomLevel.R
@@ -5,6 +5,11 @@ computePredictor.HmscRandomLevel = function(rL,Eta,Lambda,pi=c(1:nrow(Eta)),dfpi
       LRan = matrix(0,length(pi),ncol(Lambda))
       for(k in 1:rL$xDim)
          LRan = LRan + (Eta[pi,]*rL$x[as.character(dfpi),k]) %*% Lambda[,,k]
+
+      # np = length(pi)
+      # nf = ncol(Lambda)
+      # EtaArray = array(Eta[pi,], c(np,nf,rL$xDim))
+      # XArray = array(rL$x[rep(as.character(dfpi),nf),], c(np,nf,rL$xDim))
    }
    return(LRan)
 }
diff --git a/R/predict.R b/R/predict.R
@@ -70,8 +70,8 @@
 #'
 #' @export
 
-predict.Hmsc = function(object, post=poolMcmcChains(object$postList), XData=NULL,
-                        X=NULL, XRRRData=NULL, XRRR=NULL, # this has to be updated to cov-dependent associations
+predict.Hmsc = function(object, post=poolMcmcChains(object$postList), Loff=NULL,
+                        XData=NULL, X=NULL, XRRRData=NULL, XRRR=NULL, # this has to be updated to cov-dependent associations
                         studyDesign=object$studyDesign, ranLevels=object$ranLevels,
                         Gradient=NULL, Yc=NULL, mcmcStep=1, expected=FALSE,
                         predictEtaMean=FALSE, predictEtaMeanField=FALSE,
@@ -89,8 +89,7 @@ predict.Hmsc = function(object, post=poolMcmcChains(object$postList), XData=NULL
    if(!is.null(Gradient)) {
       ## don't know what to do if there is also Yc, and spatial models
       ## will trigger an error in updateEta (github issue #135)
-      if (!is.null(Yc))
-          stop("predict with arguments 'Yc' and 'Gradient' jointly is not implemented (yet)")
+      if(!is.null(Yc)) stop("predict with arguments 'Yc' and 'Gradient' jointly is not implemented (yet)")
       XData=Gradient$XDataNew
       studyDesign=Gradient$studyDesignNew
       ranLevels=Gradient$rLNew
@@ -107,14 +106,12 @@ predict.Hmsc = function(object, post=poolMcmcChains(object$postList), XData=NULL
    if(!is.null(XData)){
       switch(class(XData)[1L],
              list={
-                if (any(unlist(lapply(XData, is.na))))
-                    stop("NA values are not allowed in 'XData'")
+                if(any(unlist(lapply(XData, is.na)))) stop("NA values are not allowed in 'XData'")
                 xlev = lapply(Reduce(rbind,object$XData), levels)[unlist(lapply(Reduce(rbind,object$XData), is.factor))]
                 X = lapply(XData, function(a) model.matrix(object$XFormula, a, xlev=xlev))
              },
              data.frame={
-                if (any(is.na(XData)))
-                    stop("NA values are not allowed in 'XData'")
+                if(any(is.na(XData))) stop("NA values are not allowed in 'XData'")
                 xlev = lapply(object$XData, levels)[unlist(lapply(object$XData, is.factor))]
                 X = model.matrix(object$XFormula, XData, xlev=xlev)
              }
@@ -141,12 +138,12 @@ predict.Hmsc = function(object, post=poolMcmcChains(object$postList), XData=NULL
    )
 
    if(!is.null(Yc)){
-      if(ncol(Yc) != object$ns){
-         stop("number of columns in Yc must be equal to ns")
-      }
-      if(nrow(Yc) != nyNew){
-         stop("number of rows in Yc and X must be equal")
-      }
+      if(ncol(Yc) != object$ns) stop("number of columns in Yc must be equal to ns")
+      if(nrow(Yc) != nyNew) stop("number of rows in Yc and X must be equal")
+   }
+   if(!is.null(Loff)){
+      if(ncol(Loff) != object$ns) stop("number of columns in Loff must be equal to ns")
+      if(nrow(Loff) != nyNew) stop("number of rows in Loff and X must be equal")
    }
    if(!all(object$rLNames %in% colnames(studyDesign))){
       stop("dfPiNew does not contain all the necessary named columns")
@@ -192,7 +189,7 @@ predict.Hmsc = function(object, post=poolMcmcChains(object$postList), XData=NULL
    if (nParallel == 1) {  # non-Parallel
        pred <- lapply(seq_len(predN), function(pN, ...){
          # print(ppEta[pN,])
-         get1prediction(object, X, XRRR, Yc, rL, rLPar, post[[pN]],
+         get1prediction(object, X, XRRR, Yc, Loff, rL, rLPar, post[[pN]],
                         ppEta[pN,], PiNew, dfPiNew, nyNew, expected,
                         mcmcStep)})
 
@@ -203,14 +200,14 @@ predict.Hmsc = function(object, post=poolMcmcChains(object$postList), XData=NULL
        clusterEvalQ(cl, {
            library(Hmsc)})
        pred <- parLapply(cl, seq_len(predN), function(pN, ...)
-           get1prediction(object, X, XRRR, Yc, rL, rLPar, post[[pN]],
+           get1prediction(object, X, XRRR, Yc, Loff, rL, rLPar, post[[pN]],
                           ppEta[pN,], PiNew, dfPiNew, nyNew, expected,
                           mcmcStep, seed = seed[pN]))
        stopCluster(cl)
    } else { # fork (mac, Linux)
        seed <- sample.int(.Machine$integer.max, predN)
        pred <- mclapply(seq_len(predN), function(pN, ...)
-           get1prediction(object, X, XRRR, Yc, rL, rLPar, post[[pN]],
+           get1prediction(object, X, XRRR, Yc, Loff, rL, rLPar, post[[pN]],
                           ppEta[pN,], PiNew, dfPiNew, nyNew, expected,
                           mcmcStep, seed = seed[pN]),
            mc.cores=nParallel)
@@ -225,7 +222,7 @@ predict.Hmsc = function(object, post=poolMcmcChains(object$postList), XData=NULL
 ##  predPostEta rL rLPar
 
 get1prediction <-
-    function(object, X, XRRR, Yc, rL, rLPar, sam, predPostEta, PiNew, dfPiNew,
+    function(object, X, XRRR, Yc, Loff, rL, rLPar, sam, predPostEta, PiNew, dfPiNew,
              nyNew, expected, mcmcStep, seed = NULL)
 {
     if (!is.null(seed))
@@ -265,22 +262,23 @@ get1prediction <-
                 LRan[[r]] = LRan[[r]] + (Eta[[r]][as.character(dfPiNew[,r]),]*rL[[r]]$x[as.character(dfPiNew[,r]),k]) %*% sam$Lambda[[r]][,,k]
         }
     }
-    if(object$nr > 0){L = LFix + Reduce("+", LRan)} else L = LFix
+    L = Reduce("+", c(LFix, LRan))
+    if(!is.null(Loff)) L = L + Loff
 
     ## predict can be slow with Yc and especially with high mcmcStep
     if(!is.null(Yc) && any(!is.na(Yc))){
         Z = L
         Z = updateZ(Y=Yc, Z=Z, Beta=sam$Beta, iSigma=1/sam$sigma, Eta=Eta,
-                    Lambda=sam$Lambda, X=X, Pi=PiNew, dfPi=dfPiNew,
+                    Lambda=sam$Lambda, Loff=Loff, X=X, Pi=PiNew, dfPi=dfPiNew,
                     distr=object$distr, rL=rL)
         ## species CV from computePredictedValues runs this innermost
         ## loop nfolds * nfolds.sp * predN * mcmcStep times
         for(sN in seq_len(mcmcStep)){
             Eta = updateEta(Y=Yc, Z=Z, Beta=sam$Beta, iSigma=1/sam$sigma,
                             Eta=Eta, Lambda=sam$Lambda, Alpha=sam$Alpha,
-                            rLPar=rLPar, X=X, Pi=PiNew, dfPi=dfPiNew, rL=rL)
+                            rLPar=rLPar, Loff=Loff, X=X, Pi=PiNew, dfPi=dfPiNew, rL=rL)
             Z = updateZ(Y=Yc, Z=Z, Beta=sam$Beta, iSigma=1/sam$sigma, Eta=Eta,
-                        Lambda=sam$Lambda, X=X, Pi=PiNew, dfPi=dfPiNew,
+                        Lambda=sam$Lambda, Loff=Loff, X=X, Pi=PiNew, dfPi=dfPiNew,
                         distr=object$distr, rL=rL)
         }
         for(r in seq_len(object$nr)){

diff --git a/R/sampleMcmc.R b/R/sampleMcmc.R
@@ -255,6 +255,7 @@
       obj$hM$ranLevels = NULL
       for(r in seq_len(hM$nr)){
          obj$hM$rL[[r]]$s = NULL
+         obj$hM$rL[[r]]$sKnot = NULL
       }
     }
     obj
@@ -344,6 +345,7 @@
     adaptNf = obj$adaptNf
     Tr = hM$TrScaled
     Y = hM$YScaled
+    Loff = hM$Loff
     distr = hM$distr
     Pi = hM$Pi
     dfPi = hM$dfPi
@@ -429,7 +431,7 @@
     for(iter in seq_len(transient + samples*thin)) {
         if(!identical(updater$Gamma2, FALSE)) {
             out = try(updateGamma2(Z=Z,Gamma=Gamma,iV=iV,iSigma=iSigma,
-                                   Eta=Eta,Lambda=Lambda, X=X,Pi=Pi,
+                                   Eta=Eta,Lambda=Lambda, Loff=Loff,X=X,Pi=Pi,
                                    dfPi=dfPi,Tr=Tr,C=C,rL=hM$rL, iQg=iQg,
                                    mGamma=mGamma,iUGamma=iUGamma),
                       silent = TRUE)
@@ -444,7 +446,7 @@
                                               V=chol2inv(chol(iV)),iV=iV,
                                               id=iSigma, Eta=Eta,
                                               Lambda=Lambda,Alpha=Alpha,
-                                              X=X,Pi=Pi,dfPi=dfPi,Tr=Tr,
+                                              Loff=Loff,X=X,Pi=Pi,dfPi=dfPi,Tr=Tr,
                                               rL=hM$rL, rLPar=rLPar,
                                               Q=Qg[,,rho],iQ=iQg[,,rho],
                                               RQ=RQg[,,rho],
@@ -462,8 +464,8 @@
         if(!identical(updater$BetaLambda, FALSE)){
             BetaLambdaList = try(updateBetaLambda(Y=Y,Z=Z,Gamma=Gamma,iV=iV,
                                                   iSigma=iSigma,Eta=Eta,
-                                                  Psi=Psi,Delta=Delta,
-                                                  iQ=iQg[,,rho],X=X,Tr=Tr,
+                                                  Psi=Psi,Delta=Delta,iQ=iQg[,,rho],
+                                                  Loff=Loff,X=X,Tr=Tr,
                                                   Pi=Pi,dfPi=dfPi,C=C,
                                                   rL=hM$rL),
                                  silent = TRUE)
@@ -477,7 +479,7 @@
 
         if(!identical(updater$wRRR, FALSE) &&  hM$ncRRR>0){
             wRRRXList = try(updatewRRR(Z=Z, Beta=Beta, iSigma=iSigma,
-                                       Eta=Eta, Lambda=Lambda, X1A=X1A,
+                                       Eta=Eta, Lambda=Lambda, Loff=Loff, X1A=X1A,
                                        XRRR=hM$XRRRScaled, Pi=Pi, dfPi=dfPi,
                                        rL = hM$rL, PsiRRR=PsiRRR,
                                        DeltaRRR=DeltaRRR),
@@ -491,10 +493,10 @@
         }
 
         if(!identical(updater$BetaSel, FALSE) &&  hM$ncsel>0){
-            BetaSelXList = try(updateBetaSel(Z=Z,XSelect = hM$XSelect,
+            BetaSelXList = try(updateBetaSel(Z=Z, XSelect=hM$XSelect,
                                              BetaSel=BetaSel,Beta=Beta,
-                                             iSigma=iSigma, Lambda=Lambda,
-                                             Eta=Eta, X1=X1,Pi=Pi,dfPi=dfPi,
+                                             iSigma=iSigma, Lambda=Lambda, Eta=Eta,
+                                             Loff=Loff, X1=X1, Pi=Pi, dfPi=dfPi,
                                              rL=hM$rL),
                                silent = TRUE)
             if (!inherits(BetaSelXList, "try-error")) {
@@ -555,7 +557,7 @@
 
         if(!identical(updater$Eta, FALSE))
             out = try(updateEta(Y=Y,Z=Z,Beta=Beta,iSigma=iSigma,Eta=Eta,
-                                Lambda=Lambda,Alpha=Alpha, rLPar=rLPar, X=X,
+                                Lambda=Lambda,Alpha=Alpha, rLPar=rLPar, Loff=Loff,X=X,
                                 Pi=Pi,dfPi=dfPi,rL=hM$rL), silent = TRUE)
         if (!inherits(out, "try-error"))
             Eta <- out
@@ -572,7 +574,7 @@
 
         if(!identical(updater$InvSigma, FALSE))
             out = try(updateInvSigma(Y=Y,Z=Z,Beta=Beta,iSigma=iSigma,
-                                     Eta=Eta,Lambda=Lambda, distr=distr,X=X,
+                                     Eta=Eta,Lambda=Lambda, distr=distr,Loff=Loff,X=X,
                                      Pi=Pi,dfPi=dfPi,rL=hM$rL, aSigma=aSigma,
                                      bSigma=bSigma), silent = TRUE)
         if (!inherits(out, "try-error"))
@@ -582,7 +584,7 @@
 
         if(!identical(updater$Z, FALSE)) {
             out = try(updateZ(Y=Y,Z=Z,Beta=Beta,iSigma=iSigma,Eta=Eta,
-                              Lambda=Lambda, X=X,Pi=Pi,dfPi=dfPi,distr=distr,
+                              Lambda=Lambda, Loff=Loff,X=X,Pi=Pi,dfPi=dfPi,distr=distr,
                               rL=hM$rL))
             if (!inherits(out, "try-error"))
                 Z = out

diff --git a/R/updateBetaLambda.R b/R/updateBetaLambda.R
@@ -5,14 +5,14 @@
 #' @importFrom stats rnorm
 #' @importFrom Matrix bdiag Diagonal sparseMatrix
 #'
-updateBetaLambda = function(Y,Z,Gamma,iV,iSigma,Eta,Psi,Delta,rho, iQ, X,Tr,Pi,dfPi,C,rL){
+updateBetaLambda = function(Y,Z,Gamma,iV,iSigma,Eta,Psi,Delta,rho, iQ, Loff,X,Tr,Pi,dfPi,C,rL){
    ny = nrow(Z)
    ns = ncol(Z)
    nc = nrow(Gamma)
    nt = ncol(Tr)
    nr = ncol(Pi)
 
-   S = Z
+   if(is.null(Loff)) S = Z else S = Z - Loff
    Lambda = vector("list", nr)
    if(nr > 0){
       EtaFull = vector("list", nr)

diff --git a/R/updateBetaSel.R b/R/updateBetaSel.R
@@ -1,7 +1,7 @@
 #' @importFrom stats dnorm runif
-#' 
-updateBetaSel = function(Z=Z,XSelect, BetaSel, Beta, iSigma,
-                         Lambda, Eta, X1,Pi,dfPi,rL){
+#'
+updateBetaSel = function(Z,XSelect, BetaSel, Beta, iSigma,
+                         Lambda, Eta, Loff,X1,Pi,dfPi,rL){
 
    ny = nrow(Z)
    ns = ncol(Z)
@@ -42,11 +42,8 @@ updateBetaSel = function(Z=Z,XSelect, BetaSel, Beta, iSigma,
    LFix = matrix(NA,ny,ns)
    for(j in 1:ns)
       LFix[,j] = X[[j]]%*%Beta[,j]
-
-   if(nr > 0){
-      E = LFix + Reduce("+", LRan)
-   } else
-      E = LFix
+   E = Reduce("+", c(LFix,LRan))
+   if(!is.null(Loff)) E = E + Loff
 
    ll = matrix(NA,ny,ns)
    for (j in 1:ns){
@@ -58,7 +55,6 @@ updateBetaSel = function(Z=Z,XSelect, BetaSel, Beta, iSigma,
       XSel = XSelect[[i]]
       for (spg in 1:length(XSel$q)){
          BetaSelNew[[i]][spg] = !(BetaSel[[i]][spg])
-
          fsp = which(XSel$spGroup==spg)
 
          X2 = X0