thr3ads.net - R help - [R] Pearson corelation and p-value for matrix [Apr 2005]

If this information is useful, please help other people find it:
Share via:

Dren Scott

2005-Apr-15 21:33 UTC

[R] Pearson corelation and p-value for matrix

Hi,
 
I was trying to evaluate the pearson correlation and the p-values for an nxm
matrix, where each row represents a vector. One way to do it would be to iterate
through each row, and find its correlation value( and the p-value) with respect
to the other rows. Is there some function by which I can use the matrix as
input? Ideally, the output would be an nxn matrix, containing the p-values
between the respective vectors.
 
I have tried cor.test for the iterations, but couldn't find a function that
would take the matrix as input.
 
Thanks for the help.
 
Dren
 
 

		
---------------------------------


	[[alternative HTML version deleted]]

John Fox

2005-Apr-15 22:26 UTC

head link

[R] Pearson corelation and p-value for matrix

Dear Dren,

How about the following?

 cor.pvalues <- function(X){
    nc <- ncol(X)
    res <- matrix(0, nc, nc)
    for (i in 2:nc){
        for (j in 1:(i - 1)){
            res[i, j] <- res[j, i] <- cor.test(X[,i], X[,j])$p.value
            }
        }
    res
    }

What one then does with all of those non-independent test is another
question, I guess.

I hope this helps,
 John

