superponer - Regresión lineal rápida por grupo

superponer graficas en r

Question

Tengo 500K usuarios y necesito calcular una regresión lineal (con intercepción) para cada uno de ellos.

Cada usuario tiene alrededor de 30 registros.

Lo intenté con dplyr y lm y esto es demasiado lento. Alrededor de 2 segundos por usuario.

  df%>%                       
      group_by(user_id, add =  FALSE) %>%
      do(lm = lm(Y ~ x, data = .)) %>%
      mutate(lm_b0 = summary(lm)$coeff[1],
             lm_b1 = summary(lm)$coeff[2]) %>%
      select(user_id, lm_b0, lm_b1) %>%
      ungroup()
    )

Intenté usar lm.fit que se sabe que es más rápido, pero no parece ser compatible con dplyr .

¿Hay una manera rápida de hacer una regresión lineal por grupo?

Answer 1

Puede intentarlo usando data.table como esta. Acabo de crear algunos datos de juguetes, pero me imagino que data.table daría alguna mejora. Es bastante rápido. Pero ese es un conjunto de datos bastante grande, por lo que tal vez comparen este enfoque con una muestra más pequeña para ver si la velocidad es mucho mejor. buena suerte.

    library(data.table)

    exp <- data.table(id = rep(c("a","b","c"), each = 20), x = rnorm(60,5,1.5), y = rnorm(60,2,.2))
    # edit: it might also help to set a key on id with such a large data-set
    # with the toy example it would make no diff of course
    exp <- setkey(exp,id)
    # the nuts and bolts of the data.table part of the answer
    result <- exp[, as.list(coef(lm(y ~ x))), by=id]
    result
       id (Intercept)            x
    1:  a    2.013548 -0.008175644
    2:  b    2.084167 -0.010023549
    3:  c    1.907410  0.015823088

Answer 2

Puedes usar las fórmulas básicas para calcular la pendiente y la regresión. lm hace muchas cosas innecesarias si lo único que te importa son esos dos números. Aquí utilizo data.table para la agregación, pero también puedes hacerlo en base R (o dplyr ):

system.time(
  res <- DT[, 
    {
      ux <- mean(x)
      uy <- mean(y)
      slope <- sum((x - ux) * (y - uy)) / sum((x - ux) ^ 2)
      list(slope=slope, intercept=uy - slope * ux)
    }, by=user.id
  ]
)

Produce para 500K usuarios ~ 30 obs cada uno (en segundos):

 user  system elapsed 
 7.35    0.00    7.36 

O unos 15 microsegundos por usuario . Y para confirmar esto está funcionando como se espera:

> summary(DT[user.id==89663, lm(y ~ x)])$coefficients
             Estimate Std. Error   t value  Pr(>|t|)
(Intercept) 0.1965844  0.2927617 0.6714826 0.5065868
x           0.2021210  0.5429594 0.3722580 0.7120808
> res[user.id == 89663]
   user.id    slope intercept
1:   89663 0.202121 0.1965844

Datos:

set.seed(1)
users <- 5e5
records <- 30
x <- runif(users * records)
DT <- data.table(
  x=x, y=x + runif(users * records) * 4 - 2, 
  user.id=sample(users, users * records, replace=T)
)

Answer 3

Si todo lo que desea es coeficientes, solo usaría user_id como un factor en la regresión. Usar el código de datos simulados de @ miles2know (aunque cambiar el nombre ya que un objeto que no sea exp() compartir ese nombre me parece extraño)

dat <- data.frame(id = rep(c("a","b","c"), each = 20),
                  x = rnorm(60,5,1.5),
                  y = rnorm(60,2,.2))

mod = lm(y ~ x:id + id + 0, data = dat)

No ajustamos ninguna intercepción global ( + 0 ) de modo que la intersección para cada id sea el coeficiente de id , y no x por sí misma, de modo que las interacciones x:id sean las pendientes para cada id :

