R: Fitting Gaussian, Poisson, and Binomial Hierarchical Models
gbp
R Documentation
Fitting Gaussian, Poisson, and Binomial Hierarchical Models
Description
gbp fits Bayesian hierarchical models using the Uniform distribution on the second level variance component (variance of the prior distribution), which enables good frequentist repeated sampling properties.
a (k by 1) vector of k groups' sample means for Gaussian or of each group's number of successful trials for Binomial and Poisson data, where k is the number of groups (or units) in a dataset.
se.or.n
a (k by 1) vector composed of the standard errors of all groups for Gaussian or
of each group's total number of trials for Binomial and Poisson data.
covariates
(optional) a matrix of covariate(s) each column of which corresponds to one covariate.
mean.PriorDist
(optional) a numeric value for the second-level mean parameter, i.e. the mean of prior
distribution, if you know this value a priori.
model
a character string indicating which hierarchical model to fit. "gaussian" for Gaussian data,
"poisson" for Poisson, and "binomial" for Binomial. Default is "gaussian"
intercept
TRUE or FALSE flag indicating whether an intercept should be
included in the regression. Default is TRUE.
confidence.lvl
a float between 0 and 1 to estimate 100*confidence.lvl% intervals. Default is 0.95.
n.AR
Only for Binomial model. If n.AR = 1000, all the results will be based on 1000 posterior samples using the acceptance-rejection sampling. Default is 0.
n.AR.factor
Only for Binomial model. If n.AR = 1000 and n.AR.factor = 4, then the acceptance-rejection method will sample 4 * 1000 trial samples, and accept or reject them until the method gets n.AR posterior samples. Default is 4.
trial.scale
A scale used in the trial distribution of α. If resultant weight has too mamy 0's, scale should be smaller than before. If resultant weight has too few 0's, scale should be larger than before. If there are relatively huge weights, scale should be larger than before. Default is NA.
save.result
Only for Binomial model with the acceptance-rejection sampling. If save.result is TRUE, all the results of weights and posterior samples will be saved in the gbp object. Default is TRUE.
normal.CI
Only applicable for Gaussian data. If TRUE then a Normal
approximation is used to construct intervals for the level 1 group
means. If FALSE (default) then a Skew-Normal distribution is
used. Setting the value to TRUE may result in speed
improvements but may lead to intervals that under cover.
t
Non-negative constant to determine the hyper-prior distribution of r for the Binomial model with the acceptance-rejection method. If t is positive, then the hyper-prior distribution of r is proper, otherwise improper. dr/(t + r)^{u+1}.
u
A positive constant to determine the hyper-prior distribution of r for the Binomial model with the acceptance-rejection method. dr/(t + r)^{u+1}.
Details
gbp fits a hierarchical model whose first-level hierarchy is a distribution of observed data and second-level is a conjugate prior distribution on the first-level parameter. To be specific, for Normal data, gbp constructs a two-level Normal-Normal multilevel model. V_j (=σ^2/n_j) is assumed to be known or to be accurately estimated, and subscript j indicates j-th group
(or unit) in a dataset.
(y_j | θ_j) ~ indep N(θ_j, σ_j^2)
(θ_j | μ0_j, A) ~ indep N(μ0_j, A)
μ0_j = x_j'β
for j = 1, …, k, where k is the number of groups (units) in a dataset.
For Poisson data, gbp builds a two-level Poisson-Gamma multilevel model. A square bracket below indicates [mean, variance] of distribution, a constant multiplied to the notation representing Gamma distribution (Gam) is a scale, and y_j = z_j / n_j.
for j = 1, …, k, where k is the number of groups (units) in a dataset.
For Binomial data, gbp sets a two-level Binomial-Beta multilevel model. A square bracket below indicates [mean, variance] of distribution and y_j = z_j / n_j.
for j = 1, …, k, where k is the number of groups (units) in a dataset.
For reference, based on the above notations, the Uniform prior distribution on the second level variance component (variance of the prior distribution) is dA for Gaussian and d(1 / r)
(= dr / r^2) for Binomial and Poisson data. The second level variance component can be interpreted as variation among the first-level parameters (θ_j) or variance of ensemble information.
Under this setting, the argument y in gbp is a (k by 1) vector of k groups' sample means (y_j's in the description of Normal-Normal model above) for Gaussian or of each group's number of successful trials (z_j's) for Binomial and Poisson data, where k is the number of groups (or units) in a dataset.
The argument se.or.n is a (k by 1) vector composed of the standard errors (V_j's) of all groups for Gaussian or of each group's total number of trials (n_j's) for Binomial and Poisson data.
As for two optional arguments, covariates and mean.PriorDist, there are three feasible
combinations of them to run gbp. The first situation is when we do not have any covariate and do not
know a mean of the prior distribution (μ0) a priori. In this case, assigning none of two
optional arguments, such as "gbp(z, n, model = "binomial")", will lead to a correct model. gbp
will automatically fit a regression with only an intercept term to estimate a common mean of the prior
distribution (exchangeability).
The second situation is when we have covariate(s) and do not know a mean of the prior distribution (μ0) a priori. In this case, assigning a matrix, X, each column of which corresponds to one covariate, such as "gbp(z, n, X, model = "poisson")", will lead to a correct model. Default of gbp is to fit a regression including an intercept term to estimate a mean of the prior distribution. Double exchangeability will hold in this case.
The last case is when we know a mean of the prior distribution (μ0) a priori. Now, we do
not need to estimate regression coefficients at all because we know a true value of μ0 (strong assumption).
Designating this value into the argument of gbp like
"gbp(y, se, mean.PriorDist = 3)" is enough to account for it. For reference,
mean.PriorDist has a stronger priority than covariates, which means that when both
arguments are designated, gbp will fit a hierarchical model using the known mean of prior distribution, mean.PriorDist.
gbp returns an object of class "gbp" which provides three relevant functions plot, print, and summary.
Value
An object of class gbp comprises of:
sample.mean
sample mean of each group (or unit)
se
if Gaussian data, standard error of sample mean in each group (or unit)
n
if Binomial and Poisson data, total number of trials of each group (or unit)
prior.mean
numeric if entered, NA if not entered
prior.mean.hat
estimate of prior mean by a regression if prior mean is not assigned a priori
shrinkage
shrinkage estimate of each group (adjusted posterior mean)
sd.shrinkage
posterior standard deviation of shrinkage
post.mean
posterior mean of each group
post.sd
posterior standard deviation of each group
post.intv.low
lower bound of 100*confidence.lvl% posterior interval (quantile of posterior distribution)
post.intv.upp
upper bound of 100*confidence.lvl% posterior interval (quantile of posterior distribution)
model
"gaussian" for Gaussian, "poisson" for Poisson, and "binomial" for Binomial data
X
a covariate vector or matrix if designated. NA if not
beta.new
regression coefficient estimates
beta.var
estimated variance matrix of regression coefficient
intercept
whether TRUE or FALSE
a.new
a posterior mode of α defined as log(A) for Gaussian or log(1 / r) for Binomial and Poisson data. Practical meaning (variation of ensemble information) of estimating α will appear in summary(gbp.object).
a.var
posterior variance of α
confidence.lvl
confidence level based on which confidence interval is constructed
weight
(only for Binomial model) weights for acceptance-rejection method
p
(only for Binomial model) posterior samples of p based on the acceptance-rejection method, if this method was used. This is a (k by nsim) matrix. Each row contains posteiror samples of each random effect.
alpha
(only for Binomial model) posterior samples of alpha based on the acceptance-rejection method, if this method was used
beta
(only for Binomial model) posterior samples of beta based on the acceptance-rejection method, if this method was used
accept.rate
(only for Binomial model) the acceptance rate of the acceptance-rejection method, if this method was used
n.AR
(Only for Binomial model) If n.AR = 1000, all the results will be based on 1000 posterior samples using the acceptance-rejection sampling. Default is 0.
n.AR.factor
(only for Binomial model) If n.AR = 1000 and n.AR.factor = 4, then the acceptance-rejection method will sample 4 * 1000 trial posteior samples, and accept or reject them. Default is 4.
Author(s)
Hyungsuk Tak, Joseph Kelly, and Carl Morris
References
Morris, C. and Lysy, M. (2012). Shrinkage Estimation in Multilevel Normal Models.
Statistical Science. 27. 115-134.
Examples
# Loading datasets
data(schools)
y <- schools$y
se <- schools$se
# Arbitrary covariate for schools data
x2 <- rep(c(-1, 0, 1, 2), 2)
# baseball data where z is Hits and n is AtBats
z <- c(18, 17, 16, 15, 14, 14, 13, 12, 11, 11, 10, 10, 10, 10, 10, 9, 8, 7)
n <- c(45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45)
# One covariate: 1 if a player is an outfielder and 0 otherwise
x1 <- c(1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0)
################################################################
# Gaussian Regression Interactive Multilevel Modeling (GRIMM) #
################################################################
####################################################################################
# If we do not have any covariate and do not know a mean of the prior distribution #
####################################################################################
g <- gbp(y, se, model = "gaussian")
g
print(g, sort = FALSE)
summary(g)
plot(g)
plot(g, sort = FALSE)
### when we want to simulate pseudo datasets considering the estimated values
### as true ones.
gcv <- coverage(g, nsim = 10)
### gcv$coverageRB, gcv$coverage10, gcv$average.coverageRB, gcv$average.coverage10,
### gcv$minimum.coverageRB, gcv$minimum.coverage10, gcv$raw.resultRB, gcv$raw.result10.
### when we want to simulate pseudo datasets based on different values of A
### and a regression coefficient (intercept),
### not using estimated values as true ones.
gcv <- coverage(g, A.or.r = 9, reg.coef = 10, nsim = 10)
##################################################################################
# If we have one covariate and do not know a mean of the prior distribution yet, #
##################################################################################
g <- gbp(y, se, x2, model = "gaussian")
g
print(g, sort = FALSE)
summary(g)
plot(g)
plot(g, sort = FALSE)
### when we want to simulate pseudo datasets considering the estimated values
### as true ones.
gcv <- coverage(g, nsim = 10)
### gcv$coverageRB, gcv$coverage10, gcv$average.coverageRB, gcv$average.coverage10,
### gcv$minimum.coverageRB, gcv$minimum.coverage10, gcv$raw.resultRB, gcv$raw.result10.
### when we want to simulate pseudo datasets based on different values of A
### and regression coefficients, not using estimated values
### as true ones. Two values of reg.coef are for beta0 and beta1
gcv <- coverage(g, A.or.r = 9, reg.coef = c(10, 1), nsim = 10)
################################################
# If we know a mean of the prior distribution, #
################################################
g <- gbp(y, se, mean.PriorDist = 8, model = "gaussian")
g
print(g, sort = FALSE)
summary(g)
plot(g)
plot(g, sort = FALSE)
### when we want to simulate pseudo datasets considering the estimated values
### as true ones.
gcv <- coverage(g, nsim = 10)
### gcv$coverageRB, gcv$coverage10, gcv$average.coverageRB, gcv$average.coverage10,
### gcv$minimum.coverageRB, gcv$minimum.coverage10, gcv$raw.resultRB, gcv$raw.result10.
### when we want to simulate pseudo datasets based on different values of A and
### 2nd level mean as true ones, not using estimated values as true ones.
coverage(g, A.or.r = 9, mean.PriorDist = 5, nsim = 10)
###############################################################
# Binomial Regression Interactive Multilevel Modeling (BRIMM) #
###############################################################
####################################################################################
# If we do not have any covariate and do not know a mean of the prior distribution #
####################################################################################
b <- gbp(z, n, model = "binomial")
b
print(b, sort = FALSE)
summary(b)
plot(b)
plot(b, sort = FALSE)
### when we want to simulate pseudo datasets considering the estimated values
### as true ones.
bcv <- coverage(b, nsim = 10)
### bcv$coverageRB, bcv$coverage10, bcv$average.coverageRB, bcv$average.coverage10,
### bcv$minimum.coverageRB, bcv$minimum.coverage10, bcv$raw.resultRB, bcv$raw.result10.
### when we want to simulate pseudo datasets based on different values of r
### and a regression coefficient (intercept),
### not using estimated values as true ones.
bcv <- coverage(b, A.or.r = 60, reg.coef = -1, nsim = 10)
##################################################################################
# If we have one covariate and do not know a mean of the prior distribution yet, #
##################################################################################
b <- gbp(z, n, x1, model = "binomial")
b
print(b, sort = FALSE)
summary(b)
plot(b)
plot(b, sort = FALSE)
### when we want to simulate pseudo datasets considering the estimated values
### as true ones.
bcv <- coverage(b, nsim = 10)
### bcv$coverageRB, bcv$coverage10, bcv$average.coverageRB, bcv$average.coverage10,
### bcv$minimum.coverageRB, bcv$minimum.coverage10, bcv$raw.resultRB, bcv$raw.result10.
### when we want to simulate pseudo datasets based on different values of r
### and regression coefficients, not using estimated values
### as true ones. Two values of reg.coef are for beta0 and beta1
bcv <- coverage(b, A.or.r = 60, reg.coef = c(-1, 0), nsim = 10)
################################################
# If we know a mean of the prior distribution, #
################################################
b <- gbp(z, n, mean.PriorDist = 0.265, model = "binomial")
b
print(b, sort = FALSE)
summary(b)
plot(b)
plot(b, sort = FALSE)
### when we want to simulate pseudo datasets considering the estimated values
### as true ones.
bcv <- coverage(b, nsim = 10)
### bcv$coverageRB, bcv$coverage10, bcv$average.coverageRB, bcv$average.coverage10,
### bcv$minimum.coverageRB, bcv$minimum.coverage10, bcv$raw.resultRB, bcv$raw.result10.
### when we want to simulate pseudo datasets based on different values of r and
### 2nd level mean as true ones, not using estimated values as true ones.
bcv <- coverage(b, A.or.r = 60, mean.PriorDist = 0.3, nsim = 10)
##############################################################
# Poisson Regression Interactive Multilevel Modeling (PRIMM) #
##############################################################
################################################
# If we know a mean of the prior distribution, #
################################################
p <- gbp(z, n, mean.PriorDist = 0.265, model = "poisson")
p
print(p, sort = FALSE)
summary(p)
plot(p)
plot(p, sort = FALSE)
### when we want to simulate pseudo datasets considering the estimated values
### as true ones.
pcv <- coverage(p, nsim = 10)
### pcv$coverageRB, pcv$coverage10, pcv$average.coverageRB, pcv$average.coverage10,
### pcv$minimum.coverageRB, pcv$minimum.coverage10, pcv$raw.resultRB, pcv$raw.result10.
### when we want to simulate pseudo datasets based on different values of r and
### 2nd level mean as true ones, not using estimated values as true ones.
pcv <- coverage(p, A.or.r = 60, mean.PriorDist = 0.3, nsim = 10)
Results
R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(Rgbp)
Loading required package: sn
Loading required package: stats4
Attaching package: 'sn'
The following object is masked from 'package:stats':
sd
Loading required package: mnormt
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/Rgbp/gbp.Rd_%03d_medium.png", width=480, height=480)
> ### Name: gbp
> ### Title: Fitting Gaussian, Poisson, and Binomial Hierarchical Models
> ### Aliases: gbp
> ### Keywords: methods
>
> ### ** Examples
>
>
> # Loading datasets
> data(schools)
> y <- schools$y
> se <- schools$se
>
> # Arbitrary covariate for schools data
> x2 <- rep(c(-1, 0, 1, 2), 2)
>
> # baseball data where z is Hits and n is AtBats
> z <- c(18, 17, 16, 15, 14, 14, 13, 12, 11, 11, 10, 10, 10, 10, 10, 9, 8, 7)
> n <- c(45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45)
>
> # One covariate: 1 if a player is an outfielder and 0 otherwise
> x1 <- c(1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0)
>
> ################################################################
> # Gaussian Regression Interactive Multilevel Modeling (GRIMM) #
> ################################################################
>
> ####################################################################################
> # If we do not have any covariate and do not know a mean of the prior distribution #
> ####################################################################################
>
> g <- gbp(y, se, model = "gaussian")
> g
Summary for each group (sorted by the descending order of se):
obs.mean se prior.mean shrinkage low.intv post.mean upp.intv post.sd
8 12.00 18.0 8.17 0.734 -10.21 9.19 29.9 10.23
3 -3.00 16.0 8.17 0.685 -17.13 4.65 22.5 10.10
1 28.00 15.0 8.17 0.657 -2.32 14.98 38.8 10.56
4 7.00 11.0 8.17 0.507 -8.78 7.59 23.6 8.26
6 1.00 11.0 8.17 0.507 -13.03 4.63 20.1 8.44
2 8.00 10.0 8.17 0.459 -7.25 8.08 23.4 7.81
7 18.00 10.0 8.17 0.459 -1.29 13.48 30.8 8.18
5 -1.00 9.0 8.17 0.408 -13.30 2.74 16.7 7.63
Mean 12.5 8.17 0.552 -9.16 8.17 25.7 8.90
> print(g, sort = FALSE)
Summary for each group:
obs.mean se prior.mean shrinkage low.intv post.mean upp.intv post.sd
1 28.00 15.0 8.17 0.657 -2.32 14.98 38.8 10.56
2 8.00 10.0 8.17 0.459 -7.25 8.08 23.4 7.81
3 -3.00 16.0 8.17 0.685 -17.13 4.65 22.5 10.10
4 7.00 11.0 8.17 0.507 -8.78 7.59 23.6 8.26
5 -1.00 9.0 8.17 0.408 -13.30 2.74 16.7 7.63
6 1.00 11.0 8.17 0.507 -13.03 4.63 20.1 8.44
7 18.00 10.0 8.17 0.459 -1.29 13.48 30.8 8.18
8 12.00 18.0 8.17 0.734 -10.21 9.19 29.9 10.23
Mean 12.5 8.17 0.552 -9.16 8.17 25.7 8.90
> summary(g)
Main summary:
obs.mean se prior.mean shrinkage low.intv post.mean
Group with min(se) -1.00 9.0 8.17 0.408 -13.30 2.74
Group with median(se)1 1.00 11.0 8.17 0.507 -13.03 4.63
Group with median(se)2 7.00 11.0 8.17 0.507 -8.78 7.59
Group with max(se) 12.00 18.0 8.17 0.734 -10.21 9.19
Overall Mean 12.5 8.17 0.552 -9.16 8.17
upp.intv post.sd
Group with min(se) 16.7 7.63
Group with median(se)1 20.1 8.44
Group with median(se)2 23.6 8.26
Group with max(se) 29.9 10.23
Overall Mean 25.7 8.90
Estimation summary for the second-level variance component:
alpha = log(A) for Gaussian or alpha = log(1/r) for Binomial and Poisson data:
post.mode.alpha post.sd.alpha post.mode.A
4.77 1.14 118
Estimation summary for the regression coefficient :
estimate se z.val p.val
beta1 8.168 5.73 1.425 0.154
> plot(g)
> plot(g, sort = FALSE)
>
> ### when we want to simulate pseudo datasets considering the estimated values
> ### as true ones.
> gcv <- coverage(g, nsim = 10)
>
> ### gcv$coverageRB, gcv$coverage10, gcv$average.coverageRB, gcv$average.coverage10,
> ### gcv$minimum.coverageRB, gcv$minimum.coverage10, gcv$raw.resultRB, gcv$raw.result10.
>
> ### when we want to simulate pseudo datasets based on different values of A
> ### and a regression coefficient (intercept),
> ### not using estimated values as true ones.
> gcv <- coverage(g, A.or.r = 9, reg.coef = 10, nsim = 10)
>
> ##################################################################################
> # If we have one covariate and do not know a mean of the prior distribution yet, #
> ##################################################################################
>
> g <- gbp(y, se, x2, model = "gaussian")
> g
Summary for each group (sorted by the descending order of se):
obs.mean se X1 prior.mean shrinkage low.intv post.mean upp.intv
8 12.00 18.0 2.0 8.45 0.666 -14.83 9.64 35.2
3 -3.00 16.0 1.0 8.32 0.612 -20.05 3.93 23.9
1 28.00 15.0 -1.0 8.07 0.581 -5.14 16.43 43.4
4 7.00 11.0 2.0 8.45 0.427 -11.30 7.62 26.2
6 1.00 11.0 0.0 8.20 0.427 -14.51 4.07 20.8
2 8.00 10.0 0.0 8.20 0.381 -8.15 8.08 24.3
7 18.00 10.0 1.0 8.32 0.381 -1.59 14.31 32.3
5 -1.00 9.0 -1.0 8.07 0.333 -15.14 2.02 17.7
Mean 12.5 0.5 8.26 0.476 -11.34 8.26 28.0
post.sd
8 12.74
3 11.20
1 12.39
4 9.56
6 8.98
2 8.27
7 8.63
5 8.35
Mean 10.02
> print(g, sort = FALSE)
Summary for each group:
obs.mean se X1 prior.mean shrinkage low.intv post.mean upp.intv
1 28.00 15.0 -1.0 8.07 0.581 -5.14 16.43 43.4
2 8.00 10.0 0.0 8.20 0.381 -8.15 8.08 24.3
3 -3.00 16.0 1.0 8.32 0.612 -20.05 3.93 23.9
4 7.00 11.0 2.0 8.45 0.427 -11.30 7.62 26.2
5 -1.00 9.0 -1.0 8.07 0.333 -15.14 2.02 17.7
6 1.00 11.0 0.0 8.20 0.427 -14.51 4.07 20.8
7 18.00 10.0 1.0 8.32 0.381 -1.59 14.31 32.3
8 12.00 18.0 2.0 8.45 0.666 -14.83 9.64 35.2
Mean 12.5 0.5 8.26 0.476 -11.34 8.26 28.0
post.sd
1 12.39
2 8.27
3 11.20
4 9.56
5 8.35
6 8.98
7 8.63
8 12.74
Mean 10.02
> summary(g)
Main summary:
obs.mean se X1 prior.mean shrinkage low.intv
Group with min(se) -1.00 9.0 -1.0 8.07 0.333 -15.1
Group with median(se)1 1.00 11.0 0.0 8.20 0.427 -14.5
Group with median(se)2 7.00 11.0 2.0 8.45 0.427 -11.3
Group with max(se) 12.00 18.0 2.0 8.45 0.666 -14.8
Overall Mean 12.5 0.5 8.26 0.476 -11.3
post.mean upp.intv post.sd
Group with min(se) 2.02 17.7 8.35
Group with median(se)1 4.07 20.8 8.98
Group with median(se)2 7.62 26.2 9.56
Group with max(se) 9.64 35.2 12.74
Overall Mean 8.26 28.0 10.02
Estimation summary for the second-level variance component:
alpha = log(A) for Gaussian or alpha = log(1/r) for Binomial and Poisson data:
post.mode.alpha post.sd.alpha post.mode.A
5.09 1.17 162
Estimation summary for the regression coefficient :
estimate se z.val p.val
beta1 8.198 6.656 1.232 0.218
beta2 0.126 5.698 0.022 0.982
> plot(g)
> plot(g, sort = FALSE)
>
> ### when we want to simulate pseudo datasets considering the estimated values
> ### as true ones.
> gcv <- coverage(g, nsim = 10)
>
> ### gcv$coverageRB, gcv$coverage10, gcv$average.coverageRB, gcv$average.coverage10,
> ### gcv$minimum.coverageRB, gcv$minimum.coverage10, gcv$raw.resultRB, gcv$raw.result10.
>
> ### when we want to simulate pseudo datasets based on different values of A
> ### and regression coefficients, not using estimated values
> ### as true ones. Two values of reg.coef are for beta0 and beta1
> gcv <- coverage(g, A.or.r = 9, reg.coef = c(10, 1), nsim = 10)
>
> ################################################
> # If we know a mean of the prior distribution, #
> ################################################
>
> g <- gbp(y, se, mean.PriorDist = 8, model = "gaussian")
> g
Summary for each group (sorted by the descending order of se):
obs.mean se prior.mean shrinkage low.intv post.mean upp.intv post.sd
8 12.00 18.0 8 0.786 -6.767 8.86 26.1 8.36
3 -3.00 16.0 8 0.743 -13.448 5.18 19.3 8.38
1 28.00 15.0 8 0.718 0.595 13.64 34.6 8.95
4 7.00 11.0 8 0.578 -6.642 7.58 21.4 7.15
6 1.00 11.0 8 0.578 -10.701 5.04 18.1 7.34
2 8.00 10.0 8 0.531 -5.426 8.00 21.4 6.85
7 18.00 10.0 8 0.531 0.123 12.69 28.6 7.27
5 -1.00 9.0 8 0.478 -11.516 3.30 15.4 6.87
Mean 12.5 8 0.618 -6.723 8.04 23.1 7.65
> print(g, sort = FALSE)
Summary for each group:
obs.mean se prior.mean shrinkage low.intv post.mean upp.intv post.sd
1 28.00 15.0 8 0.718 0.595 13.64 34.6 8.95
2 8.00 10.0 8 0.531 -5.426 8.00 21.4 6.85
3 -3.00 16.0 8 0.743 -13.448 5.18 19.3 8.38
4 7.00 11.0 8 0.578 -6.642 7.58 21.4 7.15
5 -1.00 9.0 8 0.478 -11.516 3.30 15.4 6.87
6 1.00 11.0 8 0.578 -10.701 5.04 18.1 7.34
7 18.00 10.0 8 0.531 0.123 12.69 28.6 7.27
8 12.00 18.0 8 0.786 -6.767 8.86 26.1 8.36
Mean 12.5 8 0.618 -6.723 8.04 23.1 7.65
> summary(g)
Main summary:
obs.mean se prior.mean shrinkage low.intv post.mean
Group with min(se) -1.00 9.0 8 0.478 -11.52 3.30
Group with median(se)1 1.00 11.0 8 0.578 -10.70 5.04
Group with median(se)2 7.00 11.0 8 0.578 -6.64 7.58
Group with max(se) 12.00 18.0 8 0.786 -6.77 8.86
Overall Mean 12.5 8 0.618 -6.72 8.04
upp.intv post.sd
Group with min(se) 15.4 6.87
Group with median(se)1 18.1 7.34
Group with median(se)2 21.4 7.15
Group with max(se) 26.1 8.36
Overall Mean 23.1 7.65
Estimation summary for the second-level variance component:
alpha = log(A) for Gaussian or alpha = log(1/r) for Binomial and Poisson data:
post.mode.alpha post.sd.alpha post.mode.A
4.48 1.13 88.4
> plot(g)
> plot(g, sort = FALSE)
>
> ### when we want to simulate pseudo datasets considering the estimated values
> ### as true ones.
> gcv <- coverage(g, nsim = 10)
>
> ### gcv$coverageRB, gcv$coverage10, gcv$average.coverageRB, gcv$average.coverage10,
> ### gcv$minimum.coverageRB, gcv$minimum.coverage10, gcv$raw.resultRB, gcv$raw.result10.
>
> ### when we want to simulate pseudo datasets based on different values of A and
> ### 2nd level mean as true ones, not using estimated values as true ones.
> coverage(g, A.or.r = 9, mean.PriorDist = 5, nsim = 10)
$coverageRB
[1] 0.998 0.999 1.000 0.999 0.999 0.995 0.991 0.999
$coverageS
[1] 1.0 1.0 1.0 1.0 1.0 0.9 1.0 1.0
$average.coverageRB
[1] 0.998
$average.coverageS
[1] 0.988
$overall.coverageRB
[1] 0.998
$se.overal.coverageRB
[1] 8e-04
$se.coverageRB
[1] 0.0015 0.0008 0.0001 0.0004 0.0002 0.0024 0.0055 0.0006
$se.coverageS
[1] 0.0 0.0 0.0 0.0 0.0 0.1 0.0 0.0
$raw.resultRB
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] 0.9847640 0.9999998 1.0000000 0.9997908 1.0000000 0.9993732 0.9999259
[2,] 0.9997135 0.9999210 0.9999861 0.9919308 0.9999628 0.9999285 0.9997725
[3,] 1.0000000 1.0000000 0.9999972 1.0000000 0.9999992 0.9998540 0.9999694
[4,] 0.9999383 0.9988361 0.9999987 0.9999946 1.0000000 0.9997252 0.9992767
[5,] 0.9999885 0.9999871 0.9999962 0.9974014 0.9995615 0.9990831 0.9989867
[6,] 0.9999989 0.9999988 0.9824083 0.9999999 0.9807068 0.9984266 0.9998216
[7,] 0.9999048 0.9990768 0.9950878 0.9989887 0.9820918 0.9999989 0.9990731
[8,] 1.0000000 0.9950041 1.0000000 0.9999996 1.0000000 0.9999995 1.0000000
[,8] [,9] [,10]
[1,] 1.0000000 0.9991794 1.0000000
[2,] 0.9997985 0.9998286 0.9962021
[3,] 0.9999999 0.9999223 0.9984980
[4,] 0.9959410 0.9999973 0.9990709
[5,] 0.9992161 0.9998059 0.9997078
[6,] 0.9928254 0.9999976 0.9999999
[7,] 0.9443528 0.9959980 0.9996104
[8,] 0.9952561 0.9999686 1.0000000
$raw.resultS
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 1 1 1 1 1 1 1 1 1 1
[2,] 1 1 1 1 1 1 1 1 1 1
[3,] 1 1 1 1 1 1 1 1 1 1
[4,] 1 1 1 1 1 1 1 1 1 1
[5,] 1 1 1 1 1 1 1 1 1 1
[6,] 1 1 1 1 1 1 1 0 1 1
[7,] 1 1 1 1 1 1 1 1 1 1
[8,] 1 1 1 1 1 1 1 1 1 1
$confidence.lvl
[1] 0.95
$effective.n
[1] 10
$model
[1] "gr"
$case
[1] 5
$betas
[1] NA
$A.r
[1] 9
$priormeanused
[1] 5
>
> ###############################################################
> # Binomial Regression Interactive Multilevel Modeling (BRIMM) #
> ###############################################################
>
> ####################################################################################
> # If we do not have any covariate and do not know a mean of the prior distribution #
> ####################################################################################
>
> b <- gbp(z, n, model = "binomial")
> b
Summary for each group (sorted by the ascending order of n):
obs.mean n prior.mean shrinkage low.intv post.mean upp.intv post.sd
1 0.400 45 0.267 0.622 0.222 0.317 0.421 0.0507
2 0.378 45 0.267 0.622 0.218 0.309 0.408 0.0487
3 0.356 45 0.267 0.622 0.213 0.300 0.396 0.0469
4 0.333 45 0.267 0.622 0.207 0.292 0.385 0.0453
5 0.311 45 0.267 0.622 0.202 0.284 0.374 0.0440
6 0.311 45 0.267 0.622 0.202 0.284 0.374 0.0440
7 0.289 45 0.267 0.622 0.195 0.275 0.363 0.0431
8 0.267 45 0.267 0.622 0.188 0.267 0.354 0.0424
9 0.244 45 0.267 0.622 0.180 0.258 0.345 0.0421
10 0.244 45 0.267 0.622 0.180 0.258 0.345 0.0421
11 0.222 45 0.267 0.622 0.172 0.250 0.337 0.0422
12 0.222 45 0.267 0.622 0.172 0.250 0.337 0.0422
13 0.222 45 0.267 0.622 0.172 0.250 0.337 0.0422
14 0.222 45 0.267 0.622 0.172 0.250 0.337 0.0422
15 0.222 45 0.267 0.622 0.172 0.250 0.337 0.0422
16 0.200 45 0.267 0.622 0.163 0.242 0.330 0.0426
17 0.178 45 0.267 0.622 0.154 0.233 0.323 0.0433
18 0.156 45 0.267 0.622 0.144 0.225 0.318 0.0444
Mean 45 0.267 0.622 0.185 0.266 0.357 0.0439
> print(b, sort = FALSE)
Summary for each group:
obs.mean n prior.mean shrinkage low.intv post.mean upp.intv post.sd
1 0.400 45 0.267 0.622 0.222 0.317 0.421 0.0507
2 0.378 45 0.267 0.622 0.218 0.309 0.408 0.0487
3 0.356 45 0.267 0.622 0.213 0.300 0.396 0.0469
4 0.333 45 0.267 0.622 0.207 0.292 0.385 0.0453
5 0.311 45 0.267 0.622 0.202 0.284 0.374 0.0440
6 0.311 45 0.267 0.622 0.202 0.284 0.374 0.0440
7 0.289 45 0.267 0.622 0.195 0.275 0.363 0.0431
8 0.267 45 0.267 0.622 0.188 0.267 0.354 0.0424
9 0.244 45 0.267 0.622 0.180 0.258 0.345 0.0421
10 0.244 45 0.267 0.622 0.180 0.258 0.345 0.0421
11 0.222 45 0.267 0.622 0.172 0.250 0.337 0.0422
12 0.222 45 0.267 0.622 0.172 0.250 0.337 0.0422
13 0.222 45 0.267 0.622 0.172 0.250 0.337 0.0422
14 0.222 45 0.267 0.622 0.172 0.250 0.337 0.0422
15 0.222 45 0.267 0.622 0.172 0.250 0.337 0.0422
16 0.200 45 0.267 0.622 0.163 0.242 0.330 0.0426
17 0.178 45 0.267 0.622 0.154 0.233 0.323 0.0433
18 0.156 45 0.267 0.622 0.144 0.225 0.318 0.0444
Mean 45 0.267 0.622 0.185 0.266 0.357 0.0439
> summary(b)
Main summary:
obs.mean n prior.mean shrinkage low.intv
Group with min(obs.mean) 0.156 45 0.267 0.622 0.144
Group with median(obs.mean)1 0.244 45 0.267 0.622 0.180
Group with median(obs.mean)2 0.244 45 0.267 0.622 0.180
Group with max(obs.mean) 0.400 45 0.267 0.622 0.222
Overall Mean 45 0.267 0.622 0.185
post.mean upp.intv post.sd
Group with min(obs.mean) 0.225 0.318 0.0444
Group with median(obs.mean)1 0.258 0.345 0.0421
Group with median(obs.mean)2 0.258 0.345 0.0421
Group with max(obs.mean) 0.317 0.421 0.0507
Overall Mean 0.266 0.357 0.0439
Estimation summary for the second-level variance component:
alpha = log(A) for Gaussian or alpha = log(1/r) for Binomial and Poisson data:
post.mode.alpha post.sd.alpha post.mode.r
-4.31 0.82 74.1
Estimation summary for the regression coefficient :
estimate se z.val p.val
beta1 -1.012 0.1 -10.153 0
> plot(b)
> plot(b, sort = FALSE)
>
> ### when we want to simulate pseudo datasets considering the estimated values
> ### as true ones.
> bcv <- coverage(b, nsim = 10)
>
> ### bcv$coverageRB, bcv$coverage10, bcv$average.coverageRB, bcv$average.coverage10,
> ### bcv$minimum.coverageRB, bcv$minimum.coverage10, bcv$raw.resultRB, bcv$raw.result10.
>
> ### when we want to simulate pseudo datasets based on different values of r
> ### and a regression coefficient (intercept),
> ### not using estimated values as true ones.
> bcv <- coverage(b, A.or.r = 60, reg.coef = -1, nsim = 10)
>
> ##################################################################################
> # If we have one covariate and do not know a mean of the prior distribution yet, #
> ##################################################################################
>
> b <- gbp(z, n, x1, model = "binomial")
> b
Summary for each group (sorted by the ascending order of n):
obs.mean n X1 prior.mean shrinkage low.intv post.mean upp.intv post.sd
1 0.400 45 1.00 0.310 0.715 0.248 0.335 0.429 0.0462
2 0.378 45 1.00 0.310 0.715 0.244 0.329 0.420 0.0448
3 0.356 45 1.00 0.310 0.715 0.240 0.323 0.411 0.0437
4 0.333 45 1.00 0.310 0.715 0.236 0.316 0.403 0.0429
5 0.311 45 1.00 0.310 0.715 0.230 0.310 0.396 0.0424
6 0.311 45 0.00 0.233 0.715 0.179 0.256 0.341 0.0415
7 0.289 45 0.00 0.233 0.715 0.175 0.249 0.331 0.0400
8 0.267 45 0.00 0.233 0.715 0.171 0.243 0.323 0.0388
9 0.244 45 0.00 0.233 0.715 0.166 0.237 0.315 0.0380
10 0.244 45 1.00 0.310 0.715 0.210 0.291 0.379 0.0432
11 0.222 45 0.00 0.233 0.715 0.161 0.230 0.308 0.0377
12 0.222 45 0.00 0.233 0.715 0.161 0.230 0.308 0.0377
13 0.222 45 0.00 0.233 0.715 0.161 0.230 0.308 0.0377
14 0.222 45 1.00 0.310 0.715 0.202 0.285 0.375 0.0441
15 0.222 45 1.00 0.310 0.715 0.202 0.285 0.375 0.0441
16 0.200 45 0.00 0.233 0.715 0.155 0.224 0.302 0.0377
17 0.178 45 0.00 0.233 0.715 0.148 0.218 0.297 0.0381
18 0.156 45 0.00 0.233 0.715 0.140 0.211 0.292 0.0389
Mean 45 0.44 0.267 0.715 0.191 0.267 0.351 0.0410
> print(b, sort = FALSE)
Summary for each group:
obs.mean n X1 prior.mean shrinkage low.intv post.mean upp.intv post.sd
1 0.400 45 1.00 0.310 0.715 0.248 0.335 0.429 0.0462
2 0.378 45 1.00 0.310 0.715 0.244 0.329 0.420 0.0448
3 0.356 45 1.00 0.310 0.715 0.240 0.323 0.411 0.0437
4 0.333 45 1.00 0.310 0.715 0.236 0.316 0.403 0.0429
5 0.311 45 1.00 0.310 0.715 0.230 0.310 0.396 0.0424
6 0.311 45 0.00 0.233 0.715 0.179 0.256 0.341 0.0415
7 0.289 45 0.00 0.233 0.715 0.175 0.249 0.331 0.0400
8 0.267 45 0.00 0.233 0.715 0.171 0.243 0.323 0.0388
9 0.244 45 0.00 0.233 0.715 0.166 0.237 0.315 0.0380
10 0.244 45 1.00 0.310 0.715 0.210 0.291 0.379 0.0432
11 0.222 45 0.00 0.233 0.715 0.161 0.230 0.308 0.0377
12 0.222 45 0.00 0.233 0.715 0.161 0.230 0.308 0.0377
13 0.222 45 0.00 0.233 0.715 0.161 0.230 0.308 0.0377
14 0.222 45 1.00 0.310 0.715 0.202 0.285 0.375 0.0441
15 0.222 45 1.00 0.310 0.715 0.202 0.285 0.375 0.0441
16 0.200 45 0.00 0.233 0.715 0.155 0.224 0.302 0.0377
17 0.178 45 0.00 0.233 0.715 0.148 0.218 0.297 0.0381
18 0.156 45 0.00 0.233 0.715 0.140 0.211 0.292 0.0389
Mean 45 0.44 0.267 0.715 0.191 0.267 0.351 0.0410
> summary(b)
Main summary:
obs.mean n X1 prior.mean shrinkage low.intv
Group with min(obs.mean) 0.156 45 0.000 0.233 0.715 0.140
Group with median(obs.mean)1 0.244 45 0.000 0.233 0.715 0.166
Group with median(obs.mean)2 0.244 45 1.000 0.310 0.715 0.210
Group with max(obs.mean) 0.400 45 1.000 0.310 0.715 0.248
Overall Mean 45 0.444 0.267 0.715 0.191
post.mean upp.intv post.sd
Group with min(obs.mean) 0.211 0.292 0.0389
Group with median(obs.mean)1 0.237 0.315 0.0380
Group with median(obs.mean)2 0.291 0.379 0.0432
Group with max(obs.mean) 0.335 0.429 0.0462
Overall Mean 0.267 0.351 0.0410
Estimation summary for the second-level variance component:
alpha = log(A) for Gaussian or alpha = log(1/r) for Binomial and Poisson data:
post.mode.alpha post.sd.alpha post.mode.r
-4.73 0.957 113
Estimation summary for the regression coefficient :
estimate se z.val p.val
beta1 -1.194 0.131 -9.129 0.000
beta2 0.389 0.187 2.074 0.038
> plot(b)
> plot(b, sort = FALSE)
>
> ### when we want to simulate pseudo datasets considering the estimated values
> ### as true ones.
> bcv <- coverage(b, nsim = 10)
>
> ### bcv$coverageRB, bcv$coverage10, bcv$average.coverageRB, bcv$average.coverage10,
> ### bcv$minimum.coverageRB, bcv$minimum.coverage10, bcv$raw.resultRB, bcv$raw.result10.
>
> ### when we want to simulate pseudo datasets based on different values of r
> ### and regression coefficients, not using estimated values
> ### as true ones. Two values of reg.coef are for beta0 and beta1
> bcv <- coverage(b, A.or.r = 60, reg.coef = c(-1, 0), nsim = 10)
>
> ################################################
> # If we know a mean of the prior distribution, #
> ################################################
>
> b <- gbp(z, n, mean.PriorDist = 0.265, model = "binomial")
> b
Summary for each group (sorted by the ascending order of n):
obs.mean n prior.mean shrinkage low.intv post.mean upp.intv post.sd
1 0.400 45 0.265 0.65 0.230 0.312 0.401 0.0439
2 0.378 45 0.265 0.65 0.224 0.304 0.391 0.0428
3 0.356 45 0.265 0.65 0.218 0.297 0.382 0.0418
4 0.333 45 0.265 0.65 0.212 0.289 0.372 0.0409
5 0.311 45 0.265 0.65 0.206 0.281 0.363 0.0401
6 0.311 45 0.265 0.65 0.206 0.281 0.363 0.0401
7 0.289 45 0.265 0.65 0.200 0.273 0.354 0.0395
8 0.267 45 0.265 0.65 0.193 0.266 0.345 0.0390
9 0.244 45 0.265 0.65 0.186 0.258 0.337 0.0386
10 0.244 45 0.265 0.65 0.186 0.258 0.337 0.0386
11 0.222 45 0.265 0.65 0.179 0.250 0.329 0.0383
12 0.222 45 0.265 0.65 0.179 0.250 0.329 0.0383
13 0.222 45 0.265 0.65 0.179 0.250 0.329 0.0383
14 0.222 45 0.265 0.65 0.179 0.250 0.329 0.0383
15 0.222 45 0.265 0.65 0.179 0.250 0.329 0.0383
16 0.200 45 0.265 0.65 0.172 0.242 0.321 0.0381
17 0.178 45 0.265 0.65 0.164 0.234 0.313 0.0381
18 0.156 45 0.265 0.65 0.156 0.227 0.306 0.0383
Mean 45 0.265 0.65 0.191 0.265 0.346 0.0395
> print(b, sort = FALSE)
Summary for each group:
obs.mean n prior.mean shrinkage low.intv post.mean upp.intv post.sd
1 0.400 45 0.265 0.65 0.230 0.312 0.401 0.0439
2 0.378 45 0.265 0.65 0.224 0.304 0.391 0.0428
3 0.356 45 0.265 0.65 0.218 0.297 0.382 0.0418
4 0.333 45 0.265 0.65 0.212 0.289 0.372 0.0409
5 0.311 45 0.265 0.65 0.206 0.281 0.363 0.0401
6 0.311 45 0.265 0.65 0.206 0.281 0.363 0.0401
7 0.289 45 0.265 0.65 0.200 0.273 0.354 0.0395
8 0.267 45 0.265 0.65 0.193 0.266 0.345 0.0390
9 0.244 45 0.265 0.65 0.186 0.258 0.337 0.0386
10 0.244 45 0.265 0.65 0.186 0.258 0.337 0.0386
11 0.222 45 0.265 0.65 0.179 0.250 0.329 0.0383
12 0.222 45 0.265 0.65 0.179 0.250 0.329 0.0383
13 0.222 45 0.265 0.65 0.179 0.250 0.329 0.0383
14 0.222 45 0.265 0.65 0.179 0.250 0.329 0.0383
15 0.222 45 0.265 0.65 0.179 0.250 0.329 0.0383
16 0.200 45 0.265 0.65 0.172 0.242 0.321 0.0381
17 0.178 45 0.265 0.65 0.164 0.234 0.313 0.0381
18 0.156 45 0.265 0.65 0.156 0.227 0.306 0.0383
Mean 45 0.265 0.65 0.191 0.265 0.346 0.0395
> summary(b)
Main summary:
obs.mean n prior.mean shrinkage low.intv
Group with min(obs.mean) 0.156 45 0.265 0.65 0.156
Group with median(obs.mean)1 0.244 45 0.265 0.65 0.186
Group with median(obs.mean)2 0.244 45 0.265 0.65 0.186
Group with max(obs.mean) 0.400 45 0.265 0.65 0.230
Overall Mean 45 0.265 0.65 0.191
post.mean upp.intv post.sd
Group with min(obs.mean) 0.227 0.306 0.0383
Group with median(obs.mean)1 0.258 0.337 0.0386
Group with median(obs.mean)2 0.258 0.337 0.0386
Group with max(obs.mean) 0.312 0.401 0.0439
Overall Mean 0.265 0.346 0.0395
Estimation su