jCP_ar1: Compute a conditional probability of observing a set of...

In lmeNB: Compute the Personalized Activity Index Based on a Negative Binomial Model

Description

Given the parameter estimates of α,θ, δ, β[0], β[1],... of the negative binomial mixed effect AR(1) model, these functions compute the following conditional probability:

Pr(q(Y[i,new])>=q(y[i,new])|Y[i,pre]=y[i,pre]) ,

where y[i,new] and y[i,pre] are vectors of previous and new observations from subject i and q() is a function which provides a scalar summary of the new observations.

These functions are subroutines of index.batch. CP.ar1.se returns the estimate of the conditional probability and its asymptotic standard error of a subject based on AR(1) model. The evaluations for the probability is done by its subroutine jCP.ar1, which, in turn, has two subroutines CP1.ar1 and MCCP.ar1. CP1.ar1 computes the probability via the adaptive quadrature while MCCP.ar1 computes the probability via the Monte Carlo integration.

Usage

CP.ar1.se(tpar, ypre, ynew, y2m = NULL, XM, stp, 
	  RE = "G", V, MC = FALSE, qfun = "sum",i.tol=1E-75)

jCP.ar1(tpar, ypre, ynew, y2m = NULL, XM, stp, RE = "G", LG =FALSE, 
	      MC = FALSE, N.MC = 40000, qfun = "sum", oth =NULL,i.tol=1E-75)

CP1.ar1(ypre, ynew, y2m=NULL, stp, u, th, a, dt, RE = "G", oth,qfun,i.tol=1E-75)

MCCP.ar1(ypre, ynew, stp, u, th, a, dt, RE = "G", N.MC = 1000, oth, qfun = "sum") 

Arguments

tpar	A vector of length 4 + # covariates, containing the estimates of the model in the order that \log(α),\log(θ),logit(δ),β[0],β[1],.... If the semi-parametric approach is taken, then theta is a place holder and can be any value.
ypre	A vector of the length the number of previous observations, containing counts on pre-scans.
y2m	Internal use only. Set as y2m=NULL.
ynew	A vector of length the number of new observations, containing counts on new scans.
XM	See CP.se.
stp	A vector of length n[i], containing index to indicates missing scans. The first entry must be zero. For example, if there is no missing scans and there are five repeated measures, then stp=c(0,1,1,1,1). If the third scan is missing and there are four repeated measures, then stp=c(0,1,2,1).
RE	See lmeNB. Note that this option is NOT accepted in CP.ar1.se.
LG	See CP.se.
MC	If TRUE then the function MCCP.ar1 is called and the Monte carlo integration is performed to integrate out the random effect. Fast but could be unreliable; not recommended for computing the confidence intervals. If FALSE then the function CP1.ar1 is called and the adaptive quadrature is performed. Slow but reliable.
N.MC	The number of the Monte Carlo integration. Necessary if MC=TRUE.
qfun	See index.batch.
oth	See CP.se. If RE="NoN", othr must be the frequency table of the random effects, which can be obtained based on dist=obj$gi where obj is the output of fitSemiAR1.
V	See CP.se.
th	The estimated theta.
a	The estimated alpha.
dt	The estimated delta.
u	A vector of length the number of repeated measures, containing the estimated mean counts ( μ[i1],...,μ[i n[i] ]). If the mean of Y_ij is modeled linearly on beta with the log-link function, then u=exp( beta0+ XM[,1]*beta1+XM[,2]*beta2+...)).
i.tol	See lmeNB

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 internal functions of lmeNB for fitting relevant models: fitParaIND, fitParaAR1, fitSemiIND, fitSemiAR1,

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

The functions to generate simulated datasets: rNBME.R.

Examples

## Not run: 
ilgt <- function (x) 
{
    tem = exp(x)
    res = tem/(1 + tem)
    return(res)
}
lgt <- function (p) 
{
    log(p/(1 - p))
}
## the vector of a parameter estimates if log(a),log(theta),logit(delta),beta0.
tpar  <- c(log(2),log(0.5),lgt(0.5),2)
ypre <- c(0, 1)
ynew <- c(1, 0, 0)
## No covariate
XM <- NULL
## no missing visit
stp <- c(0,1,1,1,1)
RE <- "G"
## The estimate of the variance covariance matrix
V <-
matrix(
c( 0.17720309, -0.240418504,  0.093562548,  0.009141980,
  -0.24041850,  0.605132808, -0.160454773, -0.003978118,
   0.09356255, -0.160454773,  0.095101658,  0.005661923,
   0.00914198, -0.003978118,  0.005661923,  0.007574769),
nrow=4)

## the estimate of the conditional probability based on the sum summary statistics and its SE
CP.ar1.se(tpar = tpar, ypre = ypre, ynew = ynew, 
	  XM =XM, stp = stp, 
	  RE = RE, V = V, MC = FALSE, qfun = "sum")

## the estimate of the conditional probability based on the max summary statistics and its SE
CP.ar1.se(tpar = tpar, ypre = ypre, ynew = ynew, 
	  XM =XM, stp = stp, 
	  RE = RE, V = V, MC = FALSE, qfun = "max")


## CP.ar1.se calls for jCP.ar1 to compute the estimate of the conditional probability
## the estimate of the conditional probability based on the sum summary statistics
jCP.ar1(tpar = tpar, ypre = ypre, ynew = ynew,
	y2m=NULL,  XM =XM, stp = stp,
        RE = RE, LG = FALSE, MC = FALSE, N.MC = 40000, qfun = "sum", oth = NULL)

## jCP.ar1 calls for CP.ar1 to compute the estimate of the conditional probability 
## via the adaptive quadrature (MC=F)
## the estimate of the conditional probability

u <- rep(exp(tpar[4]),length(ypre)+length(ynew))

CP1.ar1(ypre = ypre, ynew =ynew, 
	stp =stp, u = u, th = exp(tpar[2]), a = exp(tpar[1]), 
	dt= ilgt(tpar[3]), RE = RE, qfun = "sum")


## jCP.ar1 calls for CP.ar1 to compute the estimate of the conditional probability 
## via the Monte Carlo method (MC=T)
## the estimate of the conditional probability
MCCP.ar1(ypre = ypre, ynew =ynew, stp = stp, 
	 u = u, th = exp(tpar[2]), a = exp(tpar[1]),  dt = ilgt(tpar[3]), 
	 RE = RE, N.MC = 1000, qfun = "sum")

## End(Not run)

lmeNB documentation built on May 2, 2019, 3:34 p.m.