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]]
use Rprof to determine where your function is spending time. What is the problem you are trying to solve? Sent from my iPhone. On Jun 23, 2010, at 5:21, li li <hannah.hlx at gmail.com> wrote:> 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 at 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.
Most of the computation time is in the functions qvnorm. You can win a
little bit by optimizing the code, but the gain is relatively small.
You can also decrease the interval used to evaluate qvnorm to win some
speed there. As you look for the upper tail, no need to evaluate the
function in negative numbers. Eg :
pair_kFWER2=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)
lpart <- 1:k # first part
hpart <- (k+1):m # second part
# make first part of the cc_z
cc_z[lpart] <- qmvnorm((k*(k-1))/(m*(m-1))*alpha, tail="upper",
interval=c(0,4),sigma=Var)$quantile
# make second part of the cc_z
cc_z[hpart] <- sapply(hpart,function(i){
qmvnorm((k*(k-1))/((m-i+k)*(m-i+k-1))*alpha,interval=c(0,4),
tail="upper", sigma=Var)$quantile
})
crit_cons <- 1-pnorm(cc_z)
# does the test whether values of crit_cons[i] are smaller than
# values pvaluessort[i] and returns a vector
# which says at which positions does not occur
# tail takes the last position. I guess that's what you wanted
out <- tail(which(!(crit_cons < pvaluessort)),1)
return(out)
}
> system.time(replicate(100,pair_kFWER(m,alpha,rho,k,pvaluessort)))
user system elapsed
5.97 0.04 6.04
> system.time(replicate(100,pair_kFWER2(m,alpha,rho,k,pvaluessort)))
user system elapsed
4.43 0.00 4.44
Check whether the function works correctly. It gives the same result
as the other one in my tests. Only in the case where your function
returns 0, the modified returns integer(0)
Cheers
Joris
On Wed, Jun 23, 2010 at 2:21 PM, li li <hannah.hlx at gmail.com>
wrote:> 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 at 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.
>
--
Joris Meys
Statistical consultant
Ghent University
Faculty of Bioscience Engineering
Department of Applied mathematics, biometrics and process control
tel : +32 9 264 59 87
Joris.Meys at Ugent.be
-------------------------------
Disclaimer : http://helpdesk.ugent.be/e-maildisclaimer.php