# R/deffCR.R In PracTools: Tools for Designing and Weighting Survey Samples

#### Documented in deffCR

```deffCR <- function(w, strvar=NULL, clvar=NULL, Wh=NULL, nest=FALSE, y){
if (any(is.null(w) | is.null(y))) stop("w and y are required for Chen-Rust deff")
n <- length(w)
sig2 <- n / (n-1) * sum(w*(y - sum(w*y)/sum(w))^2) / (sum(w)-1)

if (!is.null(strvar)){             # strata = Y
strat <- unique(strvar)
H <- length(strat)
nh <- as.vector(table(strvar))
sig2h <- deff.s <- cv2h <- vector("numeric", length=H)
if (is.null(Wh)){           # estimate Wh if it's not fed as a parameter
Wh <- by(w, INDICES=strvar, sum) / sum(w)
Wh <- as.vector(Wh)
}
for (h in strat){
wsub <- w[strvar == h]
ysub <- y[strvar == h]
cv2h[h] <- (nh[h]-1)/nh[h] * var(wsub) / mean(wsub)^2
if (is.null(clvar)){       # strata = Y, clusters = N
sig2h[h] <- nh[h] / (nh[h]-1) * sum(wsub*(ysub - sum(wsub*ysub)/sum(wsub))^2) / (sum(wsub)-1)
deff.s[h] <- Wh[h]^2/ nh[h] * n * sig2h[h] / sig2
}
}

if (is.null(clvar)) {
deff.w <- 1 + cv2h
out <- list("strata deffs" = cbind(stratum = 1:H, deff.w = deff.w, deff.c = NULL, deff.s = deff.s),
"overall deff" = sum(deff.w * deff.s)
)
}

if (!is.null(clvar)){           # strata = Y, clusters = Y
sam.dsgn <- survey::svydesign(ids = ~clvar, strata = ~strvar, data = data.frame(y), weights = w, nest = nest)
strmns <- survey::svyby(~y, by = ~strvar, design=sam.dsgn, FUN=survey::svymean, deff=TRUE)
strdeff <- strmns\$DEff.y
nh.star <- rhoh <- vector("numeric", length=H)

for (h in strat){
wsub <- w[strvar == h]
clsub <- clvar[strvar == h]
ysub <- y[strvar == h]
nh.star.num <- by(wsub, INDICES=clsub, sum)
nh.star.num <- sum(nh.star.num^2)
nh.star.den <- sum(wsub^2)
nh.star[h] <- nh.star.num / nh.star.den
rhoh[h] <- (strdeff[h] - (1 + cv2h[h])) / (1 + cv2h[h]*(nh.star[h]-1))
sig2h[h] <- nh[h] / (nh[h]-1) * sum(wsub*(ysub - sum(wsub*ysub)/sum(wsub))^2) / (sum(wsub)-1)
deff.s[h] <- Wh[h]^2/ nh[h] * n * sig2h[h] / sig2
}
deff.w = 1 + cv2h
deff.c = 1 + (nh.star-1)*rhoh
out <- list("strata components" = cbind(stratum = 1:H, deff.w = deff.w, deff.c = deff.c, deff.s = deff.s),
"overall deff" = sum(deff.w * deff.c * deff.s)
)
}
}

if (is.null(strvar)){              # strata = N
nh <- n
if (is.null(clvar)){           # strat = N, clusters = N
cv2h <- (nh-1)/nh * var(w) / mean(w)^2
deff.w = 1 + cv2h
out <- c(deff.w = deff.w)
}
else {                      # strat = N, clusters = Y
nh.star.num <- by(w, INDICES=clvar, sum)
nh.star.num <- sum(nh.star.num^2)
nh.star.den <- sum(w^2)
nh.star <- nh.star.num / nh.star.den
sam.dsgn <- survey::svydesign(ids = ~clvar, strata = NULL, data = data.frame(y), weights = w)
strmns <- survey::svymean(~y, design=sam.dsgn, deff=TRUE)
strdeff <- survey::deff(strmns)
cv2h <- (nh-1)/nh * var(w) / mean(w)^2
rhoh <- (strdeff - (1 + cv2h)) / (1 + cv2h*(nh.star-1))
deff.w = 1 + cv2h
deff.c = 1 + (nh.star-1)*rhoh
out <- list(cbind(deff.w = deff.w, deff.c = deff.c, deff.s = NULL),
"overall deff" = deff.w * deff.c
)
}
}
return(out)
}
```

## Try the PracTools package in your browser

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

PracTools documentation built on Nov. 4, 2021, 1:08 a.m.