multiple - plot en r



Regresión lineal rápida por grupo (4)

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

https://ffff65535.com

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?


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)
)

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.


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]




lm