R/design-sampling.R

In R4VN: Health Data Analysis and Publication-Ready Reporting

.r4vn_sampling_pps_pi <- function(mos,n) {
  mos<-as.numeric(mos)
  if(any(!is.finite(mos))||any(mos<=0)) stop("MOS must contain positive finite values.",call.=FALSE)
  N<-length(mos); n<-as.integer(n)
  if(n<1||n>N) stop("PPS sample size must be between 1 and the number of frame units.",call.=FALSE)
  pi<-rep(NA_real_,N); remain<-seq_len(N); nrem<-n
  while(length(remain) && nrem>0) {
    cand<-nrem*mos[remain]/sum(mos[remain])
    certainty<-which(cand>=1-1e-12)
    if(!length(certainty)) { pi[remain]<-cand; break }
    ids<-remain[certainty]
    pi[ids]<-1
    nrem<-nrem-length(ids)
    remain<-setdiff(remain,ids)
  }
  if(nrem==0 && length(remain)) pi[remain]<-0
  # numerical normalization among non-certainty units
  if(abs(sum(pi)-n)>1e-8) {
    idx<-which(pi<1 & pi>0)
    if(length(idx)) pi[idx]<-pi[idx]*(n-sum(pi[-idx]))/sum(pi[idx])
  }
  pi
}

.r4vn_sampling_pps_systematic_indices <- function(mos,n) {
  pi<-.r4vn_sampling_pps_pi(mos,n)
  cs<-cumsum(pi)
  u<-stats::runif(1,0,1)
  hits<-u+0:(n-1)
  idx<-vapply(hits,function(h) which(cs>=h-1e-12)[1],integer(1))
  list(index=idx,pi=pi,random_start=u,hits=hits)
}

.r4vn_sampling_pps_package_indices <- function(mos,n,algorithm=c("brewer","sampford","tille","pivotal","maxentropy")) {
  algorithm <- match.arg(algorithm)
  .r4vn_design_require("sampling", paste0("PPS ",algorithm," sampling"))
  pi <- .r4vn_sampling_pps_pi(mos,n)
  certainty <- which(pi >= 1-1e-10)
  remain <- which(pi < 1-1e-10 & pi > 1e-12)
  selected <- certainty
  if(length(remain)) {
    fn <- switch(algorithm,
      brewer=sampling::UPbrewer,
      sampford=sampling::UPsampford,
      tille=sampling::UPtille,
      pivotal=sampling::UPpivotal,
      maxentropy=sampling::UPmaxentropy)
    s <- fn(pi[remain])
    selected <- c(selected, remain[which(as.numeric(s)>0)])
  }
  if(length(selected) != as.integer(n)) stop("The PPS algorithm did not return the requested fixed sample size.",call.=FALSE)
  list(index=selected,pi=pi,certainty=certainty,algorithm=algorithm)
}

.r4vn_sampling_allocate_n <- function(sizes,n,method=c("proportional","equal")) {
  method<-match.arg(method); n<-as.integer(n); sizes<-as.numeric(sizes)
  if(n>sum(sizes)) stop("Requested sample exceeds available records.",call.=FALSE)
  raw<-if(method=="proportional") n*sizes/sum(sizes) else rep(n/length(sizes),length(sizes))
  alloc<-pmin(floor(raw),sizes)
  left<-n-sum(alloc)
  frac<-raw-floor(raw)
  while(left>0) {
    eligible<-which(alloc<sizes)
    if(!length(eligible)) break
    ord<-eligible[order(frac[eligible],sizes[eligible]-alloc[eligible],decreasing=TRUE)]
    for(j in ord) {
      if(left<=0) break
      alloc[j]<-alloc[j]+1; left<-left-1
    }
  }
  as.integer(alloc)
}

.r4vn_sampling_frame_hash <- function(data,idvar,extra=character()) {
  cols<-unique(c(idvar,extra)); cols<-cols[cols %in% names(data)]
  .r4vn_design_hash(list(row_order=seq_len(nrow(data)),data=data[,cols,drop=FALSE]))
}

