[6508] in s-news-athena
Re: missing values ?
daemon@ATHENA.MIT.EDU (Richard M. Heiberger)
Fri Feb 24 12:15:17 1995
Date: Fri, 24 Feb 95 11:44:27 -0500
From: rmh@astro.ocis.temple.edu (Richard M. Heiberger)
To: bodo@io.org, s-news@utstat.toronto.edu
Working with pairwise-missing correlation matrices is dangerous. You are
not guaranteed that they will be positive-semi-definite.
But it is an interesting exercise in S, and one that is better done in S
than in Fortran.
The following is the working outline of a function. It does not return
the cell-wise degrees of freedom, or other useful information.
corr.mv <- function(x, method="standard", trace=F)
{
if (method == "standard") return( cor(x) )
if (method == "listwise") return( cor(x[!apply(is.na(x), 1, any),]) )
if (method == "pairwise") {
cc <- matrix(0, ncol(x), ncol(x), dimnames=dimnames(x)[c(2,2)])
for (i in 1:ncol(x)) for (j in i:ncol(x)) {
ok.ij <- !apply(is.na(x[,c(i,j)]), 1, any)
if(trace) cat(i,j,ok.ij,"\n")
cc[i,j] <- cor(x[ok.ij,i], x[ok.ij,j])
cc[j,i] <- cc[i,j]
}
cc
}
}
> x <- matrix(rnorm(15),5,3)
> x[1,1] <- NA
> x[2,1:2] <- NA
> x
[,1] [,2] [,3]
[1,] NA 0.489922573 0.74959058
[2,] NA NA 0.18643805
[3,] -1.6732589 0.008685083 1.05607548
[4,] -0.6469167 -0.047232037 -0.18109055
[5,] -0.4464253 1.693031733 -0.08711394
> corr.mv(x)
Error in .Fortran("corvar",: subroutine corvar: Missing values in argument 1
Dumped
> corr.mv(x, method="pairwise")
[,1] [,2] [,3]
[1,] 0.9999999 0.6037866 -0.9756364
[2,] 0.6037866 1.0000000 -0.3949524
[3,] -0.9756364 -0.3949524 1.0000001
> corr.mv(x, method="pairwise", trace=T)
1 1 FALSE FALSE TRUE TRUE TRUE
1 2 FALSE FALSE TRUE TRUE TRUE
1 3 FALSE FALSE TRUE TRUE TRUE
2 2 TRUE FALSE TRUE TRUE TRUE
2 3 TRUE FALSE TRUE TRUE TRUE
3 3 TRUE TRUE TRUE TRUE TRUE
[,1] [,2] [,3]
[1,] 0.9999999 0.6037866 -0.9756364
[2,] 0.6037866 1.0000000 -0.3949524
[3,] -0.9756364 -0.3949524 1.0000001
> corr.mv(x, method="listwise")
[,1] [,2] [,3]
[1,] 0.9999999 0.6037866 -0.9756364
[2,] 0.6037866 0.9999999 -0.4141876
[3,] -0.9756364 -0.4141876 1.0000000
>
