chords (version 0.95.4)

Estimate.b.k: RDS population size estimation

Description

Estimate population size from respondent driven samples (RDS) using maximum likelihood, and several variation. The underlying idea is that the sample spreads like an epidemic in the target population as described in the reference.

Usage

Estimate.b.k(rds.object, type = "mle", jack.control = NULL)

Arguments

rds.object

A object of class rds-object as constructed by initializeRdsObject or outputted by Estimate.b.k (depending on the type used).

type

A character vector with the type of estimation. Possible values:

• mle Maximum likelihood.
• integrated Integrated maximum likelihood.
• observed Estimate with observed degrees.
• jeffreys MAP estimation with Jeffreys prior.
• parametric Assume $\beta[k]:=\beta * \theta^k$.
• rescaling Naive rescaling heuristic estimation.
• leave-d-out Leave-d-out resampling estimator.

jack.control

A object of class jack.control as constructed by makeJackControl.

Value

An rds-class object with an updated estimates slot. The estiamtes slot is list with the following components:

Details

As of version 0.95, this function is the main workhorse of the chords package. Given an rds-class object, it will return population size estimates for each degree. Note that for the rescaling and parametric estimators, the input rds-object is expected to contain some initial estimate in the estimates slot.

See the reference for a description of the likelihood problem solved. Optimization is performed by noting that likelihood is coordinate-wise convex, thus amounts to a series of line-searches.

References

[1] Berchenko, Y., Rosenblatt D.J., and S.D.W. Frost. "Modeling and Analyzing Respondent Driven Sampling as a Counting Process." arXiv:1304.3505,

See Also

initializeRdsObject, makeRdsSample, getTheta.

Examples

# Preliminaries
data(brazil)
rds.object2<- initializeRdsObject(brazil)
see <- function(x) plot(x$estimates$Nk.estimates, type='h')

# Maximum likelihood
rds.object <- Estimate.b.k(rds.object = rds.object2 )
see(rds.object)

# View estimates:
plot(rds.object$estimates$Nk.estimates, type='h')
# Population size estimate:
sum(rds.object$estimates$Nk.estimates)
plot(rds.object$estimates$log.bk.estimates, type='h')

## Recover theta assuming b.k=b_0*k^theta
getTheta(rds.object)

# How many degrees were imputed?:
table(rds.object$estimates$convergence)


# Observed degree estimation:
rds.object.4 <- Estimate.b.k(rds.object = rds.object, type='observed')
see(rds.object.4)

# Naive rescaling
rds.object.5 <- Estimate.b.k(rds.object = rds.object, type='rescaling')
see(rds.object.5)

# Parametric rates
rds.object.6 <- Estimate.b.k(rds.object = rds.object, 
                             type='parametric')
see(rds.object.6)
jack.control <- makeJackControl(3, 1e1)
rds.object.7 <- Estimate.b.k(rds.object = rds.object, 
                             type='leave-d-out', 
                             jack.control = jack.control)
see(rds.object.7)
rds.object.8 <- Estimate.b.k(rds.object = rds.object, 
                             type='integrated', 
                             jack.control = jack.control)
see(rds.object.8)
rds.object.9 <- Estimate.b.k(rds.object = rds.object, 
                             type='jeffreys')
see(rds.object.9)



## Not run: 
# ## Simulated data example:
#     dk <- c(2, 1e1) # unique degree classes
#     true.dks <- rep(0,max(dk)); true.dks[dk] <- dk
#     true.Nks <- rep(0,max(dk)); true.Nks[dk] <- 1e3
#     beta <- 1 #5e-6
#     theta <-  0.1
#     true.log.bks <- rep(-Inf, max(dk))
#     true.log.bks[dk] <- theta*log(beta*dk)
#     sample.length <- 4e2
#     nsims <- 1e2
#     
#     simlist <- list()
#     for(i in 1:nsims){
#       simlist[[i]] <- makeRdsSample(
#         N.k =true.Nks , 
#         b.k = exp(true.log.bks),
#         sample.length = sample.length)
#     }
#     
#     
#     # Estimate betas and theta with chords:
#     llvec <- rep(NA,nsims)
#     bklist <- list()
#     for(i in 1:nsims){
#       # i <- 2
#       simlist[[i]] <- Estimate.b.k(rds.object = simlist[[i]])
#       # llvec[i] <- simlist[[i]]$estimates$likelihood
#       bklist[[i]] <- simlist[[i]]$estimates$log.bk.estimates
#     }
#     b1vec <- bklist 
#     b2vec <- bklist 
#     
#     hist(b1vec)
#     abline(v=true.log.bks[2])
#     hist(b2vec)
#     abline(v=true.log.bks[10])
#     
#     beta0vec <- rep(-Inf,nsims)
#     thetavec <- rep(-Inf,nsims)
#     nvec <- rep(-Inf,nsims)
#     converged <- rep(9999,nsims)
#     
#     for(i in 1:nsims){
#       # i <- 2
#       nvec[i] <- sum(simlist[[i]]$estimates$Nk.estimates)
#       converged[i] <- sum(simlist[[i]]$estimates$convergence, na.rm=TRUE)
#       # tfit <- getTheta(simlist[[i]])
#       # beta0vec[i] <- tfit$log.beta_0
#       # thetavec[i] <- tfit$theta
#     }
#     summary(beta0vec)
#     summary(nvec)
#     # summary(thetavec)
#     # hist(thetavec)
#     # abline(v=theta)
#     hist(nvec)
#     abline(v=sum(true.Nks), col='red')
#     abline(v=median(nvec, na.rm = TRUE), lty=2)
#     table(converged)
#     
#     # Try various re-estimatinons:
#     rds.object2 <- simlist[[which(is.infinite(nvec))[1]]]
#     
#     rds.object <- Estimate.b.k(rds.object = rds.object2 )
#     see(rds.object)
#     rds.object$estimates$Nk.estimates
#     
#     rds.object.5 <- Estimate.b.k(rds.object = rds.object, type='rescaling')
#     see(rds.object.5) # will not work. less than 2 converging estimates.
#     rds.object.5$estimates$Nk.estimates
#     
#     rds.object.6 <- Estimate.b.k(rds.object = rds.object, type='parametric')
#     see(rds.object.6) # will not work. less than 2 converging estimates.
#     
#     
#     jack.control <- makeJackControl(3, 1e2)
#     rds.object.7 <- Estimate.b.k(rds.object = rds.object, 
#                              type='leave-d-out', 
#                              jack.control = jack.control)
#     see(rds.object.7)
#     rds.object.7$estimates$Nk.estimates
#     
#     
#     rds.object.8 <- Estimate.b.k(rds.object = rds.object, type='integrated')
#     see(rds.object.8)
#     rds.object.8$estimates$Nk.estimates
#     
#     
#     rds.object.9 <- Estimate.b.k(rds.object = rds.object, type='jeffreys')
#     see(rds.object.9)
#     rds.object.9$estimates$Nk.estimates
# ## End(Not run)