.r4vn_sampling_with_rng <- function(seed,code) {
  old_kind<-RNGkind()
  had_seed<-exists(".Random.seed",envir=.GlobalEnv,inherits=FALSE)
  if(had_seed) old_seed<-get(".Random.seed",envir=.GlobalEnv,inherits=FALSE)
  on.exit({
    try(do.call(RNGkind,as.list(old_kind)),silent=TRUE)
    if(had_seed) assign(".Random.seed",old_seed,envir=.GlobalEnv) else if(exists(".Random.seed",envir=.GlobalEnv,inherits=FALSE)) rm(".Random.seed",envir=.GlobalEnv)
  },add=TRUE)
  # Rejection is the post-R-3.6 sample.kind and makes the manifest explicit.
  try(RNGkind(kind="Mersenne-Twister",normal.kind="Inversion",sample.kind="Rejection"),silent=TRUE)
  set.seed(as.integer(seed))
  before<-get(".Random.seed",envir=.GlobalEnv,inherits=FALSE)
  value<-force(code)
  after<-get(".Random.seed",envir=.GlobalEnv,inherits=FALSE)
  list(value=value,rng_kind=RNGkind(),rng_state_before=before,rng_state_after=after)
}

.r4vn_sampling_block_allocate <- function(N,groups,ratio,block_sizes) {
  ratio<-as.integer(ratio); block_sizes<-sort(unique(as.integer(block_sizes)))
  if(length(groups)<2||length(groups)!=length(ratio)||any(ratio<=0)) stop("Groups and positive integer allocation ratios are required.",call.=FALSE)
  unit<-sum(ratio)
  if(N %% unit != 0L) stop(sprintf("The participant count (%d) is incompatible with the requested allocation ratio. For ratio %s, the count must be divisible by %d.",N,paste(ratio,collapse=":"),unit),call.=FALSE)
  ok<-block_sizes[is.finite(block_sizes) & block_sizes>=unit & block_sizes%%unit==0]
  if(!length(ok)) stop("Every block size must be a positive multiple of the sum of allocation-ratio units.",call.=FALSE)

  can_make <- function(total, sizes) {
    total <- as.integer(total)
    if(total==0L) return(TRUE)
    if(total<0L) return(FALSE)
    reachable <- rep(FALSE,total+1L); reachable[1L] <- TRUE
    for(i in 0:total) if(reachable[i+1L]) for(s in sizes) if(i+s<=total) reachable[i+s+1L] <- TRUE
    reachable[total+1L]
  }
  if(!can_make(N,ok)) stop(sprintf("The participant count (%d) cannot be completed exactly using the requested block sizes: %s.",N,paste(ok,collapse=", ")),call.=FALSE)

  out<-character(); blocks<-integer(); block_size<-integer(); b<-0L; remaining<-N
  while(remaining>0L) {
    feasible<-ok[ok<=remaining & vapply(ok,function(bs) can_make(remaining-bs,ok),logical(1))]
    if(!length(feasible)) stop("Unable to construct a complete permuted-block sequence with the requested settings.",call.=FALSE)
    b<-b+1L; bs<-feasible[sample.int(length(feasible),1L)]
    mult<-bs/unit
    content<-rep(groups,times=ratio*mult)
    content<-sample(content,length(content),replace=FALSE)
    out<-c(out,content)
    blocks<-c(blocks,rep(b,bs))
    block_size<-c(block_size,rep(bs,bs))
    remaining<-remaining-bs
  }
  list(group=out,block=blocks,block_size=block_size)
}

