rNBME.R: Simulate a dataset from the negative binomial mixed-effect independent/AR(1) model

Description

This function simulates a dataset based on the negative binomial mixed-effect independent/AR(1) model with two treatment groups described in Zhao et al. The group mean can be different at each time point, but no other covariates are allowed.See mle.fun, mle.ar1.fun for details of the model explanations.

Usage

rNBME.R(
	gdist = "G", n = 200, sn = 5, th = exp(1.3), 
	u1 = rep(1.5, 5), u2 = rep(1.5, 5), 
	a = exp(-0.5),d=NULL,  othrp = list(u.n = 3, s.n = 0.5, p.mx = 0.05, sh.mx = NA)
	)

Arguments

gdist

The distribution of the random effect term G_i.

If gdist="G", G_i is from the gamma distribution.

If gdist="N", G_i is from the log normal distribution.

If gdist="U", G_i (on t

The number of patients. It must be an even number.

The number of repeated measures per patient. Generated datasets are balanced design.

If gdist="G", th is a scale parameter of the gamma distribution.

If gdist="N" or gdist=="U", th is var(G_i).

If gdist="GN", see details.

If gdist=

A vector of length sn, specifying the mean of the treatment group 1 E(Y_ij) =u1[j].

A vector of length sn, specifying the mean of the treatment group 2 E(Y_ij) =u2[j].

The parameter a of the negative binomial mixed-effect independent model. See mle.fun.

If d=NULL, generate data from the independent model. If d is a scalar between 0 and 1, then d is delta in the AR(1) model, and generate datasets from the AR(1) model.

othrp

If gdist="GN", parameters for the GN option. See details. If gdist="NoN", othrp is a vector, containing a sample of G_i, which is treated as a population and G_i is resampled.

Value

idThe vector of length n*sn containing patient IDs: rep(1:n,each=sn)
vnThe vector of length n*sn containing the indicies of time points: rep(1:sn, n)
gpThe vector of length n*sn containing the indicies of the treatment groups
yThe vector of length n*sn containing generated response counts
gThe vector of length n*sn containing generated random effect terms
GparaThe record of the distribution and parameter specifications used to generate the dataset

Details

The generated datasets have equal number of scans per person.The number of patients in the two groups are the same.If gdist=="GN", datasets are generated from:

othrp$p.mx*N(mean=othrp$u.n,s.d=othrp$s.n) + (1-othrp$p.mx)*gamma(scale=th,shape),where shape of the gamma distribution is chosen to ensure E(G_i)=1.

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.

Examples

Run this code

## ================================================================================ ##
## 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 ~ log-normal(E(G_i)=1,var(G_i)=th)
set.seed(1)
sn = 5 # the number of repeated measures of each patient
n = 80 ## the number of patients
loga = - 0.5 
a = exp(loga) ## dispersion parameter 
logtheta <- 1.3
th = exp(logtheta) ## the variance of the gamma distributed random effect g_i


## No difference between the means of groups 
## The model only has an intercept term beta0 = 0.5
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

## Data 0 generated from the IND model
DT.ind= rNBME.R(gdist="N", n=n, a=a, th=th, u1=u1, u2=u2, sn=sn)
## Data 1
DT.ar= rNBME.R(gdist="N", n=n, a=a, th=th, u1=u1, u2=u2, sn=sn, d=0.4)

dt.ind=data.frame(CEL=DT.ind$y,Vcode=DT.ind$vn-2,ID=DT.ind$id)
dt.ar=data.frame(CEL=DT.ar$y,Vcode=DT.ar$vn-2,ID=DT.ar$id)
## ================================================================================ ##

#### Fit IND models
tst1=lmeNB(CEL~1, data=dt.ind, ID=dt.ind$ID, IPRT=T)
tst2=lmeNB(CEL~1, data=dt.ar, ID=dt.ar$ID, IPRT=T)
tst3=lmeNB(CEL~1, data=dt.ind, ID=dt.ind$ID, IPRT=T, RE="N")
tst4=lmeNB(CEL~1, data=dt.ar, ID=dt.ar$ID, IPRT=T, RE="N") #not printing
tst5=lmeNB(CEL~Vcode, data=dt.ind, ID=dt.ind$ID, IPRT=T, RE="N")
## conditional probability index
Psum5=index.batch(olmeNB=tst5, labelnp=dt.ind$Vcode >= 1, data=dt.ind, ID=dt.ind$ID)



## Fit the semi-parametric model
tst6=lmeNB(CEL~1, data=dt.ind, ID=dt.ind$ID, IPRT=T, RE="semipara")
tst7=lmeNB(CEL~Vcode, data=dt.ind, ID=dt.ind$ID, IPRT=T, RE="semipara")
## conditional probability index
Psum7=index.batch(olmeNB=tst7,labelnp=dt.ind$Vcode >= 1, data=dt.ind, 
                 ID=dt.ind$ID, Vcode=Vcode)


## Fit the AR1 models
tst8=lmeNB(CEL~1, data=dt.ind, ID=dt.ind$ID, IPRT=T,AR=T,Vcode=dt.ind$Vcode)
tst9=lmeNB(CEL~1, data=dt.ar, ID=dt.ar$ID, IPRT=T,AR=T,Vcode=dt.ar$Vcode)
## conditional probability index
Psum9=index.batch(olmeNB=tst9, labelnp=dt.ind$Vcode >= 1, data=dt.ind, 
                  ID=dt.ind$ID,Vcode=dt.ind$Vcode)

Run the code above in your browser using DataLab