Description Usage Arguments Value Examples
View source: R/GibbsBinomialMLB.R
This code implements a version of the collapsed Gibbs sampler from Bradley et al. (2018). The model assumes that data follow a binomial distribution with a logit link to a mixed effects model. The priors and random effects are assumed to follow a multivariate logit-beta distribution.
1 2 | GibbsBinomialMLB(B, data, nn, X, G, report = 100, rho = 0.999,
cs = c(0, 0, 0), eps = 0.01, nslice = 2, select.h = 0)
|
B |
The number of iterations of the collapsed Gibbs sampler |
data |
An n dimensional vector consisting of binomial count-valued observations. |
nn |
An n dimensional vector consisting of totals associated with binomial count-valued observations. |
X |
An nxp matrix of covariates. Each column represents a covariate. Each row corresponds to one of the n replicates. |
G |
An nxr matrix of basis functions. Each column represents a basis function. Each row corresponds to one of the n replicates. |
report |
Option to print the iteration number of the MCMC. |
rho |
A parameter between zero and one to allow for zero counts |
cs |
A three dimensional parameter vector to allow prior to have non-zero means. |
eps |
A second parameter between zero and one to allow for zero counts |
nslice |
The burnin for the slice sampler |
select.h |
Option to select pseudo-optimal values for rho and eps. The MCMC is run for only 50 iterations, and the corresponding predictions are compared to a holdout datasets. The values of rho and eps are chosen to minimize the Hellinger distance between the holdout and the naive predictions using the function optim. |
beta A pxB matrix of MCMC values for beta. The vector beta corresponds to the covariates X.
eta A rxB matrix of MCMC values for eta. The vector eta corresponds to the basis functions G.
xi A nxB matrix of MCMC values for xi. The vector xi corresponds to uncorrelated random effects.
pi A nxB matrix of MCMC values for the proportions of the binomial data, modeled with, exp(Xbeta+Geta+xi)/(1+exp(Xbeta+Geta+xi)).
pi.smooth A nxB matrix of MCMC values for the smoothed proportions associated with binomial data, modeled with, exp(Xbeta+Geta)/(1+exp(Xbeta+Geta)).
alphab A 1xB matrix of MCMC values for the shape parameter associated with beta.
alphae A 1xB matrix of MCMC values for the shape parameter associated with eta.
alphax A 1xB matrix of MCMC values for the shape parameter associated with xi.
kappab A 1xB matrix of MCMC values for the second shape parameter associated with beta.
kappae A 1xB matrix of MCMC values for the second shape parameter associated with eta.
kappax A 1xB matrix of MCMC values for the second shape parameter associated with xi
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 | library(CM)
set.seed(2)
#simulate pseudo data
n = 1000
X = cbind(1,sin(pi*seq(0,2,by = ((1 - 0)/(n - 1)))))
X=X[1:n,]
p = dim(X)[2]
output<-qr(X)
Q = qr.Q(output)
#Moran's I basis functions (Hughes and Haran, 2013)
r = 100
Psi = MoransI.Basis(Q,r,diag(n))
#large scale parameters
beta = c(0.01,-2)
mu = X%*%beta
plot(mu)
#small scale parameters
eta = sqrt(as.numeric(var(X%*%beta)/2))*rnorm(r)
rand.effects = Psi%*%eta
plot(rand.effects)
#fine scale parameters
xi = sqrt(as.numeric(var(Psi%*%eta)/2))*rnorm(n)
#simulate binomial data
nu = mu+rand.effects+xi
plot(nu)
#inverse logit
pi1 = exp(nu[1:n])/(1+exp(nu[1:n]))
#generate sample sizes
m = rpois(n,50)
data = matrix(0,n,1)
for (j in 1:n){
data[j]=rmultinom(1, m[j], prob = c(pi1[j],1-pi1[j]))[1]
}
plot(data)
#number of MCMC reps
B = 2000
burnin=1000
#Pseudo-Code 2 with updates for shape parameters
output<-GibbsBinomialMLB(B,data,m,X,Psi)
#check some trace plots
plot((as.mcmc(output$beta[1,burnin:B])))
plot((as.mcmc(output$beta[2,burnin:B])))
plot(as.mcmc(output$alphab[1,burnin:B]))
plot(as.mcmc(output$alphae[1,burnin:B]))
plot(as.mcmc(output$alphax[1,burnin:B]))
plot(as.mcmc(output$kappab[1,burnin:B]))
plot(as.mcmc(output$kappae[1,burnin:B]))
plot((as.mcmc(output$eta[5,burnin:B])))
plot((as.mcmc(output$xi[10,burnin:B])))
#hows our estimate of beta? true beta is (0.01,-2)
apply(output$beta[,burnin:B], 1, quantile,probs=c(0.025,0.5,0.975))
#hows our estimate of the proportions?
pihat = rowMeans(output$pi)
pibounds = apply(output$pi, 1, quantile, probs = c(0.025,0.975))
pilow = as.matrix(pibounds[1,],n,1)
piup = as.matrix(pibounds[2,],n,1)
#figure
plot(pi1,pihat,ylab=c("Estimate Proportion"),xlab=c("True Proportion"),main=c("Estimated Versus Truth"))
abline(0,1)
plot(1:n,pi1, type = "l", col = "black",xlab = "Arbitrary Ordering of Proportions",ylab = "Proportions")
lines(1:n,pihat, type = "l", col = "red",lty=1)
legend("topright", c("Truth","Estimated"),col=c("black", "red"), lty=c(1,1))
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.