Regresión lineal simple por pares rápida entre variables en un marco de datos

performance regression

Question

He visto muchas veces la regresión lineal simple por parejas o general en el Desbordamiento de pila. Aquí hay un conjunto de datos de juguete para este tipo de problema.

set.seed(0)
X <- matrix(runif(100), 100, 5, dimnames = list(1:100, LETTERS[1:5]))
b <- c(1, 0.7, 1.3, 2.9, -2)
dat <- X * b[col(X)] + matrix(rnorm(100 * 5, 0, 0.1), 100, 5)
dat <- as.data.frame(dat)
pairs(dat)

Así que básicamente queremos calcular 5 * 4 = 20 líneas de regresión:

----- A ~ B A ~ C A ~ D A ~ E
B ~ A ----- B ~ C B ~ D B ~ E
C ~ A C ~ B ----- C ~ D C ~ E
D ~ A D ~ B D ~ C ----- D ~ E
E ~ A E ~ B E ~ C E ~ D -----

Aquí está la estrategia de un hombre pobre :

poor <- function (dat) {
n <- nrow(dat)
p <- ncol(dat)
## all formulae
LHS <- rep(colnames(dat), p)
RHS <- rep(colnames(dat), each = p)
## function to fit and summarize a single model
fitmodel <- function (LHS, RHS) {
if (RHS == LHS) {
z <- data.frame("LHS" = LHS, "RHS" = RHS,
"alpha" = 0,
"beta" = 1,
"beta.se" = 0,
"beta.tv" = Inf,
"beta.pv" = 0,
"sig" = 0,
"R2" = 1,
"F.fv" = Inf,
"F.pv" = 0,
stringsAsFactors = FALSE)
} else {
result <- summary(lm(reformulate(RHS, LHS), data = dat))
z <- data.frame("LHS" = LHS, "RHS" = RHS,
"alpha" = result$coefficients[1, 1],
"beta" = result$coefficients[2, 1],
"beta.se" = result$coefficients[2, 2],
"beta.tv" = result$coefficients[2, 3],
"beta.pv" = result$coefficients[2, 4],
"sig" = result$sigma,
"R2" = result$r.squared,
"F.fv" = result$fstatistic[[1]],
"F.pv" = pf(result$fstatistic[[1]], 1, n - 2, lower.tail = FALSE),
stringsAsFactors = FALSE)
}
z
}
## loop through all models
do.call("rbind.data.frame", c(Map(fitmodel, LHS, RHS),
list(make.row.names = FALSE,
stringsAsFactors = FALSE)))
}

La lógica es clara: obtener todos los pares, construir la fórmula del modelo ( reformulate ), ajustar una regresión ( lm ), hacer un resumen summary , devolver todas las estadísticas y rbind para que sean un marco de datos.

OK, bien, pero ¿y si hay p variables? ¡Entonces tenemos que hacer p * (p - 1) regresiones!

Una mejora inmediata que se me ocurre, es adaptar un modelo lineal con múltiples LHS . Por ejemplo, la primera columna de esa matriz de fórmula se fusiona para

cbind(B, C, D, E) ~ A

Esto reduce el número de regresión de p * (p - 1) a p .

Pero definitivamente podemos hacerlo aún mejor sin utilizar lm y summary . Aquí está mi intento anterior: ¿Existe una estimación rápida de regresión simple (una línea de regresión con solo intercepción y pendiente)? . Es rápido porque usa la covarianza entre variables para la estimación, como resolver la ecuación normal . Pero la función simpleLM allí es bastante limitada:

1. necesita calcular vectores residuales para estimar el error estándar residual, que es un cuello de botella en el rendimiento;
2. no es compatible con múltiples LHS, por lo que debe llamarse p * (p - 1) veces en la configuración de regresión por pares).

¿Podemos generalizarlo para una rápida regresión de pares, escribiendo una función pairwise_simpleLM ?

Regresión lineal simple apareada general

