# R/gile-ss-util.R In RDS: Respondent-Driven Sampling

#### Defines functions getinclCstackedgetestCstackedprobtodist

```# turn sampling probabilities and sample counts into a synthetic population
probtodist<-function(classes,nums,prob,n,hajek=TRUE){
#classes: size value of each class
#nums: number of members of that class in the sample
#prob: Probabilities of selection for a member of that class
props<-rep(0,length(classes))
for(k in 1:length(classes)){
props[k]<-nums[k]/prob[k]
}
# Hajek is the one in Krista's thesis and it implicitly multiplies the
# HT estimate of the proportions by N as an estimate of the population size.
# Otherwise use the HT estimator which implicitly multiples by the HT
# estimate of the population size)
#
# The Hajek estimator is significantly less efficient than the HT estimator
# under PPSWOR when the outcome is approximately proportional to degree.
# When the disease status is unrelated to the degree the Hajek is slightly
# more efficient (all assuming that the population size is actually N)
#
# In this case we are estimating the counts in each class not the disease
# outcome
#
if(!hajek){
n<-round(sum(props))
}
props<-props/sum(props)
if(any(is.na(props))){
print("Error in proptodist")
}
#classes: size value of each class:
#props: proportion of the population in that class
list(classes=classes, props=props, n=n)
}

getestCstacked<-function(samp,n,nit=5, nsampsamp=1000, trace=FALSE, hajek=TRUE,SS.infinity=0.04){
# samp: sample
# n: population size
classes<-sort(unique(samp))
nums<-table(samp)
nsamp<-length(samp)
prob<-classes/sum(samp)
prob <- prob * sum(nums/prob) / n

#temp\$classes: size value of each class
#temp\$nums: number of members of that class
#temp\$prob: Probabilities of selection for a member of that class
temp<-probtodist(classes=classes,nums=nums,prob=prob,n=n,hajek=hajek)
#temp\$classes: size value of each class
#temp\$props: proportion of the population in that class
if(n*SS.infinity < nsamp){
newprobs<-getinclCstacked(classes=temp\$classes,props=temp\$props,n=n,
nsamp=nsamp,
nsampsamp=nsampsamp,
nums=nums)
#newprobs\$degvec: size value of each class
#newprobs\$pvec: Probabilities of selection for a member of that class

for(i in 1:nit){
if(!hajek){
cat(sprintf("The estimated population size is %d.\n",temp\$n))
}
temp<-probtodist(classes=newprobs\$degvec,nums=nums,prob=newprobs\$pvec,
n=temp\$n,hajek=hajek)
newprobs<-getinclCstacked(classes=temp\$classes,props=temp\$props,
n=temp\$n,
nsamp=nsamp,
nsampsamp=nsampsamp,
nums=nums)
}
}else{
if(nsamp < n){
#    VH
pvec=nsamp*(temp\$classes)/(temp\$n*sum(temp\$props*temp\$classes))
}else{
#    Arithmetic Mean
pvec=rep(nsamp,length=temp\$classes)/temp\$n
}
newprobs <- list(degvec=temp\$classes, pvec=pvec,
nbyclass=round(temp\$n*temp\$props/sum(temp\$props)))
}
temp<-probtodist(classes=newprobs\$degvec,nums=nums,prob=newprobs\$pvec,
n=temp\$n,hajek=hajek)
if(!hajek){
cat(sprintf("The estimated population size is %d.\n",temp\$n))
}
#classes: size value of each class
#props: proportion of the population in that class
#probs: Probabilities of selection for a member of that class
list(probs=newprobs\$pvec,props=temp\$props,classes=temp\$classes,n=temp\$n)
}

getinclCstacked<-function(classes,props,n,nsamp,nsampsamp,nums,hajek=TRUE){
props2<-props/sum(props)
#artificially inflate each value to at least 1 node
nbyclass<-round(n*props2)
nbyclass[nbyclass==0]<-1
offby<-n-sum(nbyclass)
if(is.na(offby)){
print("Error in getincl: offby")
}
# based on estimate of degrees, estimate true inclusion probs by class
Cret <- .C("getinclCstacked",
nbyclass=as.integer(nbyclass),
size=as.double(classes),
props2=as.double(props2),
offby=as.integer(offby),
N=as.integer(n),
K=as.integer(length(classes)),
n=as.integer(nsamp),
samplesize=as.integer(nsampsamp),
Nk=as.integer(classes),
PACKAGE="RDS")
Ninf <- (Cret\$Nk+nums)/(sum(Cret\$Nk)+nsamp)
Ninf <- nsamp*Ninf/Cret\$nbyclass
list(degvec=classes,pvec=Ninf,nbyclass=Cret\$nbyclass)
}
```

## Try the RDS package in your browser

Any scripts or data that you put into this service are public.

RDS documentation built on Dec. 2, 2017, 1:08 a.m.