BBmm | R Documentation |
BBmm
function performs beta-binomial mixed-effects models, i.e., it allows the inclusion of gaussian random effects in the linear predictor of a logistic beta-binomial regression model.
The structure of the random part of the model can be expecified by two different ways: (i) determining the random.formula
argument, or (ii) especifying the model matrix of the random effects, Z
, and determining the number of random effects in each random component, nRandComp
.
For the estimation of the parameters the joint log-likelihood can be optimized by means of two different approaches: (i) BB-Delta, specific Delta algorithm developed for beta-binomial mixed-effect regression models, and (ii) usual Newton-Raphson method.
BBmm(fixed.formula,X,y,random.formula,Z=NULL,nRandComp=NULL,m,data,
method="BB-Delta",maxiter=50,show=FALSE,nDim=1,tolerance=10^(-6))
fixed.formula |
an object of class |
X |
a matrix class object containing the covariate structure of the fixed part of the model to be fitted. If the |
y |
a vector containing the outcome variable(s). If joint analysis is expected, the outcome variables must be listed one after another in a vector. |
random.formula |
an object of class |
Z |
the design matrix of the random effects. If the |
nRandComp |
the number of random effects in each random component of the model. It must be especified as a vector, where the 'i'th value must correspond with the number of random effects of the 'i'th random component. It must be only determined when we specify the random structure of the model by the model matrix of the random effects, |
m |
maximum score number in each beta-binomial observation. |
data |
an optional data frame, list or environment (or object coercible by |
method |
the methodology for performing the estimation of the parameters. Options "BB-Delta" or "NR". Default "BB-Delta". |
maxiter |
the maximum number of iterations in the parameters estimation algorithm. Default 50. |
show |
logical parameter. If |
nDim |
number of response outcome variables involved in the multidimensional analysis. |
tolerance |
tolerance of the estimated linear predictor in the estimation process. |
BBmm
function performs beta-binomial mixed-effects regression models. It extends the beta-binomial regression model including gaussian random effects in the linear predictor. The model is defined as, conditional on some gaussian random effects u
, the response variable y
follows a beta-binomial distribution with parameters m
, p
and \phi
,
y|u \sim BB(m,p,\phi), \text{and } u \sim N(0,D)
where the regression model is defined as,
\log(p/(1-p))=X*\beta+Z*u
and D
is determined by some dispersion parameters icluded in the parameter vector \theta
.
The estimation of the regression parameters \beta
and the prediction of the random effects u
is done via the extended likelihood, where the marginal likelihood is approximated to the h-likelihood by a first order Laplace approximation,
h=f(y|\beta,u,\theta)+f(u|\theta).
The previous formula do not have a closed form and numerical methods are needed for the estimation precedure. Two approches are available in the function in order to perform the fixed and random effects estimation: (i) A special case of the Delta algorithm developed for beta-binomial mixed-effects model estimations, and (ii) the general Newton-Raphson algorithm available in R-package rootSolve
.
On the other hand, the estimation of dispersion parameters by the h-likelihood can be bias due to the previous estimation of the regression parameters. Consequenlty, a penalization of the h-likelihood must be performed to get unbiased estimations of the dispersion parameters. Lee and Nelder (1996) proposed the adjusted profile h-likelihood, where the following penalization is applied,
h(\theta)=h+(1/2)*\log[det(2\pi H^{-1})],
being H
the Hessian matrix of the model, replacing the terms \beta
and u
involved in the previous formula with their estimates.
The method iterates between both estimation processes until convergence is reached based on the pre-established tolerance of the linear predictor.
BBmm
returns an object of class "BBmm
".
The function summary
(i.e., summary.BBmm
) can be used to obtain or print a summary of the results..
fixed.coef |
estimated value of the fixed coefficients of the regression. |
fixed.vcov |
variance and covariance matrix of the estimated fixed coefficients of the regression. |
random.coef |
predicted random effects of the regression. |
sigma.coef |
estimated value of the random effects variance parameters. |
sigma.var |
variance of the estimated value of the random effects variance parameters. |
phi.coef |
estimated value of the dispersion parameter of the conditional beta-binomial distribution. |
psi.coef |
estimated value of the logrithm of the dispersion parameter of the conditional beta-binomial distribution. |
psi.var |
variance of the estimation of the logarithm of the conditional beta-binomial distribution. |
fitted.values |
the fitted mean values of the probability parameter of the conditional beta-binomial distribution. |
conv |
convergence of the methodology. If the method has converged it returns "yes", otherwise "no". |
deviance |
deviance of the model. |
df |
degrees of freedom of the model. |
null.deviance |
null-deviance, deviance of the null model. The null model will only include an intercept as the estimation of the probability parameter. |
null.df |
degrees of freedom of the null model. |
nRand |
number of random effects. |
nComp |
number of random components. |
nRandComp |
number of random effects in each random component of the model. |
namesRand |
names of the random components. |
iter |
number of iterations in the estimation method. |
nObs |
number of observations in the data. |
y |
dependent response variable in the model. |
X |
model matrix of the fixed effects. |
Z |
model matrix of the random effects. |
D |
variance and covariance matrix of the random effects. |
balanced |
if the response dependent variable is balanced it returns "yes", otherwise "no". |
m |
maximum score number in each beta-binomial observation. |
call |
the matched call. |
formula |
the formula supplied. |
J. Najera-Zuloaga
D.-J. Lee
I. Arostegui
Najera-Zuloaga, J.; Lee, D.-J.; Esteban, C. & Arostegui, I. (2023): Multidimensional beta-binomial regression model: a joint analysis of patient-reported outcomes. Statistical Modelling, 0(0).
Najera-Zuloaga J., Lee D.-J. & Arostegui, I. (2019): A beta-binomial mixed effects model approach for analysing longitudinal discrete and bounded outcomes. Biometrical Journal, 61(3), 600-615.
Najera-Zuloaga J., Lee D.-J. & Arostegui I. (2018): Comparison of beta-binomial regression model approaches to analyze health related quality of life data. Statistical Methods in Medical Research, 27(10), 2989-3009.
Lee Y. & Nelder J. A. (1996): Hierarchical generalized linear models, Journal of the Royal Statistical Society. Series B, 58, 619-678
The multiroot
and uniroot
functions of the R-package rootSolve
for the general Newton-Raphson algorithm.
set.seed(7)
# Defining the parameters
k <- 100
m <- 10
phi <- 0.5
beta <- c(1.5,-1.1)
sigma <- 0.5
# Simulating the covariate and random effects
x <- runif(k,0,10)
X <- model.matrix(~x)
z <- as.factor(rBI(k,4,0.5,2))
Z <- model.matrix(~z-1)
u <- rnorm(5,0,sigma)
# The linear predictor and simulated response variable
eta <- beta[1]+beta[2]*x+crossprod(t(Z),u)
p <- 1/(1+exp(-eta))
y <- rBB(k,m,p,phi)
dat <- data.frame(cbind(y,x,z))
dat$z <- as.factor(dat$z)
# Apply the model
model <- BBmm(fixed.formula = y~x,random.formula = ~z,m=m,data=dat)
model
#---------#
# Multidimensional regression with 2 dimensions
set.seed(5)
nId <- 25
# Simulation
m <- 10
beta1 <- c(1,-0.5)
beta2 <- c(-1,0.5)
beta <- cbind(beta1,beta2)
x1 <- rnorm(nId, 1,2)
X1 <- model.matrix(~x1)
x2 <- rnorm(nId, -1,1)
X2 <- model.matrix(~x2)
sigma <- 0.6
u <- rnorm(nId,0,sigma)
eta1 <- beta1[1]+x1*beta1[2]+u
p1 <- 1/(1+exp(-eta1))
eta2 <- beta2[1]+x2*beta2[2]+u
p2 <- 1/(1+exp(-eta2))
phi1 <- 0.3
phi2 <- 1
phi <- c(phi1,phi2)
y1 <- rBB(nId,m,p1,phi1)
y2 <- rBB(nId,m,p2,phi2)
y <- c(y1,y2)
# Define matrices
X <- Matrix::bdiag(X1,X2)
X <- as.matrix(X)
Z. <- diag(nId)
Z <- kronecker(rbind(1,1),Z.)
# Fit the model
Model.multi <- BBmm(X=X,y=y,Z=Z,nRandComp = nId,m=m,nDim=2)
Model.multi
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.