CP.ar1.se: Compute a conditional probability of observing a set of counts as extreme as the new observations of a subjectvisit given the previous observations of the same subject based on the negative binomial mixed-effect AR(1) model.

Description

Given the parameter estimates of c( alpha,theta,delta, beta0,beta1,...) 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.When input the maximum likelihood estimates of the parameters, CP.ar1.se returns the estimate of the conditional probability and its asymptotic standard error of a subject based on AR(1) model. The computation for the probability is done by its subroutine jCP.ar1, which 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, 
	  dist = "G", V, mc = FALSE, qfun = "sum")
jCP.ar1(tpar, ypre, ynew, 
	y2m=NULL, XM, stp, 
        mod = "G", LG = FALSE, MC = FALSE, N = 40000, qfun = "sum", oth = NULL)
CP1.ar1(ypre, ynew, y2m = NULL, 
	stp, u, th, a, 
	dt, mod = "G", gh, qfun = "sum")
MCCP.ar1(ypre, ynew, stp, u, th, a, dt, mod = "G", Ns = 1000, gh, qfun = "sum")

Arguments

tpar

A vector of length 4 + # covariates, containing the estimates of the model in the order that log(a),log(theta),logit(delta),beta0,beta1.... If random effects are assumed to be from the log-normal distribution, then theta is a va

ypre

A vector of counts of the length the number of previous observations.

y2m

Internal use only. Set as y2m=NULL

ynew

A vector of counts of the length the number of new observations.

A n_i by # covariates matrix containing the covariate values of subject i, whre n_i is the total number of previous and new observations. If there is no covariate, i.e., the model only has an intercept term, then set XM=NU

stp

A vector of length the total number of previous and new observations. 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

mod

If mod="G", then the conditional probability is computed by assuming the random effect is from the gamma distribution. If mod="N", then the conditional probability is computed by assuming the random effect is from the log-normal

If LG=TRUE then the logit of the conditional probability is returned.

If TRUE then the function MCCP.ar1 is called and the Monte carlo integration is performed. Fast but could be unreliable; not recommended for computing the confidence intervals. If FALSE then the function CP1.ar

The number of the Monte carlo integration. Necessary if MC=TRUE

qfun

The summary statistics q.

oth

If mod="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 mle.ar1.non3

dist

Same as mod.

The estimated variance covariance matrix of the parameters log(a),log(theta),log(delta),beta0,beta1....

The estimated theta.

The estimated alpha.

The estimated delta.

Same as N.

Same as MC.

A vector of length the number of repeated measures, containing the estimated mean ( mu_ij;j=1,...,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

Same as oth.

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

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)
dist = "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, 
	  dist = dist, 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, 
	  dist = dist, 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,
        mod = dist, LG = FALSE, MC = FALSE, N = 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, y2m = NULL, 
	stp =stp, u = u, th = exp(tpar[2]), a = exp(tpar[1]), 
	dt= ilgt(tpar[3]), mod = dist, 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]), 
	 mod =dist, Ns = 1000, qfun = "sum")

Run the code above in your browser using DataLab

Description

Usage

Arguments

References

See Also

Examples