Una variación más útil de la regresión de pares anterior es la regresión general emparejada entre un conjunto de variables de LHS y un conjunto de variables de RHS.

Ejemplo 1

Ajuste la regresión apareada entre las variables LHS A , B , C y las variables RHS D , E , es decir, ajuste 6 líneas de regresión lineal simple:

A ~ D A ~ E
B ~ D B ~ E
C ~ D C ~ E

Ejemplo 2

Ajuste una regresión lineal simple con múltiples variables de LHS a una variable de RHS en particular, por ejemplo: cbind(A, B, C, D) ~ E

Ejemplo 3

Ajuste una regresión lineal simple con una variable de LHS particular y un conjunto de variables de RHS de una en una, por ejemplo:

A ~ B A ~ C A ~ D A ~ E 

¿Podemos también tener una función rápida general_paired_simpleLM para esto?

Precaución

1. Todas las variables deben ser numéricas; los factores no están permitidos o la regresión por pares no tiene sentido.
2. La regresión ponderada no se discute, ya que el método de varianza-covarianza no está justificado en ese caso.

Answer 1

Algunos resultados estadísticos / antecedentes

(Enlace en la imagen: Función para calcular R2 (R-cuadrado) en R )

Detalles computacionales

Los cálculos involucrados aquí son básicamente el cálculo de la matriz de varianza-covarianza. Una vez que lo tenemos, los resultados para todas las regresión por pares son solo aritmética matricial de elementos.

La matriz de varianza-covarianza puede obtenerse mediante la función R cov , pero las funciones que se encuentran a continuación la calculan manualmente utilizando el crossprod . La ventaja es que, obviamente, puede beneficiarse de una biblioteca BLAS optimizada si la tiene. Tenga en cuenta que una cantidad significativa de simplificación se realiza de esta manera. La función cov R tiene use argumentos que permite el manejo de NA , pero no lo hace crossprod . ¡Estoy asumiendo que sus datos no tienen valores faltantes en absoluto! Si tiene valores faltantes, elimínelos usted mismo con na.omit(dat) .

La inicial as.matrix que convierte un marco de datos en una matriz puede ser una sobrecarga. En principio, si codifico todo en C / C ++, puedo eliminar esta coacción. Y, de hecho, muchos aritméticos matriciales de matriz de elementos pueden fusionarse en un único nido de bucle. Sin embargo, realmente me molesto en hacer esto en este momento (ya que no tengo tiempo).

Algunas personas pueden argumentar que el formato de la declaración final es inconveniente. Podría haber otro formato:

1. una lista de marcos de datos, cada uno de los cuales proporciona el resultado de la regresión para una variable LHS particular;
2. una lista de marcos de datos, cada uno de los cuales da el resultado de la regresión para una variable RHS en particular.

Esto es realmente basado en la opinión. De todos modos, siempre puede hacer un split.data.frame por la columna "LHS" o la columna "RHS" en el marco de datos que le devuelvo.

Función R pairwise_simpleLM

pairwise_simpleLM <- function (dat) {
## matrix and its dimension (n: numbeta.ser of data; p: numbeta.ser of variables)
dat <- as.matrix(dat)
n <- nrow(dat)
p <- ncol(dat)
## variable summary: mean, (unscaled) covariance and (unscaled) variance
m <- colMeans(dat)
V <- crossprod(dat) - tcrossprod(m * sqrt(n))
d <- diag(V)
## R-squared (explained variance) and its complement
R2 <- (V ^ 2) * tcrossprod(1 / d)
R2_complement <- 1 - R2
R2_complement[seq.int(from = 1, by = p + 1, length = p)] <- 0
## slope and intercept
beta <- V * rep(1 / d, each = p)
alpha <- m - beta * rep(m, each = p)
## residual sum of squares and standard error
RSS <- R2_complement * d
sig <- sqrt(RSS * (1 / (n - 2)))
## statistics for slope
beta.se <- sig * rep(1 / sqrt(d), each = p)
beta.tv <- beta / beta.se
beta.pv <- 2 * pt(abs(beta.tv), n - 2, lower.tail = FALSE)
## F-statistic and p-value
F.fv <- (n - 2) * R2 / R2_complement
F.pv <- pf(F.fv, 1, n - 2, lower.tail = FALSE)
## export
data.frame(LHS = rep(colnames(dat), times = p),
RHS = rep(colnames(dat), each = p),
alpha = c(alpha),
beta = c(beta),
beta.se = c(beta.se),
beta.tv = c(beta.tv),
beta.pv = c(beta.pv),
sig = c(sig),
R2 = c(R2),
F.fv = c(F.fv),
F.pv = c(F.pv),
stringsAsFactors = FALSE)
}