--------------------------------
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
-------------------------------- 
> -----Original Message-----
> From: r-help-bounces at stat.math.ethz.ch 
> [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Dren Scott
> Sent: Friday, April 15, 2005 4:33 PM
> To: r-help at stat.math.ethz.ch
> Subject: [R] Pearson corelation and p-value for matrix
> 
> Hi,
>  
> I was trying to evaluate the pearson correlation and the 
> p-values for an nxm matrix, where each row represents a 
> vector. One way to do it would be to iterate through each 
> row, and find its correlation value( and the p-value) with 
> respect to the other rows. Is there some function by which I 
> can use the matrix as input? Ideally, the output would be an 
> nxn matrix, containing the p-values between the respective vectors.
>  
> I have tried cor.test for the iterations, but couldn't find a 
> function that would take the matrix as input.
>  
> Thanks for the help.
>  
> Dren
>  
>  
> 
> 		
> ---------------------------------
> 
> 
> 	[[alternative HTML version deleted]]
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! 
> http://www.R-project.org/posting-guide.html

Liaw, Andy

2005-Apr-16 02:51 UTC

head link

[R] Pearson corelation and p-value for matrix

We can be a bit sneaky and `borrow' code from cor.test.default:

cor.pval <- function(x,  alternative="two-sided", ...) {
    corMat <- cor(x, ...)
    n <- nrow(x)
    df <- n - 2
    STATISTIC <- sqrt(df) * corMat / sqrt(1 - corMat^2)
    p <- pt(STATISTIC, df)
    p <- if (alternative == "less") {
        p
    } else if (alternative == "greater") {
        1 - p
    } else 2 * pmin(p, 1 - p)
    p
}

The test:
> system.time(c1 <- cor.pvals(X), gcFirst=TRUE)
[1] 13.19  0.01 13.58    NA    NA> system.time(c2 <- cor.pvalues(X), gcFirst=TRUE)
[1] 6.22 0.00 6.42   NA   NA> system.time(c3 <- cor.pval(X), gcFirst=TRUE)[1] 0.07 0.00 0.07   NA   NA

Cheers,
Andy
> From: John Fox
> 
> Dear Mark,
> 
> I think that the reflex of trying to avoid loops in R is 
> often mistaken, and
> so I decided to try to time the two approaches (on a 3GHz Windows XP
> system).
> 
> I discovered, first, that there is a bug in your function -- 
> you appear to
> have indexed rows instead of columns; fixing that:
> 
> cor.pvals <- function(mat)
> {
>   cols <- expand.grid(1:ncol(mat), 1:ncol(mat))
>   matrix(apply(cols, 1,
>                function(x) cor.test(mat[, x[1]], mat[, 
> x[2]])$p.value),
>          ncol = ncol(mat))
> }
> 
> 
> My function is cor.pvalues and yours cor.pvals. This is for a 
> data matrix
> with 1000 observations on 100 variables:
> 
> > R <- diag(100)
> > R[upper.tri(R)] <- R[lower.tri(R)] <- .5
> > library(mvtnorm)
> > X <- rmvnorm(1000, sigma=R)
> > dim(X)
> [1] 1000  100
> > 
> > system.time(cor.pvalues(X))
> [1] 5.53 0.00 5.53   NA   NA
> > 
> > system.time(cor.pvals(X))
> [1] 12.66  0.00 12.66    NA    NA
> > 
> 
> I frankly didn't expect the advantage of my approach to be 
> this large, but
> there it is.
> 
> Regards,
>  John
> 
> --------------------------------
> John Fox
> Department of Sociology
> McMaster University
> Hamilton, Ontario
> Canada L8S 4M4
> 905-525-9140x23604
> http://socserv.mcmaster.ca/jfox 
> -------------------------------- 
> 
> > -----Original Message-----
> > From: Marc Schwartz [mailto:MSchwartz at MedAnalytics.com] 
> > Sent: Friday, April 15, 2005 7:08 PM
> > To: John Fox
> > Cc: 'Dren Scott'; R-Help
> > Subject: RE: [R] Pearson corelation and p-value for matrix
> > 
> > Here is what might be a slightly more efficient way to get to
John's
> > question:
> > 
> > cor.pvals <- function(mat)
> > {
> >   rows <- expand.grid(1:nrow(mat), 1:nrow(mat))
> >   matrix(apply(rows, 1,
> >                function(x) cor.test(mat[x[1], ], mat[x[2], 
> > ])$p.value),
> >          ncol = nrow(mat))
> > }
> > 
> > HTH,
> > 
> > Marc Schwartz
> > 
> > On Fri, 2005-04-15 at 18:26 -0400, John Fox wrote:
> > > Dear Dren,
> > > 
> > > How about the following?
> > > 
> > >  cor.pvalues <- function(X){
> > >     nc <- ncol(X)
> > >     res <- matrix(0, nc, nc)
> > >     for (i in 2:nc){
> > >         for (j in 1:(i - 1)){
> > >             res[i, j] <- res[j, i] <- cor.test(X[,i], 
> X[,j])$p.value
> > >             }
> > >         }
> > >     res
> > >     }
> > > 
> > > What one then does with all of those non-independent test 
> > is another 
> > > question, I guess.
> > > 
> > > I hope this helps,
> > >  John
> > 
> > > > -----Original Message-----
> > > > From: r-help-bounces at stat.math.ethz.ch 
> > > > [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of 
> Dren Scott
> > > > Sent: Friday, April 15, 2005 4:33 PM
> > > > To: r-help at stat.math.ethz.ch
> > > > Subject: [R] Pearson corelation and p-value for matrix
> > > > 
> > > > Hi,
> > > >  
> > > > I was trying to evaluate the pearson correlation and 
> the p-values 
> > > > for an nxm matrix, where each row represents a vector. 
> > One way to do 
> > > > it would be to iterate through each row, and find its 
> correlation 
> > > > value( and the p-value) with respect to the other rows. 
> Is there 
> > > > some function by which I can use the matrix as input? 
> > Ideally, the 
> > > > output would be an nxn matrix, containing the p-values 
> > between the 
> > > > respective vectors.
> > > >  
> > > > I have tried cor.test for the iterations, but couldn't
find a
> > > > function that would take the matrix as input.
> > > >  
> > > > Thanks for the help.
> > > >  
> > > > Dren
> > 
> >
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! 
> http://www.R-project.org/posting-guide.html
> 
> 
>

Liaw, Andy

2005-Apr-16 12:51 UTC

head link

[R] Pearson corelation and p-value for matrix

> From: John Fox 
> 
> Dear Andy,
> 
> That's clearly much better -- and illustrates an effective 
> strategy for
> vectorizing (or "matricizing") a computation. I think I'll 
> add this to my
> list of programming examples. It might be a little dangerous 
> to pass ...
> through to cor(), since someone could specify 
> type="spearman", for example.
Ah, yes, the "..." isn't likely to help here!  Also, it will 
only work correctly if there are no NA's, for example (or 
else the degree of freedom would be wrong).  

Best,
Andy
 > Thanks,
>  John
> 
> --------------------------------
> John Fox
> Department of Sociology
> McMaster University
> Hamilton, Ontario
> Canada L8S 4M4
> 905-525-9140x23604
> http://socserv.mcmaster.ca/jfox 
> -------------------------------- 
> 
> > -----Original Message-----
> > From: Liaw, Andy [mailto:andy_liaw at merck.com] 
> > Sent: Friday, April 15, 2005 9:51 PM
> > To: 'John Fox'; MSchwartz at medanalytics.com
> > Cc: 'R-Help'; 'Dren Scott'
> > Subject: RE: [R] Pearson corelation and p-value for matrix
> > 
> > We can be a bit sneaky and `borrow' code from cor.test.default:
> > 
> > cor.pval <- function(x,  alternative="two-sided", ...) {
> >     corMat <- cor(x, ...)
> >     n <- nrow(x)
> >     df <- n - 2
> >     STATISTIC <- sqrt(df) * corMat / sqrt(1 - corMat^2)
> >     p <- pt(STATISTIC, df)
> >     p <- if (alternative == "less") {
> >         p
> >     } else if (alternative == "greater") {
> >         1 - p
> >     } else 2 * pmin(p, 1 - p)
> >     p
> > }
> > 
> > The test:
> > 
> > > system.time(c1 <- cor.pvals(X), gcFirst=TRUE)
> > [1] 13.19  0.01 13.58    NA    NA
> > > system.time(c2 <- cor.pvalues(X), gcFirst=TRUE)
> > [1] 6.22 0.00 6.42   NA   NA
> > > system.time(c3 <- cor.pval(X), gcFirst=TRUE)
> > [1] 0.07 0.00 0.07   NA   NA
> > 
> > Cheers,
> > Andy
> > 
> > > From: John Fox
> > > 
> > > Dear Mark,
> > > 
> > > I think that the reflex of trying to avoid loops in R is often 
> > > mistaken, and so I decided to try to time the two 
> approaches (on a 
> > > 3GHz Windows XP system).
> > > 
> > > I discovered, first, that there is a bug in your function -- you 
> > > appear to have indexed rows instead of columns; fixing that:
> > > 
> > > cor.pvals <- function(mat)
> > > {
> > >   cols <- expand.grid(1:ncol(mat), 1:ncol(mat))
> > >   matrix(apply(cols, 1,
> > >                function(x) cor.test(mat[, x[1]], mat[, 
> > > x[2]])$p.value),
> > >          ncol = ncol(mat))
> > > }
> > > 
> > > 
> > > My function is cor.pvalues and yours cor.pvals. This is 
> for a data 
> > > matrix with 1000 observations on 100 variables:
> > > 
> > > > R <- diag(100)
> > > > R[upper.tri(R)] <- R[lower.tri(R)] <- .5
> > > > library(mvtnorm)
> > > > X <- rmvnorm(1000, sigma=R)
> > > > dim(X)
> > > [1] 1000  100
> > > > 
> > > > system.time(cor.pvalues(X))
> > > [1] 5.53 0.00 5.53   NA   NA
> > > > 
> > > > system.time(cor.pvals(X))
> > > [1] 12.66  0.00 12.66    NA    NA
> > > > 
> > > 
> > > I frankly didn't expect the advantage of my approach to be 
> > this large, 
> > > but there it is.
> > > 
> > > Regards,
> > >  John
> > > 
> > > --------------------------------
> > > John Fox
> > > Department of Sociology
> > > McMaster University
> > > Hamilton, Ontario
> > > Canada L8S 4M4
> > > 905-525-9140x23604
> > > http://socserv.mcmaster.ca/jfox
> > > --------------------------------
> > > 
> > > > -----Original Message-----
> > > > From: Marc Schwartz [mailto:MSchwartz at MedAnalytics.com]
> > > > Sent: Friday, April 15, 2005 7:08 PM
> > > > To: John Fox
> > > > Cc: 'Dren Scott'; R-Help
> > > > Subject: RE: [R] Pearson corelation and p-value for matrix
> > > > 
> > > > Here is what might be a slightly more efficient way to 
> > get to John's
> > > > question:
> > > > 
> > > > cor.pvals <- function(mat)
> > > > {
> > > >   rows <- expand.grid(1:nrow(mat), 1:nrow(mat))
> > > >   matrix(apply(rows, 1,
> > > >                function(x) cor.test(mat[x[1], ], mat[x[2], 
> > > > ])$p.value),
> > > >          ncol = nrow(mat))
> > > > }
> > > > 
> > > > HTH,
> > > > 
> > > > Marc Schwartz
> > > > 
> > > > On Fri, 2005-04-15 at 18:26 -0400, John Fox wrote:
> > > > > Dear Dren,
> > > > > 
> > > > > How about the following?
> > > > > 
> > > > >  cor.pvalues <- function(X){
> > > > >     nc <- ncol(X)
> > > > >     res <- matrix(0, nc, nc)
> > > > >     for (i in 2:nc){
> > > > >         for (j in 1:(i - 1)){
> > > > >             res[i, j] <- res[j, i] <-
cor.test(X[,i],
> > > X[,j])$p.value
> > > > >             }
> > > > >         }
> > > > >     res
> > > > >     }
> > > > > 
> > > > > What one then does with all of those non-independent
test
> > > > is another
> > > > > question, I guess.
> > > > > 
> > > > > I hope this helps,
> > > > >  John
> > > > 
> > > > > > -----Original Message-----
> > > > > > From: r-help-bounces at stat.math.ethz.ch 
> > > > > > [mailto:r-help-bounces at stat.math.ethz.ch] On
Behalf Of
> > > Dren Scott
> > > > > > Sent: Friday, April 15, 2005 4:33 PM
> > > > > > To: r-help at stat.math.ethz.ch
> > > > > > Subject: [R] Pearson corelation and p-value for
matrix
> > > > > > 
> > > > > > Hi,
> > > > > >  
> > > > > > I was trying to evaluate the pearson correlation
and
> > > the p-values
> > > > > > for an nxm matrix, where each row represents a
vector.
> > > > One way to do
> > > > > > it would be to iterate through each row, and find
its
> > > correlation
> > > > > > value( and the p-value) with respect to the other
rows.
> > > Is there
> > > > > > some function by which I can use the matrix as
input?
> > > > Ideally, the
> > > > > > output would be an nxn matrix, containing the
p-values
> > > > between the
> > > > > > respective vectors.
> > > > > >  
> > > > > > I have tried cor.test for the iterations, but 
> couldn't find a 
> > > > > > function that would take the matrix as input.
> > > > > >  
> > > > > > Thanks for the help.
> > > > > >  
> > > > > > Dren
> > > > 
> > > >
> > > 
> > > ______________________________________________
> > > R-help at stat.math.ethz.ch mailing list
> > > https://stat.ethz.ch/mailman/listinfo/r-help
> > > PLEASE do read the posting guide! 
> > > http://www.R-project.org/posting-guide.html
> > > 
> > > 
> > > 
> > 
> > 
> > 
> > --------------------------------------------------------------
> > ----------------
> > Notice:  This e-mail message, together with any attachments, 
> > contains information of Merck & Co., Inc. (One Merck Drive, 
> > Whitehouse Station, New Jersey, USA 08889), and/or its 
> > affiliates (which may be known outside the United States as 
> > Merck Frosst, Merck Sharp & Dohme or MSD and in Japan, as 
> > Banyu) that may be confidential, proprietary copyrighted 
> > and/or legally privileged. It is intended solely for the use 
> > of the individual or entity named on this message.  If you 
> > are not the intended recipient, and have received this 
> > message in error, please notify us immediately by reply 
> > e-mail and then delete it from your system.
> > --------------------------------------------------------------
> > ----------------
> 
> 
>

Liaw, Andy

2005-Apr-18 19:40 UTC

head link

[R] Pearson corelation and p-value for matrix

I believe this will do:

cor.pval2 <- function(x,  alternative="two-sided") {
    corMat <- cor(x, use=if (any(is.na(x))) "pairwise.complete.obs"
else
"all")
    df <- crossprod(!is.na(x)) - 2
    STATISTIC <- sqrt(df) * corMat / sqrt(1 - corMat^2)
    p <- pt(STATISTIC, df)
    p <- if (alternative == "less") {
        p
    } else if (alternative == "greater") {
        1 - p
    } else 2 * pmin(p, 1 - p)
    p
}

Some test:
> set.seed(1)
> x <- matrix(runif(2e3 * 1e2), 2e3)
> system.time(res1 <- cor.pval(t(x)), gcFirst=TRUE)
[1] 17.28  0.77 18.16    NA    NA> system.time(res2 <- cor.pval2(t(x)), gcFirst=TRUE)
[1] 19.51  1.05 20.70    NA    NA> max(abs(res1 - res2))
[1] 0> x[c(1, 3, 6), c(2, 4)] <- NA
> x[30, 61] <- NA
> system.time(res2 <- cor.pval2(t(x)), gcFirst=TRUE)[1] 24.48  0.71 25.28    NA    NA

This is a bit slower because of the extra computation for "df".  One
can try
to save some computation by only computing with the lower (or upper)
triangular part.

Cheers,
Andy
> From: John Fox
> 
> Dear Dren,
> 
> Since cor(), on which Andy's solution is based, can compute 
> pairwise-present
> correlations, you could adapt his function -- you'll have to 
> adjust the df
> for each pair. Alternatively, you could probably save some 
> time (programming
> time + computer time) by just using my solution:
> 
> > R <- diag(100)
> > R[upper.tri(R)] <- R[lower.tri(R)] <- .5
> > library(mvtnorm)
> > X <- rmvnorm(6000, sigma=R)
> > system.time(for (i in 1:50) cor.pvalues(X), gc=TRUE)
> [1] 518.19   1.11 520.23     NA     NA
> 
> I know that time is money, but nine minutes (on my machine) 
> probably won't
> bankrupt anyone.
> 
> Regards,
>  John
> 
> --------------------------------
> John Fox
> Department of Sociology
> McMaster University
> Hamilton, Ontario
> Canada L8S 4M4
> 905-525-9140x23604
> http://socserv.mcmaster.ca/jfox 
> -------------------------------- 
> 
> > -----Original Message-----
> > From: r-help-bounces at stat.math.ethz.ch 
> > [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Dren Scott
> > Sent: Monday, April 18, 2005 12:41 PM
> > To: 'R-Help'
> > Subject: RE: [R] Pearson corelation and p-value for matrix
> > 
> > Hi all,
> >  
> > Thanks Andy, Mark and John for all the help. I really 
> > appreciate it. I'm new to both R and statistics, so please 
> > excuse any gaffes on my part. 
> >  
> > Essentially what I'm trying to do, is to evaluate for each 
> > row, how many other rows would have a p-value < 0.05. So, 
> > after I get my N x N p-value matrix, I'll just filter out 
> > values that are > 0.05.
> >  
> > Each of my datasets (6000 rows x 100 columns) would consist 
> > of some NA's. The iterative procedure (cor.pvalues) proposed 
> > by John would yield the values, but it would take an 
> > inordinately long time (I have 50 of these datasets to 
> > process). The solution proposed by Andy, although fast, would 
> > not be able to incorporate the NA's.
> >  
> > Is there any workaround for the NA's? Or possibly do you 
> > think I could try something else?
> >  
> > Thanks very much. 
> >  
> > Dren
> > 
> > 
> > "Liaw, Andy" <andy_liaw at merck.com> wrote:
> > > From: John Fox
> > > 
> > > Dear Andy,
> > > 
> > > That's clearly much better -- and illustrates an 
> effective strategy 
> > > for vectorizing (or "matricizing") a computation. I
think
> I'll add 
> > > this to my list of programming examples. It might be a little 
> > > dangerous to pass ...
> > > through to cor(), since someone could specify 
> type="spearman", for 
> > > example.
> > 
> > Ah, yes, the "..." isn't likely to help here! Also, it
will
> > only work correctly if there are no NA's, for example (or 
> > else the degree of freedom would be wrong). 
> > 
> > Best,
> > Andy
> > 
> > > Thanks,
> > > John
> > > 
> > > --------------------------------
> > > John Fox
> > > Department of Sociology
> > > McMaster University
> > > Hamilton, Ontario
> > > Canada L8S 4M4
> > > 905-525-9140x23604
> > > http://socserv.mcmaster.ca/jfox
> > > --------------------------------
> > > 
> > > > -----Original Message-----
> > > > From: Liaw, Andy [mailto:andy_liaw at merck.com]
> > > > Sent: Friday, April 15, 2005 9:51 PM
> > > > To: 'John Fox'; MSchwartz at medanalytics.com
> > > > Cc: 'R-Help'; 'Dren Scott'
> > > > Subject: RE: [R] Pearson corelation and p-value for matrix
> > > > 
> > > > We can be a bit sneaky and `borrow' code from
cor.test.default:
> > > > 
> > > > cor.pval <- function(x,
alternative="two-sided", ...) {
> corMat <- 
> > > > cor(x, ...) n <- nrow(x) df <- n - 2 STATISTIC <- 
> > sqrt(df) * corMat 
> > > > / sqrt(1 - corMat^2) p <- pt(STATISTIC, df) p <- if 
> > (alternative == 
> > > > "less") { p } else if (alternative ==
"greater") {
> > > > 1 - p
> > > > } else 2 * pmin(p, 1 - p)
> > > > p
> > > > }
> > > > 
> > > > The test:
> > > > 
> > > > > system.time(c1 <- cor.pvals(X), gcFirst=TRUE)
> > > > [1] 13.19 0.01 13.58 NA NA
> > > > > system.time(c2 <- cor.pvalues(X), gcFirst=TRUE)
> > > > [1] 6.22 0.00 6.42 NA NA
> > > > > system.time(c3 <- cor.pval(X), gcFirst=TRUE)
> > > > [1] 0.07 0.00 0.07 NA NA
> > > > 
> > > > Cheers,
> > > > Andy
> > > > 
> > > > > From: John Fox
> > > > > 
> > > > > Dear Mark,
> > > > > 
> > > > > I think that the reflex of trying to avoid loops in R 
> is often 
> > > > > mistaken, and so I decided to try to time the two
> > > approaches (on a
> > > > > 3GHz Windows XP system).
> > > > > 
> > > > > I discovered, first, that there is a bug in your 
> > function -- you 
> > > > > appear to have indexed rows instead of columns; fixing
that:
> > > > > 
> > > > > cor.pvals <- function(mat)
> > > > > {
> > > > > cols <- expand.grid(1:ncol(mat), 1:ncol(mat)) 
> > matrix(apply(cols, 
> > > > > 1,
> > > > > function(x) cor.test(mat[, x[1]], mat[, 
> x[2]])$p.value), ncol = 
> > > > > ncol(mat)) }
> > > > > 
> > > > > 
> > > > > My function is cor.pvalues and yours cor.pvals. This is
> > > for a data
> > > > > matrix with 1000 observations on 100 variables:
> > > > > 
> > > > > > R <- diag(100)
> > > > > > R[upper.tri(R)] <- R[lower.tri(R)] <- .5
> > > > > > library(mvtnorm)
> > > > > > X <- rmvnorm(1000, sigma=R)
> > > > > > dim(X)
> > > > > [1] 1000 100
> > > > > > 
> > > > > > system.time(cor.pvalues(X))
> > > > > [1] 5.53 0.00 5.53 NA NA
> > > > > > 
> > > > > > system.time(cor.pvals(X))
> > > > > [1] 12.66 0.00 12.66 NA NA
> > > > > > 
> > > > > 
> > > > > I frankly didn't expect the advantage of my
approach to be
> > > > this large,
> > > > > but there it is.
> > > > > 
> > > > > Regards,
> > > > > John
> > > > > 
> > > > > --------------------------------
> > > > > John Fox
> > > > > Department of Sociology
> > > > > McMaster University
> > > > > Hamilton, Ontario
> > > > > Canada L8S 4M4
> > > > > 905-525-9140x23604
> > > > > http://socserv.mcmaster.ca/jfox
> > > > > --------------------------------
> > > > > 
> > > > > > -----Original Message-----
> > > > > > From: Marc Schwartz [mailto:MSchwartz at
MedAnalytics.com]
> > > > > > Sent: Friday, April 15, 2005 7:08 PM
> > > > > > To: John Fox
> > > > > > Cc: 'Dren Scott'; R-Help
> > > > > > Subject: RE: [R] Pearson corelation and p-value
for matrix
> > > > > > 
> > > > > > Here is what might be a slightly more efficient
way to
> > > > get to John's
> > > > > > question:
> > > > > > 
> > > > > > cor.pvals <- function(mat)
> > > > > > {
> > > > > > rows <- expand.grid(1:nrow(mat), 1:nrow(mat)) 
> > matrix(apply(rows, 
> > > > > > 1,
> > > > > > function(x) cor.test(mat[x[1], ], mat[x[2], 
> > ])$p.value), ncol = 
> > > > > > nrow(mat)) }
> > > > > > 
> > > > > > HTH,
> > > > > > 
> > > > > > Marc Schwartz
> > > > > > 
> > > > > > On Fri, 2005-04-15 at 18:26 -0400, John Fox wrote:
> > > > > > > Dear Dren,
> > > > > > > 
> > > > > > > How about the following?
> > > > > > > 
> > > > > > > cor.pvalues <- function(X){
> > > > > > > nc <- ncol(X)
> > > > > > > res <- matrix(0, nc, nc)
> > > > > > > for (i in 2:nc){
> > > > > > > for (j in 1:(i - 1)){
> > > > > > > res[i, j] <- res[j, i] <-
cor.test(X[,i],
> > > > > X[,j])$p.value
> > > > > > > }
> > > > > > > }
> > > > > > > res
> > > > > > > }
> > > > > > > 
> > > > > > > What one then does with all of those
non-independent test
> > > > > > is another
> > > > > > > question, I guess.
> > > > > > > 
> > > > > > > I hope this helps,
> > > > > > > John
> > > > > > 
> > > > > > > > -----Original Message-----
> > > > > > > > From: r-help-bounces at
stat.math.ethz.ch
> > > > > > > > [mailto:r-help-bounces at
stat.math.ethz.ch] On Behalf Of
> > > > > Dren Scott
> > > > > > > > Sent: Friday, April 15, 2005 4:33 PM
> > > > > > > > To: r-help at stat.math.ethz.ch
> > > > > > > > Subject: [R] Pearson corelation and
p-value for matrix
> > > > > > > > 
> > > > > > > > Hi,
> > > > > > > > 
> > > > > > > > I was trying to evaluate the pearson
correlation and
> > > > > the p-values
> > > > > > > > for an nxm matrix, where each row
represents a vector.
> > > > > > One way to do
> > > > > > > > it would be to iterate through each row,
and find its
> > > > > correlation
> > > > > > > > value( and the p-value) with respect to
the other rows.
> > > > > Is there
> > > > > > > > some function by which I can use the
matrix as input?
> > > > > > Ideally, the
> > > > > > > > output would be an nxn matrix,
containing the p-values
> > > > > > between the
> > > > > > > > respective vectors.
> > > > > > > > 
> > > > > > > > I have tried cor.test for the
iterations, but
> > > couldn't find a
> > > > > > > > function that would take the matrix as
input.
> > > > > > > > 
> > > > > > > > Thanks for the help.
> > > > > > > > 
> > > > > > > > Dren
> > > > > > 
> > > > > >
> > > > > 
> > > > > ______________________________________________
> > > > > R-help at stat.math.ethz.ch mailing list 
> > > > > https://stat.ethz.ch/mailman/listinfo/r-help
> > > > > PLEASE do read the posting guide! 
> > > > > http://www.R-project.org/posting-guide.html
> > > > > 
> > > > > 
> > > > > 
> > > > 
> > > > 
> > > > 
> > > >
--------------------------------------------------------------
> > > > ----------------
> > > > Notice: This e-mail message, together with any 
> > attachments, contains 
> > > > information of Merck & Co., Inc. (One Merck Drive,
Whitehouse
> > > > Station, New Jersey, USA 08889), and/or its affiliates 
> > (which may be 
> > > > known outside the United States as Merck Frosst, Merck 
> > Sharp & Dohme 
> > > > or MSD and in Japan, as
> > > > Banyu) that may be confidential, proprietary copyrighted
and/or
> > > > legally privileged. It is intended solely for the use of the
> > > > individual or entity named on this message. If you are not
the
> > > > intended recipient, and have received this message in 
> > error, please 
> > > > notify us immediately by reply e-mail and then delete it 
> > from your 
> > > > system.
> > > >
--------------------------------------------------------------
> > > > ----------------
> > > 
> > > 
> > > 
> > 
> > 
> > 
> > --------------------------------------------------------------
> > ----------------
> > 
> > --------------------------------------------------------------
> > ----------------
> > 
> > 		
> > ---------------------------------
> > 
> > 
> > 	[[alternative HTML version deleted]]
> > 
> > ______________________________________________
> > R-help at stat.math.ethz.ch mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide! 
> > http://www.R-project.org/posting-guide.html
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! 
> http://www.R-project.org/posting-guide.html
> 
> 
>

Liaw, Andy

2005-Apr-19 00:33 UTC

head link

[R] Pearson corelation and p-value for matrix

> From: John Fox 
> 
> Dear Andy,
> 
> Very nice! (My point was that if this is a one-time thing, for Dren to
> puzzle over it is probably more time-consuming than simply doing it
> inefficiently.)
There's a memory/speed tradeoff.  The loop avoids creating multiple
6000x6000 matrices that my version would require, and one of those would
take up 274MB!  I could not run my function on a 6000x6000 case on my laptop
(nor rcorr).  The example I showed was 2000x2000.  Also, as Frank pointed
out, it's not flexible.  As John said, it's good for one-time use.  For
more
general consumption, one would want to write more flexible and robust code,
in which case one might very well re-invent rcorr...

Cheers,
Andy
 > Regards,
>  John
> 
> --------------------------------
> John Fox
> Department of Sociology
> McMaster University
> Hamilton, Ontario
> Canada L8S 4M4
> 905-525-9140x23604
> http://socserv.mcmaster.ca/jfox 
> -------------------------------- 
> 
> > -----Original Message-----
> > From: Liaw, Andy [mailto:andy_liaw at merck.com] 
> > Sent: Monday, April 18, 2005 2:40 PM
> > To: 'John Fox'; 'Dren Scott'
> > Cc: 'R-Help'
> > Subject: RE: [R] Pearson corelation and p-value for matrix
> > 
> > I believe this will do:
> > 
> > cor.pval2 <- function(x,  alternative="two-sided") {
> >     corMat <- cor(x, use=if (any(is.na(x))) 
> > "pairwise.complete.obs" else
> > "all")
> >     df <- crossprod(!is.na(x)) - 2
> >     STATISTIC <- sqrt(df) * corMat / sqrt(1 - corMat^2)
> >     p <- pt(STATISTIC, df)
> >     p <- if (alternative == "less") {
> >         p
> >     } else if (alternative == "greater") {
> >         1 - p
> >     } else 2 * pmin(p, 1 - p)
> >     p
> > }
> > 
> > Some test:
> > 
> > > set.seed(1)
> > > x <- matrix(runif(2e3 * 1e2), 2e3)
> > > system.time(res1 <- cor.pval(t(x)), gcFirst=TRUE)
> > [1] 17.28  0.77 18.16    NA    NA
> > > system.time(res2 <- cor.pval2(t(x)), gcFirst=TRUE)
> > [1] 19.51  1.05 20.70    NA    NA
> > > max(abs(res1 - res2))
> > [1] 0
> > > x[c(1, 3, 6), c(2, 4)] <- NA
> > > x[30, 61] <- NA
> > > system.time(res2 <- cor.pval2(t(x)), gcFirst=TRUE)
> > [1] 24.48  0.71 25.28    NA    NA
> > 
> > This is a bit slower because of the extra computation for 
> > "df".  One can try to save some computation by only
computing
> > with the lower (or upper) triangular part.
> > 
> > Cheers,
> > Andy
> > 
> > > From: John Fox
> > > 
> > > Dear Dren,
> > > 
> > > Since cor(), on which Andy's solution is based, can compute 
> > > pairwise-present correlations, you could adapt his function 
> > -- you'll 
> > > have to adjust the df for each pair. Alternatively, you 
> > could probably 
> > > save some time (programming time + computer time) by just 
> using my 
> > > solution:
> > > 
> > > > R <- diag(100)
> > > > R[upper.tri(R)] <- R[lower.tri(R)] <- .5
> > > > library(mvtnorm)
> > > > X <- rmvnorm(6000, sigma=R)
> > > > system.time(for (i in 1:50) cor.pvalues(X), gc=TRUE)
> > > [1] 518.19   1.11 520.23     NA     NA
> > > 
> > > I know that time is money, but nine minutes (on my machine) 
> > probably 
> > > won't bankrupt anyone.
> > > 
> > > Regards,
> > >  John
> > > 
> > > --------------------------------
> > > John Fox
> > > Department of Sociology
> > > McMaster University
> > > Hamilton, Ontario
> > > Canada L8S 4M4
> > > 905-525-9140x23604
> > > http://socserv.mcmaster.ca/jfox
> > > --------------------------------
> > > 
> > > > -----Original Message-----
> > > > From: r-help-bounces at stat.math.ethz.ch 
> > > > [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of 
> Dren Scott
> > > > Sent: Monday, April 18, 2005 12:41 PM
> > > > To: 'R-Help'
> > > > Subject: RE: [R] Pearson corelation and p-value for matrix
> > > > 
> > > > Hi all,
> > > >  
> > > > Thanks Andy, Mark and John for all the help. I really 
> > appreciate it. 
> > > > I'm new to both R and statistics, so please excuse any 
> > gaffes on my 
> > > > part.
> > > >  
> > > > Essentially what I'm trying to do, is to evaluate for 
> > each row, how 
> > > > many other rows would have a p-value < 0.05. So, after I 
> > get my N x 
> > > > N p-value matrix, I'll just filter out values that are
> 0.05.
> > > >  
> > > > Each of my datasets (6000 rows x 100 columns) would 
> > consist of some 
> > > > NA's. The iterative procedure (cor.pvalues) proposed by 
> > John would 
> > > > yield the values, but it would take an inordinately 
> long time (I 
> > > > have 50 of these datasets to process). The solution proposed
by
> > > > Andy, although fast, would not be able to incorporate the
NA's.
> > > >  
> > > > Is there any workaround for the NA's? Or possibly do 
> you think I 
> > > > could try something else?
> > > >  
> > > > Thanks very much. 
> > > >  
> > > > Dren
> > > > 
> > > > 
> > > > "Liaw, Andy" <andy_liaw at merck.com> wrote:
> > > > > From: John Fox
> > > > > 
> > > > > Dear Andy,
> > > > > 
> > > > > That's clearly much better -- and illustrates an
> > > effective strategy
> > > > > for vectorizing (or "matricizing") a
computation. I think
> > > I'll add
> > > > > this to my list of programming examples. It might be a
little
> > > > > dangerous to pass ...
> > > > > through to cor(), since someone could specify
> > > type="spearman", for
> > > > > example.
> > > > 
> > > > Ah, yes, the "..." isn't likely to help here!
Also, it
> will only 
> > > > work correctly if there are no NA's, for example (or
else
> > the degree 
> > > > of freedom would be wrong).
> > > > 
> > > > Best,
> > > > Andy
> > > > 
> > > > > Thanks,
> > > > > John
> > > > > 
> > > > > --------------------------------
> > > > > John Fox
> > > > > Department of Sociology
> > > > > McMaster University
> > > > > Hamilton, Ontario
> > > > > Canada L8S 4M4
> > > > > 905-525-9140x23604
> > > > > http://socserv.mcmaster.ca/jfox
> > > > > --------------------------------
> > > > > 
> > > > > > -----Original Message-----
> > > > > > From: Liaw, Andy [mailto:andy_liaw at merck.com]
> > > > > > Sent: Friday, April 15, 2005 9:51 PM
> > > > > > To: 'John Fox'; MSchwartz at
medanalytics.com
> > > > > > Cc: 'R-Help'; 'Dren Scott'
> > > > > > Subject: RE: [R] Pearson corelation and p-value
for matrix
> > > > > > 
> > > > > > We can be a bit sneaky and `borrow' code from 
> > cor.test.default:
> > > > > > 
> > > > > > cor.pval <- function(x,
alternative="two-sided", ...) {
> > > corMat <-
> > > > > > cor(x, ...) n <- nrow(x) df <- n - 2
STATISTIC <-
> > > > sqrt(df) * corMat
> > > > > > / sqrt(1 - corMat^2) p <- pt(STATISTIC, df) p
<- if
> > > > (alternative => > > > > >
"less") { p } else if (alternative == "greater") {
> > > > > > 1 - p
> > > > > > } else 2 * pmin(p, 1 - p)
> > > > > > p
> > > > > > }
> > > > > > 
> > > > > > The test:
> > > > > > 
> > > > > > > system.time(c1 <- cor.pvals(X),
gcFirst=TRUE)
> > > > > > [1] 13.19 0.01 13.58 NA NA
> > > > > > > system.time(c2 <- cor.pvalues(X),
gcFirst=TRUE)
> > > > > > [1] 6.22 0.00 6.42 NA NA
> > > > > > > system.time(c3 <- cor.pval(X),
gcFirst=TRUE)
> > > > > > [1] 0.07 0.00 0.07 NA NA
> > > > > > 
> > > > > > Cheers,
> > > > > > Andy
> > > > > > 
> > > > > > > From: John Fox
> > > > > > > 
> > > > > > > Dear Mark,
> > > > > > > 
> > > > > > > I think that the reflex of trying to avoid
loops in R
> > > is often
> > > > > > > mistaken, and so I decided to try to time the
two
> > > > > approaches (on a
> > > > > > > 3GHz Windows XP system).
> > > > > > > 
> > > > > > > I discovered, first, that there is a bug in
your
> > > > function -- you
> > > > > > > appear to have indexed rows instead of
columns;
> fixing that:
> > > > > > > 
> > > > > > > cor.pvals <- function(mat)
> > > > > > > {
> > > > > > > cols <- expand.grid(1:ncol(mat),
1:ncol(mat))
> > > > matrix(apply(cols,
> > > > > > > 1,
> > > > > > > function(x) cor.test(mat[, x[1]], mat[,
> > > x[2]])$p.value), ncol > > > > > > >
ncol(mat)) }
> > > > > > > 
> > > > > > > 
> > > > > > > My function is cor.pvalues and yours
cor.pvals. This is
> > > > > for a data
> > > > > > > matrix with 1000 observations on 100
variables:
> > > > > > > 
> > > > > > > > R <- diag(100)
> > > > > > > > R[upper.tri(R)] <- R[lower.tri(R)]
<- .5
> > > > > > > > library(mvtnorm)
> > > > > > > > X <- rmvnorm(1000, sigma=R)
> > > > > > > > dim(X)
> > > > > > > [1] 1000 100
> > > > > > > > 
> > > > > > > > system.time(cor.pvalues(X))
> > > > > > > [1] 5.53 0.00 5.53 NA NA
> > > > > > > > 
> > > > > > > > system.time(cor.pvals(X))
> > > > > > > [1] 12.66 0.00 12.66 NA NA
> > > > > > > > 
> > > > > > > 
> > > > > > > I frankly didn't expect the advantage of
my approach to be
> > > > > > this large,
> > > > > > > but there it is.
> > > > > > > 
> > > > > > > Regards,
> > > > > > > John
> > > > > > > 
> > > > > > > -------------------------------- John Fox
Department of
> > > > > > > Sociology McMaster University Hamilton,
Ontario
> > Canada L8S 4M4
> > > > > > > 905-525-9140x23604
> > > > > > > http://socserv.mcmaster.ca/jfox
> > > > > > > --------------------------------
> > > > > > > 
> > > > > > > > -----Original Message-----
> > > > > > > > From: Marc Schwartz [mailto:MSchwartz at
MedAnalytics.com]
> > > > > > > > Sent: Friday, April 15, 2005 7:08 PM
> > > > > > > > To: John Fox
> > > > > > > > Cc: 'Dren Scott'; R-Help
> > > > > > > > Subject: RE: [R] Pearson corelation and
p-value
> for matrix
> > > > > > > > 
> > > > > > > > Here is what might be a slightly more
efficient way to
> > > > > > get to John's
> > > > > > > > question:
> > > > > > > > 
> > > > > > > > cor.pvals <- function(mat)
> > > > > > > > {
> > > > > > > > rows <- expand.grid(1:nrow(mat),
1:nrow(mat))
> > > > matrix(apply(rows, 
> > > > > > > > 1,
> > > > > > > > function(x) cor.test(mat[x[1], ],
mat[x[2],
> > > > ])$p.value), ncol = 
> > > > > > > > nrow(mat)) }
> > > > > > > > 
> > > > > > > > HTH,
> > > > > > > > 
> > > > > > > > Marc Schwartz
> > > > > > > > 
> > > > > > > > On Fri, 2005-04-15 at 18:26 -0400, John
Fox wrote:
> > > > > > > > > Dear Dren,
> > > > > > > > > 
> > > > > > > > > How about the following?
> > > > > > > > > 
> > > > > > > > > cor.pvalues <- function(X){
> > > > > > > > > nc <- ncol(X)
> > > > > > > > > res <- matrix(0, nc, nc)
> > > > > > > > > for (i in 2:nc){
> > > > > > > > > for (j in 1:(i - 1)){
> > > > > > > > > res[i, j] <- res[j, i] <-
cor.test(X[,i],
> > > > > > > X[,j])$p.value
> > > > > > > > > }
> > > > > > > > > }
> > > > > > > > > res
> > > > > > > > > }
> > > > > > > > > 
> > > > > > > > > What one then does with all of
those
> > non-independent test
> > > > > > > > is another
> > > > > > > > > question, I guess.
> > > > > > > > > 
> > > > > > > > > I hope this helps,
> > > > > > > > > John
> > > > > > > > 
> > > > > > > > > > -----Original Message-----
> > > > > > > > > > From: r-help-bounces at
stat.math.ethz.ch
> > > > > > > > > > [mailto:r-help-bounces at
stat.math.ethz.ch]
> On Behalf Of
> > > > > > > Dren Scott
> > > > > > > > > > Sent: Friday, April 15, 2005
4:33 PM
> > > > > > > > > > To: r-help at
stat.math.ethz.ch
> > > > > > > > > > Subject: [R] Pearson
corelation and p-value
> for matrix
> > > > > > > > > > 
> > > > > > > > > > Hi,
> > > > > > > > > > 
> > > > > > > > > > I was trying to evaluate the
pearson correlation and
> > > > > > > the p-values
> > > > > > > > > > for an nxm matrix, where each
row represents
> > a vector. 
> > > > > > > > One way to do
> > > > > > > > > > it would be to iterate through
each row,
> and find its
> > > > > > > correlation
> > > > > > > > > > value( and the p-value) with
respect to the
> > other rows. 
> > > > > > > Is there
> > > > > > > > > > some function by which I can
use the matrix
> as input? 
> > > > > > > > Ideally, the
> > > > > > > > > > output would be an nxn matrix,
containing
> the p-values
> > > > > > > > between the
> > > > > > > > > > respective vectors.
> > > > > > > > > > 
> > > > > > > > > > I have tried cor.test for the
iterations, but
> > > > > couldn't find a
> > > > > > > > > > function that would take the
matrix as input.
> > > > > > > > > > 
> > > > > > > > > > Thanks for the help.
> > > > > > > > > > 
> > > > > > > > > > Dren
> > > > > > > > 
> > > > > > > >
> > > > > > > 
> > > > > > >
______________________________________________
> > > > > > > R-help at stat.math.ethz.ch mailing list 
> > > > > > > https://stat.ethz.ch/mailman/listinfo/r-help
> > > > > > > PLEASE do read the posting guide! 
> > > > > > > http://www.R-project.org/posting-guide.html
> > > > > > > 
> > > > > > > 
> > > > > > > 
> > > > > > 
> > > > > > 
> > > > > > 
> > > > > > 
> --------------------------------------------------------------
> > > > > > ----------------
> > > > > > Notice: This e-mail message, together with any 
> > > > attachments, contains 
> > > > > > information of Merck & Co., Inc. (One Merck
Drive,
> Whitehouse 
> > > > > > Station, New Jersey, USA 08889), and/or its
affiliates
> > > > (which may be 
> > > > > > known outside the United States as Merck Frosst,
Merck
> > > > Sharp & Dohme 
> > > > > > or MSD and in Japan, as
> > > > > > Banyu) that may be confidential, proprietary 
> > copyrighted and/or 
> > > > > > legally privileged. It is intended solely for the 
> use of the 
> > > > > > individual or entity named on this message. If you
> > are not the 
> > > > > > intended recipient, and have received this message
in
> > > > error, please 
> > > > > > notify us immediately by reply e-mail and then
delete it
> > > > from your 
> > > > > > system.
> > > > > > 
> --------------------------------------------------------------
> > > > > > ----------------
> > > > > 
> > > > > 
> > > > > 
> > > > 
> > > > 
> > > > 
> > > >
--------------------------------------------------------------
> > > > ----------------
> > > > 
> > > >
--------------------------------------------------------------
> > > > ----------------
> > > > 
> > > > 		
> > > > ---------------------------------
> > > > 
> > > > 
> > > > 	[[alternative HTML version deleted]]
> > > > 
> > > > ______________________________________________
> > > > R-help at stat.math.ethz.ch mailing list
> > > > https://stat.ethz.ch/mailman/listinfo/r-help
> > > > PLEASE do read the posting guide! 
> > > > http://www.R-project.org/posting-guide.html
> > > 
> > > ______________________________________________
> > > R-help at stat.math.ethz.ch mailing list
> > > https://stat.ethz.ch/mailman/listinfo/r-help
> > > PLEASE do read the posting guide! 
> > > http://www.R-project.org/posting-guide.html
> > > 
> > > 
> > > 
> > 
> > 
> > 
> > --------------------------------------------------------------
> > ----------------
> > Notice:  This e-mail message, together with any attachments, 
> > contains information of Merck & Co., Inc. (One Merck Drive, 
> > Whitehouse Station, New Jersey, USA 08889), and/or its 
> > affiliates (which may be known outside the United States as 
> > Merck Frosst, Merck Sharp & Dohme or MSD and in Japan, as 
> > Banyu) that may be confidential, proprietary copyrighted 
> > and/or legally privileged. It is intended solely for the use 
> > of the individual or entity named on this message.  If you 
> > are not the intended recipient, and have received this 
> > message in error, please notify us immediately by reply 
> > e-mail and then delete it from your system.
> > --------------------------------------------------------------
> > ----------------
> 
> 
> 
>

Apparently Analagous Threads

Search for more reasonably related threads

R help - Apr 2005 - Pearson corelation and p-value for matrix

[R] Pearson corelation and p-value for matrix

[R] Pearson corelation and p-value for matrix

[R] Pearson corelation and p-value for matrix

[R] Pearson corelation and p-value for matrix

[R] Pearson corelation and p-value for matrix

[R] Pearson corelation and p-value for matrix

Apparently Analagous Threads