受访者驱动抽样(RDS)的相关R程序

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

受访者驱动抽样（RDS）是一种专门针对隐匿人群的抽样方法，常用于跨性别者、暗娼、MSM等因耻辱感和法律制度的约束难以识别和接触的人群，并逐渐应用于一般人群。这里给出国外学者写的一个相关R程序：
基本参数设定：
n <- 20000 #pop size
ns <- 500 # rds sample size
seeds <- 7 # #of seeds
#d <- round(rlnorm(n,meanlog = 1)) + 1
d <- round(rexp(n,1/5)) + 1 #heavy tailed degrees
#d <- rpois(n,lambda = 5) + 1 #light tail
#d <- rep(10, n)
g <- rep(1,n)

首先构造邻接矩阵的函数：
makeGraph <- function(d){
n <- length(d)
edges <- list()
d1 <- d
while(sum(d1)>0.5){
  l <- sample.int(n,1,prob=d1)
  dt <- d1
  dt[l] <- d1[l] - 1
  if(sum(dt) > 0.5){
   k <- sample.int(n,1,prob=dt)
   v <- TRUE
   d1[k] <- d1[k] - 1
  }else{
   k <- sample.int(n,1,prob=d)
   v <- FALSE
  }
  d1[l] <- d1[l] - 1
  edges[[length(edges)+1]] <- c(l,k,v)
}
el <- do.call(rbind,edges)
el
}
el <- makeGraph(d)   #el即为相关联的矩阵

构造RDS抽样函数：
sampRDS <- function(el, d, ns, g, ss, biased=TRUE,pr = c(0,.1,.9)){
maxR <- length(pr) -1
n <- length(g)
ml <- if(is.factor(g)) max(levels(g)) else max(g)
if(biased){
  seeds <- sample.int(n,ns,prob=as.numeric(as.factor(g))-1)
}else{
  seeds <- sample.int(n,ns, prob = d)
}
samp <- seeds
recr <- rep(-1,ns)
time <- 0 + (1:ns)/10000000
rcTime <- rexp(ns)
v <- rep(-1,ns)
t1 <- time
while(length(samp) < ss){
  subjIndex <- which.min(t1 + rcTime)
  if(length(subjIndex)==0){
   print("redraw")
   subj <- sample( (1:n)[-samp],1)
   samp <- c(samp,subj)
   recr <- c(recr,-1)
   time <- c(time,max(time+1))
   rcTime <- c(rcTime,rexp(1))
   t1 <- c(t1,max(time+1))
   v <- c(v,-1)
  }else{
   t <- t1[subjIndex] + rcTime[subjIndex]
   t1[subjIndex] <- NA
   subj <- samp[subjIndex]
   nr <- sample(0:maxR,1,replace=FALSE,prob=pr)
   nbrs <- rbind(el[el[,1]==subj,2:3,drop=FALSE],
                 el[el[,2]==subj & el[,3]>0.5,c(1,3),drop=FALSE]
   )
   nbrs <- nbrs[!(nbrs[,1] %in% samp) & nbrs[,1]!=subj,,drop=FALSE]
   nr <- min(nr,nrow(nbrs))
   if(nr>0){
    s <- sample.int(nrow(nbrs),nr,replace=FALSE)
    s <- s[!duplicated(nbrs[s,1])]
    nr <- length(s)
    samp <- c(samp,nbrs[s,1])
    recr <- c(recr,rep(subj,nr))
    tm <- t + t + (0:(nr-1)) / 1000000
    time <- c(time,tm)
    t1 <- c(t1,tm)
    rcTime <- c(rcTime,rexp(nr))
    v <- c(v,nbrs[s,2])
   }
  }
}
data.frame(subject=samp,recruiter=recr,time=time,v=v)
}


rds <- sampRDS(el, d, seeds,g,ns,FALSE, pr = c(0,.1,.9))

基本函数如上。

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏0 回帖

关键词：RDS 受访者 R程序 法律制度 参数设定

虽然可以跑出结果，但是说实话这两个函数中的一些命令还是没太看明白。

欢迎有兴趣的朋友一起参与讨论。

感谢分享，有原链接或者参考文献吗

关于RDS检验可以看一下康奈尔Heckathorn教授的文章。

从邻接矩阵el中进行抽取的关键函数部分：while(length(samp) < ss){
  subjIndex <- which.min(t1 + rcTime)
  if(length(subjIndex)==0){
   print("redraw")
   subj <- sample( (1:n)[-samp],1)
   samp <- c(samp,subj)
   recr <- c(recr,-1)
   time <- c(time,max(time+1))
   rcTime <- c(rcTime,rexp(1))
   t1 <- c(t1,max(time+1))
   v <- c(v,-1)
  }else{
   t <- t1[subjIndex] + rcTime[subjIndex]
   t1[subjIndex] <- NA
   subj <- samp[subjIndex]
   nr <- sample(0:maxR,1,replace=FALSE,prob=pr)
   nbrs <- rbind(el[el[,1]==subj,2:3,drop=FALSE],
                 el[el[,2]==subj & el[,3]>0.5,c(1,3),drop=FALSE]
   )
   nbrs <- nbrs[!(nbrs[,1] %in% samp) & nbrs[,1]!=subj,,drop=FALSE]
   nr <- min(nr,nrow(nbrs))
   if(nr>0){
    s <- sample.int(nrow(nbrs),nr,replace=FALSE)
    s <- s[!duplicated(nbrs[s,1])]
    nr <- length(s)
    samp <- c(samp,nbrs[s,1])
    recr <- c(recr,rep(subj,nr))
    tm <- t + t + (0:(nr-1)) / 1000000
    time <- c(time,tm)
    t1 <- c(t1,tm)
    rcTime <- c(rcTime,rexp(nr))
    v <- c(v,nbrs[s,2])
   }
  }
这个地方始终没研究明白，不知有没有高人可以解释一下这其中的步骤。

如果从RDS函数中把上面的程序单独抽取出来运行的话，结果就会报错。

还有一点不懂的：在function（RDS）中写的这个循环：while(length(samp) < ss) {}中为什么没有终止命令？不知哪位高人可以指点一二，谢谢。

在function（RDS）中还有这一部分：
if(nr>0){
    s <- sample.int(nrow(nbrs),nr,replace=FALSE)
    s <- s[!duplicated(nbrs[s,1])]
    nr <- length(s)
    samp <- c(samp,nbrs[s,1])
    recr <- c(recr,rep(subj,nr))
    tm <- t + t + (0:(nr-1)) / 1000000
    time <- c(time,tm)
    t1 <- c(t1,tm)
    rcTime <- c(rcTime,rexp(nr))
    v <- c(v,nbrs[s,2])
   }
然后没有对应的else{}