.huberpsiw <-
function(x, k, w=FALSE, asymmetric=FALSE){
# test whether k is valid
if (k <= 0) stop("k must not be <= 0!\n")
# weights
if (w){
if (asymmetric){
res <- pmin(1, k/x)
}else{
res <- pmin(1, k / abs(x))
}
}else{
# psi
if (asymmetric){
res <- pmin(k, x)
}else{
res <- pmin(k, pmax(-k, x))
}
}
return(res)
}
.inflmad <-
function(x, weight){
est <- .svyquantiles(x, weight, probs=0.5)
S <- .svyquantiles(abs(x - est), weight, probs=0.5) * 1.4826
# density estimate at the median
f <- .svydensity(x, weight, quant=0.5, ngrid=401)
# quantile of T+S
at.plus <- sum( weight[1:length(which(x <= est+S))] ) / sum(weight)
# quantile of T-S
at.minus <- sum( weight[1:length(which(x <= est-S))] ) / sum(weight)
#density estimate at the quantile(T+S)
f.plus <- .svydensity(x, weight, quant=at.plus, ngrid=401)
#density estimate at the quantile(TSS)
f.minus <- .svydensity(x, weight, quant=at.minus, ngrid=401)
numerator <- sign(abs(x - est) - S) - ((f.plus - f.minus) / f) * sign(x - est)
denominator <- 2 * (f.plus + f.minus)
res <- as.vector(numerator / denominator)
}
.inflmest <-
function(x, weight, est, k, S, asymmetric, exact){
u <- (x - est) / S
# psi
psi.value <- .huberpsiw(u, k, w=FALSE, asymmetric=asymmetric)
# psi prime
if (asymmetric){
psiprime.value <- 1 * (u <= k)
}else{
psiprime.value <- 1 * (abs(u) <= k)
}
# exact or not?
if(exact){
numerator <- psi.value * S - .inflmad(x, weight) * (sum(weight * u * psiprime.value)) / (sum(weight))
}
else{
numerator <- psi.value
}
denominator <- sum(weight * psiprime.value) / sum(weight)
res <- numerator / denominator
return(res)
}
.infltrim <-
function(x, weight, lo=0, hi=1){
# x <- sort(x)
if(is.null(weight)){
weight <- rep(1, length(x))
}
xilo <- .svyquantiles(x, weight, probs=lo)
xihi <- .svyquantiles(x, weight, probs=hi)
below <- floor(lo * length(x))
above <- ceiling(hi * length(x))
mat <- c(rep((1 - lo)*xilo - (1 - hi)*xihi, below),
rep(-lo*xilo - (1 - hi)*xihi, (above - below)),
rep(hi*xihi - lo*xilo, (length(x) - above)))
if(below != 0){
x[1:below] <- 0
}
if(above != length(x)){
x[(above + 1):length(x)] <- 0
}
phi <- (x + mat)*(1/(hi - lo)) - .trimmedmean(x, weight, lo, hi)$tm
return(phi)
}
.inflwinsor <-
function(x, weight, lo=0, hi=1){
# x <- sort(x)
if(is.null(weight)){
weight <- rep(1, length(x))
}
xilo <- .svyquantiles(x, weight, probs=lo)
xihi <- .svyquantiles(x, weight, probs=hi)
below <- floor(lo * length(x))
above <- ceiling(hi * length(x))
dens.lo <- .svydensity(x, weight, quant=lo)
dens.hi <- .svydensity(x, weight, quant=hi)
winlo <- lo - c(rep(1, below),rep(0, (length(x) - below)))
winhi <- hi - c(rep(1, above),rep(0, (length(x) - above)))
mat <- c(rep((1 - lo)*xilo - (1 - hi)*xihi, below),
rep(-lo*xilo - (1 - hi)*xihi, (above - below)),
rep(hi*xihi - lo*xilo, (length(x) - above)))
if(below != 0){
x[1:below] <- 0
}
if(above != length(x)){
x[(above + 1):length(x)] <- 0
}
phi1 <- (x + mat) - .winsorizedmean(x, weight, lo, hi)*(hi - lo)
if(dens.lo==0){
phi <- phi1 + (1-hi)*(winhi/dens.hi)
}
else{
phi <- phi1 + lo*(winlo/dens.lo) + (1-hi)*(winhi/dens.hi)
}
return(phi)
}
.irls <-
function(k, mat, acc, steps, design, beta.med, mer=FALSE, true.mean=NULL, asymmetric){
niter <- 0
beta.0 <- beta.med + 2 * acc
beta.r <- beta.med
x.sqrt <- sqrt(mat[, 3])
while ( (abs(beta.r - beta.0) > acc ) & ( niter < steps) ){
beta.0 <- beta.r
stand.res <- (mat[, 1] - beta.0 * mat[, 3]) / x.sqrt
scale <- .svyquantiles(abs(stand.res), mat[, 2]) * 1.4829
if (scale == 0) {
scale <- .svyquantiles(abs(stand.res), mat[, 2], probs=0.8) * 0.7803
cat("MAD is 0, trying quantile 0.8 of absolute deviations\n")
}
u <- .huberpsiw(stand.res / scale, k, w=TRUE, asymmetric=asymmetric)
beta.r <- sum(mat[, 2] * u * mat[, 1]) / sum(mat[, 2] * u * mat[, 3])
niter <- niter + 1
}
ures <- u * (mat[, 1] - beta.r * mat[, 3])
design1 <- update(design, ures)
se <- as.vector(sum(mat[, 2]) * sqrt(vcov(svymean(~ures, design1))) / sum(mat[, 2] * u * mat[, 3]))
if(mer){
mse <- as.vector(se^2 + (beta.r - true.mean)^2)
return(mse)
}
else{
result <- list(m=beta.r, se=se, ures=ures, u=u, niter=niter, scale=scale)
return(result)
}
}
.rsvycheck <-
function(x, design, na.rm=TRUE){
if(!inherits(design, "survey.design2")){
stop("Design must be an object of the 'survey.design (2)' class\n")
}
method <- ""
if(design$call[[1]] == "subset"){
method <- paste(method, " (for domains)", sep="")
}
original.design <- design
if(!inherits(x, "formula")){
stop("y must be a formula object\n")
}
else{
# ordering the data acording to the study variable
# and generate a new svydesign with the ordered values
ordereddata <- design$variables[order(design$variables[paste(x[2])]), ]
design$call$data <- substitute(ordereddata)
design <- eval(design$call)
# extract the variables from the svydesign -> xmat
mf <- model.frame(x, design$variables, na.action = na.pass)
xx <- lapply(attr(terms(x), "variables")[-1], function(tt) model.matrix(eval(bquote(~0 + .(tt))), mf))
cols <- sapply(xx, NCOL)
xmat <- matrix(nrow = NROW(xx[[1]]), ncol = sum(cols))
scols <- c(0, cumsum(cols))
for (i in 1:length(xx)) {
xmat[, scols[i] + 1:cols[i]] <- xx[[i]]
}
colnames(xmat) <- do.call("c", lapply(xx, colnames))
xmat <- cbind(xmat, as.matrix(weights(design)))
colnames(xmat)[NCOL(xmat)] <- "w"
}
# check if xmat has NAs
if(na.rm){
nas <- rowSums(is.na(xmat))
design <- design[nas == 0]
xmat <- xmat[nas == 0, , drop=F]
na <- length(which(nas!=0))
}
else{
na <- NA
}
if(any(xmat[,NCOL(xmat)] == 0)) cat("There are weights=0\n")
res <- list(mat=xmat, method=method, na=na, design=design)
class(res) <- "rsvycheck"
return(res)
}
.status <-
function(at, total){
if(total >= 30){
p.indicator <- diff(floor(30 * seq(1, total+1) / total))
if(at==1){
cat(" Progress \n")
cat(" +------+------+------+------+ \n")
cat(" 0% 25% 50% 75% 100% \n")
cat(" ")
}
else{
if(p.indicator[at]==1){
cat("*")
}
}
}
else{
repeating <- diff(floor(seq(1, 30, length.out=total+1)))
if(at==1){
cat(" Progress \n")
cat(" +------+------+------+------+ \n")
cat(" 0% 25% 50% 75% 100% \n")
cat(" ")
for(j in 1:repeating[at]){
cat("*")
}
}
else{
for(j in 1:repeating[at]){
cat("*")
}
}
}
if(at == total){
cat("\n ...done \n")
}
}
.svydensity <-
function(x, weight, quant=0.5, ngrid=401){
bwd <- dpik(x, scalest="minim", level=2, kernel="normal",
canonical=FALSE, gridsize=ngrid, range.x=range(x),
truncate=TRUE)
at <- seq(min(x), max(x), length=ngrid)
nx <- rowsum(c(rep(0,ngrid),weight), c(1:ngrid, findInterval(x,at)))
dens <-locpoly(rep(1,ngrid),nx*ngrid/(diff(range(x))*sum(weight)),
binned=TRUE, bandwidth=bwd, range.x=range(x))
res <- dens$y[floor(length(dens$y)*quant)]
if(length(res)==0) res <- 0
return(res)
}
.svyquantiles <-
function(x, w, ties=NULL, probs = 0.5, method="constant", f=0){
ord <- order(x)
x <- x[ord]
w <- w[ord]
if(is.null(ties)){
w.ord <- cumsum(w) / sum(w)
}
else{
w <- rowsum(w, ties, reorder = F)
x <- sort(aggregate(x, by=list(ties), unique)$x)
w.ord <- cumsum(w) / sum(w)
}
cdf <- approxfun(w.ord, x, method = method, f = f,
yleft = min(x), yright = max(x), rule=2)
cdf(probs)
}
.trimmedmean <-
function(x, weight, lo=0, hi=1){
xilo <- .svyquantiles(x, weight, probs=lo)
xihi <- .svyquantiles(x, weight, probs=hi)
w <- weight[x >= xilo & x <= xihi]
sw <- sum(w)
tm <- sum(w * x[x>=xilo & x<= xihi]) / sw
u <- rep(1, length(weight))
u[x < xilo] <- 0
u[x > xihi] <- 0
u <- u * (sum(weight) / sw)
res <- list(tm=tm, u=u)
return(res)
}
.updte <-
function(y, design, multiply){
variable.string <- paste(y[[2]])
myvar <- design$variables[, variable.string]
if (length(multiply) != length(myvar)) stop("dimension of 'multiply' is not correct!\n")
design$variables[, variable.string] <- myvar * multiply
return(design)
}
.winsorizedmean <-
function(x, weight, lo=0, hi=1){
xilo <- .svyquantiles(x, weight, probs=lo)
xihi <- .svyquantiles(x, weight, probs=hi)
x[x < xilo] <- xilo
x[x > xihi] <- xihi
res <- sum(weight*x) / sum(weight)
return(res)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.