.r4vn_sampling_run <- function(data,method,idvar,seed,params=list()) {
  if(!is.data.frame(data)||!nrow(data)) stop("A non-empty sampling frame is required.",call.=FALSE)
  if(!idvar %in% names(data)) stop("ID variable was not found in the sampling frame.",call.=FALSE)
  if(anyNA(data[[idvar]])||anyDuplicated(data[[idvar]])) stop("The ID variable must be complete and unique.",call.=FALSE)
  seed<-as.integer(seed)
  if(!is.finite(seed)) stop("Seed must be an integer.",call.=FALSE)

  extra<-unique(unlist(params[c("strata","cluster","mos")],use.names=FALSE))
  frame_hash<-.r4vn_sampling_frame_hash(data,idvar,extra)
  run_id<-.r4vn_design_id("SP")

  rr<-.r4vn_sampling_with_rng(seed,{
    N<-nrow(data)
    full<-data
    full$sample_selected<-FALSE
    full$sample_order<-NA_integer_
    full$sampling_probability<-NA_real_
    full$sampling_weight<-NA_real_
    full$sampling_run<-run_id
    selected_rows<-integer()
    details<-list()

    if(method=="simple") {
      n<-as.integer(params$n)
      if(n<1||n>N) stop("n must be between 1 and the frame size.",call.=FALSE)
      selected_rows<-sample.int(N,n,replace=FALSE)
      full$sampling_probability<-n/N; full$sampling_weight<-N/n
      details<-list(n=n,population_n=N)
    }

    if(method=="systematic") {
      n<-as.integer(params$n)
      if(n<1||n>N) stop("n must be between 1 and the frame size.",call.=FALSE)
      k<-N/n; start<-stats::runif(1,0,k)
      pos<-floor(start+(0:(n-1))*k)+1L
      pos<-pmin(pos,N)
      if(anyDuplicated(pos)) stop("Systematic selection produced duplicate positions; check n and frame size.",call.=FALSE)
      selected_rows<-pos
      full$sampling_probability<-n/N; full$sampling_weight<-N/n
      details<-list(interval=k,random_start=start,first_position=pos[1])
    }

    if(method=="stratified") {
      strata<-params$strata
      if(!strata %in% names(full)) stop("Stratification variable not found.",call.=FALSE)
      n<-as.integer(params$n); allocation<-params$allocation %||% "proportional"
      lev<-unique(as.character(full[[strata]])); sizes<-vapply(lev,function(z) sum(as.character(full[[strata]])==z),integer(1))
      if(identical(allocation,"custom")) {
        custom<-.r4vn_design_parse_named_counts(params$custom_counts %||% "")
        alloc<-as.integer(custom[lev]); alloc[is.na(alloc)]<-0L
        if(sum(alloc)!=n) stop("Custom stratum counts must sum to n and use stratum labels exactly.",call.=FALSE)
        if(any(alloc>sizes)) stop("A custom stratum count exceeds the available stratum size.",call.=FALSE)
      } else alloc<-.r4vn_sampling_allocate_n(sizes,n,allocation)
      for(j in seq_along(lev)) {
        rows<-which(as.character(full[[strata]])==lev[j]); if(alloc[j]) selected_rows<-c(selected_rows,sample(rows,alloc[j],replace=FALSE))
        full$sampling_probability[rows]<-alloc[j]/length(rows)
        full$sampling_weight[rows]<-if(alloc[j]>0) length(rows)/alloc[j] else NA_real_
      }
      details<-list(strata=strata,stratum_sizes=setNames(sizes,lev),stratum_n=setNames(alloc,lev),allocation=allocation)
    }

    if(method=="pps_systematic") {
      mos<-params$mos; n<-as.integer(params$n)
      if(!mos %in% names(full)) stop("MOS variable not found.",call.=FALSE)
      z<-.r4vn_sampling_pps_systematic_indices(full[[mos]],n)
      selected_rows<-z$index
      full$sampling_probability<-z$pi
      full$sampling_weight<-ifelse(z$pi>0,1/z$pi,NA_real_)
      details<-list(mos=mos,random_start=z$random_start,hits=z$hits,certainty_ids=full[[idvar]][z$pi>=1-1e-12])
    }

    if(method %in% c("pps_brewer","pps_sampford","pps_tille","pps_pivotal","pps_maxentropy")) {
      mos<-params$mos; n<-as.integer(params$n)
      if(!mos %in% names(full)) stop("MOS variable not found.",call.=FALSE)
      alg<-sub("^pps_","",method)
      z<-.r4vn_sampling_pps_package_indices(full[[mos]],n,alg)
      selected_rows<-z$index
      full$sampling_probability<-z$pi
      full$sampling_weight<-ifelse(z$pi>0,1/z$pi,NA_real_)
      details<-list(mos=mos,algorithm=alg,certainty_ids=full[[idvar]][z$certainty])
    }

    if(method=="pps_wr") {
      mos<-params$mos; n<-as.integer(params$n)
      if(!mos %in% names(full)) stop("MOS variable not found.",call.=FALSE)
      m<-as.numeric(full[[mos]]); if(any(!is.finite(m))||any(m<=0)) stop("MOS must be positive and complete.",call.=FALSE)
      p<-m/sum(m); draws<-sample.int(N,n,replace=TRUE,prob=p)
      counts<-tabulate(draws,nbins=N)
      selected_rows<-which(counts>0)
      full$selection_count<-counts
      full$draw_probability<-p
      full$sampling_probability<-1-(1-p)^n
      full$sampling_weight<-ifelse(full$sampling_probability>0,1/full$sampling_probability,NA_real_)
      details<-list(mos=mos,draws=draws)
    }

    if(method=="cluster") {
      cluster<-params$cluster; ncl<-as.integer(params$n_clusters)
      if(!cluster %in% names(full)) stop("Cluster variable not found.",call.=FALSE)
      clu<-unique(as.character(full[[cluster]])); if(ncl<1||ncl>length(clu)) stop("Invalid number of clusters.",call.=FALSE)
      selcl<-sample(clu,ncl,replace=FALSE); selected_rows<-which(as.character(full[[cluster]]) %in% selcl)
      full$sampling_probability<-ncl/length(clu); full$sampling_weight<-length(clu)/ncl
      details<-list(cluster=cluster,selected_clusters=selcl)
    }

    if(method=="multistage_pps") {
      cluster<-params$cluster; mos<-params$mos; ncl<-as.integer(params$n_clusters); within<-as.integer(params$n_within)
      if(!all(c(cluster,mos)%in%names(full))) stop("Cluster and MOS variables are required.",call.=FALSE)
      spl<-split(seq_len(N),as.character(full[[cluster]]))
      clu<-names(spl)
      cm<-vapply(spl,function(ix){z<-unique(full[[mos]][ix]); if(length(z)!=1||!is.finite(z)||z<=0) stop("MOS must be positive and constant within each cluster.",call.=FALSE); as.numeric(z)},numeric(1))
      pps<-.r4vn_sampling_pps_systematic_indices(cm,ncl)
      selcl<-clu[pps$index]
      selected_rows<-integer()
      pi_psu<-setNames(pps$pi,clu)
      full$psu_inclusion_probability<-pi_psu[as.character(full[[cluster]])]
      full$within_cluster_probability<-NA_real_
      for(cl in selcl) {
        rows<-spl[[cl]]; take<-min(within,length(rows)); rr2<-sample(rows,take,replace=FALSE); selected_rows<-c(selected_rows,rr2)
        full$within_cluster_probability[rows]<-take/length(rows)
      }
      full$sampling_probability<-full$psu_inclusion_probability*full$within_cluster_probability
      full$sampling_weight<-ifelse(full$sampling_probability>0,1/full$sampling_probability,NA_real_)
      details<-list(cluster=cluster,mos=mos,selected_clusters=selcl,psu_random_start=pps$random_start,n_within=within)
    }

    if(method %in% c("rct_complete","rct_simple","rct_block","rct_stratified_block")) {
      groups<-trimws(unlist(strsplit(params$groups %||% "Control,Intervention",",",fixed=TRUE)))
      ratio<-.r4vn_design_parse_numlist(params$ratio %||% "1,1")
      if(length(groups)!=length(ratio)||length(groups)<2||any(ratio<=0)||any(abs(ratio-round(ratio))>1e-8)) stop("RCT group names and allocation ratios are invalid. Ratios must be positive integers such as 1,1 or 2,1.",call.=FALSE)
      ratio<-as.integer(round(ratio)); unit<-sum(ratio)
      full$random_group<-NA_character_; full$random_order<-NA_integer_; full$random_block<-NA_integer_; full$random_block_size<-NA_integer_
      allocate_rows<-function(rows) {
        nr<-length(rows)
        if(nr %% unit != 0L) stop(sprintf("Exact allocation is required. A stratum/participant set of %d cannot be allocated exactly in the ratio %s. The count must be divisible by %d.",nr,paste(ratio,collapse=":"),unit),call.=FALSE)
        ord<-sample(rows,nr,replace=FALSE)
        if(method %in% c("rct_simple","rct_complete")) {
          cnt<-(nr/unit)*ratio
          g<-sample(rep(groups,times=cnt),nr,replace=FALSE); b<-rep(NA_integer_,nr); bsize<-rep(NA_integer_,nr)
        } else {
          bs<-.r4vn_design_parse_numlist(params$block_sizes %||% "2,4,6")
          al<-.r4vn_sampling_block_allocate(nr,groups,ratio,bs); g<-al$group; b<-al$block; bsize<-al$block_size
        }
        full$random_group[ord]<<-g; full$random_order[ord]<<-seq_len(nr); full$random_block[ord]<<-b; full$random_block_size[ord]<<-bsize
      }
      if(method=="rct_stratified_block") {
        strata<-params$strata; if(!strata%in%names(full)) stop("RCT stratification variable not found.",call.=FALSE)
        for(z in unique(as.character(full[[strata]]))) allocate_rows(which(as.character(full[[strata]])==z))
      } else allocate_rows(seq_len(N))
      selected_rows<-seq_len(N)
      details<-list(groups=groups,ratio=ratio,block_sizes=if(method%in%c("rct_block","rct_stratified_block")) .r4vn_design_parse_numlist(params$block_sizes %||% "2,4,6") else NULL,strata=if(method=="rct_stratified_block") params$strata else NULL)
      counts<-table(full$random_group)
      expected<-(N/unit)*ratio; names(expected)<-groups
      if(!identical(as.integer(counts[groups]),as.integer(expected))) stop("Internal allocation check failed: the final group totals do not match the requested allocation ratio exactly.",call.=FALSE)
    }

    if(!length(selected_rows)) stop("Sampling method did not select any records.",call.=FALSE)
    # sample_order refers to unique selected records; PPS-WR retains draw order separately.
    urows<-unique(selected_rows)
    full$sample_selected[urows]<-TRUE
    full$sample_order[urows]<-seq_along(urows)
    list(full=full,selected=full[urows,,drop=FALSE],selected_rows=urows,details=details)
  })

  val<-rr$value
  result<-list(
    run_id=run_id,method=method,seed=seed,rng_kind=rr$rng_kind,
    rng_state_before=rr$rng_state_before,rng_state_after=rr$rng_state_after,
    frame_hash=frame_hash,idvar=idvar,params=params,n_frame=nrow(data),
    n_selected=nrow(val$selected),selected=val$selected,full=val$full,
    details=val$details,created=as.character(Sys.time()),r_version=R.version.string,
    r4vn_version=.r4vn_design_r4vn_version()
  )
  result$summary_text<-sprintf("R4VN used %s on a frame of %d records and returned %d selected records. The run is identified by %s and was generated with seed %d.",method,nrow(data),result$n_selected,run_id,seed)
  class(result)<-"r4vn_sampling"
  result
}

