BBmm function performs beta-binomial mixed-effects models, i.e., it allows the inclusion of gaussian random effects in the linear predictor of a logitic 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. The estimation of the parameters can also be done by means of two approaches: (i) BB-NR, special Newton-Raphson algorithm developed for beta-binomial mixed-effect models, and (ii) using the rootSolve R-package.BBmm(fixed.formula,random.formula,Z=NULL,nRandComp=NULL,m,data,
method="BB-NR",maxiter=100)"formula" (or one that can be coerced to that class): a symbolic description of the fixed part of the model to be fitted."formula" (or one that can be coerced to that class): a symbolic description of the random part of the model to be fitted. It must be specified in cases where the model matrix of the random effects Z is not determined.random.formula argument is specified this argument should not be specified, as we will be specifying twice the random structure of the model.Z.as.data.frame to a data frame) containing the variables in the model. If not found in data, the variables are taken from environment(formula).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.. BBmm function performs beta-binomial mixed effects models. It extends the beta-binomial logistic regression to the inclusion of random effects in the linear predictor of the model. The model is defined as, conditional on some gaussian random effects \(u\) the response variable \(y\) follows a beta-binomial distribution of parameters \(m\), \(p\) and \(phi\),
$$y|u \sim BB(m,p,phi), u \sim N(0,D)$$
where
$$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 paramters \(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 a Newton-Raphson 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 estimation of the regression parameters. Consequenlty, a penalization of the h-likelihood must be performed to get an unbiased profile h-likelihood 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})]$$
where \(H\) is the Hessian matrix of the model, and the terms \(beta\) and \(u\) involved in the previous formula are replaced by their estimates. The method iterates between both estimation processes until convergence is reached.multiroot and uniroot functions of the R-package rootSolve for the general Newton-Raphson algorithm.set.seed(14)
# 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
Run the code above in your browser using DataLab