coef(mod)
#      ida      idb      idc    x:ida    x:idb    x:idc 
# 1.779686 1.893582 1.946069 0.039625 0.033318 0.000353 

Entonces, para a nivel de id , el coeficiente de ida , 1.78, es la intersección y el coeficiente x:ida , 0.0396, es la pendiente.

Dejaré la recopilación de estos coeficientes en las columnas apropiadas de un marco de datos para usted ...

Esta solución debería ser muy rápida porque no tiene que lidiar con subconjuntos de marcos de datos. Probablemente podría acelerarse aún más con fastLm o tal.

Nota sobre la escalabilidad:

Acabo de probar esto en los datos de tamaño completo simulados de @nrussell y tuve problemas de asignación de memoria. Dependiendo de la cantidad de memoria que tenga, es posible que no funcione de una sola vez, pero probablemente podría hacerlo en lotes de ID de usuario. Alguna combinación de su respuesta y mi respuesta podría ser la más rápida en general, o nrussell podría ser más rápida, expandir el factor de identificación de usuario en miles de variables ficticias podría no ser computacionalmente eficiente, ya que he estado esperando más de un Un par de minutos ahora para una ejecución en solo 5000 ID de usuario.

Answer 4

Un ejemplo utilizando Rfast.

Suponiendo una respuesta única y variables predictoras 500K.

y <- rnorm(30)
x <- matrnorm(500*1000,30)
system.time( Rfast::univglms(y, x,"normal") )  ## 0.70 seconds

Suponiendo variables de respuesta 500K y una variable de predictor de singl.

system.time( Rfast::mvbetas(x,y) )  ## 0.60 seconds

Nota: Los tiempos anteriores disminuirán en el futuro cercano.

Answer 5

Actualización: Como lo señaló Dirk, mi enfoque original puede mejorarse en gran medida especificando x e Y directamente en lugar de usar la interfaz basada en fastLm de fastLm , lo que fastLm una sobrecarga de procesamiento (bastante significativa). Para la comparación, utilizando el conjunto de datos de tamaño completo original,

R> system.time({
  dt[,c("lm_b0", "lm_b1") := as.list(
    unname(fastLm(x, Y)$coefficients))
    ,by = "user_id"]
})
#  user  system elapsed 
#55.364   0.014  55.401 
##
R> system.time({
  dt[,c("lm_b0","lm_b1") := as.list(
    unname(fastLm(Y ~ x, data=.SD)$coefficients))
    ,by = "user_id"]
})
#   user  system elapsed 
#356.604   0.047 356.820

este simple cambio produce aproximadamente una aceleración de 6.5x .

[Enfoque original]

Probablemente haya algún margen de mejora, pero lo siguiente tomó aproximadamente 25 minutos en una máquina virtual de Linux (procesador de 2.6 GHz), ejecutando R de 64 bits:

library(data.table)
library(RcppArmadillo)
##
dt[
  ,c("lm_b0","lm_b1") := as.list(
    unname(fastLm(Y ~ x, data=.SD)$coefficients)),
  by=user_id]
##
R> dt[c(1:2, 31:32, 61:62),]
   user_id   x         Y     lm_b0    lm_b1
1:       1 1.0 1674.8316 -202.0066 744.6252
2:       1 1.5  369.8608 -202.0066 744.6252
3:       2 1.0  463.7460 -144.2961 374.1995
4:       2 1.5  412.7422 -144.2961 374.1995
5:       3 1.0  513.0996  217.6442 261.0022
6:       3 1.5 1140.2766  217.6442 261.0022

Datos:

dt <- data.table(
  user_id = rep(1:500000,each=30))
##
dt[, x := seq(1, by=.5, length.out=30), by = user_id]
dt[, Y := 1000*runif(1)*x, by = user_id]
dt[, Y := Y + rnorm(
  30, 
  mean = sample(c(-.05,0,0.5)*mean(Y),1), 
  sd = mean(Y)*.25), 
  by = user_id]