.r4vn_sampling_verify_frame <- function(run,data) {
  extra<-unique(unlist(run$params[c("strata","cluster","mos")],use.names=FALSE))
  current<-.r4vn_sampling_frame_hash(data,run$idvar,extra)
  list(ok=identical(current,run$frame_hash),original_hash=run$frame_hash,current_hash=current,original_n=run$n_frame,current_n=nrow(data))
}

.r4vn_sampling_reproduce <- function(run,data) {
  if(!inherits(run,"r4vn_sampling")) stop("run must be an r4vn_sampling object.",call.=FALSE)
  chk<-.r4vn_sampling_verify_frame(run,data)
  if(!isTRUE(chk$ok)) stop(sprintf("Sampling frame does not match the original run (original N=%d, current N=%d). Exact reproduction was blocked.",chk$original_n,chk$current_n),call.=FALSE)
  .r4vn_sampling_run(data,run$method,run$idvar,run$seed,run$params)
}

.r4vn_sampling_write_template <- function(file,type=c("simple","pps","multistage","rct")) {
  type<-match.arg(type)
  .r4vn_design_require("openxlsx","Excel sampling-frame templates")
  wb<-openxlsx::createWorkbook()
  openxlsx::addWorksheet(wb,"Instructions")
  instructions<-data.frame(
    Field=c("R4VN template","Purpose","Important"),
    Description=c(type,"Replace example rows with your own sampling frame. Column names may be changed because the Studio allows column mapping.","ID values must be unique and complete. PPS MOS values must be positive."),
    stringsAsFactors=FALSE)
  openxlsx::writeData(wb,"Instructions",instructions)
  if(type=="simple") dat<-data.frame(id=sprintf("ID%03d",1:6),stratum=c("A","A","A","B","B","B"),eligible=1)
  if(type=="pps") dat<-data.frame(cluster_id=sprintf("C%03d",1:6),stratum=c("Urban","Urban","Rural","Rural","Rural","Urban"),MOS=c(12540,8320,4210,7130,9500,6000),eligible=1)
  if(type=="multistage") dat<-data.frame(person_id=sprintf("P%03d",1:12),cluster_id=rep(sprintf("C%02d",1:4),each=3),MOS=rep(c(1200,850,2300,1400),each=3),stratum=rep(c("A","A","B","B"),each=3))
  if(type=="rct") dat<-data.frame(participant_id=sprintf("P%03d",1:12),site=rep(c("H1","H2"),each=6),sex=rep(c("Female","Male"),6))
  openxlsx::addWorksheet(wb,"Sampling_Frame"); openxlsx::writeData(wb,"Sampling_Frame",dat)
  openxlsx::addWorksheet(wb,"_R4VN_INFO"); openxlsx::writeData(wb,"_R4VN_INFO",data.frame(key=c("template_type","template_version","generated_by"),value=c(type,"1.0","R4VN design()")))
  openxlsx::saveWorkbook(wb,file,overwrite=TRUE)
  invisible(file)
}

