lmeNB (version 1.2)

lmeNB: Performs the maximum likelihood estimation for the negative binomial mixed-effect model. This function is a wrapper for mle.fun, mle.ar1.fun, mle.a3.fun and mle.ar1.non3.

Description

The model assumes that given the random effect G_i=g_i, the count responses Y_ijs of subject i, (i = 1, ..., N), at time points j (=1,...,n_i) follow the negative binomial distribution: Y_ij | G_i =g_i ~ NB(r_ij,p_i), where r_{ij} = exp(X_ij^T beta)/a and beta is fixed effect coefficients. The failure probability p_i of subject i, is parametrized as: p_i= 1/(g_i*a+1), and a > 0. The model assumes E(G_i) = 1 so that E(Y_ij|G_i=g_i)=r_ij*g_i*a and E(Y_ij)= r_ij*g_i = exp(X_ij^T beta). This assumption allows the interpretation of the latent random variable G_i as the subject i's activity rate relative to the overall cohort. Regarding the dependence structures of Y_ij and Y_ij' conditional on the random effect, we consider two scenarios. [1]: Y_ij and Y_ij' are independent conditionally on Gi. This assumption leads us to cov(Y_ij,Y_ij'|G_i=g_i)=0 cov(Y_ij,Y_ij')=E(Y_ij)^2*Var(Y_ij)+E(Y_ij)*(1+(Var(Y_ij)+1)*a) [2]: AR(1) structures for Y_ij and Y_ij' conditionally on Gi. Given emph{ G_i=g_i}, Y_ij depends on Y_i(j-1) through the beta binomial thinning and is conditionally independent on Y_ij' given Y_i(j-1) for all j' < j-1. The beta binomial thinning operator depends on a parameter d which indicates the strength of the positive AR(1) association between repeated measures of a subject: the larger d, the stronger the positive correlations between the repeated measures of the same subject are. This interpretation depends on the result: cov(Y_ij,Y_ij'|G_i=g_i)=d^{j-j'} E(Y_ij')*(a*g_i^2+g_i). Regarding the random effect G_is distribution, lmNB allows three scenarios: (1) The log-normal distribution with E(G_i)=1 and Var(G_i)=theta (2) The gamma distribution with E(G_i)=1 and Var(G_i)=theta. (3) No distributional assumption and the random effect distribution is approximated by the estimated values of the quantity: gamma_i = w_i (y_i+/ mu_i+) + (1-w_i) , i=1,...,N, where y_i+ = sum_j=1^n_i y_ij, mu_i+ = sum_j=1^n_i mu_ij and, w_i = sqrt( Var(G_i)/Var(Y_i+/mu_i+) ). See Zhao et al. for more details.

Usage

lmeNB(formula,data,ID,p.ini=NULL,IPRT=FALSE,AR=FALSE,RE=c("G","N","semipara"),
      deps=0.001,Vcode,i.tol=1e-75,o.tol=1.e-3,maxit=100)

Arguments

formula

An object of class "formula" (or one that can be coerced to that class): a symbolic description of the model to be fitted. The formula must contain an intercept term.

data

A data frame, list or environment (or object coercible by as.data.frame to a data frame) containing the variables in the model. Each row must contain the data corresponding to the repeated measure j of a subject and the rows (i,j)

ID

A vector of length sum_i^N n_i, containing the patient IDs. i.e., c(rep(ID_1,n_1),rep(ID_2,n_2),...,rep(ID_N,n_N)). The length must be the same as the number of rows of data. Missing values are NOT accepted.

p.ini

The initial values of the parameters: If AR=0 and (RE="G" or RE="N"),p.ini=(log(a), log(var(G)), beta0, beta1, ...). If AR=0 and RE="semipara",

IPRT

A logical, passed to Iprt of function optim. If TRUE then print iterations.

AR

A logical, if TRUE, then the AR(1) structure is assumed among the responses.

RE

The distribution of random effects G_i. If model="G" then the random effects are assumed to be from the gamma distribution. If model="N" then they are assumed to be form the log-normal distribution.

deps

Passed to mle.a3.fun and mle.ar1.non3.

Vcode

A vector of length the total number of repeated measures, containing the indices of time point. For example, there are three subjects and the first two subjects do not have missing visits and completed five visits while the last subject missed the third

i.tol

A real number to determine the tolerance for integrate. Necessary only for semiparametric random effect model.

o.tol

A real number to determine the tolerance for optim. Necessary only for semiparametric random effect model.

maxit

The maximum number of iterations. Necessary only for semiparametric random effect model.

Value

• opt
  The values returned by optim.
• nlk
  The value of the negative log-likelihood corresponding to opt$par
• V
  The approximated asymptotic covariance matrix of the maximum likelihood estimators. V=solve(opt$hessian)
• est
  A (3 + # covariates) by 2 matrix. The first column contains the estimates of the model parameters, log(a), log(theta),beta0,beta1,...)The second column contains the approximated standard deviations of the estimators, i.e., sqrt(diag(V))
• mod
  If model="G" then mod="G". If model="N" then mod="N".
• idat
  A dataframe, containing ID, CEL, x.1, x.2, ...The column labeled as CEL contains the response counts.
• cor
  "ind", indicating that the model assumes independent structure of the count responses given the random effects.

Details

mle.fun calls optim to minimize the negative log-likelihood of the negative binomial model with respect to the model parameters: c(log(a), log(theta), beta0, beta1, ...).The Nelder-Mead algorithm is employed.The log-likelihood is obtained by marginalizing out the random effects.The numerical integration is carried out using adaptive quadrature.

The missing count responses, if assumed to be missing at random, can be ignored.Other types of missing data are currently not accepted.

When the estimated over-dispersion parameter (a) is close to zero, the negative binomial model reduces to the poisson model, suggesting that the negative binomial mixed-effect model might not be appropriate.When AR=1 and the estimated auto-correlation parameter (d) is close to zero, the model is suggesting that there is no AR(1) structure among the sequentially collected responses. Hence user might use AR=0 setting which assume no AR(1) structure.We note that the results could be sensitive to initial values.

References

Zhao, Y., Li, D.K.B., Petkau, J.A., Riddehough, A. & Traboulsee, A. Detection of unusual increases in MRI lesion counts in multiple sclerosis patients.

See Also

The subroutines of this function is: mle.fun, mle.ar1.fun, mle.a3.fun, mle.ar1.non3,

The subroutines of index.batch to compute the conditional probability index: jCP.ar1, CP1.ar1, MCCP.ar1, CP.ar1.se, CP.se, jCP,

The functions to generate simulated datasets: rNBME.R.

Examples

## See the examples in help files of rNBME.R.