R: избегая summary.plm

Question

Я использую R для запуска симуляции Монте-Карло, изучая производительность панельных оценщиков данных. Поскольку я буду выполнять большое количество испытаний, мне нужно получить хотя бы приличную производительность из своего кода.

Использование Rprof в 10 испытаниях моего моделирования показывает, что значительная часть времени затрачивается на вызовы summary.plm. Первые несколько строк Rprofsummary приведены ниже:

$by.total
                            total.time total.pct self.time self.pct
"trial"                          54.48     100.0      0.00      0.0
"coefs"                          53.90      98.9      0.06      0.1
"model.matrix"                   36.72      67.4      0.10      0.2
"model.matrix.pFormula"          35.98      66.0      0.06      0.1
"summary"                        33.82      62.1      0.00      0.0
"summary.plm"                    33.80      62.0      0.08      0.1
"r.squared"                      29.00      53.2      0.02      0.0
"FUN"                            24.84      45.6      7.52     13.8

Я звоню summary в своем коде, потому что мне нужно получить стандартные ошибки оценок коэффициентов, а также сами коэффициенты (которые я мог бы получить только из объекта plm). Мой звонок выглядит как

regression <- plm(g ~ y0 + Xit, data=panel_data, model=model, index=c("country","period"))

coefficients_estimated <- summary(regression)$coefficients[,"Estimate"]
ses_estimated <- summary(regression)$coefficients[,"Std. Error"]

У меня есть ноющее чувство, что это огромная трата процессорного времени, но я не знаю достаточно о том, как R делает вещи, чтобы избежать вызова итогов. Я был бы признателен за любую информацию о том, что происходит за кулисами здесь, или о каком-либо способе сокращения времени, необходимого для того, чтобы это оправдало себя.

Answer 1

Вам просто нужно заглянуть внутрь plm:::summary.plm, чтобы увидеть, что он делает. Когда вы это сделаете, вы увидите, что две строки, вызывающие summary() в вашей модели, могут быть заменены на:

coefficients_estimated <- coef(regression)
ses_estimated <- sqrt(diag(vcov(regression)))

Например:

require(plm)
data("Produc", package = "plm")
zz <- plm(log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp, 
          data = Produc, index = c("state","year"))

summary(zz) дает:

> summary(zz)
Oneway (individual) effect Within Model

....

Coefficients :
             Estimate  Std. Error t-value  Pr(>|t|)    
log(pcap) -0.02614965  0.02900158 -0.9017    0.3675    
log(pc)    0.29200693  0.02511967 11.6246 < 2.2e-16 ***
log(emp)   0.76815947  0.03009174 25.5273 < 2.2e-16 ***
unemp     -0.00529774  0.00098873 -5.3582 1.114e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
....

и две строки, которые я показал, возвращаются для zz:

> coef(zz)
   log(pcap)      log(pc)     log(emp)        unemp 
-0.026149654  0.292006925  0.768159473 -0.005297741 
> sqrt(diag(vcov(zz)))
   log(pcap)      log(pc)     log(emp)        unemp 
0.0290015755 0.0251196728 0.0300917394 0.0009887257

Вы на самом деле не предоставляете достаточно информации (например, свой код симуляции или полный вывод из Rprof()), чтобы сказать, поможет ли это - это, конечно, не похоже, что в * 1018 тратится огромное количество времени. ; FUN намного дороже, чем все, что вы показываете, и из элементов, которые вы показываете, r.squared() - единственный, который появляется в plm:::summary.plm(), и, похоже, совсем не занимает времени.

Итак, заметно ли ускорение вышеизложенного, еще неизвестно.

Answer 2

Просто используйте coeftest(zz).coeftest находится в пакете lmtest;это даст вам коэффициенты и стандартные ошибки от plm объектов намного быстрее, чем summary.plm.

Answer 3

Если вы хотите продвинуться дальше, взгляните на фактический код функции plm:::plm. Вы заметите, что перед окончательным вызовом plm:::plm.fit есть много проверок аргументов.), перейдите прямо к plm.fit.

Один последний момент.Вы упоминаете, что ваша проблема - симуляция Монте-Карло.Можете ли вы использовать параллельные вычисления для увеличения скорости?