.r4vn_sampling_manifest <- function(run, include_rng_state = TRUE) {
  if(!inherits(run,"r4vn_sampling")) stop("run must be an r4vn_sampling object.",call.=FALSE)
  x <- run

  # Keep only the minimum result fingerprint needed to verify an exact rerun.
  # The full sampling frame and selected records are deliberately NOT embedded
  # in a project file because they may contain sensitive research data.
  if(!is.null(run$selected) && run$idvar %in% names(run$selected)) {
    x$selected_ids <- as.character(run$selected[[run$idvar]])
  } else {
    x$selected_ids <- NULL
  }

  if(!is.null(run$full) && run$idvar %in% names(run$full) && "random_group" %in% names(run$full)) {
    keep0 <- c(run$idvar,"random_group","random_order","random_block","random_block_size")
    if(!is.null(run$params$strata)) keep0 <- c(keep0, run$params$strata)
    keep <- intersect(unique(keep0),names(run$full))
    x$randomization_fingerprint <- run$full[,keep,drop=FALSE]
    x$randomization_fingerprint[[run$idvar]] <- as.character(x$randomization_fingerprint[[run$idvar]])
  } else {
    x$randomization_fingerprint <- NULL
  }

  # Keep explicit NULL fields so `$selected` cannot partially match
  # `selected_ids` (and similarly for full data) after serialization.
  x["selected"] <- list(NULL)
  x["full"] <- list(NULL)
  if(!include_rng_state) {
    x$rng_state_before <- NULL
    x$rng_state_after <- NULL
  }
  x$manifest_only <- TRUE
  class(x) <- c("r4vn_sampling_manifest","r4vn_sampling")
  x
}


