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