Dear all,
   I have the following program for a multiple comparison procedure.
There are two functions for the two steps. First step is to calculate the
critical values,
while the second step is the actual procedure [see below: program with two
functions].
   This  work fine. However, However I want to put them into one function
for the convenience
of later use [see below: program with one function]. Some how the big
function works extremely
slow.  For example I chose
m <- 10
rho <- 0.1
k <- 2
alpha <- 0.05
pvaluessort <- sort(1-pnorm(rnorm(10))

Is there anything that I did wrong?
 Thank you!
                    Hannah


######Program with two functions############
## first step
library(mvtnorm)
cc_f <- function(m, rho, k, alpha) {

cc_z <- numeric(m)
var <- matrix(c(1,rho,rho,1), nrow=2, ncol=2, byrow=T)
for (i in 1:m){
   if (i <= k) {cc_z[i] <- qmvnorm((k*(k-1))/(m*(m-1))*alpha, tail="upper",
sigma=var)$quantile} else
               {cc_z[i] <- qmvnorm((k*(k-1))/((m-i+k)*(m-i+k-1))*alpha,
tail="upper", sigma=var)$quantile}
               }
cc <- 1-pnorm(cc_z)
return(cc)
                             }
## second step
pair_kFWER=function(m,crit_cons,pvaluessort){
k=0
while((k <m)&&(crit_cons[m-k] < pvaluessort[m-k])){
k=k+1
}
return(m-k)
}
###########################################
##############Program with one function ##############

pair_kFWER=function(m,alpha,rho,k,pvaluessort){
library(mvtnorm)
cc_z <- numeric(m)
var <- matrix(c(1,rho,rho,1), nrow=2, ncol=2, byrow=T)
for (i in 1:m){
   if (i <= k) {cc_z[i] <- qmvnorm((k*(k-1))/(m*(m-1))*alpha, tail="upper",
sigma=var)$quantile} else
               {cc_z[i] <- qmvnorm((k*(k-1))/((m-i+k)*(m-i+k-1))*alpha,
tail="upper", sigma=var)$quantile}
               }
crit_cons <- 1-pnorm(cc_z)

k=0
while((k <m)&&(crit_cons[m-k] < pvaluessort[m-k])){
k=k+1
}
return(m-k)
}
#####################################################

        [[alternative HTML version deleted]]

______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Reply via email to