broom::augment trop lent avec le modèle mixte du panel

Question

Je voudrais demander,

lors de l'élaboration de modèles mixtes de données de panel, puis de récupérer le résultat complet de la régression en utilisant broom::augment

il est très lent par rapport à d'autres modèles comme le modèle aléatoire, le modèle à effet fixe, etc.

exemple :

#Load packages
library(foreign)
library(plm)
library(stargazer)
library(dplyr)
library(tidyverse)
library(quantreg)

#Load in Wooldridge data on crime
crime <- read.dta("http://fmwww.bc.edu/ec-p/data/wooldridge/crime4.dta")

#Declare our data to be a panel data set
crime.p <- pdata.frame(crime,index=c("county","year"))
crime.p %>% head

panel1 <- plm(crmrte ~ polpc + prbconv + avgsen + density, data = crime.p, model = "within")
panel2 <- plm(crmrte ~ polpc + prbconv + avgsen + density, data = crime.p, model = "random")
panel3 <- plm(crmrte ~ polpc + prbconv + avgsen + density, data = crime.p, model = "pool")
panel4 <- rq(crmrte ~ polpc + prbconv + avgsen + density + factor(county)-1, data = crime.p)
panel5 <- lmer(crmrte ~ polpc + prbconv + avgsen + density + (1 | county), data = crime.p)

broom::augment(panel1)
broom::augment(panel2)
broom::augment(panel4)

# This one is very slow
broom::augment(panel5)

En cela, il est lent, mais avec de plus grands ensembles de données, il peut prendre jusqu'à 20 minutes, Existe-t-il un moyen de retrouver le résultat qui serait donné par broom::augment(mixed) mais beaucoup plus rapidement ?

Answer 1

panel5 n'est pas particulièrement lent sur ma machine (~0.01secondes). Vous pouvez voir où le plus de temps est pris en calculant chaque colonne séparément et en chronométrant :

J'ai créé une fonction personnalisée pour calculer chacun des éléments de l'élément augment résultat :

fe <- function(panel){
  f <- fixef(panel)
  int <- f[1]
  coef <- f[-1]
  d <- panel5@frame
  d1 <- d[names(f) %in% names(d)]

  int + rowSums(mapply(`*`, d1, coef))
}

augment2 <- function(panel) {
  d <- panel@frame
  p <- predict(panel)
  r <- resid(panel)
  hat <- hatvalues(panel)
  cook <- cooks.distance(panel)
  fe <- fe(panel)
  mu <- panel@resp$mu
  wr <- panel@resp$wtres
  cbind(d,fitted=p,redidual = r,hat,cooksd=cook,fixed=fe,mu,wtres=wr)
}

Si nous chronométrons les deux, ils sont assez similaires puisque tout ce que j'ai fait, c'est calculer les mêmes choses qui augment calcule :

microbenchmark::microbenchmark(
  augment = broom::augment(panel5),
  cust = augment2(panel5)
)
#Unit: milliseconds
#    expr       min       lq     mean   median       uq      max neval
# augment 12.332313 12.85186 16.83977 13.31599 15.69942 59.84045   100
#    cust  9.976812 10.38016 11.94285 10.66314 11.32606 43.50708   100

Si nous enlevons chacune des colonnes à tour de rôle de notre personnalisé augment2 nous pouvons voir que la plupart du temps est pris dans le calcul de la fonction hatvalues y cooks.distance . Si vous n'avez pas vraiment besoin de l'une ou des deux colonnes, cela accélérera considérablement votre calcul. En excluant les deux, nous pouvons constater une augmentation de 25x de la vitesse :

augment3 <- function(panel) {
  d <- panel@frame
  p <- predict(panel)
  r <- resid(panel)
  #hat <- hatvalues(panel)
  #cook <- cooks.distance(panel)
  fe <- fe(panel)
  mu <- panel@resp$mu
  wr <- panel@resp$wtres
  cbind(d,fitted=p,redidual = r,fixed=fe,mu,wtres=wr)
}

microbenchmark::microbenchmark(
  augment = broom::augment(panel5),
  cust = augment3(panel5)
)

#Unit: microseconds
#    expr       min        lq      mean    median         uq       max neval
# augment 12549.486 12924.666 17760.754 13341.712 15480.9555 87276.359   100
#    cust   406.818   474.119   694.698   532.005   586.5685  4702.398   100

En fait, le temps pris est dû à cooks.distance y hat.values qui, je l'imagine, ne peuvent être calculés de manière beaucoup plus efficace qu'ils ne le sont déjà.