library(groupTesting)
## To illustrate 'glm.gt', we use data simulated
## by the functions 'hier.gt.simulation' and 'array.gt.simulation'.
## Note: The simulated data-structures are consistent
## with the data-structure required for 'gtData'.
## Example 1: MLE from 3-stage hierarchical group testing data.
## The data used is simulated by 'hier.gt.simulation'.
N <- 200 # Sample size
S <- 3 # 3-stage hierarchical testing
psz <- c(6,2,1) # Pool sizes used in stages 1-3
Se <- c(.95,.95,.98) # Sensitivities in stages 1-3
Sp <- c(.95,.98,.96) # Specificities in stages 1-3
assayID <- c(1,2,3) # Assays used in stages 1-3
param.t <- c(-3,2,1) # The TRUE parameter to be estimated
# Simulating covariates:
set.seed(123)
x1 <- rnorm(N, mean=0, sd=0.75)
x2 <- rbinom(N, size=1, prob=0.5)
X <- cbind(1, x1, x2)
colnames( X ) <- c("Intercept", "Predictor 1", "Predictor 2")
# Note: Because the 1st column of X is 1, intercept model will be fit.
# Specifying logit inverse link:
g <- function(t){exp(t)/(1+exp(t))}
pReg <- g(X%*%param.t)
# Simulating test responses:
gtOut <- hier.gt.simulation(N,pReg,S,psz,Se,Sp,assayID)$gtData
# Fitting the model (with intercept):
param0 <- param.t + 0.2 # Initial value
res <- glm.gt(beta0=param0,gtData=gtOut,X=X,
g=g,dg=NULL,d2g=NULL,
grdMethod="central",covariance=TRUE,
nburn=2000,ngit=5000,maxit=200,
tol=1e-03,tracing=TRUE,conf.level=0.95)
# Note: Because dg and d2g are NULL (i.e., the exact derivatives
# are not given), numerical derivatives are used.
# Estimation results:
# > res
# $param
# [1] -2.840802 1.992916 0.677176
# $covariance
# [,1] [,2] [,3]
# [1,] 0.2134439 -0.10147555 -0.16693776
# [2,] -0.1014756 0.16855122 0.02997113
# [3,] -0.1669378 0.02997113 0.26324589
# $iterUsed
# [1] 10
# $convergence
# [1] 0
# $summary
# Estimate Std.Err 95%lower 95%upper
# Intercept -2.841 0.462 -3.746 -1.935
# Predictor 1 1.993 0.411 1.188 2.798
# Predictor 2 0.677 0.513 -0.328 1.683
## Example 2: MLE from two-dimensional array testing data.
## The data used is simulated by 'array.gt.simulation'.
N <- 200 # Sample size
protocol <- "A2" # 2-stage array without testing initial master pool
n <- 5 # Row/column size
Se <- c(0.95, 0.95) # Sensitivities
Sp <- c(0.98, 0.98) # Specificities
assayID <- c(1, 1) # The same assay in both stages
param <- c(-4,1,1) # The TRUE parameter to be estimated
# Simulating data:
set.seed(123)
x1 <- runif(N)
x2 <- rnorm(N, mean=0, sd=0.5)
x3 <- rbinom(N, size=1, prob=0.5)
X <- cbind(x1, x2, x3)
# Note: Because the 1st column of X is not 1,
# the model without intercept will be fit.
# Finding g, dg, and d2g from the function 'glmLink':
res0 <- glmLink(fn.name="logit")
g <- res0$g # Logit inverse link g()
dg <- res0$dg # The exact first derivate of g
d2g <- res0$d2g # The exact second derivate of g
pReg <- g(X%*%param) # Individual probabilities
gtOut <- array.gt.simulation(N,pReg,protocol,n,Se,Sp,assayID)$gtData
# Fitting the model (without intercept):
param0 <- param + 0.2
res <- glm.gt(beta0=param0,gtData=gtOut,X=X,g=g,
dg=dg,d2g=d2g,covariance=TRUE,
nburn=2000,ngit=5000,maxit=200,
tol=1e-03,tracing=TRUE,conf.level=0.95)
print(res)
# \donttest{
## Example 3: MLE from non-overlapping initial pooled responses.
## The data used is simulated by 'hier.gt.simulation'.
## Note: With initial pooled responses, our MLE is equivalent
## to the MLE in Vansteelandt et al. (2000).
N <- 1000 # Sample size
psz <- 5 # Pool size
S <- 1 # 1-stage testing
Se <- 0.95 # Sensitivity
Sp <- 0.99 # Specificity
assayID <- 1 # Assay used for all pools
param <- c(-3,2,1) # The TRUE parameter to be estimated
# Simulating data:
set.seed(123)
x1 <- rnorm(N, mean=0, sd=0.75)
x2 <- rbinom(N, size=1, prob=0.5)
X <- cbind(1, x1, x2)
# Finding g, dg, and d2g by the function 'glmLink':
res0 <- glmLink(fn.name="probit") # Probit link
g <- res0$g
dg <- res0$dg
d2g <- res0$d2g
pReg <- g(X%*%param)
gtOut <- hier.gt.simulation(N,pReg,S,psz,Se,Sp,assayID)$gtData
# Fitting the model:
param0 <- param + 0.2
res <- glm.gt(beta0=param0,gtData=gtOut,X=X,g=g,
dg=dg,d2g=d2g,covariance=TRUE,
nburn=2000,ngit=5000,maxit=200,
tol=1e-03,tracing=TRUE,conf.level=0.95)
print(res)
## Example 4: MLE from individual (one-by-one) testing data.
## The data used is simulated by 'hier.gt.simulation'.
N <- 1000 # Sample size
psz <- 1 # Pool size 1 (i.e., individual testing)
S <- 1 # 1-stage testing
Se <- 0.95 # Sensitivity
Sp <- 0.99 # Specificity
assayID <- 1 # Assay used for all pools
param <- c(-3,2,1) # The TRUE parameter to be estimated
# Simulating data:
set.seed(123)
x1 <- rnorm(N, mean=0, sd=0.75)
x2 <- rbinom(N, size=1, prob=0.5)
X <- cbind(1, x1, x2)
g <- function(t){exp(t)/(1+exp(t))} # Inverse logit
pReg <- g(X%*%param)
gtOut <- hier.gt.simulation(N,pReg,S,psz,Se,Sp,assayID)$gtData
# Fitting the model:
param0 <- param + 0.2
res <- glm.gt(beta0=param0,gtData=gtOut,
X=X,g=g,dg=NULL,d2g=NULL,
grdMethod="central",covariance=TRUE,
nburn=2000,ngit=5000,maxit=200,
tol=1e-03,tracing=TRUE,conf.level=0.95)
print(res)
## Example 5: Using pooled testing data.
# Pooled test outcomes:
Z <- c(1, 0, 1, 0, 1, 0, 1, 0, 0)
# Design matrix, X:
x1 <- c(0.8,1.2,0.4,1.5,1.8,1.8,0.1,1.6,0.2,0.2,1.8,0.2)
x2 <- c(31,56,45,64,26,47,22,60,35,41,32,41)
X <- cbind(x1, x2)
# Pool sizes used:
psz <- c(6, 6, 2, 2, 2, 1, 1, 1, 1)
# Pool-specific Se & Sp:
Se <- c(.90, .90, .95, .95, .95, .92, .92, .92, .92)
Sp <- c(.92, .92, .96, .96, .96, .90, .90, .90, .90)
# Assays used:
Assay <- c(1, 1, 2, 2, 2, 3, 3, 3, 3)
# Pool members:
Memb <- rbind(
c(1, 2, 3, 4, 5, 6),
c(7, 8, 9, 10, 11, 12),
c(1, 2, -9, -9, -9, -9),
c(3, 4, -9, -9, -9, -9),
c(5, 6, -9, -9, -9, -9),
c(1,-9, -9, -9, -9, -9),
c(2,-9, -9, -9, -9, -9),
c(5,-9, -9, -9, -9, -9),
c(6,-9, -9, -9, -9, -9)
)
# The data-structure suited for 'gtData':
gtOut <- cbind(Z, psz, Se, Sp, Assay, Memb)
# Fitting the model with logit link:
g <- function(t){exp(t)/(1+exp(t))}
param0 <- c(0, 0)
res <- glm.gt(beta0=param0,gtData=gtOut,X=X,
g=g,dg=NULL,d2g=NULL,
grdMethod="central",covariance=TRUE,
nburn=2000,ngit=5000,maxit=200,
tol=1e-03,tracing=TRUE,conf.level=0.95)
print(res)
# }
Run the code above in your browser using DataLab