1svysmooth<-function(formula,design,...) UseMethod("svysmooth", design) 2svysmooth.default<-function(formula, design,method=c("locpoly","quantreg"),bandwidth=NULL,quantile,df=4,...){ 3 switch(match.arg(method), 4 locpoly=svylocpoly(formula,design,bandwidth=bandwidth,...), 5 quantreg=svyrqss(formula,design,quantile=quantile,df=df,...) 6 ) 7} 8 9fitted.rq<-function(object,...) object$x%*% object$coefficients/object$weights 10 11svyrqss<-function(formula,design,quantile=0.5,df=4,...){ 12 mf<-model.frame(formula, model.frame(design), na.action=na.omit) 13 naa<-attr(mf,"na.action") 14 15 tt<-attr(terms(formula),"term.labels") 16 df<-rep(df, length=length(tt)) 17 quantile<-rep(quantile, length=length(tt)) 18 19 if (length(formula)==3){ 20 density<-FALSE 21 } else { 22 density<-TRUE 23 stop("type='quantreg' does not do densities") 24 } 25 26 w<-weights(design,type="sampling") 27 if (length(naa)) w<-w[-naa] 28 environment(formula)<-environment() 29 30 ll<-vector("list", length(tt)) 31 for(i in 1:length(tt)){ 32 termi<-as.name(tt[i]) 33 ff<-eval(bquote(update(formula,.~splines::bs(.(termi),df=.(df[i]))))) 34 rqfit<-quantreg::rq(ff, tau=quantile[i],weights=w,data=mf,...) 35 xx<-mf[,i+1] 36 oo<-order(xx) 37 ll[[i]]<-list(x=xx[oo],y=fitted.rq(rqfit)[oo]) 38 } 39 names(ll)<-attr(terms(formula),"term.labels") 40 41 attr(ll,"call")<-sys.call(-2) 42 attr(ll,"density")<-density 43 if(density) 44 attr(ll,"ylab")<-"Density" 45 else 46 attr(ll,"ylab")<-paste(deparse(formula[[2]]),collapse="") 47 48 class(ll)<-"svysmooth" 49 50 ll 51 52} 53 54svylocpoly<-function(formula, design, ngrid=401, xlim=NULL, 55 ylim=NULL, bandwidth=NULL,...){ 56 57 mf<-model.frame(formula,model.frame(design)) 58 mm<-model.matrix(terms(formula),mf) 59 if(attr(terms(formula),"intercept")) 60 mm<-mm[,-1,drop=FALSE] 61 62 naa<-attr(mf,"na.action") 63 64 65 if (length(formula)==3){ 66 Y<-model.response(mf) 67 density<-FALSE 68 } else density<-TRUE 69 70 71 if (is.null(xlim)){ 72 xlim<-apply(mm,2,range) 73 } 74 if (!is.matrix(xlim)) 75 xlim<-matrix(xlim,nrow=2) 76 77 78 if (is.null(bandwidth)){ 79 bandwidth<-numeric(ncol(mm)) 80 for(i in 1:ncol(mm)){ 81 bandwidth[i]<-if(density) KernSmooth::dpik(mm[,i],gridsize=ngrid) else KernSmooth::dpill(mm[,i],Y,gridsize=ngrid) 82 } 83 } else { 84 bandwidth<-rep(bandwidth, length=ncol(mm)) 85 } 86 87 w<-weights(design,type="sampling") 88 if (length(naa)) w<-w[-naa] 89 90 ll<-vector("list", ncol(mm)) 91 for(i in 1:NCOL(mm)){ 92 gx<-seq(min(xlim[,i]), max(xlim[,i]), length=ngrid) 93 nx<-rowsum(c(rep(0,ngrid),w), c(1:ngrid, findInterval(mm[,i],gx))) 94 if (density){ 95 ll[[i]]<-KernSmooth::locpoly(rep(1,ngrid),nx*ngrid/(diff(xlim[,i])*sum(w)), 96 binned=TRUE, bandwidth=bandwidth[i], range.x=xlim[,i]) 97 }else{ 98 ny<-rowsum(c(rep(0,ngrid), Y*w), c(1:ngrid, findInterval(mm[,i],gx))) 99 ll[[i]]<-KernSmooth::locpoly(nx, ny, binned=TRUE, bandwidth=bandwidth[i], range.x=xlim[,i]) 100 } 101 names(ll)<-attr(terms(formula),"term.labels") 102 } 103 attr(ll,"call")<-sys.call(-2) 104 attr(ll,"density")<-density 105 if(density) 106 attr(ll,"ylab")<-"Density" 107 else 108 attr(ll,"ylab")<-paste(deparse(formula[[2]]),collapse="") 109 110 class(ll)<-"svysmooth" 111 112 ll 113 114} 115 116print.svysmooth<-function(x,...){ 117 if(attr(x,"density")) 118 cat("Density estimate: :") 119 else 120 cat("Scatterplot smoother :") 121 print(attr(x,"call")) 122 invisible(x) 123} 124 125 126 127plot.svysmooth<-function(x, which=NULL,type="l",xlabs=NULL,ylab=NULL,...){ 128 if (is.null(which)) 129 which<-seq(length=length(x)) 130 if (is.character(which)) 131 which<-match(which,names(x)) 132 133 if(is.null(xlabs)) xlabs<-names(x)[which] 134 if(is.null(ylab)) ylab<-attr(x,"ylab") 135 136 for(i in seq(length=length(which))) 137 plot(x[[which[i]]], type=type, xlab=xlabs[i], ylab=ylab, ...) 138 invisible(NULL) 139} 140 141lines.svysmooth<-function(x,which=NULL,...){ 142 for(i in names(x)) lines(x[[i]],...) 143} 144