1
2
3svyranktest<-function(formula,design,test=c('wilcoxon','vanderWaerden','median',"KruskalWallis"),...){
4        UseMethod("svyranktest", design)
5}
6
7svyranktest.survey.design<-svyranktest.svyrep.design<-function(formula, design,
8                                                               test=c('wilcoxon','vanderWaerden','median',"KruskalWallis"),...)
9{
10  mf<-model.frame(formula,model.frame(design),na.action=na.omit)
11  if (!is.null(naa<-attr(mf,"na.action"))){
12    design<-design[-naa,]
13    mf<-model.frame(formula,model.frame(design))
14  }
15  y<-mf[,1]
16  g<-mf[,2]
17
18  if (length(unique(g))!=2) {
19      return(multiranktest(formula,design, test,...))
20  }
21
22  if (is.character(test)) {
23    test<-match.arg(test)
24    testf<-switch(test, wilcoxon=,KruskalWallis=function(r,N) r/N,
25		  vanderWaerden=function(r,N) qnorm(r/N),
26		  median=function(r,N) as.numeric(r>N/2))
27  } else{
28    testf<-test
29  }
30  if (identical(test,"wilcoxon")) test<-"KruskalWallis"
31
32  ii<-order(y)
33  n<-length(y)
34  rankhat<-numeric(n)
35
36  w<-weights(design,"sampling")
37  ## calibrated designs have weights set to zero, not removed
38  na.fixup<-FALSE
39  if (is.calibrated(design) && length(w)>n && !is.null(naa)){
40      w<-w[-naa]
41      na.fixup<-TRUE
42  }
43
44
45  N<-sum(w)
46  rankhat[ii]<-ave(cumsum(w[ii])-w[ii]/2,factor(y[ii]))
47  rankscore<-testf(rankhat,N)
48  m <- lm(rankscore~g, weights=w)
49  delta<-coef(m)[2]
50  xmat<-model.matrix(m)
51
52  if (na.fixup){
53      infn<-matrix(0,nrow=nrow(xmat)+length(naa),ncol=ncol(xmat))
54      infn[-naa,]<-(xmat*(rankscore-fitted(m)))%*%summary(m)$cov.unscaled
55  } else {
56      infn<- (xmat*(rankscore-fitted(m)))%*%summary(m)$cov.unscaled
57  }
58  tot.infn<-svytotal(infn,design)
59  if (is.character(test))
60    method<-paste("Design-based",test,"test")
61  else if (!is.null(attr(test,"name")))
62    method<-paste("Design-based",attr(test,"name"),"test")
63  else method<-"Design-based rank test"
64
65  rval <- list(statistic = coef(m)[2]/SE(tot.infn)[2], parameter = degf(design) -
66               1, estimate = coef(m)[2], null.value = 0, alternative = "two.sided",
67               method = method, data.name = deparse(formula))
68  rval$p.value <- 2 * pt(-abs(rval$statistic), df = rval$parameter)
69  names(rval$statistic) <- "t"
70  names(rval$parameter) <- "df"
71  names(rval$estimate) <- "difference in mean rank score"
72  names(rval$null.value) <- "difference in mean rank score"
73  class(rval) <- "htest"
74  rval
75}
76
77multiranktest<-function(formula,design,test=c('wilcoxon','vanderWaerden','median','KruskalWallis'),...){
78    mf<-model.frame(formula,model.frame(design),na.action=na.omit)
79    if (!is.null(naa<-attr(mf,"na.action"))){
80        design<-design[-naa,]
81        #mf<-model.frame(formula,model.frame(design),na.action=na.fail)
82    }
83    y<-mf[,1]
84    g<-mf[,2]
85
86    if (is.character(test)) {
87        test<-match.arg(test)
88        testf<-switch(test, wilcoxon=,KruskalWallis=function(r,N) r/N,
89                      vanderWaerden=function(r,N) qnorm(r/N),
90                      median=function(r,N) as.numeric(r>N/2))
91    } else{
92        testf<-test
93    }
94    if (identical(test,"wilcoxon")) test<-"KruskalWallis"
95
96    ii<-order(y)
97    n<-length(y)
98    rankhat<-numeric(n)
99    w<-weights(design,"sampling")
100     ## calibrated designs have weights set to zero, not removed
101    if (is.calibrated(design) && length(w)>n && !is.null(naa)){
102      w<-w[-naa]
103     }
104
105    N<-sum(w)
106    rankhat[ii]<-ave(cumsum(w[ii])-w[ii]/2,factor(y[ii]))
107    rankscore<-testf(rankhat,N)
108    m <- glm(rankscore~factor(g),weights=w)
109    m$na.action<-naa
110    V<-svy.varcoef(m,design)
111    ndf<-length(unique(g))-1
112    beta<-coef(m)[-1]
113    V<-V[-1,-1]
114    chisq<-beta%*%solve(V,beta)
115    ddf<-degf(design)-ndf
116    if (is.character(test))
117        method<-paste("Design-based",test,"test")
118    else if (!is.null(attr(test,"name")))
119        method<-paste("Design-based",attr(test,"name"),"test")
120    else method<-"Design-based rank test"
121    names(chisq)<-"Chisq"
122    names(ndf)<-"df"
123    rval<-list(parameter=chisq,statistic=ndf,ddf=ddf,p.value=pf(chisq/ndf,ndf,ddf,lower.tail=FALSE),
124               method=method, data.name = deparse(formula))
125    class(rval)<-"htest"
126    rval
127}
128