.r4vn_sampling_method_paragraph <- function(run) {
  if(!inherits(run, "r4vn_sampling")) return("")
  method <- run$method
  n <- run$n_selected
  N <- run$n_frame
  seed <- run$seed
  idvar <- run$idvar %||% "ID"
  d <- run$details %||% list()
  p <- run$params %||% list()

  if(method == "simple") {
    return(sprintf(
      "A simple random sample of %d units was selected without replacement from %d eligible units. Each unit had an equal inclusion probability of %.4f. The random sample was generated in R4VN using seed %d.",
      n, N, n / N, seed
    ))
  }

  if(method == "systematic") {
    k <- d$interval %||% (N / n)
    st <- d$random_start %||% NA_real_
    fp <- d$first_position %||% NA_integer_
    return(sprintf(
      "Systematic random sampling was used to select %d units from an ordered frame of %d units. The sampling interval was k = N/n = %d/%d = %.4f. A random start of %.4f on the interval scale was generated, giving first selected frame position %s, after which every k-th position was selected. The procedure was generated in R4VN using seed %d.",
      n, N, N, n, k, st, ifelse(is.finite(fp), as.character(fp), "NA"), seed
    ))
  }

  if(method == "stratified") {
    strata <- d$strata %||% p$strata %||% "the selected stratum variable"
    alloc <- d$allocation %||% p$allocation %||% "proportional"
    sn <- d$stratum_n
    detail <- if(!is.null(sn) && length(sn)) paste(paste0(names(sn), " = ", as.integer(sn)), collapse = "; ") else ""
    return(sprintf(
      "Stratified random sampling was performed using %s as the stratification variable. A total sample of %d was allocated using %s allocation%s, and simple random sampling without replacement was then performed within each stratum. The procedure was generated in R4VN using seed %d.",
      strata, n, alloc, if(nzchar(detail)) paste0(" (", detail, ")") else "", seed
    ))
  }

  if(method == "pps_systematic") {
    mos <- d$mos %||% p$mos %||% "the measure-of-size variable"
    rs <- d$random_start %||% NA_real_
    return(sprintf(
      "Systematic probability-proportional-to-size (PPS) sampling without replacement was used to select %d units from %d units, with %s as the measure of size. Inclusion probabilities were proportional to the measure of size, with certainty units handled automatically where applicable. The systematic PPS random start was %.4f. The sample was generated in R4VN using seed %d.",
      n, N, mos, rs, seed
    ))
  }

  if(method %in% c("pps_brewer","pps_sampford","pps_tille","pps_pivotal","pps_maxentropy")) {
    alg <- d$algorithm %||% sub("^pps_", "", method)
    mos <- d$mos %||% p$mos %||% "the measure-of-size variable"
    return(sprintf(
      "%s PPS sampling without replacement was used to select %d units from %d units, with %s as the measure of size. First-order inclusion probabilities and sampling weights were retained in the output. The sample was generated in R4VN using seed %d.",
      tools::toTitleCase(alg), n, N, mos, seed
    ))
  }

  if(method == "pps_wr") {
    mos <- d$mos %||% p$mos %||% "the measure-of-size variable"
    return(sprintf(
      "Probability-proportional-to-size sampling with replacement was performed using %s as the measure of size. A total of %d draws were made from %d source units. Draw probabilities, selection counts, inclusion probabilities, and sampling weights were retained in the output. The sample was generated in R4VN using seed %d.",
      mos, as.integer(p$n %||% n), N, seed
    ))
  }

  if(method == "cluster") {
    cluster <- d$cluster %||% p$cluster %||% "the cluster variable"
    nc <- length(d$selected_clusters %||% character())
    return(sprintf(
      "One-stage cluster sampling was performed using %s as the cluster identifier. %d clusters were selected by simple random sampling, and all eligible units within selected clusters were included. The sample was generated in R4VN using seed %d.",
      cluster, nc, seed
    ))
  }

  if(method == "multistage_pps") {
    cluster <- d$cluster %||% p$cluster %||% "the cluster variable"
    mos <- d$mos %||% p$mos %||% "the measure-of-size variable"
    nc <- length(d$selected_clusters %||% character())
    nw <- d$n_within %||% p$n_within %||% NA_integer_
    return(sprintf(
      "A multistage probability sample was selected. At the first stage, %d clusters defined by %s were selected using systematic PPS with %s as the measure of size. At the second stage, up to %s units were selected by simple random sampling within each selected cluster. Overall selection probabilities and sampling weights were retained. The procedure was generated in R4VN using seed %d.",
      nc, cluster, mos, as.character(nw), seed
    ))
  }

  sprintf("A probability sample of %d units was selected from %d source units using %s in R4VN with random seed %d.", n, N, method, seed)
}