Vamos a comparar el resultado en el conjunto de datos de juguete en la pregunta.

oo <- poor(dat)
rr <- pairwise_simpleLM(dat)
all.equal(oo, rr)
#[1] TRUE

Veamos su salida:

rr[1:3, ]
# LHS RHS alpha beta beta.se beta.tv beta.pv sig
#1 A A 0.00000000 1.0000000 0.00000000 Inf 0.000000e+00 0.0000000
#2 B A 0.05550367 0.6206434 0.04456744 13.92594 5.796437e-25 0.1252402
#3 C A 0.05809455 1.2215173 0.04790027 25.50126 4.731618e-45 0.1346059
# R2 F.fv F.pv
#1 1.0000000 Inf 0.000000e+00
#2 0.6643051 193.9317 5.796437e-25
#3 0.8690390 650.3142 4.731618e-45

Cuando tenemos el mismo LHS y RHS, la regresión no tiene sentido, por lo tanto, la intersección es 0, la pendiente es 1, etc.

¿Qué pasa con la velocidad? Todavía usando este ejemplo de juguete:

library(microbenchmark)
microbenchmark("poor_man''s" = poor(dat), "fast" = pairwise_simpleLM(dat))
#Unit: milliseconds
# expr min lq mean median uq max
# poor_man''s 127.270928 129.060515 137.813875 133.390722 139.029912 216.24995
# fast 2.732184 3.025217 3.381613 3.134832 3.313079 10.48108

La brecha se va haciendo cada vez más amplia ya que tenemos más variables. Por ejemplo, con 10 variables tenemos:

set.seed(0)
X <- matrix(runif(100), 100, 10, dimnames = list(1:100, LETTERS[1:10]))
b <- runif(10)
DAT <- X * b[col(X)] + matrix(rnorm(100 * 10, 0, 0.1), 100, 10)
DAT <- as.data.frame(DAT)
microbenchmark("poor_man''s" = poor(DAT), "fast" = pairwise_simpleLM(DAT))
#Unit: milliseconds
# expr min lq mean median uq max
# poor_man''s 548.949161 551.746631 573.009665 556.307448 564.28355 801.645501
# fast 3.365772 3.578448 3.721131 3.621229 3.77749 6.791786

Función R general_paired_simpleLM

