1svyquantile<-function(x,design,quantiles,...) UseMethod("svyquantile", design)
2
3svyquantile.survey.design <- function (x, design, quantiles, alpha = 0.05,
4                                       interval.type = c("mean", "beta","xlogit", "asin","score"),
5                                       na.rm = FALSE,
6                                       ci=TRUE, se = ci,
7                                       qrule=c("math","school","shahvaish","hf1","hf2","hf3","hf4","hf5","hf6","hf7","hf8","hf9"),
8                                       df = NULL, ...) {
9
10   if (inherits(x, "formula"))
11        x <- model.frame(x, model.frame(design), na.action = na.pass)
12    else if (typeof(x) %in% c("expression", "symbol"))
13        x <- eval(x, model.frame(design, na.action = na.pass))
14    if (na.rm) {
15        nas <- rowSums(is.na(x))
16        design <- design[nas == 0, ]
17        if (length(nas) > length(design$prob))
18            x <- x[nas == 0, , drop = FALSE]
19        else x[nas > 0, ] <- 0
20    }
21
22    if(is.null(df))
23        df<-degf(design)
24
25    if (is.character(qrule)){
26        qrulename<-paste("qrule",match.arg(qrule),sep="_")
27        qrule<-get(qrulename, mode="function")
28    }
29
30    qcrit<-if(df==Inf) qnorm else function(...) qt(...,df=df)
31
32
33    interval.type<-match.arg(interval.type)
34    if (interval.type=="score"){
35        ci_fun<-ffullerCI
36    } else {
37        ci_fun<-function(...) woodruffCI(...,method=interval.type)
38    }
39
40    w<-weights(design)
41
42    rvals<-lapply(x,
43                  function(xi){
44                      r<-t(sapply(quantiles,
45                                  function(p){
46                                      qhat<-qrule(xi,w,p)
47                                      if(ci){
48                                          ci<-ci_fun(xi,qhat,p,design,qrule,alpha,df)
49                                          names(ci)<-c(round(100*alpha/2,2),round(100-100*alpha/2,2))
50                                          c(quantile=qhat,ci=ci,
51                                            se=as.numeric(diff(ci)/(2*qcrit(1-alpha/2)))
52                                            )
53                                      } else {
54                                          c(quantile=qhat)
55                                      }
56                                  }
57                                  ))
58                      if(!ci)
59                          colnames(r)<-quantiles
60                      else
61                          rownames(r)<-quantiles
62                      r
63                  }
64                  )
65    attr(rvals,"hasci")<-ci
66    class(rvals)<-"newsvyquantile"
67    rvals
68}
69
70svyquantile.svyrep.design <- function (x, design, quantiles, alpha = 0.05,
71                                       interval.type = c("mean", "beta","xlogit", "asin","quantile"),
72                                       na.rm = FALSE,
73                                       ci=TRUE, se = ci,
74                                       qrule=c("math","school","shahvaish","hf1","hf2","hf3","hf4","hf5","hf6","hf7","hf8","hf9"),
75                                       df = NULL, return.replicates=FALSE,...) {
76    interval.type <- match.arg(interval.type)
77
78    if (design$type %in% c("JK1", "JKn") && interval.type == "quantile")
79        warning("Jackknife replicate weights may not give valid standard errors for quantiles")
80    if (design$type %in% "other" && interval.type == "quantile")
81        message("Not all replicate weight designs give valid standard errors for quantiles.")
82
83    if (inherits(x, "formula"))
84        x <- model.frame(x, design$variables, na.action = if (na.rm)
85            na.pass
86        else na.fail)
87    else if (typeof(x) %in% c("expression", "symbol"))
88        x <- eval(x, design$variables)
89    if (na.rm) {
90        nas <- rowSums(is.na(x))
91        design <- design[nas == 0, ]
92        if (length(nas) > length(design$prob))
93            x <- x[nas == 0, , drop = FALSE]
94        else x[nas > 0, ] <- 0
95    }
96
97    if (is.character(qrule)){
98        qrulename<-paste("qrule",match.arg(qrule),sep="_")
99        qrule<-get(qrulename, mode="function")
100    }
101
102    if (is.null(df)) df <- degf(design)
103    if (df == Inf)  qcrit <- qnorm else qcrit <- function(...) qt(..., df = df)
104
105    w<-weights(design,"sampling")
106
107    if (interval.type=="quantile"){
108        ci_fun<-function(...) repCI(...,return.replicates=return.replicates)
109    } else {
110        if(return.replicates) warning("return.replicates=TRUE only implemented for interval.type='quantile'")
111        ci_fun<-woodruffCI
112    }
113
114    if ((interval.type=="quantile") && return.replicates){
115              rvals<-lapply(x,
116                      function(xi){ lapply(quantiles,
117                                           function(p){
118                                               qhat<-qrule(xi,w,p)
119                                               ci<-ci_fun(xi,qhat,p,design,qrule,alpha,df)
120                                               list(quantile=qhat, replicates=attr(ci,"replicates"))
121                                           }
122                                           )
123                      }
124                      )
125              ests<-sapply(rvals, function(v) sapply(v, function(qi) qi$qhat))
126              attr(ests, "scale") <- design$scale
127              attr(ests, "rscales") <- design$rscales
128              attr(ests, "mse") <- design$mse
129              reps<-sapply(rvals, function(v) t(sapply(v, function(qi) qi$replicates)))
130              rval<-list(quantile=ests,replicates=reps)
131
132              attr(rval,"var")<-svrVar(reps, design$scale, design$rscales, mse=design$mse, coef=ests)
133              class(rval)<-"svrepstat"
134              return(rval)
135    } else {
136        rvals<-lapply(x,
137                      function(xi){
138                          r<- t(sapply(quantiles,
139                                       function(p){
140                                           qhat<-qrule(xi,w,p)
141                                           if(ci){
142                                               ci<-ci_fun(xi,qhat,p,design,qrule,alpha,df)
143                                               names(ci)<-c(round(100*alpha/2,2),round(100-100*alpha/2,2))
144                                               c(quantile=qhat,ci=ci,
145                                                 se=diff(ci)/(2*qcrit(1-alpha/2))
146                                                 )} else{
147                                                      c(quantile=qhat)
148                                                  }
149
150                                       }
151                                       ))
152                          if(!ci)
153                              colnames(r)<-quantiles
154                          else
155                              rownames(r)<-quantiles
156                          r
157                      }
158                      )
159        attr(rvals,"hasci")<-ci | se
160        class(rvals)<-"newsvyquantile"
161    }
162
163    rvals
164}
165
166SE.newsvyquantile<-function(object,...) {
167    if(!attr(object,"hasci")) stop("object does not have uncertainty estimates")
168    do.call(c,lapply(object,function(ai) ai[,4]))
169}
170
171vcov.newsvyquantile<-function(object,...) {
172    if(!attr(object,"hasci")) stop("object does not have uncertainty estimates")
173    r<-do.call(c,lapply(object,function(ai) ai[,4]))^2
174    v<-matrix(NA,nrow=length(r),ncol=length(r))
175    diag(v)<-r
176    v
177}
178
179coef.newsvyquantile<-function(object,...){
180    if(attr(object,"hasci"))
181        do.call(c,lapply(object,function(ai) ai[,1]))
182    else
183        do.call(c,lapply(object,function(ai) ai[1,]))
184}
185
186
187
188
189confint.newsvyquantile<-function(object,...){
190    if(!attr(object,"hasci")) stop("object does not have uncertainty estimates")
191    l<-do.call(c,lapply(object,function(ai) ai[,2]))
192    u<-do.call(c,lapply(object,function(ai) ai[,3]))
193    cbind(l,u)
194}
195
196
197woodruffCI<-function(x, qhat,p, design, qrule,alpha,df,method=c("mean","beta","xlogit","asin")){
198
199    method<-match.arg(method)
200    m<-svymean(x<=qhat, design)
201    names(m)<-"pmed"
202
203    pconfint<-switch(method,
204                    mean=as.vector(confint(m, 1, level = 1-alpha, df = df)),
205                    xlogit= {xform <- svycontrast(m, quote(log(`pmed`/(1 - `pmed`))));
206                        expit(as.vector(confint(xform, 1, level = 1-alpha,  df = df)))},
207                    beta={n.eff <- coef(m) * (1 - coef(m))/vcov(m);
208                        rval <- coef(m)[1]
209                        n.eff <- n.eff * (qt(alpha/2, nrow(design) - 1)/qt(alpha/2, degf(design)))^2
210                        c(qbeta(alpha/2, n.eff * rval, n.eff * (1 - rval) +  1), qbeta(1 - alpha/2, n.eff * rval + 1, n.eff *  (1 - rval)))
211                    },
212                    asin={xform <- svycontrast(m, quote(asin(sqrt(`pmed`))))
213                        sin(as.vector(confint(xform, 1, level = 1-alpha,  df = df)))^2
214                    }
215                    )
216    lower<-if(is.nan(pconfint[1]) || (pconfint[1]<0)) NaN else qrule(x,weights(design,"sampling"), pconfint[1])
217    upper<-if(is.nan(pconfint[2])|| (pconfint[2]>1)) NaN else qrule(x,weights(design,"sampling"), pconfint[2])
218    rval<-c(lower, upper)
219    names(rval)<-c(round(100*alpha/2,1),round(100*(1-alpha/2),1))
220    rval
221
222}
223
224
225
226ffullerCI<-function(x, qhat, p, design, qrule, alpha, df){
227        qcrit<-if(df==Inf) qnorm else function(...) qt(...,df=df)
228        U <- function(theta) {
229            ((x > theta) - (1 - p))
230        }
231        scoretest <- function(theta, qlimit) {
232            umean <- svymean(U(theta), design)
233            umean/SE(umean) - qlimit
234        }
235        iqr <- IQR(x)
236        lowerT <- min(x) + iqr/100
237        upperT <- max(x) - iqr/100
238        tol <- 1/(100 * sqrt(nrow(design)))
239
240        qlow<- uniroot(scoretest, interval = c(lowerT, upperT), qlimit = qcrit(alpha/2, lower.tail = FALSE), tol = tol)$root
241        qup<-uniroot(scoretest, interval = c(lowerT, upperT), qlimit = qcrit(alpha/2, lower.tail = TRUE), tol = tol)$root
242
243        w<-weights(design)
244
245        c(qrule(x, w, mean(x<=qlow)), qrule(x, w, mean(x<=qup)))
246
247    }
248
249
250repCI<-function(x, qhat, p, design, qrule, alpha, df,return.replicates){
251    qcrit<-if(df==Inf) qnorm else function(...) qt(...,df=df)
252
253    wrep<-weights(design,"analysis")
254    reps<-apply(wrep,2, function(wi) qrule(x,wi,p))
255    v<-with(design, svrVar(reps, scale=scale, rscales=rscales, mse=mse,coef=qhat))
256
257    ci<- qhat+ c(-1,1)*sqrt(v)*qcrit(1-alpha/2)
258
259    if (return.replicates) attr(ci,"replicates")<-reps
260
261    ci
262}
263