# RDS population size estimation 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. ```r 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: * mleMaximum likelihood. * integratedIntegrated maximum likelihood. * observedEstimate with observed degrees. * jeffreysMAP estimation with Jeffreys prior. * parametricAssume $\beta[k]:=\beta * \theta^k$. * rescalingNaive rescaling heuristic estimation. * leave-d-outLeave-d-out resampling estimator. - `jack.control`: A object of class `jack.control` as constructed by `makeJackControl.` ## 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. ## Returns An `rds-class` object with an updated `estimates` slot. The `estiamtes` slot is `list` with the following components: - **call**: The function call. - **Nk.estimates**: The estimated degree frequencies. - **log.bk.estimates**: The estimated sampling rates for each degree. In log scale. - **convergence**: 0 if estimation of $N[k]$'s converged. Otherwise, 1 or -1, depending on the sign of the score function at the MLE. ## 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 ```r # 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) ```

Estimate.b.k function