general_paired_simpleLM <- function (dat_LHS, dat_RHS) {
## matrix and its dimension (n: numbeta.ser of data; p: numbeta.ser of variables)
dat_LHS <- as.matrix(dat_LHS)
dat_RHS <- as.matrix(dat_RHS)
if (nrow(dat_LHS) != nrow(dat_RHS)) stop("''dat_LHS'' and ''dat_RHS'' don''t have same number of rows!")
n <- nrow(dat_LHS)
pl <- ncol(dat_LHS)
pr <- ncol(dat_RHS)
## variable summary: mean, (unscaled) covariance and (unscaled) variance
ml <- colMeans(dat_LHS)
mr <- colMeans(dat_RHS)
vl <- colSums(dat_LHS ^ 2) - ml * ml * n
vr <- colSums(dat_RHS ^ 2) - mr * mr * n
##V <- crossprod(dat - rep(m, each = n)) ## cov(u, v) = E[(u - E[u])(v - E[v])]
V <- crossprod(dat_LHS, dat_RHS) - tcrossprod(ml * sqrt(n), mr * sqrt(n)) ## cov(u, v) = E[uv] - E{u]E[v]
## R-squared (explained variance) and its complement
R2 <- (V ^ 2) * tcrossprod(1 / vl, 1 / vr)
R2_complement <- 1 - R2
## slope and intercept
beta <- V * rep(1 / vr, each = pl)
alpha <- ml - beta * rep(mr, each = pl)
## residual sum of squares and standard error
RSS <- R2_complement * vl
sig <- sqrt(RSS * (1 / (n - 2)))
## statistics for slope
beta.se <- sig * rep(1 / sqrt(vr), each = pl)
beta.tv <- beta / beta.se
beta.pv <- 2 * pt(abs(beta.tv), n - 2, lower.tail = FALSE)
## F-statistic and p-value
F.fv <- (n - 2) * R2 / R2_complement
F.pv <- pf(F.fv, 1, n - 2, lower.tail = FALSE)
## export
data.frame(LHS = rep(colnames(dat_LHS), times = pr),
RHS = rep(colnames(dat_RHS), each = pl),
alpha = c(alpha),
beta = c(beta),
beta.se = c(beta.se),
beta.tv = c(beta.tv),
beta.pv = c(beta.pv),
sig = c(sig),
R2 = c(R2),
F.fv = c(F.fv),
F.pv = c(F.pv),
stringsAsFactors = FALSE)
}

Aplique esto al Ejemplo 1 en la pregunta.

general_paired_simpleLM(dat[1:3], dat[4:5])
# LHS RHS alpha beta beta.se beta.tv beta.pv sig
#1 A D -0.009212582 0.3450939 0.01171768 29.45071 1.772671e-50 0.09044509
#2 B D 0.012474593 0.2389177 0.01420516 16.81908 1.201421e-30 0.10964516
#3 C D -0.005958236 0.4565443 0.01397619 32.66585 1.749650e-54 0.10787785
#4 A E 0.008650812 -0.4798639 0.01963404 -24.44040 1.738263e-43 0.10656866
#5 B E 0.012738403 -0.3437776 0.01949488 -17.63426 3.636655e-32 0.10581331
#6 C E 0.009068106 -0.6430553 0.02183128 -29.45569 1.746439e-50 0.11849472
# R2 F.fv F.pv
#1 0.8984818 867.3441 1.772671e-50
#2 0.7427021 282.8815 1.201421e-30
#3 0.9158840 1067.0579 1.749650e-54
#4 0.8590604 597.3333 1.738263e-43
#5 0.7603718 310.9670 3.636655e-32
#6 0.8985126 867.6375 1.746439e-50

Aplique esto al Ejemplo 2 en la pregunta.

general_paired_simpleLM(dat[1:4], dat[5])
# LHS RHS alpha beta beta.se beta.tv beta.pv sig
#1 A E 0.008650812 -0.4798639 0.01963404 -24.44040 1.738263e-43 0.1065687
#2 B E 0.012738403 -0.3437776 0.01949488 -17.63426 3.636655e-32 0.1058133
#3 C E 0.009068106 -0.6430553 0.02183128 -29.45569 1.746439e-50 0.1184947
#4 D E 0.066190196 -1.3767586 0.03597657 -38.26820 9.828853e-61 0.1952718
# R2 F.fv F.pv
#1 0.8590604 597.3333 1.738263e-43
#2 0.7603718 310.9670 3.636655e-32
#3 0.8985126 867.6375 1.746439e-50
#4 0.9372782 1464.4551 9.828853e-61

Aplique esto al Ejemplo 3 en la pregunta.

