Unidentifiable parameters in multiple sharing per-contact infectivity model
Erik Erhardt & Yuridia Leyva
October 06, 2018
This whitepaper provides a simulation study to understand why only a single-parameter model is reliable in an infectivity study where surveys ask about multiple sharing events in a previous time period with a response of negative/positive for the infection.
1 Executive summary
The one-parameter maximum likelihood estimate (MLE) for per-event probability of infection given below as eq. 1 with the special case of eq. 2 works very well (Jewell and Shiboski 1990, Boelen et al. (2014)). However, the two-parameter model as in (Hughes et al. 2012, Magaret et al. (2015)) in those same equations has extreme bias due to unidentifiability issues. In general, the multiparameter model gives unreliable results due to parameter identifiability issues when multiple sharing events occured over the survey period and the response changed from negative to positive. Descriptively, because we don’t know which mode of sharing caused the response change and we don’t know which (or how many) sharing events could have individually been responsible for the change, accurate inference is not possible and this is reflected in the shape of the likelihood function and ultimately reflected in the wildly incorrect results.
2 Data description
The UFO (yoU Find Out) study was used to obtain an estimate of the per-contact infectivity of HCV infection associated with young IDUs engaging in RNS and/or RES. The UFO study is a prospective study of viral Hepatitis C among non-infected young adults generally under the age of 30 who engaged in the use of injection drugs within the past month. Subjects were sampled from San Francisco, CA by outreach workers and by word-of-mouth between 2003-2008 and 2010-2014 . A possible limitation of this sampling strategy is the uncertainty of how representative the sample is of the young injection drug user population in San Francisco or elsewhere. Participating young IDUs required an HCV-negative status at baseline screening; they were tested using a viremia test (HCV RNA) or tested for HCV antibodies (anti-HCV). Every three months, they received follow-up testing for HCV and were questioned by an interviewer regarding their exposures to HCV via injection drug use. Subjects received $10 USD for the screening visit and $20 upon return for their HCV test results. Descriptions of the study design and methods for the UFO cohort have been published in detail .
2.1 Data Description and Data Cleaning
Some of the participants in the UFO study were tested using both anti-HCV and HCV RNA tests, but had serological results that were non-coincidental. Individuals who were HCV RNA or anti-HCV positive were placed in the seroconverted category (status = 1). The remaining individuals were placed in the non-seroconverted category (status = 0), as they had results where either both types of HCV tests were negative, or one of the HCV tests was negative and the other was unknown.
The variables used for this analysis consisted strictly of receptive exposures to HCV. Receptive exposures consist of the young IDU having used a needle or injecting paraphernalia after it was used by someone else. The variables were placed into one of two exposure categories: needle or equipment exposures for subjects who engaged in RNS or RES, respectively. Table~ lists the questions corresponding to each of the variables used, as well as its classification of exposure type.
Total number of injection exposures in the last 30 days (both receptive and non-receptive) were calculated as the product of the number of days injected in the last 30 days and the number of times injected per day.
Variables with Likert scale responses were expressed as probabilities: (1) “Always” 1.00, (2) “Usually” 0.75, (3) “Sometimes” 0.50, (4) “Rarely” 0.25, (5) “Never” 0.00.
UFO staff interviewed and surveyed subjects every three months, but some subjects had long time gaps between surveys. Within 180 days, a small proportion of UFO subjects have cleared and become reinfected with HCV . Series of surveys were thus separated into multiple monitoring windows, as shown in Figure~. For example, suppose a person has been interviewed a total of 8 times, but there is a gap of 180 days or more between the 4th and 5th survey. To account for this long time period between surveys, the first 4 surveys are placed in monitoring window 1, and the last 4 surveys are placed in monitoring window 2.
For each subject, a decision was then made as to which observations to keep in order to conduct the analysis. Subjects required a serological status of 0 (HCV-negative) at baseline in order for probabilities to be calculated. Only the first observation with serological status of 1 (HCV-positive) was useful to help quantify the number of exposures it took to seroconvert an HCV-negative person to HCV-positive. Therefore, leading 1s, as well as any values after the first HCV positive status of 1 following any zeros, were removed. An example of this is illustrated in Table~.
It may be the case that IDUs who naturally clear HCV (have a positive serostatus before a negative in their survey timeline) have better capability to “deal with” the virus than other IDUs. If this is the case, then including these IDUs in the sample may underestimate the per-contact infectivity rate. Therefore, a second analysis with the additional selection criterion that excludes IDUs who were positive at their first survey was also conducted. (We call this L0s for “leading 0s, only” in the raw data.)
We adjusted the number of exposures relative to the actual intersurvey times. First, the number of days between surveys was divided by 90 days (the length of time about which the surveys asked) to obtain a “stretch” multiplier. Then the reported 90-day exposure values were multiplied by this stretch factor to account for intersurvey times more or less than 90 days. An example of this is in Figure~.
3 Statistical model and maximum likelihood estimate for per-contact infectivity rates
Our goal is to estimate the per-contact infectivity rates, \(\mathbf{\beta_{\textrm{N}}}\) and \(\mathbf{\beta_{\textrm{E}}}\), of HCV transmission for receptive needle sharing and receptive ancillary equipment sharing, respectively. We estimate \(\mathbf{\beta_{\textrm{N}}}\) and \(\mathbf{\beta_{\textrm{E}}}\) via a maximum likelihood estimate (MLE), which is the value of the parameter that makes the data most likely under the model. We use the following likelihood function, \(L\), for a sample size of \(N\) participants and \(S_i\) surveys): \[ L(\mathbf{\beta_{\textrm{N}}},\mathbf{\beta_{\textrm{E}}} \,|\,n_{\textrm{N}ij}, n_{\textrm{E}ij}, y_{ij}) = \prod_{i=1}^{N}\prod_{j=1}^{S_i}f_{ij}^{y_{ij}}(1-f_{ij})^{(1-y_{ij})}. \qquad(1)\]
To model the UFO data, the probabilities of the data must be considered a function of the \(\beta\) parameters in the model. The following probability mass function (pmf) was used for each subject \(i=1,\ldots,N\) and all subjects’ surveys \(j=1,\ldots,S_i\): \[ f_{ij} = f_{ij}(\mathbf{\beta_{\textrm{N}}},\mathbf{\beta_{\textrm{E}}} \,|\,n_{\textrm{N}ij},n_{\textrm{E}ij}) = 1-(1-\mathbf{\beta_{\textrm{N}}})^{n_{\textrm{N}ij}}(1-\mathbf{\beta_{\textrm{E}}})^{n_{\textrm{E}ij}}, \qquad(2)\] where \(n_{\textrm{N}ij}\) is the number of receptive exposures to needles and \(n_{\textrm{E}ij}\) is the number of exposures associated with receptive sharing of ancillary injecting equipment (such as cookers or cottons). \(y_{ij}\) is the status of seroconversion; \(y_{ij}=0\) for HCV-negative participants, and \(y_{ij}=1\) for HCV-positive participants. This probability model assumes exposure probabilities are independent.
To obtain the MLEs, we use the log of \(L\) rather than the likelihood function itself in order to avoid numerical overflow issues. To maximize the log-likelihood, we minimize the negative log-likelihood using the optim() procedure in R. A uni-parameter analysis of per-incident seroconversion associated with needles alone assumes \(\mathbf{\beta_{\textrm{E}}} \equiv 0\). There is no general closed-form solution for \(\mathbf{\beta_{\textrm{N}}}\), even in the uni-parameter case, thus a numerical optimization procedure is necessary.
The per-contact infectivities of Hepatitis C Virus associated with RNS and RES were estimated using bootstrap intervals at the \(95\%\) confidence level. To construct a bootstrap estimate of the MLE sampling distribution, we obtained 1000 bootstrap resamples by sampling with replacement from IDUs including all their “clean” surveys. Then we performed uni-parameter and bi-parameter maximum likelihood estimates on each resample.
4 Simulation study
To understand the behavior of the (negative log) likelihood function and the MLE, a series of small cases are followed by a larger data-realistic case.
Define functions for simulating data and evaluating the negative log-likelihood function for optimization.
# likelihood function (log-likelihood is often more computationally stable)
f.log.like.beta <- function(b, n, y) {
# b is a vector, n is a data frame, and y is a vector
n <- as.matrix(n)
#library(Rmpfr)
#n <- mpfr(n, 200) # 120 bit numbers, 53 is "double"
#org# # pmf
#org# f.i <- 1 - exp(apply( t(t(n) * log(1 - b)), 1, sum))
#org# # log likelihood
#org# l.i <- y * log(f.i) + (1 - y) * log(1 - f.i)
#org# # compute sum, removing NaN and Inf values
#org# sum.log.likelihood <- sum( l.i[is.finite(l.i)])
## Improve numerical stability for y=0 (can't do for y=1)
# pmf
log.f.i.stab <- apply( t(t(n) * log(1 - b)), 1, sum)
#log.f.i <- 1 - exp(apply( t(t(n) * log(1 - b)), 1, sum))
# log likelihood
#l.i.stab <- y * log(1 - exp(log.f.i.stab)) + (1 - y) * log.f.i.stab
l.i.stab.1_y <- (1 - y) * log.f.i.stab # fixes y=0 loglike for larger b
library(Rmpfr) # fixes y=1 loglike for larger b
l.i.stab.y <- y * log(1 - exp(mpfr(log.f.i.stab, 200)))
l.i.stab <- l.i.stab.1_y + l.i.stab.y
# compute sum, removing NaN and Inf values
sum.log.likelihood <- sum( l.i.stab[is.finite(l.i.stab)])
sum.log.likelihood <- as.double(sum.log.likelihood)
#detach(package:Rmpfr)
# return the negative log-likelihood value to minimize
return( - sum.log.likelihood )
}
##### improved above by fixing overflow issues
### # likelihood function (log-likelihood is often more computationally stable)
### f.log.like.beta <- function(b, n, y) {
### # b is a vector, n is a data frame, and y is a vector
### n <- as.matrix(n)
###
### # pmf
### f.i <- 1 - exp(apply( t(t(n) * log(1 - b)), 1, sum))
### # log likelihood
### l.i <- y * log(f.i) + (1 - y) * log(1 - f.i)
### # compute sum, removing NaN and Inf values
### sum.log.likelihood <- sum( l.i[is.finite(l.i)])
###
### #print(b)
### #print( - sum.log.likelihood)
###
### # return the negative log-likelihood value to minimize
### return( - sum.log.likelihood )
### }
bound.lower <- 1e-10
# function to test small simulations
f.est.betas.sim <- function(
N.types
, probs.seroconversion
, dat.clean
, bound.lower = 1e-10) {
log.likelihood.true.betas <- f.log.like.beta( b = probs.seroconversion
, n = dat.clean[,grep("^n.", colnames(dat.clean))]
, y = dat.clean$y.seroconversion
)
#print("Optim probs.seroconversion")
# Find best beta using optim to maximize the log-likelihood
# (by minimizing the negative log-likelihood)
# beta optim
optim.log.likelihood.beta.dat <-
optim(
# 1 2 3
par = rep(0.1, length(grep("^n.", colnames(dat.clean)))) #probs.seroconversion #
, n = dat.clean[,grep("^n.", colnames(dat.clean))]
, lower = rep(bound.lower, length(grep("^n.", colnames(dat.clean))))
, upper = rep(0.999, length(grep("^n.", colnames(dat.clean))))
, fn = f.log.like.beta
, y = dat.clean$y.seroconversion
, method = "L-BFGS-B"
, control = list(factr=1e8, pgtol=1e-6)
#, control = list(trace=TRUE) #, maxit = 1e3) #, factr=1e11, pgtol=1e-8)
)
# print("Convergence Status")
# print(optim.log.likelihood.beta.dat$message)
# print("True -loglike")
optim.log.likelihood.beta.dat$log.likelihood.true.betas <- log.likelihood.true.betas
# print(log.likelihood.true.betas)
# print("Est -loglike")
# print(optim.log.likelihood.beta.dat$value)
# print("")
# print("True probs.seroconversion")
# print(probs.seroconversion)
# print("Est probs.seroconversion")
# print(optim.log.likelihood.beta.dat$par)
return(optim.log.likelihood.beta.dat)
} # f.est.betas.sim
## Plot 2D surface
f.plot.2D <- function(
N.types
, probs.seroconversion
, dat.clean
, optim.results
, f.name = NULL
, seq.length = 51
, seq.limits = c(1e-20, 1 - 1e-20)) {
# name used in plot below
#f.name <- "dat.clean, 1 and 2"
# plot the function
# define ranges of x to plot over and put into matrix
#b1 <- seq(1e-20, 3 * probs.seroconversion[1], length = seq.length)
#b2 <- seq(1e-20, 3 * probs.seroconversion[2], length = seq.length)
b1 <- seq(seq.limits[1], seq.limits[2], length = seq.length)
b2 <- seq(seq.limits[1], seq.limits[2], length = seq.length)
b.grid <- as.matrix(expand.grid(b1, b2))
colnames(b.grid) <- c("b1", "b2")
# evaluate function
ff <- rep(NA, nrow(b.grid))
for (i in 1:nrow(b.grid)) {
ff[i] <- f.log.like.beta( b = b.grid[i,]
, n = dat.clean[,grep("^n.", colnames(dat.clean))]
, y = dat.clean$y.seroconversion
)
}
# put X and y values in a data.frame for plotting
df <- data.frame(ff, b.grid)
# # plot the function
# library(lattice) # use the lattice package
# wireframe(ff ~ b1 * b2 # y, x1, and x2 axes to plot
# , data = df # data.frame with values to plot
# #, main = f.name # name the plot
# , shade = TRUE # make it pretty
# , scales = list(arrows = FALSE) # include axis ticks
# , screen = list(z = 30, x = -70) # view position
# )
image(sort(unique(df$b1)), sort(unique(df$b2)), matrix(df$ff, nrow = length(b1))
, main = f.name, sub = "+ = true, x = est"
, xlab = "beta1", ylab = "beta2"
, col = cm.colors(12))
contour(sort(unique(df$b1)), sort(unique(df$b2)), matrix(df$ff, nrow = length(b1)), nlevels = 25, add=TRUE)
points(probs.seroconversion[1], probs.seroconversion[2], pch = 3, cex = 3, col = "red")
points(optim.results$par[1], optim.results$par[2] , pch = 4, cex = 3, col = "green4")
#library(rgl)
#plot3d(df$b1, df$b2, df$ff)
##surface3d(unique(df$b1), unique(df$b2), df$ff)
##persp3d(unique(df$b1), unique(df$b2), df$ff)
}## Bootstrap
f.est.betas.bs <- function(
dat.clean
, R = 1000) {
# Resample id with replacement
id.unique <- sort(unique(dat.clean$id))
dat.resample <- rbind(dat.clean, dat.clean, dat.clean)
bound.lower <- 1e-10
optim.beta <- matrix(NA, nrow=R, ncol=length(grep("^n.", colnames(dat.clean)))
, dimnames = list(1:R, paste("beta", 1:length(grep("^n.", colnames(dat.clean))), sep="")))
optim.status <- data.frame(convergence = rep(NA, R)
, message = rep(NA, R))
r <- 0 # init
r.every <- 0
ptm <- proc.time()
while( (r < R) & (r.every < 10 * R)) {
## for (r in 1:R) {
r <- r + 1
r.every <- r.every + 1
cat(" ", r)
if ((r %% 10) == 0) {cat("\n")}
id.resample <- sort(sample(id.unique, replace = TRUE))
dat.resample[,] <- NA
ind <- 0
for (i.id in id.resample) {
one.resample <- subset(dat.clean, id == i.id)
len.resample <- dim(one.resample)[1]
dat.resample[ind + 1:len.resample, ] <- one.resample
ind <- ind + len.resample
}
dat.use <- na.omit(dat.resample)
bound.lower = 1e-10
optim.log.likelihood.beta.dat <-
optim(
# 1 2 3
par = rep(0.1, length(grep("^n.", colnames(dat.use)))) #probs.seroconversion #
, n = dat.use[,grep("^n.", colnames(dat.use))]
, lower = rep(bound.lower, length(grep("^n.", colnames(dat.use))))
, upper = rep(0.999, length(grep("^n.", colnames(dat.use))))
, fn = f.log.like.beta
, y = dat.use$y.seroconversion
, method = "L-BFGS-B"
, control = list(factr=1e8, pgtol=1e-6)
#, control = list(trace=TRUE) #, maxit = 1e3) #, factr=1e11, pgtol=1e-8)
)
optim.beta[r,] <- optim.log.likelihood.beta.dat$par
optim.status$convergence[r] <- optim.log.likelihood.beta.dat$convergence
optim.status$message [r] <- optim.log.likelihood.beta.dat$message
# if (any(optim.log.likelihood.beta$par == bound.lower)) {
# # if a value is at the lower bound, retry
# r <- r - 1
# cat("r")
# } else {
# optim.beta[r,] <- optim.log.likelihood.beta$par
# }
}
time.bs <- proc.time() - ptm
# # optim success
# r
# r.every
# r/r.every
return(data.frame(optim.beta, optim.status))
} # f.est.betas.bs4.1 Simulations summary
Estimation of per-contact infectivty is often accurate. For any number of types of exposures, when there is a unique exposure of one type at seroconversion, the estimation is unbiased. When there are multiple exposures of one type at seroconversion, it is not possible to know how many of the exposures would have resulted in seroconversion – thus the per-contact probability is estimated between the minimum and maximum probabilities, closer to the minimum.
With two exposure types when there are multiple exposures at seroconversion the MLE may locate on the boundary of the parameter space. In this case, the loglikelihood function favors the larger \(\beta\) and sends the other to zero. Also, when the number of exposures becomes very large, and the computer code for the (negative log) likelihood function is written in a direct way, the function becomes jagged due to numerical underflow (a probability raised to a high power goes to zero). Therefore, the R package for arbitrarily precise numbers was implemented to increase bit precision from 53 (“double”) to 200, and this relieved the large-\(n\) underflow issue.
For the notation in the following tables, let there be two exposure types, and let n.1 and n.2 be the number of exposures of each type for four survey observations (4 rows of data). Let y.1 and y.2 indicate the unobserved true exposure set that caused a seroconversion for each type. Let y.seroconversion be the observed HCV status (0=negative, 1=positive).
The data-realistic simulation estimates the parameter on the boundary, which is consistent with the simpler simulation cases. It favors the larger \(\beta\) and sends the other to zero.
4.2 Simulations: one type with nonoverlapping/separate exposures
4.2.1 Simulation: Univariate, one seroconverted obs: \(\beta_1 = 1/21\)
Univariate optimization is exact when there’s one observation at seroconversion.
| y.seroconversion | n.1 | y.1 |
|---|---|---|
| 0 | 20 | 0 |
| 1 | 1 | 1 |
# Find best beta over a grid of possible beta values
beta.range <- seq(1e-20, 1 - 1e-20, length = 1000)
log.likelihood.range <- rep(NA, length(beta.range))
for (i.beta in 1:length(beta.range)) {
log.likelihood.range[i.beta] <- f.log.like.beta(b = beta.range[i.beta]
, n = dat.clean[,grep("^n.1", colnames(dat.clean))]
, y = dat.clean$y.seroconversion
)
}
log.likelihood.true.betas <- f.log.like.beta(b = probs.seroconversion
, n = dat.clean[,grep("^n.1", colnames(dat.clean))]
, y = dat.clean$y.seroconversion
)
grid.log.likelihood.beta <- beta.range[which(log.likelihood.range == min(log.likelihood.range))]
#grid.log.likelihood.beta
########################
optim.log.likelihood.beta.dat <-
optim(
# 1 2 3
par = rep(0.1, length(grep("^n.", colnames(dat.clean)))) #probs.seroconversion #
, n = dat.clean[,grep("^n.", colnames(dat.clean))]
, lower = rep(bound.lower, length(grep("^n.", colnames(dat.clean))))
, upper = rep(0.999, length(grep("^n.", colnames(dat.clean))))
, fn = f.log.like.beta
, y = dat.clean$y.seroconversion
, method = "L-BFGS-B"
, control = list(factr=1e11, pgtol=1e-16)
#, control = list(trace=TRUE) #, maxit = 1e3) #, factr=1e11, pgtol=1e-8)
)
optim.log.likelihood.beta.dat$log.likelihood.true.betas <- log.likelihood.true.betas
# plot log-likelihood with true and estimated betas
plot(beta.range, log.likelihood.range, type = "l"
, main = "negative log-likelihood by beta"
, sub = "red=true, blue=grid, green=optim")
abline(v = probs.seroconversion, col = "red") # line at true value
abline(v = grid.log.likelihood.beta, col = "blue") # line at best value via grid
abline(v = optim.log.likelihood.beta.dat$par, col = "green") # line at best value via optim
Model results: CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH
| Type | \(-\log(L)\) | \(\beta_1\) |
|---|---|---|
| True | 4.02033 | 0.047619 |
| Est | 4.02033 | 0.0476348 |
4.2.2 Simulation: Univariate, four seroconverted obs: \(\beta_1 = 4/21\)
When there’s multiple observations at seroconversion, convergence is inside the range of plausible correct values. In this case, the estimated \(\beta_1\) is between \(1/21\) and \(4/21\).
| y.seroconversion | n.1 | y.1 |
|---|---|---|
| 0 | 17 | 0 |
| 1 | 4 | 1 |
# Find best beta over a grid of possible beta values
beta.range <- seq(1e-20, 1 - 1e-20, length = 1000)
log.likelihood.range <- rep(NA, length(beta.range))
for (i.beta in 1:length(beta.range)) {
log.likelihood.range[i.beta] <- f.log.like.beta(b = beta.range[i.beta]
, n = dat.clean[,grep("^n.1", colnames(dat.clean))]
, y = dat.clean$y.seroconversion
)
}
log.likelihood.true.betas <- f.log.like.beta(b = probs.seroconversion
, n = dat.clean[,grep("^n.1", colnames(dat.clean))]
, y = dat.clean$y.seroconversion
)
grid.log.likelihood.beta <- beta.range[which(log.likelihood.range == min(log.likelihood.range))]
#grid.log.likelihood.beta
########################
optim.log.likelihood.beta.dat <-
optim(
# 1 2 3
par = rep(0.1, length(grep("^n.", colnames(dat.clean)))) #probs.seroconversion #
, n = dat.clean[,grep("^n.", colnames(dat.clean))]
, lower = rep(bound.lower, length(grep("^n.", colnames(dat.clean))))
, upper = rep(0.999, length(grep("^n.", colnames(dat.clean))))
, fn = f.log.like.beta
, y = dat.clean$y.seroconversion
, method = "L-BFGS-B"
, control = list(factr=1e11, pgtol=1e-16)
#, control = list(trace=TRUE) #, maxit = 1e3) #, factr=1e11, pgtol=1e-8)
)
optim.log.likelihood.beta.dat$log.likelihood.true.betas <- log.likelihood.true.betas
# plot log-likelihood with true and estimated betas
plot(beta.range, log.likelihood.range, type = "l"
, main = "negative log-likelihood by beta"
, sub = "red=true, blue=grid, green=optim")
abline(v = probs.seroconversion, col = "red") # line at true value
abline(v = grid.log.likelihood.beta, col = "blue") # line at best value via grid
abline(v = optim.log.likelihood.beta.dat$par, col = "green") # line at best value via optim
probs.seroconversion2 <- c(1/21)
log.likelihood.true.betas2 <- f.log.like.beta(b = probs.seroconversion2
, n = dat.clean[,grep("^n.1", colnames(dat.clean))]
, y = dat.clean$y.seroconversion
)
abline(v = probs.seroconversion2, col = "red") # line at true value
Model results: CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH
| Type | \(-\log(L)\) | \(\beta_1\) |
|---|---|---|
| True | 4.15342 | 0.190476 |
| Est | 2.55629 | 0.051463 |
| True 2 | 2.55936 | 0.047619 |
4.3 Simulations: two types with nonoverlapping exposures
In the following examples, when there is a single exposure at seroconversion the \(\beta\) MLE estimates the true parameters without bias. When there are multiple exposures at seroconversion for a single exposure type, then there is censoring, that is, the specific exposure(s) responsible for seroconversion is unknown. In this case, the MLE is between the lower and upper bounds of the minimum and maximum number of exposures responsible for seroconversion.
4.3.1 Separate, one seroconverted obs: \(\beta_1 = 1/3\), \(\beta_2 = 1/5\)
Convergence to the correct value.
| y.seroconversion | n.1 | y.1 | n.2 | y.2 |
|---|---|---|---|---|
| 0 | 0 | 0 | 4 | 0 |
| 1 | 0 | 0 | 1 | 1 |
| 0 | 2 | 0 | 0 | 0 |
| 1 | 1 | 1 | 0 | 0 |
optim.log.likelihood.beta.dat <- f.est.betas.sim(N.types, probs.seroconversion, dat.clean, bound.lower = 1e-10)
f.plot.2D(N.types, probs.seroconversion, dat.clean, optim.log.likelihood.beta.dat, f.name)
Model results: CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH
| Type | \(-\log(L)\) | \(\beta_1\) | \(\beta_2\) | plot symbol |
|---|---|---|---|---|
| True | 4.41155 | 0.333333 | 0.2 | + |
| Est | 4.41155 | 0.333333 | 0.200001 | x |
4.3.2 Separate, one seroconverted obs: \(\beta_1 = 1/4\), \(\beta_2 = 1/6\)
Convergence to the correct value.
| y.seroconversion | n.1 | y.1 | n.2 | y.2 |
|---|---|---|---|---|
| 0 | 0 | 0 | 5 | 0 |
| 1 | 0 | 0 | 1 | 1 |
| 0 | 3 | 0 | 0 | 0 |
| 1 | 1 | 1 | 0 | 0 |
optim.log.likelihood.beta.dat <- f.est.betas.sim(N.types, probs.seroconversion, dat.clean, bound.lower = 1e-10)
f.plot.2D(N.types, probs.seroconversion, dat.clean, optim.log.likelihood.beta.dat, f.name)
Model results: CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH
| Type | \(-\log(L)\) | \(\beta_1\) | \(\beta_2\) | plot symbol |
|---|---|---|---|---|
| True | 4.95271 | 0.25 | 0.166667 | + |
| Est | 4.95271 | 0.250001 | 0.166668 | x |
4.3.3 Separate, two seroconverted obs: \(\beta_1 = 2/4\), \(\beta_2 = 2/6\)
Convergence to inside the range of plausible correct values.
| y.seroconversion | n.1 | y.1 | n.2 | y.2 |
|---|---|---|---|---|
| 0 | 0 | 0 | 4 | 0 |
| 1 | 0 | 0 | 2 | 1 |
| 0 | 2 | 0 | 0 | 0 |
| 1 | 2 | 1 | 0 | 0 |
optim.log.likelihood.beta.dat <- f.est.betas.sim(N.types, probs.seroconversion, dat.clean, bound.lower = 1e-10)
f.plot.2D(N.types, probs.seroconversion, dat.clean, optim.log.likelihood.beta.dat, f.name)
f.name <- "Separate, two obs: 2/4, 2/6"
probs.seroconversion2 <- c(2/4, 2/6) # upper bound
optim.log.likelihood.beta.dat2 <- f.est.betas.sim(N.types, probs.seroconversion2, dat.clean, bound.lower = 1e-10)
points(probs.seroconversion2[1], probs.seroconversion2[2], pch = 3, cex = 3, col = "red")
#f.plot.2D(N.types, probs.seroconversion2, dat.clean, optim.log.likelihood.beta.dat2, f.name)Model results: CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH
| Type | \(-\log(L)\) | \(\beta_1\) | \(\beta_2\) | plot symbol |
|---|---|---|---|---|
| True | 3.31695 | 0.25 | 0.166667 | + |
| Est | 3.29584 | 0.292894 | 0.183505 | x |
| True 2 | 3.88362 | 0.5 | 0.333333 | + |
4.3.4 Separate, four seroconverted obs: \(\beta_1 = 4/6\), \(\beta_2 = 4/8\)
Convergence to inside the range of plausible correct values.
| y.seroconversion | n.1 | y.1 | n.2 | y.2 |
|---|---|---|---|---|
| 0 | 0 | 0 | 4 | 0 |
| 1 | 0 | 0 | 4 | 1 |
| 0 | 2 | 0 | 0 | 0 |
| 1 | 4 | 1 | 0 | 0 |
optim.log.likelihood.beta.dat <- f.est.betas.sim(N.types, probs.seroconversion, dat.clean, bound.lower = 1e-10)
f.plot.2D(N.types, probs.seroconversion, dat.clean, optim.log.likelihood.beta.dat, f.name)
f.name <- "Separate, two obs: 4/6, 4/8"
probs.seroconversion2 <- c(4/6, 4/8) # upper bound
optim.log.likelihood.beta.dat2 <- f.est.betas.sim(N.types, probs.seroconversion2, dat.clean, bound.lower = 1e-10)
points(probs.seroconversion2[1], probs.seroconversion2[2], pch = 3, cex = 3, col = "red")
#f.plot.2D(N.types, probs.seroconversion2, dat.clean, optim.log.likelihood.beta.dat2, f.name)Model results: CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL
| Type | \(-\log(L)\) | \(\beta_1\) | \(\beta_2\) | plot symbol |
|---|---|---|---|---|
| True | 2.43937 | 0.166667 | 0.125 | + |
| Est | 2.34107 | 0.240165 | 0.159105 | x |
| True 2 | 5.04677 | 0.666667 | 0.5 | + |
4.4 Simulations: two types with overlapping exposures
In the following examples, when there is are overlapping exposures at seroconversion the \(\beta\) MLE estimates the true parameters {} bias. There are conditions (not shown) where the sample sizes are very large and the estimates are brought back off the boundary.
A future area of work is to use a Bayesian method to include prior information that may help this boundary issue.
4.4.1 Overlapping, one seroconverted obs: \(\beta_1 = 1/3\), \(\beta_2 = 1/5\)
If both exposures have only 1 exposure on seroconversion, then the negative log likelihood function has a minimum on the boundary, sending the exposure with more observations to probability equal to 0.
| y.seroconversion | n.1 | y.1 | n.2 | y.2 |
|---|---|---|---|---|
| 0 | 0 | 0 | 3 | 0 |
| 1 | 1 | 0 | 1 | 1 |
| 0 | 1 | 0 | 0 | 0 |
| 1 | 1 | 1 | 1 | 0 |
optim.log.likelihood.beta.dat <- f.est.betas.sim(N.types, probs.seroconversion, dat.clean, bound.lower = 1e-10)
f.plot.2D(N.types, probs.seroconversion, dat.clean, optim.log.likelihood.beta.dat, f.name)
Model results: CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL
| Type | \(-\log(L)\) | \(\beta_1\) | \(\beta_2\) | plot symbol |
|---|---|---|---|---|
| True | 2.59918 | 0.333333 | 0.2 | + |
| Est | 1.90954 | 0.666666 | 10^{-10} | x |
4.4.2 Overlapping on nonserconversion, one seroconverted obs: \(\beta_1 = 1/5\), \(\beta_2 = 1/7\)
If both exposures have only 1 exposure on seroconversion, but overlapping events on nonseroconversion, then there is no additional issue. The estimates are similar to as when there was no overlap.
| y.seroconversion | n.1 | y.1 | n.2 | y.2 |
|---|---|---|---|---|
| 0 | 1 | 0 | 4 | 0 |
| 1 | 0 | 0 | 2 | 1 |
| 0 | 2 | 0 | 1 | 0 |
| 1 | 2 | 1 | 0 | 0 |
optim.log.likelihood.beta.dat <- f.est.betas.sim(N.types, probs.seroconversion, dat.clean, bound.lower = 1e-10)
f.plot.2D(N.types, probs.seroconversion, dat.clean, optim.log.likelihood.beta.dat, f.name)
f.name <- "Overlapping both on non-seroconversion, two obs: 2/5, 2/7"
f.name <- "2/5, 2/7"
probs.seroconversion2 <- c(2/5, 2/7) # upper bound
optim.log.likelihood.beta.dat2 <- f.est.betas.sim(N.types, probs.seroconversion2, dat.clean, bound.lower = 1e-10)
points(probs.seroconversion2[1], probs.seroconversion2[2], pch = 3, cex = 3, col = "red")
#f.plot.2D(N.types, probs.seroconversion2, dat.clean, optim.log.likelihood.beta.dat2, f.name)Model results: CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL
| Type | \(-\log(L)\) | \(\beta_1\) | \(\beta_2\) | plot symbol |
|---|---|---|---|---|
| True | 3.78871 | 0.2 | 0.142857 | + |
| Est | 3.77647 | 0.225404 | 0.154848 | x |
| True 2 | 4.37489 | 0.4 | 0.285714 | + |
4.4.3 Overlapping nonseroconversion with both seroconvertion obs: \(\beta_1 = 1/5\), \(\beta_2 = 1/7\)
Similar boundary issue when multiple exposure types overlap with one or more exposures.
| y.seroconversion | n.1 | y.1 | n.2 | y.2 |
|---|---|---|---|---|
| 0 | 0 | 0 | 4 | 0 |
| 1 | 1 | 0 | 2 | 1 |
| 0 | 2 | 0 | 0 | 0 |
| 1 | 2 | 1 | 1 | 0 |
optim.log.likelihood.beta.dat <- f.est.betas.sim(N.types, probs.seroconversion, dat.clean, bound.lower = 1e-10)
f.plot.2D(N.types, probs.seroconversion, dat.clean, optim.log.likelihood.beta.dat, f.name)
#f.name <- "Overlapping both on non-seroconversion, two obs: 2/5, 2/7"
f.name <- "2/5, 2/7"
probs.seroconversion2 <- c(2/5, 2/7) # upper bound
optim.log.likelihood.beta.dat2 <- f.est.betas.sim(N.types, probs.seroconversion2, dat.clean, bound.lower = 1e-10)
points(probs.seroconversion2[1], probs.seroconversion2[2], pch = 3, cex = 3, col = "red")
#f.plot.2D(N.types, probs.seroconversion2, dat.clean, optim.log.likelihood.beta.dat2, f.name)Model results: CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL
| Type | \(-\log(L)\) | \(\beta_1\) | \(\beta_2\) | plot symbol |
|---|---|---|---|---|
| True | 2.74437 | 0.2 | 0.142857 | + |
| Est | 2.35365 | 0.459688 | 10^{-10} | x |
| True 2 | 3.03025 | 0.4 | 0.285714 | + |
4.4.4 Overlapping n.1 on n.2 on seroconvertion obs: \(\beta_1 = 1/5\), \(\beta_2 = 1/6\)
n.1 predicted when n.2 caused, but the model is blind to this.
| y.seroconversion | n.1 | y.1 | n.2 | y.2 |
|---|---|---|---|---|
| 0 | 0 | 0 | 4 | 0 |
| 1 | 1 | 0 | 2 | 1 |
| 0 | 2 | 0 | 0 | 0 |
| 1 | 2 | 1 | 0 | 0 |
optim.log.likelihood.beta.dat <- f.est.betas.sim(N.types, probs.seroconversion, dat.clean, bound.lower = 1e-10)
f.plot.2D(N.types, probs.seroconversion, dat.clean, optim.log.likelihood.beta.dat, f.name)
f.name <- "Overlapping n1 on n2 on seroconversion, two obs: 2/5, 2/6"
probs.seroconversion2 <- c(2/5, 2/6) # upper bound
optim.log.likelihood.beta.dat2 <- f.est.betas.sim(N.types, probs.seroconversion2, dat.clean, bound.lower = 1e-10)
points(probs.seroconversion2[1], probs.seroconversion2[2], pch = 3, cex = 3, col = "red")
#f.plot.2D(N.types, probs.seroconversion2, dat.clean, optim.log.likelihood.beta.dat2, f.name)Model results: CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL
| Type | \(-\log(L)\) | \(\beta_1\) | \(\beta_2\) | plot symbol |
|---|---|---|---|---|
| True | 3.00815 | 0.2 | 0.166667 | + |
| Est | 2.35365 | 0.459688 | 10^{-10} | x |
| True 2 | 3.39995 | 0.4 | 0.333333 | + |
4.5 Simulation: data-realistic case
The data-realistic simulation strategy is to choose \(N\) subjects to study and survey each subject, \(i=1,\ldots,N\), for \(S_i\) surveys, where each subject has a common per-contact rate of seroconversion. Assume that each subject is surveyed at perfect 90-day intervals, so “stretching” is not necessary for the simulations, and that subjects report perfectly on their number of exposures. Let the maximum number of surveys for each subject be 15 and choose a rate of observations per subject \(\lambda\) (e.g., 3). Simulate the number of surveys iid for each subject as \(S_i \sim \textrm{Poisson}(\lambda)\), and let \(s=1,\ldots,S_i\). For all subjects, choose a number of exposure types, \(E\), and let \(j=1,\ldots,E\), to study (e.g., \(E=1\) for only needles, \(E=2\) for needles and equipment). For each exposure, choose a 30-day exposure rate \(r\) (e.g., 20 and 10). Simulate the number of exposures per survey for each subject as \(n_{isj} \sim \textrm{Poisson}(r_j)\). For each exposure, choose a per-contact seroconversion probability \(\beta\) (e.g., 0.004 and 0.0004). Simulate the seroconversion due to contact for each subject as \(y_{isj} \sim \textrm{Binomial}(n_{isj}, \beta_j)\). Finally, “clean” the data by excluding leading positive seroconversion surveys and any surveys after the first seroconversion.
# function to simulate and clean data
f.sim.dat <- function(N.users
, N.types
, probs.seroconversion
, rate.exposure.30.days
, rate.number.of.observations.per.subject
, n.max.surveys) {
dat <- data.frame(id = sort(rep(1:N.users, n.max.surveys))
, survey.number = rep(1:n.max.surveys, N.users)
, y.seroconversion = NA
)
for (i.N.types in 1:N.types) {
# exposures
dat[,paste("n.", i.N.types, sep="")] <- rpois(N.users * n.max.surveys, rate.exposure.30.days[i.N.types])
# sero for this exposure
dat[,paste("y.", i.N.types, sep="")] <- NA
}
#str(dat)
# increase to 90-day exposures by multiplying by 3
dat[,paste("n.", i.N.types, sep="")] <- 3 * dat[,paste("n.", i.N.types, sep="")]
# determine serostatus
for (i.obs in 1:nrow(dat)) {
for (i.N.types in 1:N.types) {
dat[i.obs, paste("y.", i.N.types, sep="")] <- rbinom(1, size = dat[i.obs, paste("n.", i.N.types, sep="")], probs.seroconversion[i.N.types])
}
}
# seroconversion by any exposure type
dat$y.seroconversion <- as.numeric(apply( as.matrix(dat[, grep("^y.", colnames(dat))[-1]]), 1, sum) > 0)
#str(dat)
# reduce number of observations per subject
library(plyr)
dat <- ddply(dat
, .(id)
, function (.X) {
n.obs <- min(n.max.surveys, rpois(1, rate.number.of.observations.per.subject) + 1)
.X <- .X[1:n.obs, ]
return(.X)
}
)
#str(dat)
# remove leading 1s and values after first 1
# there are 3 possible methods for doing this
dat.clean <- ddply(dat, "id",
function(.X){
#print(.X$id)
# Step 0, remove NAs
.X <- na.omit(.X)
#out <- data.frame(id=.X$id[1], intdate=.X$intdate[1], n.exposures.syringe=NA, n.exposures.equip=NA, n.exposures.back=NA, y.seroconversion=NA)
ind.1to0 <- which(diff(.X$y.seroconversion) == -1)
# Step 1, remove leading sero = 1
if (.X$y.seroconversion[1] == 1) {
if (length(ind.1to0) == 0) {
last.y1 <- nrow(.X) # then all sero = 1
} else {
last.y1 <- ind.1to0[1] # last sero = 1 before going to 0
}
# remove leading 1s
.X <- .X[-c(1:last.y1), ]
}
# no observations left (they were all 1s)
if (nrow(.X) == 0) {
return(NULL)
}
# Step 2, remove all values after 1s sero = 1
ind.0to1 <- which(diff(.X$y.seroconversion) == +1)
if (length(ind.0to1) > 0) {
.X <- .X[1:(ind.0to1[1] + 1), ]
}
## Step 3, return exposures
# select which aggregation method
e.method <- 2
## ## method 1 (from Erik's 2/15/2015 notes)
## # aggregate over all measurement periods
## if (e.method == 1) {
## out <- .X[nrow(.X), ]
## out[, grep("^n.", colnames(dat))] <- colSums(.X[, grep("^n.", colnames(dat))])
## return(out)
## }
## method 2 (from Erik's 2/15/2015 notes)
# separate each measurement period
if (e.method == 2) {
return(.X)
}
## ## method 3 (from Erik's 2/15/2015 notes)
## # partition by status: aggregate 0s, separate last 1 (if exists)
## if (e.method == 3) {
## ind.0to1 <- which(diff(.X$y.seroconversion) == +1)
## if (length(ind.0to1) == 0) {
## out <- .X[nrow(.X), ]
## out[, grep("^n.", colnames(dat))] <- colSums(.X[, grep("^n.", colnames(dat))])
## return(out)
## }
## if (length(ind.0to1) > 0) {
## out0 <- .X[1, ]
## out0[, grep("^n.", colnames(dat))] <- colSums(.X[1:(nrow(.X)-1), grep("^n.", colnames(dat))])
## out1 <- .X[nrow(.X), ]
## out <- rbind(out0, out1)
## return(out)
## }
##
## return(.X)
## }
}
)
#str(dat.clean)
} # f.sim.dat4.5.1 Simulation: Real-data: \(\beta_1 = 0.004\), \(\beta_2 = 0.0004\)
################################################################################
# 12/8/2014
# 1/15/2015
# 2/15/2015, multiple observation windows (surveys) per subject
# path <- "C:/Dropbox/UNM/advise/2015_YuridiaLeyva_MS_KimPage/ErikYuridia_shared/sim"
# setwd(path)
## Simulation settings
# number of users
N.users <- 500
# number of exposure types
N.types <- 2
# true probability of seroconversion (to estimate this quantity)
probs.seroconversion <- c(0.004, 0.0004, 0.01, 0.01)[1:N.types]
# average number of exposures in 30 days (arbitrary time period, not used in calculation)
rate.exposure.30.days <- c(20, 10, 4, 4)[1:N.types] # 4 is roughly once a week
rate.number.of.observations.per.subject <- 3 # poisson rv
## Simulate data
# create an empty data.frame to hold data
# Allow each subject to have up to n.max.surveys surveys
# then determine a number of surveys for each subject and remove the extras
n.max.surveys <- 10 # maximum number of surveys for each subject
dat.clean <- f.sim.dat(N.users
, N.types
, probs.seroconversion
, rate.exposure.30.days
, rate.number.of.observations.per.subject
, n.max.surveys)Data were simulated with the following parameters.
| Name | Parameter | value |
|---|---|---|
| Users | \(N\) | 500 |
| Exposure types | \(E\) | 2 |
| Seroconversion probabilities | \(\beta\) | 0.004, 410^{-4} |
| Number of surveys for each subject | \(\lambda\) | 3 |
| 30-day Poisson exposure rate | \(r\) | 20, 10 |
| Max number of surveys | 10 | |
| Resulting data: | ||
| Number of Users | 500 | |
| Number of surveys | 2000 |
optim.log.likelihood.beta.dat <- f.est.betas.sim(N.types, probs.seroconversion, dat.clean, bound.lower = 1e-10)
write.csv(optim.log.likelihood.beta.dat$par, "sim1_optim_beta.csv", row.names = FALSE)R <- 100
optim.log.likelihood.beta.bs <- f.est.betas.bs(dat.clean, R = R) 1 2 3 4 5 6 7 8 9 10
11 12 13 14 15 16 17 18 19 20
21 22 23 24 25 26 27 28 29 30
31 32 33 34 35 36 37 38 39 40
41 42 43 44 45 46 47 48 49 50
51 52 53 54 55 56 57 58 59 60
61 62 63 64 65 66 67 68 69 70
71 72 73 74 75 76 77 78 79 80
81 82 83 84 85 86 87 88 89 90
91 92 93 94 95 96 97 98 99 100write.csv(optim.log.likelihood.beta.bs, "sim1_optim_beta_CI_bs.csv", row.names = FALSE)
# sort for CI
optim.log.likelihood.beta.bs.sorted <- optim.log.likelihood.beta.bs[,grep("^beta", colnames(optim.log.likelihood.beta.bs))]
for (i.col in 1:ncol(optim.log.likelihood.beta.bs.sorted)) {
optim.log.likelihood.beta.bs.sorted[,i.col] <- sort(optim.log.likelihood.beta.bs.sorted[,i.col])
}
CI.bs <- optim.log.likelihood.beta.bs.sorted[c(max(1, round(0.025*R)), min(R, round(0.975*R+1))), ]f.name <- "Full data sim: 0.004, 0.0004"
f.plot.2D(N.types, probs.seroconversion, dat.clean, optim.log.likelihood.beta.dat, f.name, seq.length = 21)
f.plot.2D(N.types, probs.seroconversion, dat.clean, optim.log.likelihood.beta.dat, f.name, seq.length = 21, seq.limits = c(1e-20, 0.01))
points(optim.log.likelihood.beta.bs[,1], optim.log.likelihood.beta.bs[,2], pch = ".")
Model results: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
| Type | \(-\log(L)\) | \(\beta_1\) | \(\beta_2\) | plot symbol |
|---|---|---|---|---|
| True | 432.135 | 0.004 | 410^{-4} | + |
| Est | 425.94 | 0.0035451 | 10^{-10} | x |
| 95% CI | lower | 0.0013963 | 10^{-10} | . |
| 95% CI | upper | 0.0041275 | 9.8899410^{-4} | . |
4.5.2 Simulation: Real-data: \(\beta_1 = 0.02\), \(\beta_2 = 0.005\)
################################################################################
# 12/8/2014
# 1/15/2015
# 2/15/2015, multiple observation windows (surveys) per subject
# path <- "C:/Dropbox/UNM/advise/2015_YuridiaLeyva_MS_KimPage/ErikYuridia_shared/sim"
# setwd(path)
## Simulation settings
# number of users
N.users <- 500
# number of exposure types
N.types <- 2
# true probability of seroconversion (to estimate this quantity)
probs.seroconversion <- c(0.02, 0.005, 0.01, 0.01)[1:N.types]
# average number of exposures in 30 days (arbitrary time period, not used in calculation)
rate.exposure.30.days <- c(10, 5, 4, 4)[1:N.types] # 4 is roughly once a week
rate.number.of.observations.per.subject <- 4 # poisson rv
## Simulate data
# create an empty data.frame to hold data
# Allow each subject to have up to n.max.surveys surveys
# then determine a number of surveys for each subject and remove the extras
n.max.surveys <- 15 # maximum number of surveys for each subject
dat.clean <- f.sim.dat(N.users
, N.types
, probs.seroconversion
, rate.exposure.30.days
, rate.number.of.observations.per.subject
, n.max.surveys)Data were simulated with the following parameters.
| Name | Parameter | value |
|---|---|---|
| Users | \(N\) | 500 |
| Exposure types | \(E\) | 2 |
| Seroconversion probabilities | \(\beta\) | 0.02, 0.005 |
| Number of surveys for each subject | \(\lambda\) | 4 |
| 30-day Poisson exposure rate | \(r\) | 10, 5 |
| Max number of surveys | 20 | |
| Resulting data: | ||
| Number of Users | 500 | |
| Number of surveys | 2000 |
optim.log.likelihood.beta.dat <- f.est.betas.sim(N.types, probs.seroconversion, dat.clean, bound.lower = 1e-10)
write.csv(optim.log.likelihood.beta.dat$par, "sim2_optim_beta.csv", row.names = FALSE)R <- 100
optim.log.likelihood.beta.bs <- f.est.betas.bs(dat.clean, R = R) 1 2 3 4 5 6 7 8 9 10
11 12 13 14 15 16 17 18 19 20
21 22 23 24 25 26 27 28 29 30
31 32 33 34 35 36 37 38 39 40
41 42 43 44 45 46 47 48 49 50
51 52 53 54 55 56 57 58 59 60
61 62 63 64 65 66 67 68 69 70
71 72 73 74 75 76 77 78 79 80
81 82 83 84 85 86 87 88 89 90
91 92 93 94 95 96 97 98 99 100write.csv(optim.log.likelihood.beta.bs, "sim2_optim_beta_CI_bs.csv", row.names = FALSE)
# sort for CI
optim.log.likelihood.beta.bs.sorted <- optim.log.likelihood.beta.bs[,grep("^beta", colnames(optim.log.likelihood.beta.bs))]
for (i.col in 1:ncol(optim.log.likelihood.beta.bs.sorted)) {
optim.log.likelihood.beta.bs.sorted[,i.col] <- sort(optim.log.likelihood.beta.bs.sorted[,i.col])
}
CI.bs <- optim.log.likelihood.beta.bs.sorted[c(max(1, round(0.025*R)), min(R, round(0.975*R+1))), ]f.name <- "Full data sim: 0.004, 0.0004"
f.plot.2D(N.types, probs.seroconversion, dat.clean, optim.log.likelihood.beta.dat, f.name, seq.length = 21)
f.plot.2D(N.types, probs.seroconversion, dat.clean, optim.log.likelihood.beta.dat, f.name, seq.length = 21, seq.limits = c(1e-20, 0.01))
points(optim.log.likelihood.beta.bs[,1], optim.log.likelihood.beta.bs[,2], pch = ".")
Model results: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
| Type | \(-\log(L)\) | \(\beta_1\) | \(\beta_2\) | plot symbol |
|---|---|---|---|---|
| True | 787.014 | 0.02 | 0.005 | + |
| Est | 772.004 | 0.0164325 | 0.0026726 | x |
| 95% CI | lower | 0.011127 | 10^{-10} | . |
| 95% CI | upper | 0.0206002 | 0.005463 | . |
Bibliography
Boelen, Lies, Suzy Teutsch, David P Wilson, Kate Dolan, Greg J Dore, Andrew R Lloyd, Fabio Luciani, HITS investigators, and others. 2014. “Per-Event Probability of Hepatitis C Infection During Sharing of Injecting Equipment.” PloS One 9 (7). Public Library of Science: e100749.
Hughes, James P, Jared M Baeten, Jairam R Lingappa, Amalia S Magaret, Anna Wald, Guy De Bruyn, James Kiarie, et al. 2012. “Determinants of Per-Coital-Act Hiv-1 Infectivity Among African Hiv-1–serodiscordant Couples.” Journal of Infectious Diseases 205 (3). Oxford University Press: 358–65.
Jewell, Nicholas P, and Stephen C Shiboski. 1990. “Statistical Analysis of Hiv Infectivity Based on Partner Studies.” Biometrics. JSTOR, 1133–50.
Magaret, Amalia S, Andrew Mujugira, James P Hughes, Jairam Lingappa, Elizabeth A Bukusi, Guy DeBruyn, Sinead Delany-Moretlwe, et al. 2015. “Effect of Condom Use on Per-Act Hsv-2 Transmission Risk in Hiv-1, Hsv-2-Discordant Couples.” Clinical Infectious Diseases 62 (4). Oxford University Press: 456–61.