¿Correlación ponderada de Pearson?

Question

Tengo un 2396x34 double matrix llamado y donde cada fila (2396) representa una situación separada que consta de 34 segmentos de tiempo consecutivos.

También tengo un numeric[34] llamado x que representa una única situación de 34 segmentos de tiempo consecutivos.

Actualmente estoy calculando la correlación entre cada fila de y y x así:

crs[,2] <- cor(t(y),x)

Lo que necesito ahora es reemplazar la función cor en la declaración anterior con un correlación ponderada . El vector de peso xy.wt tiene 34 elementos de longitud, por lo que se puede asignar un peso diferente a cada uno de los 34 segmentos de tiempo consecutivos.

Encontré la función Weighted Covariance Matrixcov.wt y pensé que si primero scale los datos debería funcionar como la función cor. De hecho, puede especificar que la función también devuelva una matriz de correlación. Lamentablemente, no parece que pueda usarlo de la misma manera porque no puedo suministrar mis dos variables (x y y) por separado.

¿Alguien sabe de una manera en que puedo obtener una correlación ponderada de la manera que describí sin sacrificar mucha velocidad?

Editar: Quizá alguna función matemática se podría aplicar a y antes de la función cor con el fin de obtener los mismos resultados que yo estoy buscando. ¿Tal vez si multiplico cada elemento por xy.wt/sum(xy.wt)?

Edición # 2 me encontré con otra función corr en el paquete boot.

corr(d, w = rep(1, nrow(d))/nrow(d)) 

d 
A matrix with two columns corresponding to the two variables whose correlation we wish to calculate. 

w 
A vector of weights to be applied to each pair of observations. The default is equal weights for each pair. Normalization takes place within the function so sum(w) need not equal 1. 

Esto tampoco es lo que necesito, pero está más cerca.

Edición # 3 Aquí hay un código para generar el tipo de datos que estoy trabajando con:

x<-cumsum(rnorm(34)) 
y<- t(sapply(1:2396,function(u) cumsum(rnorm(34)))) 
xy.wt<-1/(34:1) 

crs<-cor(t(y),x) #this works but I want to use xy.wt as weight 

Answer 1

Usted puede volver a la definición de la correlación.

f <- function(x, y, w = rep(1,length(x))) { 
    stopifnot(length(x) == dim(y)[2]) 
    w <- w/sum(w) 
    # Center x and y, using the weighted means 
    x <- x - sum(x*w) 
    y <- y - apply(t(y) * w, 2, sum) 
    # Compute the variance 
    vx <- sum(w * x * x) 
    vy <- rowSums(w * y * y) # Incorrect: see Heather's remark, in the other answer 
    # Compute the covariance 
    vxy <- colSums(t(y) * x * w) 
    # Compute the correlation 
    vxy/sqrt(vx * vy) 
} 
f(x,y)[1] 
cor(x,y[1,]) # Identical 
f(x, y, xy.wt) 

Answer 2

Por desgracia, la respuesta aceptada es mal cuando y es una matriz de más de una fila. El error está en la línea

vy <- rowSums(w * y * y) 

Queremos multiplicar las columnas de y por w, pero esto va a multiplicar las filas de los elementos de w, reciclados según sea necesario.Por lo tanto

> f(x, y[1, , drop = FALSE], xy.wt) 
[1] 0.103021 

es correcto, porque en este caso la multiplicación se realiza elemento a elemento, lo que equivale a la multiplicación por columnas aquí, pero

> f(x, y, xy.wt)[1] 
[1] 0.05463575 

da una respuesta incorrecta debido a la fila- multiplicación sabia.

Podemos corregir la función de la siguiente manera

f2 <- function(x, y, w = rep(1,length(x))) { 
    stopifnot(length(x) == dim(y)[2]) 
    w <- w/sum(w) 
    # Center x and y, using the weighted means 
    x <- x - sum(x * w) 
    ty <- t(y - colSums(t(y) * w)) 
    # Compute the variance 
    vx <- sum(w * x * x) 
    vy <- colSums(w * ty * ty) 
    # Compute the covariance 
    vxy <- colSums(ty * x * w) 
    # Compute the correlation 
    vxy/sqrt(vx * vy) 
} 

y comprobar los resultados con los producidos por corr del boot paquete:

> res1 <- f2(x, y, xy.wt) 
> res2 <- sapply(1:nrow(y), 
+    function(i, x, y, w) corr(cbind(x, y[i,]), w = w), 
+    x = x, y = y, w = xy.wt) 
> all.equal(res1, res2) 
[1] TRUE 

que en sí mismo proporciona otra forma de que este problema podría ser resuelto

Answer 3

Aquí es una generalización para calcular la correlación de Pearson ponderada entre dos matrices (en lugar de un vector y una matriz, como en la pregunta original):

matrix.corr <- function (a, b, w = rep(1, nrow(a))/nrow(a)) 
{ 
    # normalize weights 
    w <- w/sum(w) 

    # center matrices 
    a <- sweep(a, 2, colSums(a * w)) 
    b <- sweep(b, 2, colSums(b * w)) 

    # compute weighted correlation 
    t(w*a) %*% b/sqrt(colSums(w * a**2) %*% t(colSums(w * b**2))) 
} 

Utilizando el ejemplo anterior y la función de correlación de Heather , podemos verificar que:

> sum(matrix.corr(as.matrix(x, nrow=34),t(y),xy.wt) - f2(x,y,xy.wt)) 
[1] 1.537507e-15 

En términos de llamar a la sintaxis, esto se asemeja a la no ponderada cor:

> a <- matrix(c(1,2,3,1,3,2), nrow=3) 
> b <- matrix(c(2,3,1,1,7,3,5,2,8,1,10,12), nrow=3) 
> matrix.corr(a,b) 
    [,1]  [,2] [,3]  [,4] 
[1,] -0.5 0.3273268 0.5 0.9386522 
[2,] 0.5 0.9819805 -0.5 0.7679882 
> cor(a, b) 
    [,1]  [,2] [,3]  [,4] 
[1,] -0.5 0.3273268 0.5 0.9386522 
[2,] 0.5 0.9819805 -0.5 0.7679882