general_paired_simpleLM(dat[1], dat[2:5])
# LHS RHS alpha beta beta.se beta.tv beta.pv sig
#1 A B 0.112229318 1.0703491 0.07686011 13.92594 5.796437e-25 0.16446951
#2 A C 0.025628210 0.7114422 0.02789832 25.50126 4.731618e-45 0.10272687
#3 A D -0.009212582 0.3450939 0.01171768 29.45071 1.772671e-50 0.09044509
#4 A E 0.008650812 -0.4798639 0.01963404 -24.44040 1.738263e-43 0.10656866
# R2 F.fv F.pv
#1 0.6643051 193.9317 5.796437e-25
#2 0.8690390 650.3142 4.731618e-45
#3 0.8984818 867.3441 1.772671e-50
#4 0.8590604 597.3333 1.738263e-43

Incluso podemos simplemente hacer una regresión lineal simple entre dos variables:

general_paired_simpleLM(dat[1], dat[2])
# LHS RHS alpha beta beta.se beta.tv beta.pv sig
#1 A B 0.1122293 1.070349 0.07686011 13.92594 5.796437e-25 0.1644695
# R2 F.fv F.pv
#1 0.6643051 193.9317 5.796437e-25

Esto significa que la función simpleLM ahora está obsoleta.

Apéndice: Markdown (necesita compatibilidad con MathJax) para la imagen

Denote our variables by $x_1$, $x_2$, etc, a pairwise simple linear regression takes the form $$x_i = /alpha_{ij} + /beta_{ij}x_j$$ where $/alpha_{ij}$ and $/beta_{ij}$ is the intercept and the slope of $x_i /sim x_j$, respectively. We also denote $m_i$ and $v_i$ as the sample mean and **unscaled** sample variance of $x_i$. Here, the unscaled variance is just the sum of squares without dividing by sample size, that is $v_i = /sum_{k = 1}^n(x_{ik} - m_i)^2 = (/sum_{k = 1}^nx_{ik}^2) - n m_i^2$. We also denote $V_{ij}$ as the **unscaled** covariance between $x_i$ and $x_j$: $V_{ij} = /sum_{k = 1}^n(x_{ik} - m_i)(x_{jk} - m_j)$ = $(/sum_{k = 1}^nx_{ik}x_{jk}) - nm_im_j$.
Using the results for a simple linear regression given in [Function to calculate R2 (R-squared) in R](https://stackoverflow.com/a/40901487/4891738), we have $$/beta_{ij} = V_{ij} / / / v_j,/quad /alpha_{ij} = m_i - /beta_{ij}m_j,/quad r_{ij}^2 = V_{ij}^2 / / / (v_iv_j),$$ where $r_{ij}^2$ is the R-squared. Knowing $r_{ij}^2 = RSS_{ij} / / / TSS_{ij}$ where $RSS_{ij}$ and $TSS_{ij} = v_i$ are residual sum of squares and total sum of squares of $x_i /sim x_j$, we can derive $RSS_{ij}$ and residual standard error $/sigma_{ij}$ **without actually computing residuals**: $$RSS_{ij} = (1 - r_{ij}^2)v_i,/quad /sigma_{ij} = /sqrt{RSS_{ij} / / / (n - 2)}.$$
F-statistic $F_{ij}$ and associated p-value $p_{ij}^F$ can also be obtained from sum of squares: $$F_{ij} = /tfrac{(TSS_{ij} - RSS_{ij}) / / / 1}{RSS_{ij} / / / (n - 2)} = (n - 2) r_{ij}^2 / / / (1 - r_{ij}^2),/quad p_{ij}^F = 1 - /texttt{CDF_F}(F_{ij};/ 1,/ n - 2),$$ where $/texttt{CDF_F}$ denotes the CDF of F-distribution.
The only thing left is the standard error $e_{ij}$, t-statistic $t_{ij}$ and associated p-value $p_{ij}^t$ for $/beta_{ij}$, which are $$e_{ij} = /sigma_{ij} / / / /sqrt{v_i},/quad t_{ij} = /beta_{ij} / / / e_{ij},/quad p_{ij}^t = 2 * /texttt{CDF_t}(-|t_{ij}|; / n - 2),$$ where $/texttt{CDF_t}$ denotes the CDF of t-distribution.