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