######################################################################
# survey mean function for non-replicated designs that also prints design effects
# G Raab may 04
#
 svymean.deff<-function(model,design){
#
# calculates a design effect 
# from a survey estimate of the population variance
#
 res<-svymean(model,design)
 varsrs<-svyvar(model,design)/dim(design$variables)[1]
print(sqrt(varsrs))
 deff<-attr(res,"var")/varsrs
 print(cat("Design effect",deff,"\n"))
  res
}
######################################################################
# survey mean function for replicated designs that also prints design effects
# G Raab may 04
#
svrepmean.deff<-function(model,design){
#
# calculates a design effect for a replicated design
# using an estimate of the population variance
# from a non-replicated survey
#
 res<-svrepmean(model,design)
 varsrs<-svrepvar(model,design)/dim(design$variables)[1]
 deff<-attr(res,"var")/varsrs
 print(cat("Design effect",deff,"\n"))
  res
}
##############################################################################
#
#  alternative survey quantile function
#  G Raab may 04
#
my.svyquantile<-function (formula = ~ANNETINC, design = hbai.des1, quantiles = c(0.1, 
    0.5, 0.95), size = 0.95) 
{
    if (size <= 0 || size >= 1) 
        stop("Size needs to be between 0 and 1 to give a valid conf interval")
    if (any(quantiles <= 0) || any(quantiles >= 1)) 
        stop("percentile points need to be between 0 and 1")
    x <- formula
    if (inherits(x, "formula")) 
        x <- model.frame(x, design$variables)
    else if (typeof(x) %in% c("expression", "symbol")) 
        x <- eval(x, design$variables)
    w <- weights(design)
    computeQuantiles <- function(xx, qq) {
        oo <- order(xx)
        cum.w <- cumsum(w[oo])/sum(w)
        cdf <- approxfun(cum.w, xx[oo], method = "linear", f = 1, 
            yleft = min(xx), yright = max(xx))
        cdf(qq)
    }
    if (!is.null(dim(x)) & dim(x)[2] > 1) 
        stop("Only does s.e.s for one variable at a time")
    if (!is.null(dim(x))) 
        x <- x[, 1]
    Quantile <- computeQuantiles(x, quantiles)
    getpse <- function(rv) {
        desn <- update(design, pct = 1 * (x < rv))
        if (sum(x < rv) < 10 || sum(x > rv) < 10) 
            cat("WARNING!!!!! Sample numbers too small for approximation to be reliable at value", 
                rv, "\n")
        se <- sqrt(attr(svymean(~pct, desn), "var"))
    }
    sep <- sapply(Quantile, getpse)
    q.minus.1se <- computeQuantiles(x, quantiles - sep)
    q.plus.1se <- computeQuantiles(x, quantiles + sep)
    se <- (q.plus.1se - q.minus.1se)/2
    l.limit <- computeQuantiles(x, quantiles + qnorm((1 - size)/2) * 
        sep)
    u.limit <- computeQuantiles(x, quantiles + qnorm(1 - (1 - 
        size)/2) * sep)
    cbind(quantiles, Quantile, se, l.limit, u.limit)
}


######################################################
# alternative survey quantile function
# may 04
#
my.svrepquantile<-function (formula = ~ANNETINC, design = hbai.des1, quantiles = c(0.1, 
    0.5, 0.95), size = 0.95) 
{
  if (!inherits(design, "svyrep.design")) 
        stop("Not a survey replicates object")    
if (size <= 0 || size >= 1) 
        stop("Size needs to be between 0 and 1 to give a valid conf interval")
    if (any(quantiles <= 0) || any(quantiles >= 1)) 
        stop("percentile points need to be between 0 and 1")
    x <- formula
    if (inherits(x, "formula")) 
        x <- model.frame(x, design$variables)
    else if (typeof(x) %in% c("expression", "symbol")) 
        x <- eval(x, design$variables)
    w <- design$pweights

    computeQuantiles <- function(xx, qq) {
        oo <- order(xx)
        cum.w <- cumsum(w[oo])/sum(w)
        cdf <- approxfun(cum.w, xx[oo], method = "linear", f = 1, 
            yleft = min(xx), yright = max(xx))
        cdf(qq)
    }
    if (!is.null(dim(x)) & dim(x)[2] > 1) 
        stop("Only does s.e.s for one variable at a time")
    if (!is.null(dim(x))) 
        x <- x[, 1]

    Quantile <- computeQuantiles(x, quantiles)
    getpse <- function(rv) {
        desn <- update(design, pct = 1 * (x < rv))
        if (sum(x < rv) < 10 || sum(x > rv) < 10) 
            cat("WARNING!!!!! Sample numbers too small for approximation to be reliable at value", 
                rv, "\n")
        se <- sqrt(attr(svrepmean(~pct, desn), "var"))
    }
    sep <- sapply(Quantile, getpse)
    q.minus.1se <- computeQuantiles(x, quantiles - sep)
    q.plus.1se <- computeQuantiles(x, quantiles + sep)
    se <- (q.plus.1se - q.minus.1se)/2
    l.limit <- computeQuantiles(x, quantiles + qnorm((1 - size)/2) * 
        sep)
    u.limit <- computeQuantiles(x, quantiles + qnorm(1 - (1 - 
        size)/2) * sep)
    cbind(quantiles, Quantile, se, l.limit, u.limit)
}