R/PoiSimMixModSMFVB.r In densEstBayes: Density Estimation via Bayesian Inference Engines

Defines functions PoiSimMixModSMFVB

```########## R function: PoiSimMixModSMFVB ##########

# For fitting and inference in a simple (single variance
# parameter) general design Poisson mixed model via
# mean field variational Bayes.

# Last changed: 14 SEp 2020

PoiSimMixModSMFVB <- function(y,X,Z,sigmabeta,ssigma,convToler,maxIter)
{
# Set required dimension variables:

ncX <- ncol(X)  ;  ncZ <- ncol(Z)  ;  ncC <- ncX + ncZ

# Set constant variational parameter values:

kappa.q.sigsq <- 0.5*(ncZ + 1)
kappa.q.a <- 1

# Set initial values of non-constant variational parameters:

groupVec <- as.factor(rep(1,length(y)))
fitPQL <- MASS::glmmPQL(y~-1+X,random=list(groupVec=nlme::pdIdent(~-1+Z)),family=poisson,
verbose=FALSE)

mu.q.betau <- c(fitPQL\$coef\$fixed,unlist(fitPQL\$coef\$random))
Sigma.q.betau <- diag(ncC)
mu.q.recip.sigsq <- as.numeric(exp(-2*unlist(fitPQL\$modelStruct)))

# Form the C matrix:

Cmat <- cbind(X,Z)

# Perform iterations:

converged <- FALSE
itNum <- 0
while (!converged)
{
itNum <- itNum + 1

# Update q(beta,u) parameters:

w.q.betau <- as.vector(exp(crossprod(t(Cmat),mu.q.betau) + 0.5*rowSums((crossprod(t(Cmat),Sigma.q.betau)*Cmat))))
M.q.recip.sigsq <- diag(c(rep((1/sigmabeta^2),ncX),rep(mu.q.recip.sigsq,ncZ)))
Sigma.q.betau <- try(solve(crossprod(Cmat,(w.q.betau*Cmat)) + M.q.recip.sigsq),silent=TRUE)

# Check whether solve() successfuly or not:

failedSolveCall <-  is.null(dim(Sigma.q.betau))
if (failedSolveCall)
{
errMsg1 <- "\n\n Call to densEstBayes() with method=\"SMFVB\" has failed due\n"
errMsg2 <- "to a singular matrix and/or non-convergence problem. Use\n"
errMsg3 <- "method=\"slice\", method=\"NUTS\" or method=\"HMC\" instead.\n\n"
stop(paste(errMsg1,errMsg2,errMsg3),call.=FALSE)
}

if(!failedSolveCall)
{
mu.q.betau <- mu.q.betau + crossprod(Sigma.q.betau,(crossprod(Cmat,(y-w.q.betau)) - M.q.recip.sigsq%*%mu.q.betau))

# Update q(a) parameters:

lambda.q.a <- mu.q.recip.sigsq + (1/ssigma^2)
mu.q.recip.a <- kappa.q.a/lambda.q.a

# Update q(sigsq) parameters:

lambda.q.sigsq <- mu.q.recip.a + 0.5*(sum(mu.q.betau[-(1:ncX)]^2) + sum(diag(Sigma.q.betau)[-(1:ncX)]))

mu.q.recip.sigsq <- kappa.q.sigsq/lambda.q.sigsq

if (itNum>1)
relErr <- abs((mu.q.recip.sigsq - mu.q.recip.sigsq.OLD)/mu.q.recip.sigsq.OLD)

mu.q.recip.sigsq.OLD <- mu.q.recip.sigsq

# Do convergence check:

if (itNum>1)
if (relErr<convToler)
converged <- TRUE

if (itNum>=maxIter)
{
converged <- TRUE
warning("Semiparametric mean field variational Bayes did not converge after maximum number of iterations.\n")
}
}
}

# Return q-density parameters:

return(list(mu.q.betau=mu.q.betau,Sigma.q.betau=Sigma.q.betau,
kappa.q.sigsq=kappa.q.sigsq,lambda.q.sigsq=lambda.q.sigsq))
}

############ End of PoiSimMixModSMFVB ############
```

Try the densEstBayes package in your browser

Any scripts or data that you put into this service are public.

densEstBayes documentation built on Aug. 19, 2021, 9:06 a.m.