R: Fit the semi-parametric negative binomial mixed-effect...
fitSemiIND
R Documentation
Fit the semi-parametric negative binomial mixed-effect independent model.
Description
This function fits the semi-parametric negative binomial mixed-effect independent model to repeated count responses (Zhao et al. 2013).
The conditional distribution of response count given random effect is modelled by Negative Binomial as described in description of lmeNB.
The conditional dependence among the response counts of a subject is assumed independent.
The semiparametric procedure is employed for random effects.
See descriptions of lmeNB.
A vector of length 3 + # covariates, containing the initial values for
the parameters (log(alpha), log(Var(G[i])), beta[0], beta[1],...).
NULL is accepted.
IPRT
See lmeNB.
deps
See lmeNB.
maxit
See lmeNB.
u.low
See lmeNB.
Details
The algorithm repeats the following four steps until a stoping criterion is satisfied:
Step 1) Given α, Estimate the coefficients of covariates by the method of generalized Least Squares.
That is, this step solves for:
argmin_{β} ∑[i=1]^N (Y[i]-E(Y[i];β))^T W[i] (Y[i]-E(Y[i];β))
where the weight matrix for each patient W[i] is selected to Var(Y[i])^{-1}
(which is a function of α)
if it exists, else it is set to be an identity matrix.
Step 2) Approximate the distribution of the random effect G[i] by γ.
Step 3) Estimate α by minimizing the negative psudo-profile likelihood.
The numerical minimization is carried out using optimize and the numerical integration is carried out using adaptive quadrature.
Step 4) Estimate Var(G[i]) by the medhod of moment and update the weights.
All the computations are done in R.
Value
opt
See lmeNB.
diffPara
The largest absolute difference of parameter vectors between the current and previous iterations.
V
NULL
est
See lmeNB.
gtb
The relative frequency table of G[i], (i=1,...,N).
gh1 (the second column) contains the unique values in ascending order and ghw (the first column) contains the corresponding relative frequencies.
counter
The number of iterations before the algorithm was terminated
gi
A vector of length N, containing the approximated random effect G[i],i=1,...,N.
RE
"NoN", denoting that the fitted model is a semi-parametric mixed-effect model.
AR
FALSE
paraAll
Record estimated parameters at every iteration.
Author(s)
Zhao, Y. and Kondo, Y.
References
Detection of unusual increases in MRI lesion counts in individual multiple sclerosis patients. (2013) Zhao, Y., Li, D.K.B., Petkau, A.J., Riddehough, A., Traboulsee, A., Journal of the American Statistical Association.
See Also
The main function to fit the Negative Binomial mixed-effect model:
lmeNB,
The functions to fit the other models:
fitParaIND,
fitParaAR1,
fitSemiAR1,
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
## Not run:
## generate a simulated dataset from the negative binomial
## mixed-effect independent model:
## Y_ij | G_i = g_i ~ NB(r_ij,p_i) where r_ij = exp(X^T beta)/a , p_i =1/(a*g_i+1)
## with G_i is from unknown distribution
## For the simulation purpose, G_i's are from the mixture of
## the gamma and the log-normal distributions.
sn <- 5 ## the number of repeated measures of each subject
n <- 80 ## the number of subjects
logtheta <- 1.3
th <- exp(logtheta) ## the parameter for the gamma distribution of g_i
loga <- -0.5
## the parameter for the failure probability of the negative binomial distribution
a <- exp(loga)
b0 <- 0.5
u1 <- rep(exp(b0),sn) ## the mean vector for group 1 at time point 1,...,sn
u2 <- rep(exp(b0),sn) ## the mean vector for group 2 at time point 1,...,sn
DT3 <- rNBME.R(gdist="GN", n=n, a=a, th=th, u1=u1, u2=u2, sn=sn,
othrp=list(p.mx=0.1,u.n=3,s.n=1,sh.mx = NA) ## 0 < p.mx < 1
)
ID <- DT3$id
dt3 <- data.frame(CEL=DT3$y)
Vcode <- rep(-1:(sn-2),n) # scan number -1, 0, 1, 2, 3
new <-Vcode>0 # new scans: 1,2,3
## 1) Fit the negative binomial mixed-effect AR(1) model
## where random effects is from the gamma distribution
logitd <- -0.2
re.gamma.ar1 <- fitParaAR1(formula=CEL~1,data=dt3,ID=ID,
Vcode=Vcode,
p.ini=c(loga,logtheta,logitd,b0),
## log(a), log(theta), logit(d), b0
RE="G",
IPRT=TRUE)
## compute the estimates of the conditional probabilities
## with sum of the new repeated measure as a summary statistics
## i.se=FALSE,C=TRUE options for speed up!
Psum <- index.batch(olmeNB=re.gamma.ar1,data=dt3,ID=ID, Vcode=Vcode,
labelnp=new,qfun="sum", IPRT=TRUE,i.se=FALSE,C=TRUE,i.tol=1E-3)
## 2) Fit the negative binomial mixed-effect AR(1) model
## where random effects is from the log-normal distribution
re.logn.ar1 <- fitParaAR1(formula=CEL~1,data=dt3,ID=ID,
Vcode=Vcode, RE="N", IPRT=TRUE)
## REQUIRES SOME TIME..
Psum <- index.batch(olmeNB=re.logn.ar1, data=dt3,ID=ID,Vcode=Vcode,
labelnp=new,qfun="sum", IPRT=TRUE,i.se=FALSE,C=TRUE,i.tol=1E-3)
## 3) Fit the negative binomial independent model
## where random effects is from the gamma distribution
re.gamma.ind <- fitParaIND(formula=CEL~1,data=dt3,ID=ID,
RE="G",IPRT=TRUE)
Psum <- index.batch(olmeNB=re.gamma.ind, data=dt3,ID=ID,Vcode=Vcode,
labelnp=new,qfun="sum", IPRT=TRUE,i.se=TRUE)
## 4) Fit the negative binomial independent model
## where random effects is from the lognormal distribution
re.logn.ind <- fitParaIND(formula=CEL~1,data=dt3,ID=ID,
RE="N",
p.ini=c(loga,logtheta,b0),
IPRT=TRUE)
Psum <- index.batch(olmeNB=re.logn.ind,data=dt3,ID=ID,labelnp=new,qfun="sum", IPRT=TRUE)
## 5) Fit the semi-parametric negative binomial AR(1) model
logvarG <- -0.4
re.semi.ar1 <- fitSemiAR1(formula=CEL~1,data=dt3,ID=ID,Vcode=Vcode)
Psum <- index.batch(olmeNB=re.semi.ar1,data=dt3,ID=ID,Vcode=Vcode,
labelnp=new,qfun="sum", IPRT=TRUE,MC=TRUE,i.se=FALSE)
## 6) Fit the semi-parametric negative binomial independent model
## This is closest to the true model
re.semi.ind <- fitSemiIND(formula=CEL~1,data=dt3,ID=ID, p.ini=c(loga, logvarG, b0))
## compute the estimates of the conditional probabilities
## with sum of the new repeated measure as a summary statistics
Psum <- index.batch(olmeNB=re.semi.ind,data=dt3,ID=ID, labelnp=new,
qfun="sum", IPRT=TRUE,i.se=FALSE)
## compute the estimates of the conditional probabilities
## with max of the new repeated measure as a summary statistics
Pmax <- index.batch(olmeNB=re.semi.ind, data=dt3,ID=ID,labelnp=new, qfun="max",
IPRT=TRUE,i.se=FALSE)
## Which patient's estimated probabilities based on the sum and max
## statistics disagrees the most?
( IDBigDif <- which(rank(abs(Pmax$condProbSummary[,1]-Psum$condProbSummary[,1]))==80) )
## Show the patient's CEL counts
dt3$CEL[ID==IDBigDif]
## Show the estimated conditional probabilities based on the sum summary statistics
Psum$condProbSummary[IDBigDif,]
## Show the estimated conditional probabilities based on the max summary statistics
Pmax$condProbSummary[IDBigDif,]
## ======================================================================== ##
## == Which model performed the best in terms of the estimation of beta0 == ##
## ======================================================================== ##
getpoints <- function(y,estb0,sdb0=NULL,crit=qnorm(0.975))
{
points(estb0,y,col="blue",pch=16)
if (!is.null(sdb0))
{
points(c(estb0-crit*sdb0,estb0+crit*sdb0),rep(y,2),col="red",type="l")
}
}
ordermethod <- c("gamma.ar1","logn.ar1","gamma.ind","logn.ind","semi.ar1","semi.ind")
estb0s <- c(
re.gamma.ar1$est[4,1],
re.logn.ar1$est[4,1],
re.gamma.ind$est[3,1],
re.logn.ind$est[3,1],
re.semi.ar1$est[4],
re.semi.ind$est[3]
)
## The true beta0 is:
b0
c <- 1.1
plot(0,0,type="n",xlim=c(min(estb0s)-0.5,max(estb0s)*c),ylim=c(0,7),yaxt="n",
main="Simulated from the independent model \n with random effect ~ mixture of normal and gamma")
legend("topright",
legend=ordermethod)
abline(v=b0,lty=3)
## 1) gamma.ar1
sdb0 <- re.gamma.ar1$est[4,2]
getpoints(6,estb0s[1],sdb0)
## 2)logn.ar1
sdb0 <- re.logn.ar1$est[4,2]
getpoints(5,estb0s[2],sdb0)
## 3) gamma.ind
sdb0 <- re.gamma.ind$est[3,2]
getpoints(4,estb0s[3],sdb0)
## 4) logn.ind
sdb0 <- re.logn.ind$est[3,2]
getpoints(3,estb0s[4],sdb0)
## 5) semi.ar1
getpoints(2,estb0s[5])
## 6) semi.ind
getpoints(1,estb0s[6])
## End(Not run)