CovYiYj.approx: Approximates the conditional covariance between all pairs of...

In gregorkb/aenetgt: Adaptive Elastic Net for Group Testing Data

Approximates the conditional covariance between all pairs of disease statuses with Gibbs sampling.

Description

Approximates the conditional covariance between all pairs of disease statuses with Gibbs sampling.

Usage

CovYiYj.approx(Z, Y, X, b, Se, Sp, EY, GI = 50000)

Arguments

Z	Group testing output from one of the functions individual.assay.gen, masterpool.assay.gen, dorfman.assay.gen, or array.assay.gen.
Y	Group testing output from one of the functions individual.assay.gen, masterpool.assay.gen, dorfman.assay.gen, or array.assay.gen.
X	Design matrix with first column a column of 1s.
b	Parameter values at which to compute the conditional expectations.
Se	A vector of testing sensitivities of length max(Z[,3]).
Sp	A vector of testing specificities of length max(Z[,3]).
EY	The vector of conditional expectations of individual disease statuses.
GI	The length of the Gibbs sampling Markov chain.

Value

An approximation to the matrix of conditional covariances between all pairs of disease statuses.

This function uses a Gibbs sampler to approximate the conditional covariance of all pairs of disease statuses, conditional on the observed assay data and the disease statuses of all other individuals. This is used to implement Louis's method for computing the observed information matrix associated with the likelihood function, from which Wald-type standard errors can be computed for the maximum likelihood estimators. For more details see Gregory et al. (2018+).

Examples

# generate individual covariate values and disease statuses
N <- 160
data <- model0(N)
X <- data$X
Y.true <- data$Yi
Se <- c(.95,.92) # set master pool and individual assay sensitivity
Sp <- c(.97,.98) # set master pool and individual assay specificity
cj <- 4 # set size of master pools
# subject individuals to dorfman testing
assay.data <- dorfman.assay.gen(Y.true,Se,Sp,cj)
Z <- assay.data$Z
Y <- assay.data$Y
b.mle <- mlegt(X, Y, Z, Se, Sp, tol = .01)$b.mle # compute mle
eta.mle <- X %*% b.mle
EY <- EYexact(Z,Y,eta.mle,Se,Sp)
px <- as.numeric(logit(b.mle,X))
CovYiYj <- CovYiYj.approx(Z,Y,X,b.mle,Se,Sp,EY)
# use Louis' method to get the observed information matrix
Hess <- t(X) %*% ( - diag(px * (1 - px))  +  CovYiYj  ) %*% X	
b.cov.est <-  solve( - Hess ) 
b.mle.se <- sqrt(diag(b.cov.est)) # estimated standard errors of mles	

gregorkb/aenetgt documentation built on Oct. 12, 2022, 11:51 a.m.