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