.r4vn_sampling_generated_frame <- function(n, id_name = "ID") {
  n <- as.integer(n)
  if (!is.finite(n) || n < 1L) stop("Population/participant count must be at least 1.", call. = FALSE)
  if (n > 10000000L) stop("Generated lists are limited to 10,000,000 IDs in the Studio.", call. = FALSE)
  out <- data.frame(seq_len(n), stringsAsFactors = FALSE)
  names(out) <- id_name
  out
}

.r4vn_randomization_print_table <- function(run) {
  if (!inherits(run, "r4vn_sampling")) stop("run must be an r4vn_sampling object.", call. = FALSE)
  idvar <- run$idvar
  if (!is.null(run$full) && "random_group" %in% names(run$full)) {
    z <- run$full
  } else if (!is.null(run$randomization_fingerprint) && "random_group" %in% names(run$randomization_fingerprint)) {
    z <- run$randomization_fingerprint
  } else {
    stop("The run does not contain a randomization allocation table.", call. = FALSE)
  }

  strata_var <- run$params$strata %||% NULL
  keep <- c(idvar, "random_order", "random_group", "random_block", "random_block_size")
  if (!is.null(strata_var) && strata_var %in% names(z)) keep <- c(keep, strata_var)
  keep <- unique(intersect(keep, names(z)))
  z <- z[, keep, drop = FALSE]

  stratified <- !is.null(strata_var) && strata_var %in% names(z)
  if (stratified && "random_order" %in% names(z)) {
    z <- z[order(as.character(z[[strata_var]]), z$random_order, na.last = TRUE), , drop = FALSE]
  } else if ("random_order" %in% names(z)) {
    z <- z[order(z$random_order, na.last = TRUE), , drop = FALSE]
  }

  sequence_value <- if (!stratified && "random_order" %in% names(z)) z$random_order else seq_len(nrow(z))
  out <- data.frame(
    Sequence = sequence_value,
    Group = as.character(z$random_group),
    stringsAsFactors = FALSE
  )
  if (stratified) {
    out$Stratum <- as.character(z[[strata_var]])
    if ("random_order" %in% names(z)) out[["Order within stratum"]] <- z$random_order
  }
  if ("random_block" %in% names(z) && any(!is.na(z$random_block))) out$Block <- z$random_block
  if ("random_block_size" %in% names(z) && any(!is.na(z$random_block_size))) out[["Block size"]] <- z$random_block_size
  rownames(out) <- NULL
  out
}


