From f26855d7c34b5d98d294873300e2db7be66f22c5 Mon Sep 17 00:00:00 2001 From: Gleb Tikhonov Date: Fri, 22 Mar 2024 11:26:12 +0200 Subject: [PATCH] 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. --- R/Hmsc.R | 14 ++++++++-- R/biPlot.R | 4 +-- R/c.Hmsc.R | 2 +- R/computeInitialParameters.R | 2 +- R/computePredictedValues.R | 7 +++-- R/computePredictor.HmscRandomLevel.R | 5 ++++ R/predict.R | 42 +++++++++++++--------------- R/sampleMcmc.R | 24 ++++++++-------- R/updateBetaLambda.R | 4 +-- R/updateBetaSel.R | 14 ++++------ R/updateEta.R | 10 +++---- R/updateGamma2.R | 8 ++---- R/updateGammaEta.R | 10 +++---- R/updateInvSigma.R | 8 ++---- R/updateZ.R | 8 ++---- R/updatewRRR.R | 11 ++------ tests/testthat/test-sampling.R | 14 +++++----- 17 files changed, 92 insertions(+), 95 deletions(-) diff --git a/R/Hmsc.R b/R/Hmsc.R index a41ec54..049f813 100644 --- 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 index 3fe39e4..1d30279 100644 --- a/R/biPlot.R +++ b/R/biPlot.R @@ -22,7 +22,7 @@ #' # 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 @@ -30,7 +30,7 @@ #' #' @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 index 453dc74..201bea0 100644 --- 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 index 3c83aa3..8adfec0 100644 --- 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 index 19134ec..2d1b620 100644 --- 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 index 053119b..c08aed6 100644 --- 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 index c9e56f6..3cba9a7 100644 --- 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 index ca34b2d..57e0da8 100644 --- 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 index 97fcc9c..7a6a2af 100644 --- 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 index 7b74fb2..3bdf086 100644 --- 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 diff --git a/R/updateEta.R b/R/updateEta.R index af9aad3..63de714 100644 --- a/R/updateEta.R +++ b/R/updateEta.R @@ -1,7 +1,7 @@ #' @importFrom stats rnorm #' @importFrom Matrix bdiag Diagonal sparseMatrix t Matrix #' -updateEta = function(Y,Z,Beta,iSigma,Eta,Lambda,Alpha, rLPar, X,Pi,dfPi,rL){ +updateEta = function(Y,Z,Beta,iSigma,Eta,Lambda,Alpha, rLPar, Loff,X,Pi,dfPi,rL){ ny = nrow(Z) ns = ncol(Z) nr = ncol(Pi) @@ -30,11 +30,9 @@ updateEta = function(Y,Z,Beta,iSigma,Eta,Lambda,Alpha, rLPar, X,Pi,dfPi,rL){ } for(r in seq_len(nr)){ rnames=rownames(Eta[[r]]) - if(nr > 1){ - S = Z - (LFix + Reduce("+", LRan[setdiff(1:nr, r)])) - } else{ - S = Z - LFix - } + S = Z - Reduce("+", c(LFix,LRan[setdiff(1:nr, r)])) + if(!is.null(Loff)) S = S - Loff + lambda = Lambda[[r]] nf = dim(lambda)[1] lPi = Pi[,r] diff --git a/R/updateGamma2.R b/R/updateGamma2.R index 929daca..f4a35c3 100644 --- a/R/updateGamma2.R +++ b/R/updateGamma2.R @@ -3,7 +3,7 @@ #' @importFrom stats rnorm #' -updateGamma2 = function(Z,Gamma=Gamma,iV,iSigma,Eta,Lambda, X,Pi,dfPi,Tr,C,rL, iQg, mGamma,iUGamma){ +updateGamma2 = function(Z,Gamma=Gamma,iV,iSigma,Eta,Lambda, Loff,X,Pi,dfPi,Tr,C,rL, iQg, mGamma,iUGamma){ ns = ncol(Z) ny = nrow(Z) nc = nrow(iV) @@ -27,10 +27,8 @@ updateGamma2 = function(Z,Gamma=Gamma,iV,iSigma,Eta,Lambda, X,Pi,dfPi,Tr,C,rL, i LRan[[r]] = LRan[[r]] + (Eta[[r]][Pi[,r],]*rL[[r]]$x[as.character(dfPi[,r]),k]) %*% Lambda[[r]][,,k] } } - if(nr > 0){ - S = Z - Reduce("+", LRan) - } else - S = Z + if(nr > 0) S = Z - Reduce("+", LRan) else S = Z + if(!is.null(Loff)) S = S - Loff if(is.null(C)){ # if(all(iSigma==1)){ diff --git a/R/updateGammaEta.R b/R/updateGammaEta.R index c5debdc..8c92ec6 100644 --- a/R/updateGammaEta.R +++ b/R/updateGammaEta.R @@ -4,7 +4,7 @@ #' @importFrom stats rnorm #' @importFrom Matrix Diagonal sparseMatrix bdiag #' -updateGammaEta = function(Z,Gamma,V,iV,id,Eta,Lambda,Alpha, X,Tr,Pi,dfPi,rL, rLPar,Q,iQ,RQ,mGamma,U,iU){ +updateGammaEta = function(Z,Gamma,V,iV,id,Eta,Lambda,Alpha, Loff,X,Tr,Pi,dfPi,rL, rLPar,Q,iQ,RQ,mGamma,U,iU){ ny = nrow(Z) ns = ncol(Z) nr = ncol(Pi) @@ -44,11 +44,9 @@ updateGammaEta = function(Z,Gamma,V,iV,id,Eta,Lambda,Alpha, X,Tr,Pi,dfPi,rL, rLP for(r in seq_len(nr)){ if(rL[[r]]$xDim == 0){ - if(nr > 1){ - S = Z - Reduce("+", LRan[setdiff(1:nr,r)]) - } else{ - S = Z - } + if(nr > 1) S = Z - Reduce("+", LRan[setdiff(1:nr,r)]) else S = Z + if(!is.null(Loff)) S = S - Loff + Lam = Lambda[[r]] nf = nrow(Lam) lPi = Pi[,r] diff --git a/R/updateInvSigma.R b/R/updateInvSigma.R index 516eb0e..5e2894d 100644 --- a/R/updateInvSigma.R +++ b/R/updateInvSigma.R @@ -1,6 +1,6 @@ #' @importFrom stats rgamma #' -updateInvSigma = function(Y,Z,Beta,iSigma,Eta,Lambda, distr,X,Pi,dfPi,rL, aSigma,bSigma){ +updateInvSigma = function(Y,Z,Beta,iSigma,Eta,Lambda, distr,Loff,X,Pi,dfPi,rL, aSigma,bSigma){ indVarSigma = (distr[,2]==1) if(any(indVarSigma)){ ny = nrow(Z) @@ -27,10 +27,8 @@ updateInvSigma = function(Y,Z,Beta,iSigma,Eta,Lambda, distr,X,Pi,dfPi,rL, aSigma LRan[[r]] = LRan[[r]] + (Eta[[r]][Pi[,r],]*rL[[r]]$x[as.character(dfPi[,r]),k]) %*% Lambda[[r]][,,k] } } - if(nr > 0){ - Eps = Z - (LFix + Reduce("+", LRan)) - } else - Eps = Z - LFix + Eps = Z - Reduce("+", c(LFix,LRan)) + if(!is.null(Loff)) Eps = Eps - Loff Yx = !is.na(Y) nyx = colSums(Yx) diff --git a/R/updateZ.R b/R/updateZ.R index a827910..61351f6 100644 --- a/R/updateZ.R +++ b/R/updateZ.R @@ -1,7 +1,7 @@ #' @importFrom stats dnorm pnorm rnorm #' @importFrom truncnorm rtruncnorm #' @importFrom BayesLogit rpg -updateZ = function(Y,Z,Beta,iSigma,Eta,Lambda, X,Pi,dfPi,distr,rL, ind){ +updateZ = function(Y,Z,Beta,iSigma,Eta,Lambda, Loff,X,Pi,dfPi,distr,rL, ind){ ZPrev = Z ny = nrow(Y) ns = ncol(Y) @@ -28,10 +28,8 @@ updateZ = function(Y,Z,Beta,iSigma,Eta,Lambda, X,Pi,dfPi,distr,rL, ind){ LRan[[r]] = LRan[[r]] + (Eta[[r]][Pi[,r],]*rL[[r]]$x[as.character(dfPi[,r]),k]) %*% Lambda[[r]][,,k] } } - if(nr > 0){ - E = LFix + Reduce("+", LRan) - } else - E = LFix + E = Reduce("+", c(LFix,LRan)) + if(!is.null(Loff)) E = E + Loff Z = matrix(NA,ny,ns) indNA = is.na(Y) diff --git a/R/updatewRRR.R b/R/updatewRRR.R index d55d12b..808465a 100644 --- a/R/updatewRRR.R +++ b/R/updatewRRR.R @@ -4,7 +4,7 @@ # #' @importFrom stats rnorm # -updatewRRR = function(Z,Beta,iSigma,Eta,Lambda,X1A,XRRR,Pi,dfPi,rL,PsiRRR,DeltaRRR){ +updatewRRR = function(Z,Beta,iSigma,Eta,Lambda, Loff,X1A,XRRR,Pi,dfPi,rL,PsiRRR,DeltaRRR){ ny = nrow(Z) ns = ncol(Z) nr = ncol(Pi) @@ -40,11 +40,8 @@ updatewRRR = function(Z,Beta,iSigma,Eta,Lambda,X1A,XRRR,Pi,dfPi,rL,PsiRRR,DeltaR LRan[[r]] = LRan[[r]] + (Eta[[r]][Pi[,r],]*rL[[r]]$x[as.character(dfPi[,r]),r]) %*% Lambda[[r]][,,k] } } - if(nr > 1){ - S = Z - (LFix + Reduce("+", LRan)) - } else{ - S = Z - LFix - } + S = Z - Reduce("+", c(LFix,LRan)) + if(!is.null(Loff)) S = S - Loff A1 = BetaRRR%*%diag(iSigma,nrow = length(iSigma))%*%t(BetaRRR) A2 = t(XRRR)%*%XRRR @@ -60,7 +57,6 @@ updatewRRR = function(Z,Beta,iSigma,Eta,Lambda,X1A,XRRR,Pi,dfPi,rL,PsiRRR,DeltaR wRRR = matrix(we,nrow = ncRRR) X = X1A - if(ncRRR>0){ XB=XRRR%*%t(wRRR) if(is.matrix(X)){ @@ -75,6 +71,5 @@ updatewRRR = function(Z,Beta,iSigma,Eta,Lambda,X1A,XRRR,Pi,dfPi,rL,PsiRRR,DeltaR wRRRXList=list() wRRRXList$wRRR = wRRR wRRRXList$X = X - return(wRRRXList) } diff --git a/tests/testthat/test-sampling.R b/tests/testthat/test-sampling.R index 2005be1..f64d267 100644 --- a/tests/testthat/test-sampling.R +++ b/tests/testthat/test-sampling.R @@ -5,7 +5,7 @@ test_that("updateGamma2 is correct", { parList = computeInitialParameters(TD$m,initPar=NULL) dataParList = computeDataParameters(TD$m) Gamma = updateGamma2(Z=parList$Z,Gamma=parList$Gamma,iV=chol2inv(chol(parList$V)),iSigma=sqrt(parList$sigma), - Eta=parList$Eta,Lambda=parList$Lambda, X=TD$m$X,Pi=TD$m$Pi,dfPi=TD$m$dfPi,Tr=TD$m$Tr, + Eta=parList$Eta,Lambda=parList$Lambda, Loff=TD$m$Loff,X=TD$m$X,Pi=TD$m$Pi,dfPi=TD$m$dfPi,Tr=TD$m$Tr, C=TD$m$C,rL=TD$m$rL, iQg=dataParList$iQg, mGamma=TD$m$mGamma,iUGamma=chol2inv(chol(TD$m$UGamma))) expect_equal(ncol(Gamma),3) expect_equal(nrow(Gamma),3) @@ -17,7 +17,7 @@ test_that("UpdateGammaEta is correct",{ parList = computeInitialParameters(TD$m,initPar=NULL) dataParList = computeDataParameters(TD$m) GammaEtaList = updateGammaEta(Z=parList$Z,Gamma=parList$Gamma,V=parList$V,iV=chol2inv(chol(parList$V)),id=sqrt(parList$sigma), - Eta=parList$Eta,Lambda=parList$Lambda,Alpha=parList$Alpha, X=TD$m$X,Pi=TD$m$Pi, + Eta=parList$Eta,Lambda=parList$Lambda,Alpha=parList$Alpha, Loff=TD$m$Loff,X=TD$m$X,Pi=TD$m$Pi, dfPi=TD$m$dfPi,Tr=TD$m$Tr,rL=TD$m$rL, rLPar=dataParList$rLPar,Q=dataParList$Qg[,,parList$rho], iQ=dataParList$iQg[,,parList$rho],RQ=dataParList$RQg[,,parList$rho], mGamma=TD$m$mGamma,U=TD$m$UGamma,iU=chol2inv(chol(TD$m$UGamma))) @@ -38,7 +38,7 @@ test_that("updateBetaLambda is correct",{ dataParList = computeDataParameters(TD$m) BetaLambdaList = updateBetaLambda(Y=TD$Y,Z=parList$Z,Gamma=parList$Gamma,iV=chol2inv(chol(parList$V)), iSigma=sqrt(parList$sigma),Eta=parList$Eta,Psi=parList$Psi,Delta=parList$Delta, - iQ=dataParList$iQg[,,parList$rho],X=TD$m$X,Tr=TD$m$Tr,Pi=TD$m$Pi,dfPi=TD$m$dfPi,C=TD$m$C,rL=TD$m$rL) + iQ=dataParList$iQg[,,parList$rho],Loff=TD$m$Loff,X=TD$m$X,Tr=TD$m$Tr,Pi=TD$m$Pi,dfPi=TD$m$dfPi,C=TD$m$C,rL=TD$m$rL) Beta = BetaLambdaList$Beta Lambda = BetaLambdaList$Lambda expect_equal(length(Lambda),2) @@ -97,9 +97,9 @@ test_that("updateEta is correct", { set.seed(200) parList = computeInitialParameters(TD$m,initPar=NULL) dataParList = computeDataParameters(TD$m) - eta = updateEta(Y = TD$m$Y, Z=parList$Z,Beta=parList$Beta,iSigma=sqrt(parList$sigma),Eta=parList$Eta, - Lambda=parList$Lambda, Alpha = parList$Alpha, rLPar = dataParList$rLPar, X = TD$m$X, Pi = TD$m$Pi, - dfPi = TD$m$dfPi,rL=TD$m$rL) + eta = updateEta(Y=TD$m$Y, Z=parList$Z, Beta=parList$Beta,iSigma=sqrt(parList$sigma),Eta=parList$Eta, + Lambda=parList$Lambda, Alpha=parList$Alpha, rLPar=dataParList$rLPar, Loff=TD$m$Loff, X=TD$m$X, Pi=TD$m$Pi, + dfPi=TD$m$dfPi,rL=TD$m$rL) expect_equal(length(eta),2) expect_equal(length(eta[[1]]),100) expect_equal(length(eta[[2]]),20) @@ -126,7 +126,7 @@ test_that("updateInvSigma is correct",{ parList = computeInitialParameters(m,initPar='fixed effects') dataParList = computeDataParameters(m) iSigma = updateInvSigma(Y=m$Y,Z=parList$Z,Beta=parList$Beta,iSigma=sqrt(parList$sigma), - Eta=parList$Eta,Lambda=parList$Lambda, distr=m$distr,X=m$X,Pi=m$Pi, + Eta=parList$Eta,Lambda=parList$Lambda, distr=m$distr,Loff=m$Loff,X=m$X,Pi=m$Pi, dfPi=m$dfPi,rL=m$rL, aSigma=m$aSigma,bSigma=m$bSigma) expect_equal(round(iSigma),c(0,1,0,7)) # 1 and 4 are arbitrary & depend on set.seed() })