.r4vn_randomization_method_paragraph <- function(run) {
  if(!inherits(run,"r4vn_sampling")) stop("run must be an r4vn_sampling object.",call.=FALSE)
  groups <- run$details$groups %||% character()
  ratio <- run$details$ratio %||% numeric()
  ratio_txt <- if(length(ratio)) paste(ratio,collapse=":") else "the prespecified"
  group_txt <- if(length(groups)==2L) paste0(groups[1]," and ",groups[2]) else paste(groups,collapse=", ")
  method_txt <- switch(run$method,
    rct_simple = "random allocation with fixed group totals",
    rct_complete = "complete randomization with fixed group totals",
    rct_block = "permuted-block randomization",
    rct_stratified_block = "stratified permuted-block randomization",
    run$method_label %||% run$method
  )
  extra <- character()
  if(run$method %in% c("rct_block","rct_stratified_block") && !is.null(run$details$block_sizes)) {
    extra <- c(extra,paste0("randomly varying block sizes of ",paste(run$details$block_sizes,collapse=", ")))
  }
  if(run$method=="rct_stratified_block" && !is.null(run$details$strata)) {
    extra <- c(extra,paste0("stratification by ",run$details$strata))
  }
  extra_txt <- if(length(extra)) paste0(" using ",paste(extra,collapse=" and ")) else ""
  paste0(
    "Participants were allocated in a ",ratio_txt," ratio to ",group_txt,
    " using ",method_txt,extra_txt,". Exact final group totals were enforced to match the prespecified allocation ratio. ",
    "The allocation sequence was generated in R4VN using random seed ",run$seed,"."
  )
}